用于语音和音乐的基频(音高)估计算法：YIN

周期信号的基频估计（fundamental frequency estimation）在许多应用中都起着重要的作用，例如：在汉语语音识别中，由于不同的声调对应不同的字词，所以音高是非常重要的特征，另外在识别说话人情感的时候，说话的语调也含有显著的情绪倾向；还有在数字音乐应用中，包括乐谱转录、哼唱识别、K歌软件等中都有用武之地。

音高识别算法非常之多，包括在近些年发展起来的深度学习领域中也出现了音高估计模型的身影，下面介绍一下比较通用的YIN算法，并罗列出用于阅读代码的主要几大步骤，并增加一些个人演算；更详细的见大佬的原论文《YIN, a fundamental frequency estimator for speech and music》。进一步，工程中可以使用更加稳定的probabilisitc YIN算法，我们将在以后的文章进行介绍。

第一步：自相关函数和差函数

若信号有周期，则将信号在时间轴上平移一个周期将与自身重合。因而自相关函数（见下式）（autocorrelation function）将在周期的整倍数时刻取最大值，因而自相关函数显然可以用来做基频估计。

$\begin {aligned}&\boldsymbol r_t[\tau] = \sum_{n=t+1}^{t+W} x[n]\cdot x[n+\tau] \\\end{aligned}$

其中： $W$ 是执行一次自相关运算的窗口大小； $x[n]$ 是时域信号。

自相关函数用来估计基频的方法是，穷举一定范围内所有的 $\tau$ 值，并且挑选最高的非零 $\tau$ 作为基本周期。换言之，自相关方法倾向于挑选较大的周期，即倾向于较小的基频。当 $\tau$ 的取值上限较大时，很容易发生低估基频的错误。

还有一种自相关函数的定义：

$\begin {aligned}&\boldsymbol r_t^{\prime}[\tau] = \sum_{n=t+1}^{t+W-\tau} x[n]\cdot x[n+\tau] \end{aligned}$

该定义中，随着 $\tau$ 的增加，被加项越来越少，将会形成振幅逐渐衰减的自相关函数形状。因而在新的定义中，尽管 $\tau$ 有较大的取值上界，理论上也不会容易出现低估基频的错误。但是当 $\tau$ 的取值下限足够接近0时，由于 $\boldsymbol r_t^{\prime}[\tau]$ 振幅逐渐衰减的性质，会倾向选择接近0的基本周期，因而容易造成高估基频的错误。

该论文首先就提出了若干论点，说明自相关函数在实际应用中不如新提出来的差函数（difference function）好用，其为：

$\begin {aligned}&\boldsymbol d_t[\tau] = \sum_{n=t+1}^{t+W} (x[n] - x[n+\tau])^2 \end{aligned}$

换成差函数后，基本周期是使得该函数取最小值的那个 $\tau$ 值。并且差函数和自相关函数之间有如下关系：（其中，对中间的那个式子做个简单的变量变换就可以了）

$\begin {aligned}\boldsymbol d_t[\tau] &= \sum_{n=t+1}^{t+W} (x[n] - x[n+\tau])^2 \\&= \sum_{n=t+1}^{t+W} x[n]^2 + \sum_{n=t+1}^{t+W} x[n+\tau]^2 -2 \sum_{n=t+1}^{t+W} x[n]\cdot x[n+\tau] \\&=\boldsymbol r_{t}[0] + \boldsymbol r_{t+\tau}[0] - 2\boldsymbol r_{t}[\tau]\end{aligned}$

此时，若 $\boldsymbol r_{t}[0] + \boldsymbol r_{t+\tau}[0]$ 为常量，则最大化 $\boldsymbol r_{t}[\tau]$ 就等价于最小化 $\boldsymbol d_{t}[\tau]$ 。但是其中 $\boldsymbol r_{t+\tau}[0]$ 这一项是和 $\tau$ 有关的项，所以这两个准则的确会导致不同的估计。仅仅由于这一个改变就在测试数据上，将误差从10%降低到了1.95%。作者给的解释是自相关函数对信号的幅度变化非常敏感；当信号振幅随着时间增大时，自相关函数在基本周期的一定范围内的倍数时刻 $kT_0,(0<k<k_0)$ 的值会呈现增加的现象，因此由于自相关函数选择最大周期对应的那个频率作为识别的基础频率，因此会高估基本周期，从而出现低估基础频率的错误。

而作者提出的差函数并没有这类问题，因为可以想象一下一个正弦函数的周期保持不变，但是幅度变化，但是不管是逐渐变小还是逐渐变大，都是在平移一个基本周期时出现的差函数低洼比其他更大周期出现的差函数低洼更低，甚至当振幅忽高忽低时，由于信号的规律性，差函数会在所有的 $\tau$ 上都变大，所以对变化的振幅会不那么敏感。

第二步：累计均值规范化差函数

直接使用差函数会出现两个问题：1）和自相关函数类似，由于差函数仍然在 $\tau$ 为0的时候给出最小值0，并且由于实际音频数据不是完美周期性，导致差函数在真实的基本周期处的取值仍然小于在0处的取值，我们可以通过设置一个基本周期的取值下限还改善该问题。2）当在第一共振峰（first formant）处出现一个较强的共振，仍然会在2倍基频处（即1/2倍基本周期处）出现一些次级的低洼，这些次级低洼可能比基本周期对应的低洼还要更低；而对于这第二个问题，并不能通过增加一个取值下限来改善，因为基本周期 $F_0$ 和第一共振峰 $F_1$ 的取值范围是重合的。所以作者进一步提出了“累计均值规范化差函数”（cumulative mean normalised difference function），如下式所示。

$\begin{aligned}\boldsymbol d_t^{\prime}[\tau] = \begin {cases}1, &\text{if } \tau = 0 \\ \boldsymbol d[\tau]/(\frac{1}{\tau} \sum_{j=1}^{\tau} \boldsymbol d_t[j]), & \text{otherwise} \end{cases} \end{aligned}$

直接设置差函数在0处取1，并且用直到 $\tau$ 时刻的差函数相应取值的均值（累积均值）来规范化；可以想象，当信号原始的差函数在低于基本周期的地方出现更深的次级低洼时，累积均值同样也会很小，因而除以很小的数相当于放大原始的差函数的值，因而更好地降低次级共振峰引起的错误率。并且直到 $\tau$ 达到真正的基本周期处出现很深的低洼。见下图，图a）是原始的差函数形状，图b）是均值规范化后的。新的差函数显然可以抑制低估周期（高估基频）的现象，并且新的差函数并不需要设置显式的基频上限。

YIN论文插图

第三步：绝对阈值

实际应用中很容易出现差函数在更大的 $\tau$ 上出现更深的低洼，例如上图b），若图中第二个很深的低洼处于候选区域内，则很容易选择第二个低洼对应的基频，而出现低估基频的错误。为了应对这个问题，我们可以考虑设置一个绝对的阈值（如图中虚线），找到候选区域内所有低于虚线的那些低洼所对应的频率，将最小的 $\tau$ 对应的频率（更高的频率）作为基频返回；而我们之所以可以放心地这么做，是因为在上一步我们用“均值规范化差函数”大大降低了高估基频的错误。但是因为阈值是固定的，一旦在阈值下方没有找到仍和候选值，则算法直接返回全局的最小值。

通过阈值，作者将测试数据上的错误率从1.69%降低到0.78%；大大降低了低估错误率，仅仅轻微地增加了高估错误率，而这是很容易理解并接受的。那么如何选择该阈值呢？作者给了阈值的意义一个解释：“一个近似周期信号功率中的非周期性功率所占的比例”。作者给出了下面的恒等式：

$\frac{1}{2W} \sum_{j=t+1}^{t+W}(x[j]^2 + x[j+T_0]^2) = \frac{1}{4W}\sum_{j=t+1}^{t+W}(x[j] + x[j+T_0])^2 + \frac{1}{4W}\sum_{j=t+1}^{t+W}(x[j] - x[j+T_0])^2$

等式左边可以近似视为信号的功率；而右边作为信号功率的两部分，第二部分可以用来度量信号的非周期性：显然若信号是周期的，并且基本周期为 $T_0$ ，那么该项为0，表示信号是完全周期的；当其非零时，度量了非周期性强弱；并且若再叠加一个周期为 $T_0$ 的完美周期信号，则显然不会改变第二项的取值；因而第二项的确是度量周期信号中非周期性的部分功率。为了看清“绝对阈值”对应物理上的非周期性功率和总功率之比率的解释，在 $\tau=T_0$ 时我们做一个简单的推理：“累积均值规范化差函数”的分子项为： $\boldsymbol d_t[T_0] = \sum_{j=t+1}^{t+W}(x[j] - x[j+T_0])^2$ ，而分母项为：

$\begin {aligned}\frac{1}{T_0} \sum_{j=1}^{T_0} \boldsymbol d_t[j] &= \frac{1}{T_0} \sum_{j=1}^{T_0} \sum_{i=t+1}^{t+W}(x[i] - x[i+j])^2 \\&=\sum_{i=t+1}^{t+W} \left ( \frac{1}{T_0} \sum_{j=1}^{T_0}(x[i]^2 + x[i+j]^2) - 2\frac{1}{T_0}\sum_{j=1}^{T_0}x[i]x[i+j] \right) \\&=\sum_{i=t+1}^{t+W} \left (x[i]^2 + \frac{1}{T_0} \sum_{j=1}^{T_0}x[i+j]^2 - 2\frac{x[i]}{T_0}\sum_{j=1}^{T_0}x[i+j] \right)\end {aligned}$

其中第三项中的 $\sum_{j=1}^{T_0}x[i+j]$ 表示近似周期函数在一个周期内积分，其结果可以近似视为0，所以第三项可以忽略。所以只剩下 $\sum_{i=t+1}^{t+W} x[i]^2 + \sum_{i=t+1}^{t+W}\frac{1}{T_0} \sum_{j=1}^{T_0}x[i+j]^2$ 两项，设 ${\rm P}_x$ 为信号的平均功率，则可将第2项近似为： $\sum_{i=t+1}^{t+W}\frac{1}{T_0} \sum_{j=1}^{T_0}x[i+j]^2 = \sum_{i=t+1}^{t+W} {\rm P}_x = W{\rm P}_x$ ，第1项近似为： $\sum_{i=t+1}^{t+W} x[i]^2 = W \frac{1}{W}\sum_{i=t+1}^{t+W} x[i]^2 = W{\rm P}_x$ ；所以：

$\boldsymbol d_t^{\prime}[T_0] = \frac{\boldsymbol d_t[T_0]}{2W{\rm p}_x} = \frac{\frac{1}{4W}\boldsymbol d_t[T_0]}{\frac{1}{2}{\rm p}_x}$

这是一个信号的非周期性功率与信号的平均功率的比率形式，但是我这里推导的分母是1/2的平均功率，而不是论文里说的2倍的平均功率，所以是作者的1/4。我暂时没有继续推敲了，不过这的确是作者所述的功率比率形式，这说明了绝对阈值的物理意义。

所以对这部分总结一下就是：“累计均值归一化差函数”的阈值对应的是信号的非周期性程度，我们的目的是估计周期信号的基频，所以当我们将其设置为很小的值时，我们就是期望信号是周期信号的置信度很高，因而我们对我们估计出来的基频抱有较大的信心。

第四步：抛物线插值

若信号的基本周期是刚好是采样周期的整数倍，则我们使用前面的方法找到正确的差函数低洼，对应的就是正确的基频。但是一般不会满足这个条件，所以使用样本序列在时域进行基频估计显然是不准确的，因此作者提出了抛物线插值：在函数 $\boldsymbol d_t^{\prime}[\tau]$ 的每个极小值点，使用左右两个相邻点加上极小值点本身三个点进行抛物线插值，将插值得到的抛物线的最小值对应的横坐标作为候选者放到最后的基频选择列表里。

事实上，由于使用了“累积均值规范化差函数”，而规范化后的差函数和原始的差函数对极小值点横坐标的估计有一定的偏差。所以，实际会使用原始的差函数对应的三点进行插值。

设 $x[n]$ 是无周期的离散时间信号， $X(\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}$ 为其傅立叶变换，且 $X(\omega)$ 为周期 $2\pi$ 的周期函数。所以信号 $x[n]$ 的功率谱的傅立叶变换为：

$\begin{aligned}&\frac{1}{2\pi}\int_{-\pi}^{\pi} |X(\omega)|^2 e^{-j\omega n} {\rm d}\omega \\&= \frac{1}{2\pi}\int_{-\pi}^{\pi} X(\omega) X(\omega)^\ast e^{-j\omega n} {\rm d}\omega \\&= \frac{1}{2\pi}\int_{-\pi}^{\pi} X(\omega) [\sum_{m=-\infty}^{\infty} x[m] e^{j\omega m}] e^{-j\omega n} {\rm d}\omega \\&= \sum_{m=-\infty}^{\infty} x[m] \frac{1}{2\pi}\int_{-\pi}^{\pi} X(\omega) e^{j\omega (m-n)} \\&= \sum_{m=-\infty}^{\infty} x[m] x[m-n] \\&=r_{xx}[n]\end{aligned}$

其中，利用了 $x[m]$ 为实数序列，并交换了积分和求和次序。

所以我们得到了功率谱密度的傅立叶变换是信号 $x[n]$ 的自相关函数 $r_{xx}[n]$ ，并且再次利用实信号性质，得到 $|X(-\omega)|^2 = |X(\omega)^{\ast}|^2 = |X(\omega)|^2$ ，所以功率谱密度函数是偶函数，因而其傅立叶变换中正弦变换部分消失了，又假设 $X(\omega)$ 是带限的，所以功率谱密度的傅立叶变换仅由一组余弦函数就能表达。

又每个余弦函数的在0点的泰勒展开式仅仅包括偶次项，即常数项、二次项、4次项等等，所以总体上 $r_{xx}[n]$ 在0点可以由一组偶次项的多项式函数级数来表示。进一步，而4次以上的多项式多半来自于信号中的高频部分，所以若假设信号中的高频部分很弱，则可以近似用0次项和2次项就可以近似 $r_{xx}[n]$ 对应的实数函数在0点的形状，即可以用抛物线来进行插值。

一些证明

下面证明几个阅读源码时使用傅立叶变换进行加速的证明

1. 设 $x[n] \in \mathbb R^N$ 是实数信号， $\mathcal F [x(-n)](k) = \overline {X[k]}$ ，（其中 $\mathcal F[\cdot]$ 表示傅立叶变换， $x[-n]$ 表示逆序信号）

$\sum_{n=0}^{N-1}x[N-n] e^{-j2\pi kn/N} = \sum_{n=0}^{N-1}x[N-n] e^{j2\pi k(N-n)/N} = \overline { \sum_{n=0}^{N-1}x[m] e^{-j2\pi km/N} }= \overline {X[k]}$

2. 对任意N维信号 $x[n] \in \mathbb C^N$ ，则 $\mathcal F [r_{xx}(\tau)](k) = X[k]\cdot X[-k]$ 。

$\begin{aligned}\mathcal F[r_{xx}(\tau)](k) &= \sum_{n=0}^{N-1} r_{xx}(\tau) e^{-j2\pi kn/N} \\&=\sum_{n=0}^{N-1} \left (\sum_{m=0}^{N-1} x(m)x(m+n)\right)e^{-j2\pi kn/N} \\&=\sum_{m=0}^{N-1}x(m) \left (\sum_{n=0}^{N-1} x(m+n)e^{-j2\pi kn/N} \right) \\&=\sum_{m=0}^{N-1}x(m) \left (\sum_{l=m}^{m+N-1} x(l)e^{-j2\pi k(l-m)/N} \right) \\&=\sum_{m=0}^{N-1}x(m) e^{-j2\pi (-k)m/N} \left (\sum_{l=0}^{N-1} x(l)e^{-j2\pi kl/N} \right) \\&=X[-k]\cdot X[k] = S_{xx}[k]\end{aligned}$

因而可知 $S_{xx}[k]$ 的傅立叶逆变换是 $r_{xx}(\tau)$ ，所以可以分别先用傅立叶正变换求出 $X[-k]$ 和 $X[k]$ ，计算复杂度为 $O(n \ {\rm log}n)$ ，再执行线性复杂度的复数乘法，最后再执行复杂度为 $O(n \ {\rm log}n)$ 的傅立叶逆变换，总体上计算复杂度为 $O(n \ {\rm log}n)$ ，当n足够大时，速度将远远快于直接计算相关函数 $r_{xx}(\tau)$ 的 $O(n^2)$ 的复杂度。

最后编辑于：2022.10.23 16:27:44

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 217,406评论 6赞 503
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 92,732评论 3赞 393
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 163,711评论 0赞 353
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 58,380评论 1赞 293
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 67,432评论 6赞 392
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 51,301评论 1赞 301
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 40,145评论 3赞 418
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 39,008评论 0赞 276
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 45,443评论 1赞 314
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 37,649评论 3赞 334
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 39,795评论 1赞 347
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 35,501评论 5赞 345
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 41,119评论 3赞 328
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 31,731评论 0赞 22
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,865评论 1赞 269
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 47,899评论 2赞 370
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 44,724评论 2赞 354