用于語音和音樂的基頻(音高)估計算法：YIN

????????周期信號的基頻估計（fundamental frequency estimation）在許多應用中都起著重要的作用辛孵，例如：在漢語語音識別中议谷，由于不同的聲調(diào)對應不同的字詞，所以音高是非常重要的特征弯囊，另外在識別說話人情感的時候氢拥，說話的語調(diào)也含有顯著的情緒傾向；還有在數(shù)字音樂應用中暂刘，包括樂譜轉(zhuǎn)錄饺谬、哼唱識別、K歌軟件等中都有用武之地鸳惯。

????????音高識別算法非常之多商蕴，包括在近些年發(fā)展起來的深度學習領域中也出現(xiàn)了音高估計模型的身影，下面介紹一下比較通用的YIN算法芝发，并羅列出用于閱讀代碼的主要幾大步驟绪商，并增加一些個人演算；更詳細的見大佬的原論文《YIN, a fundamental frequency estimator for speech and music》辅鲸。進一步格郁，工程中可以使用更加穩(wěn)定的probabilisitc YIN算法，我們將在以后的文章進行介紹独悴。

第一步：自相關函數(shù)和差函數(shù)

????????若信號有周期例书，則將信號在時間軸上平移一個周期將與自身重合。因而自相關函數(shù)（見下式）（autocorrelation function）將在周期的整倍數(shù)時刻取最大值刻炒，因而自相關函數(shù)顯然可以用來做基頻估計决采。

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

其中： $W$ 是執(zhí)行一次自相關運算的窗口大小坟奥； $x[n]$ 是時域信號树瞭。

? ? ? ? 自相關函數(shù)用來估計基頻的方法是，窮舉一定范圍內(nèi)所有的 $\tau$ 值爱谁，并且挑選最高的非零 $\tau$ 作為基本周期晒喷。換言之，自相關方法傾向于挑選較大的周期访敌，即傾向于較小的基頻凉敲。當 $\tau$ 的取值上限較大時，很容易發(fā)生低估基頻的錯誤寺旺。

????????還有一種自相關函數(shù)的定義：

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

? ? ? ? 該定義中爷抓，隨著 $\tau$ 的增加，被加項越來越少阻塑，將會形成振幅逐漸衰減的自相關函數(shù)形狀蓝撇。因而在新的定義中，盡管 $\tau$ 有較大的取值上界叮姑，理論上也不會容易出現(xiàn)低估基頻的錯誤唉地。但是當 $\tau$ 的取值下限足夠接近0時，由于 $\boldsymbol r_t^{\prime}[\tau]$ 振幅逐漸衰減的性質(zhì)传透，會傾向選擇接近0的基本周期耘沼，因而容易造成高估基頻的錯誤。

? ? ? ? 該論文首先就提出了若干論點朱盐，說明自相關函數(shù)在實際應用中不如新提出來的差函數(shù)（difference function）好用群嗤，其為：

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

? ? ? ? 換成差函數(shù)后，基本周期是使得該函數(shù)取最小值的那個 $\tau$ 值兵琳。并且差函數(shù)和自相關函數(shù)之間有如下關系：（其中狂秘，對中間的那個式子做個簡單的變量變換就可以了）

$\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$ 有關的項，所以這兩個準則的確會導致不同的估計清女。僅僅由于這一個改變就在測試數(shù)據(jù)上钱烟，將誤差從10%降低到了1.95%。作者給的解釋是自相關函數(shù)對信號的幅度變化非常敏感嫡丙；當信號振幅隨著時間增大時拴袭，自相關函數(shù)在基本周期的一定范圍內(nèi)的倍數(shù)時刻 $kT_0,(0<k<k_0)$ 的值會呈現(xiàn)增加的現(xiàn)象，因此由于自相關函數(shù)選擇最大周期對應的那個頻率作為識別的基礎頻率曙博，因此會高估基本周期拥刻，從而出現(xiàn)低估基礎頻率的錯誤。

? ? ? ? 而作者提出的差函數(shù)并沒有這類問題父泳，因為可以想象一下一個正弦函數(shù)的周期保持不變般哼，但是幅度變化，但是不管是逐漸變小還是逐漸變大尘吗，都是在平移一個基本周期時出現(xiàn)的差函數(shù)低洼比其他更大周期出現(xiàn)的差函數(shù)低洼更低逝她，甚至當振幅忽高忽低時，由于信號的規(guī)律性睬捶，差函數(shù)會在所有的 $\tau$ 上都變大黔宛，所以對變化的振幅會不那么敏感。

第二步：累計均值規(guī)范化差函數(shù)

? ? ? ? 直接使用差函數(shù)會出現(xiàn)兩個問題：1）和自相關函數(shù)類似擒贸，由于差函數(shù)仍然在 $\tau$ 為0的時候給出最小值0臀晃，并且由于實際音頻數(shù)據(jù)不是完美周期性，導致差函數(shù)在真實的基本周期處的取值仍然小于在0處的取值介劫，我們可以通過設置一個基本周期的取值下限還改善該問題徽惋。2）當在第一共振峰（first formant）處出現(xiàn)一個較強的共振，仍然會在2倍基頻處（即1/2倍基本周期處）出現(xiàn)一些次級的低洼座韵，這些次級低洼可能比基本周期對應的低洼還要更低险绘；而對于這第二個問題踢京，并不能通過增加一個取值下限來改善，因為基本周期 $F_0$ 和第一共振峰 $F_1$ 的取值范圍是重合的宦棺。所以作者進一步提出了“累計均值規(guī)范化差函數(shù)”（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}$

? ? ? ? 直接設置差函數(shù)在0處取1代咸，并且用直到 $\tau$ 時刻的差函數(shù)相應取值的均值（累積均值）來規(guī)范化蹈丸；可以想象，當信號原始的差函數(shù)在低于基本周期的地方出現(xiàn)更深的次級低洼時呐芥，累積均值同樣也會很小逻杖，因而除以很小的數(shù)相當于放大原始的差函數(shù)的值，因而更好地降低次級共振峰引起的錯誤率思瘟。并且直到 $\tau$ 達到真正的基本周期處出現(xiàn)很深的低洼荸百。見下圖，圖a）是原始的差函數(shù)形狀滨攻，圖b）是均值規(guī)范化后的管搪。新的差函數(shù)顯然可以抑制低估周期（高估基頻）的現(xiàn)象，并且新的差函數(shù)并不需要設置顯式的基頻上限铡买。

YIN論文插圖

第三步：絕對閾值

? ? ? ? 實際應用中很容易出現(xiàn)差函數(shù)在更大的 $\tau$ 上出現(xiàn)更深的低洼更鲁，例如上圖b），若圖中第二個很深的低洼處于候選區(qū)域內(nèi)奇钞，則很容易選擇第二個低洼對應的基頻澡为，而出現(xiàn)低估基頻的錯誤。為了應對這個問題景埃，我們可以考慮設置一個絕對的閾值（如圖中虛線）媒至，找到候選區(qū)域內(nèi)所有低于虛線的那些低洼所對應的頻率，將最小的 $\tau$ 對應的頻率（更高的頻率）作為基頻返回谷徙；而我們之所以可以放心地這么做拒啰，是因為在上一步我們用“均值規(guī)范化差函數(shù)”大大降低了高估基頻的錯誤。但是因為閾值是固定的完慧，一旦在閾值下方?jīng)]有找到仍和候選值谋旦，則算法直接返回全局的最小值。

? ? ? ? 通過閾值屈尼，作者將測試數(shù)據(jù)上的錯誤率從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$ 時我們做一個簡單的推理：“累積均值規(guī)范化差函數(shù)”的分子項為： $\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]$ 表示近似周期函數(shù)在一個周期內(nèi)積分，其結(jié)果可以近似視為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。我暫時沒有繼續(xù)推敲了不恭，不過這的確是作者所述的功率比率形式叶雹，這說明了絕對閾值的物理意義。

? ? ? ? 所以對這部分總結(jié)一下就是：“累計均值歸一化差函數(shù)”的閾值對應的是信號的非周期性程度换吧，我們的目的是估計周期信號的基頻折晦，所以當我們將其設置為很小的值時，我們就是期望信號是周期信號的置信度很高沾瓦，因而我們對我們估計出來的基頻抱有較大的信心满着。

第四步：拋物線插值

? ? ? ? 若信號的基本周期是剛好是采樣周期的整數(shù)倍，則我們使用前面的方法找到正確的差函數(shù)低洼贯莺，對應的就是正確的基頻风喇。但是一般不會滿足這個條件，所以使用樣本序列在時域進行基頻估計顯然是不準確的缕探，因此作者提出了拋物線插值：在函數(shù) $\boldsymbol d_t^{\prime}[\tau]$ 的每個極小值點魂莫，使用左右兩個相鄰點加上極小值點本身三個點進行拋物線插值，將插值得到的拋物線的最小值對應的橫坐標作為候選者放到最后的基頻選擇列表里爹耗。

????????事實上耙考，由于使用了“累積均值規(guī)范化差函數(shù)”谜喊，而規(guī)范化后的差函數(shù)和原始的差函數(shù)對極小值點橫坐標的估計有一定的偏差。所以倦始，實際會使用原始的差函數(shù)對應的三點進行插值斗遏。

? ? ? ? 設 $x[n]$ 是無周期的離散時間信號， $X(\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}$ 為其傅立葉變換鞋邑，且 $X(\omega)$ 為周期 $2\pi$ 的周期函數(shù)诵次。所以信號 $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]$ 為實數(shù)序列枚碗，并交換了積分和求和次序逾一。

????????所以我們得到了功率譜密度的傅立葉變換是信號 $x[n]$ 的自相關函數(shù) $r_{xx}[n]$ ，并且再次利用實信號性質(zhì)视译，得到 $|X(-\omega)|^2 = |X(\omega)^{\ast}|^2 = |X(\omega)|^2$ ，所以功率譜密度函數(shù)是偶函數(shù)归敬，因而其傅立葉變換中正弦變換部分消失了酷含，又假設 $X(\omega)$ 是帶限的，所以功率譜密度的傅立葉變換僅由一組余弦函數(shù)就能表達汪茧。

????????又每個余弦函數(shù)的在0點的泰勒展開式僅僅包括偶次項椅亚，即常數(shù)項、二次項舱污、4次項等等呀舔，所以總體上 $r_{xx}[n]$ 在0點可以由一組偶次項的多項式函數(shù)級數(shù)來表示。進一步扩灯，而4次以上的多項式多半來自于信號中的高頻部分媚赖，所以若假設信號中的高頻部分很弱，則可以近似用0次項和2次項就可以近似 $r_{xx}[n]$ 對應的實數(shù)函數(shù)在0點的形狀珠插，即可以用拋物線來進行插值惧磺。

一些證明

????????下面證明幾個閱讀源碼時使用傅立葉變換進行加速的證明

1. 設 $x[n] \in \mathbb R^N$ 是實數(shù)信號， $\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)$ 江解，再執(zhí)行線性復雜度的復數(shù)乘法设预，最后再執(zhí)行復雜度為 $O(n \ {\rm log}n)$ 的傅立葉逆變換，總體上計算復雜度為 $O(n \ {\rm log}n)$ 犁河，當n足夠大時絮缅，速度將遠遠快于直接計算相關函數(shù) $r_{xx}(\tau)$ 的 $O(n^2)$ 的復雜度鲁沥。

最后編輯于：2022.10.23 16:27:44

?著作權歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末，一起剝皮案震驚了整個濱河市耕魄，隨后出現(xiàn)的幾起案子画恰，更是在濱河造成了極大的恐慌，老刑警劉巖吸奴，帶你破解...
沈念sama閱讀 218,546評論 6贊 507
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件允扇，死亡現(xiàn)場離奇詭異，居然都是意外死亡则奥，警方通過查閱死者的電腦和手機考润，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 93,224評論 3贊 395
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來读处，“玉大人糊治，你說我怎么就攤上這事》２眨” “怎么了井辜？”我有些...
開封第一講書人閱讀 164,911評論 0贊 354
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵，是天一觀的道長管闷。經(jīng)常有香客問我粥脚，道長，這世上最難降的妖魔是什么包个？我笑而不...
開封第一講書人閱讀 58,737評論 1贊 294
?港島之戀（遺憾婚禮）
正文為了忘掉前任刷允，我火速辦了婚禮，結(jié)果婚禮上碧囊，老公的妹妹穿的比我還像新娘树灶。我一直安慰自己，他們只是感情好糯而，可當我...
茶點故事閱讀 67,753評論 6贊 392
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布破托。她就那樣靜靜地躺著，像睡著了一般歧蒋。火紅的嫁衣襯著肌膚如雪土砂。梳的紋絲不亂的頭發(fā)上，一...
開封第一講書人閱讀 51,598評論 1贊 305
城市分裂傳說
那天谜洽，我揣著相機與錄音萝映，去河邊找鬼。笑死阐虚，一個胖子當著我的面吹牛序臂，可吹牛的內(nèi)容都是我干的。我是一名探鬼主播，決...
沈念sama閱讀 40,338評論 3贊 418
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼奥秆，長吁一口氣：“原來是場噩夢啊……” “哼逊彭！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起构订，我...
開封第一講書人閱讀 39,249評論 0贊 276
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤侮叮，失蹤者是張志新（化名）和其女友劉穎，沒想到半個月后悼瘾，有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體囊榜，經(jīng)...
沈念sama閱讀 45,696評論 1贊 314
?護林員之死
正文獨居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 37,888評論 3贊 336
?白月光啟示錄
正文我和宋清朗相戀三年亥宿，在試婚紗的時候發(fā)現(xiàn)自己被綠了卸勺。大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
茶點故事閱讀 40,013評論 1贊 348
活死人
序言：一個原本活蹦亂跳的男人離奇死亡烫扼，死狀恐怖曙求，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情映企，我是刑警寧澤悟狱，帶...
沈念sama閱讀 35,731評論 5贊 346
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布，位于F島的核電站卑吭，受9級特大地震影響芽淡，放射性物質(zhì)發(fā)生泄漏马绝。R本人自食惡果不足惜豆赏，卻給世界環(huán)境...
茶點故事閱讀 41,348評論 3贊 330
男人毒藥：我在死后第九天來索命
文/蒙蒙一、第九天我趴在偏房一處隱蔽的房頂上張望富稻。院中可真熱鬧掷邦，春花似錦、人聲如沸椭赋。這莊子的主人今日做“春日...
開封第一講書人閱讀 31,929評論 0贊 22
一樁弒父案，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽哪怔。三九已至宣蔚，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間认境，已是汗流浹背胚委。一陣腳步聲響...
開封第一講書人閱讀 33,048評論 1贊 270
情欲美人皮
我被黑心中介騙來泰國打工，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留叉信，地道東北人亩冬。一個月前我還...
沈念sama閱讀 48,203評論 3贊 370
代替公主和親
正文我出身青樓，卻偏偏與公主長得像硼身，于是被迫代替她去往敵國和親硅急。傳聞我的和親對象是個殘疾皇子覆享，可洞房花燭夜當晚...
茶點故事閱讀 44,960評論 2贊 355