【基频提取算法-YIN】

浩哥依然

已于 2024-04-24 10:20:14 修改

阅读量1k

点赞数 7

分类专栏：基频提取文章标签：音频信号处理基频检测算法 YIN YIN 算法代码学习

于 2024-02-27 19:59:23 首次发布

本文链接：https://blog.csdn.net/wjrenxinlei/article/details/136330690

版权

基频提取专栏收录该内容

4 篇文章 1 订阅

订阅专栏

本文对基频提取算法 YIN 做以介绍。如有表述不当之处欢迎批评指正。欢迎任何形式的转载，但请务必注明出处。

文章目录

1. 引言
2. YIN 各模块代码讲解
3. 总结

1. 引言

之前的好几个项目都用到了基频提取算法 PYIN，实验结果显示其准确度较高，能满足项目需求。PYIN算法是在 YIN 算法的基础上改进而来，因此，有必要深入研究学习下 YIN 算法。

本文结合开源代码详细介绍了 YIN 算法各个模块的实现细节，如果不了解该算法，可以先阅读笔者另一篇关于该算法论文的博客【论文笔记之 YIN】。开源代码的地址为：https://code.soundsoftware.ac.uk/projects/pyin/files。该开源代码并不是 YIN 的原作者实现的，而是 PYIN 的作者实现的，网页上共开源了三个版本，笔者以最早的那个版本为主进行介绍。

2. YIN 各模块代码讲解

本章将针对 YIN 算法的各个模块，结合开源代码来深入学习。

2.1. 差分函数的实现

代码实现的差分函数是论文中的公式 $(7)$ ，如下：
$\begin{align} d_t(\tau) = r_t(0) + r_{t+\tau}(0) - 2r_t(\tau)\tag{7} \end{align}$

该差分函数由三项构成，都可以由公式 $(1)$ 计算得到：
$\begin{align} r_t(\tau) = \sum_{j=t+1}^{t+W} x_j x_{j+\tau} \tag{1} \end{align}$

其中，第一项可以根据公式 $(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(0)$ 和 $r_{t+1}(0)$ 的 $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$ 的取值。窗长应该满足 $\leq$ 音频帧长，这块取最大值 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]; // real
    yinStyleACFImag[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 part
    yinBuffer[j] = powerTerms[0] + powerTerms[j] - 2 * audioTransformedReal[j+yinBufferSize-1];
}

2.2. 累积均值归一化差分函数的实现

这块直接根据公式 $(8)$ 的定义来实现:
$\begin{align} d_t^{'}(\tau) = \left\{ \begin{array}{lc} 1, & \text{if} \; \tau = 0, \\ d_t(\tau)/[(1/\tau)\sum_{j=1}^{\tau}{d_t(j)}], & \text{otherwise}. \\ \end{array} \right. \end{align} \tag{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 TarsosDSP
    tau = 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 implementation
    if (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 -1
        betterTau = tau + (s2 - s0) / (2 * (2 * s1 - s2 - s0));
        
        // std::cerr << tau << " --> " << betterTau << std::endl;
        
    }
    return betterTau;
}

代码的前半部分都是在处理边界情况，核心是在最后的 else 部分，其对应的公式为：
$\begin{align} x_{peak} = x_{0} + \frac{1}{2} \frac{y_{+1}-y_{-1}}{2y_{0} - y_{-1} - y_{+1}} \notag \end{align}$