本文对基频提取算法 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,\cdots ,$ 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. 总结

结合公式和开源代码的实现可以发现，差分函数的快速实现以及抛物线插值的简化版本这块是值得学习和吸收的，其余部分基本按照公式实现，相对而言比较简单。