本文對基頻提取算法 YIN 做以介紹。如有表述不當之處歡迎批評指正。歡迎任何形式的轉載，但請務必注明出處。

1. 引言

之前的好幾個項目都用到了基頻提取算法 Pyin，實驗結果顯示其準確度較高，能滿足項目需求。Pyin算法是在 YIN 算法的基礎上改進而來，因此，有必要深入研究學習下 YIN 算法。

本文結合開源代碼詳細介紹了 YIN 算法各個模塊的實現細節，如果不了解該算法，可以先閱讀筆者另一篇關于該算法論文的博客 https://blog.csdn.net/wjrenxinlei/article/details/136306282。開源代碼的地址為：https://code.soundsoftware.ac.uk/projects/pyin/files

2. YIN 各模塊代碼講解

本章將針對 YIN 算法的各個模塊，結合開源代碼來深入學習。

2.1. 差分函數的實現

代碼實現的差分函數是論文中的公式 $(7)$，如下：
$$\begin{array}{r}{d}_t(\tau )=r_t(0)+r_{t+\tau }(0)?2r_t(\tau )\end{array}\text{\hspace{1em}(7)?}$$

該差分函數由三項構成，都可以由公式 $(1)$ 計算得到：
$$\begin{array}{r}{r}_t(\tau )=\sum _{j=t+1}^{t+W}{x}_j{x}_{j+\tau }\end{array}\text{\hspace{1em}(1)?}$$

其中，第一項可以根據公式 $(1)$ 直接計算，開源代碼實現如下：

// POWER TERM CALCULATION
// ... for the power terms in equation (7) in the Yin paper
powerTerms[0] = 0.0;
for (size_t j = 0; j < yinBufferSize; ++j) {powerTerms[0] += in[j] * in[j];
}

第二項可以遞歸計算，因為 $r_t(\tau )$ 和 $r_t(\tau +1)$ 的 $W$ 個求和項中，只有一項是不一樣的，其余 $W?1$ 項是完全相同的，開源代碼實現如下：

// now iteratively calculate all others (saves a few multiplications)
for (size_t tau = 1; tau < yinBufferSize; ++tau) {powerTerms[tau] = powerTerms[tau-1] - in[tau-1] * in[tau-1] + in[tau+yinBufferSize] * in[tau+yinBufferSize];  
}
// 筆者注：源碼實現有誤，上述等號右邊最后一項應該是 in[tau+yinBufferSize-1] * in[tau+yinBufferSize-1];

這兒需要注意的是，窗長 $W$ 和 延遲 $\tau$ 的取值。窗長應該滿足 $2W\le$ 音頻幀長，這塊取最大值 yinBufferSize（該值是音頻幀長的一半，具體可見源碼）。窗長確定好之后，$\tau$ 能取到的最大值也就隨之確定了，其能取到的最大值應該是 幀長 - 窗長 = yinBufferSize，這塊可以看到源碼實現中 $\tau$ 的最大取值為 yinBuferSize - 1。即在源碼實現中 $\tau =0,1,?,$ yinBuferSize - 1。

第三項是序列的自（互）相關，這塊的計算用到了數字信號處理中的一個基本知識點，即可以使用 FFT 來加速計算自（互）相關。下面詳細描述下該計算過程。
首先，得到參與相關運算的兩個序列。第一個序列是輸入給算法的音頻幀。第二個序列是由該音頻幀的前 yinBufferSize 個采樣點（這是由于有窗長的限制）得到的。首先，對這 yinBufferSize 個采樣點先做逆序操作（這是相關和卷積的主要區別，如果不做逆序操作，那么接下來的步驟算出來的將會是卷積），然后對其補零（使其和輸入給算法的音頻幀長度一樣）得到第二個序列。
接著，對上述得到的兩個序列分別做 FFT 運算。
然后，將上述得到的兩個 FFT 序列相乘，并對相乘得到的結果做 IFFT 運算。
最后，對上述得到的結果取下標（從 0 開始）為 [ yinBufferSize - 1, 2 * yinBufferSize - 2] 的結果。

開源代碼實現如下：

// YIN-STYLE AUTOCORRELATION via FFT
// 1. data
esfft(uint(frameSize), false, in, nullImag, audioTransformedReal, audioTransformedImag);// 2. half of the data, disguised as a convolution kernel
for (size_t j = 0; j < yinBufferSize; ++j) {kernel[j] = in[yinBufferSize-1-j];kernel[j+yinBufferSize] = 0;
}
esfft(uint(frameSize), false, kernel, nullImag, kernelTransformedReal, kernelTransformedImag);// 3. convolution via complex multiplication -- written into
for (size_t j = 0; j < frameSize; ++j) {yinStyleACFReal[j] = audioTransformedReal[j]*kernelTransformedReal[j] - audioTransformedImag[j]*kernelTransformedImag[j]; // realyinStyleACFImag[j] = audioTransformedReal[j]*kernelTransformedImag[j] + audioTransformedImag[j]*kernelTransformedReal[j]; // imaginary
}
esfft(uint(frameSize), true, yinStyleACFReal, yinStyleACFImag, audioTransformedReal, audioTransformedImag);

其中，frameSize = yinBufferSize * 2 是音頻幀長。in 是輸入給算法的音頻幀。

最后，將上面三項相加，得到差分函數公式 $(7)$ 的最終結果。開源代碼實現如下：

// CALCULATION OF difference function
// ... according to (7) in the Yin paper.
for (size_t j = 0; j < yinBufferSize; ++j) {// taking only the real partyinBuffer[j] = powerTerms[0] + powerTerms[j] - 2 * audioTransformedReal[j+yinBufferSize-1];
}

2.2. 累積均值歸一化差分函數的實現

這塊直接根據公式 $(8)$ 的定義來實現:
$$\begin{array}{r}{d}_t^{{}^{\prime }}(\tau )=\{\begin{array}{lc}1,& \text{if}\tau =0,\\ d_t(\tau )/[(1/\tau ){\sum }_{j=1}^{\tau }{d}_t(j)],& \text{otherwise}.\end{array}\end{array}\text{\hspace{1em}(8)}$$

開源代碼實現如下：

void cumulativeDifference(double *yinBuffer, const size_t yinBufferSize)
{    size_t tau;yinBuffer[0] = 1;double runningSum = 0;for (tau = 1; tau < yinBufferSize; ++tau) {runningSum += yinBuffer[tau];if (runningSum == 0){yinBuffer[tau] = 1;} else {yinBuffer[tau] *= tau / runningSum;}}    
}

2.3. 絕對閾值

這一步的目的是，找到第一個小于閾值的谷值所對應的 $\tau$，如果都大于閾值，那么找到最小的谷值所對應的 $\tau$。開源代碼實現如下：

int absoluteThreshold(const double *yinBuffer, const size_t yinBufferSize, const double thresh)
{size_t tau;size_t minTau = 0;double minVal = 1000.;// using Joren Six's "loop construct" from TarsosDSPtau = 2;while (tau < yinBufferSize){if (yinBuffer[tau] < thresh){while (tau+1 < yinBufferSize && yinBuffer[tau+1] < yinBuffer[tau]){++tau;}return tau;} else {if (yinBuffer[tau] < minVal){minVal = yinBuffer[tau];minTau = tau;}}++tau;}if (minTau > 0){return -minTau;}return 0;
}

其中，第二個 while 循環是為了找到谷值，else 部分是為了找到最小的谷值所對應的 $\tau$.

2.4. 拋物線插值

插值是為了找到更精確的周期估計。開源代碼中是在相鄰的三組樣本上，使用簡化的拋物線插值公式來尋找，實現如下：

double parabolicInterpolation(const double *yinBuffer, const size_t tau, const size_t yinBufferSize) 
{// this is taken almost literally from Joren Six's Java implementationif (tau == yinBufferSize) // not valid anyway.{return static_cast<double>(tau);}double betterTau = 0.0;size_t x0;size_t x2;if (tau < 1) {x0 = tau;} else {x0 = tau - 1;}if (tau + 1 < yinBufferSize) {x2 = tau + 1;} else {x2 = tau;}if (x0 == tau) {if (yinBuffer[tau] <= yinBuffer[x2]) {betterTau = tau;} else {betterTau = x2;}} else if (x2 == tau) {if (yinBuffer[tau] <= yinBuffer[x0]) {betterTau = tau;} else {betterTau = x0;}} else {float s0, s1, s2;s0 = yinBuffer[x0];s1 = yinBuffer[tau];s2 = yinBuffer[x2];// fixed AUBIO implementation, thanks to Karl Helgason:// (2.0f * s1 - s2 - s0) was incorrectly multiplied with -1betterTau = tau + (s2 - s0) / (2 * (2 * s1 - s2 - s0));// std::cerr << tau << " --> " << betterTau << std::endl;}return betterTau;
}

代碼的前半部分都是在處理邊界情況，核心是在最后的 else 部分，其對應的公式為：
$$\begin{array}{r}{x}_{peak}=x_0+\frac{1}{2}\frac{y_{+1}?y_{?1}}{2y_0?y_{?1}?y_{+1}}\end{array}$$

2.5. 最優局部估計

開源代碼中暫且沒有該步驟的實現。

3. 總結

結合公式和開源代碼的實現可以發現，差分函數的快速實現以及拋物線插值的簡化版本這塊是值得學習和吸收的，其余部分基本按照公式實現，相對而言比較簡單。