长春火车站防疫政策,asp化妆品网站 后台,自己怎么设计网页,电销卡代理加盟本文对基频提取算法 YIN 做以介绍。如有表述不当之处欢迎批评指正。欢迎任何形式的转载#xff0c;但请务必注明出处。 文章目录 1. 引言2. YIN 各模块代码讲解2.1. 差分函数的实现2.2. 累积均值归一化差分函数的实现2.3. 绝对阈值2.4. 抛物线插值2.5. 最优局部估计 3. 总结 1…本文对基频提取算法 YIN 做以介绍。如有表述不当之处欢迎批评指正。欢迎任何形式的转载但请务必注明出处。 文章目录 1. 引言2. YIN 各模块代码讲解2.1. 差分函数的实现2.2. 累积均值归一化差分函数的实现2.3. 绝对阈值2.4. 抛物线插值2.5. 最优局部估计 3. 总结 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 ) (7) (7)如下 d t ( τ ) r t ( 0 ) r t τ ( 0 ) − 2 r t ( τ ) \begin{align} d_t(\tau) r_t(0) r_{t\tau}(0) - 2r_t(\tau)\tag{7} \end{align} dt(τ)rt(0)rtτ(0)−2rt(τ)(7)
该差分函数由三项构成都可以由公式 ( 1 ) (1) (1) 计算得到 r t ( τ ) ∑ j t 1 t W x j x j τ \begin{align} r_t(\tau) \sum_{jt1}^{tW} x_j x_{j\tau} \tag{1} \end{align} rt(τ)jt1∑tWxjxjτ(1)
其中第一项可以根据公式 ( 1 ) (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 ( τ ) r_t(\tau) rt(τ) 和 r t ( τ 1 ) r_t(\tau1) rt(τ1) 的 W W W 个求和项中只有一项是不一样的其余 W − 1 W-1 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[tauyinBufferSize] * in[tauyinBufferSize];
}
// 笔者注源码实现有误上述等号右边最后一项应该是 in[tauyinBufferSize-1] * in[tauyinBufferSize-1];这儿需要注意的是窗长 W W W 和 延迟 τ \tau τ 的取值。窗长应该满足 2 W ≤ 2W \leq 2W≤ 音频帧长这块取最大值 yinBufferSize该值是音频帧长的一半具体可见源码。窗长确定好之后 τ \tau τ 能取到的最大值也就随之确定了其能取到的最大值应该是 帧长 - 窗长 yinBufferSize这块可以看到源码实现中 τ \tau τ 的最大取值为 yinBuferSize - 1。即在源码实现中 τ 0 , 1 , ⋯ , \tau 0, 1, \cdots, τ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[jyinBufferSize] 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 ) (7) (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[jyinBufferSize-1];
}2.2. 累积均值归一化差分函数的实现
这块直接根据公式 ( 8 ) (8) (8) 的定义来实现: d t ′ ( τ ) { 1 , if τ 0 , d t ( τ ) / [ ( 1 / τ ) ∑ j 1 τ d t ( j ) ] , otherwise . (8) \begin{align} d_t^{}(\tau) \left\{ \begin{array}{lc} 1, \text{if} \; \tau 0, \\ d_t(\tau)/[(1/\tau)\sum_{j1}^{\tau}{d_t(j)}], \text{otherwise}. \\ \end{array} \right. \end{align} \tag{8} dt′(τ){1,dt(τ)/[(1/τ)∑j1τdt(j)],ifτ0,otherwise.(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 Sixs loop construct from TarsosDSPtau 2;while (tau yinBufferSize){if (yinBuffer[tau] thresh){while (tau1 yinBufferSize yinBuffer[tau1] 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 Sixs Java implementationif (tau yinBufferSize) // not valid anyway.{return static_castdouble(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 部分其对应的公式为 x p e a k x 0 1 2 y 1 − y − 1 2 y 0 − y − 1 − y 1 \begin{align} x_{peak} x_{0} \frac{1}{2} \frac{y_{1}-y_{-1}}{2y_{0} - y_{-1} - y_{1}} \notag \end{align} xpeakx0212y0−y−1−y1y1−y−1
2.5. 最优局部估计
开源代码中暂且没有该步骤的实现。
3. 总结
结合公式和开源代码的实现可以发现差分函数的快速实现以及抛物线插值的简化版本这块是值得学习和吸收的其余部分基本按照公式实现相对而言比较简单。
