AI(I)語音(I):MFCC特徵參數提取

網上很多關於MFCC提取的文章，但本文純粹我自己手碼，本來不想寫的，但這東西忘記的快，所以記錄我自己看一個python demo並且自己本地debug的過程，在此把這個demo的步驟記下來，所以文章主要傾向說怎麼做，而不是道理論述。由於python的matplotlib.pyplot庫沒有下載成功不會畫圖，文中大部分圖片是我網上找的。

必備基礎知知識:

1. 對數指數自然對數正弦餘弦求和求積向量相乘導數一階二階微分

2. 矩陣及其乘法，A[i][j],i行j列，A[m][n]*B[x][y]必須滿足n=x，此處為了提高用戶體驗度，我一般用二維數組來描述矩陣，不到萬不得已我不會搬出矩陣來嚇人。

3. 傅利葉變換，容易忘記，忘記就參考[傅里葉變換 掐死教程]

4. 格式為.wav的音頻文件，注意注意不是.avi，.wav文件簡單理解就是無損音頻文件，該種格式文件可以通過c#快速編寫一個錄音程序得到。個人理解這種文件其實存儲的是錄音設備按固定頻率採取的真實聲波的某個點，並記錄下該刻度值，從而做到從真實信號到模擬信號的轉換。比如聲波在實際中是一個連續的波形，但弱水三千只取一瓢飲，錄音設備每隔很短的時間就取一個點，把這些點記錄下來，成為離散的值，如果間隔時間足夠短的話可以基本反映聲音的特徵。

5.採樣定律，第4步說的每隔很短的時間就採取一個點，那究竟是多短？答案是大概1秒鐘採取8000或16000個點，這8000和16000也叫做採樣頻率，採樣頻率需要滿足真實信號最大頻率處的2倍，這樣可以基本真實還原聲音。比如音效卡採樣率是8000HZ，那麼可以認為原始信號的最大頻率處是4000HZ，注意，原始信號的源頭是聲帶一張一合把肺部氣體排出，這一張一合的頻率叫做基音頻率（聲音波由3個階段，聲帶發出的基音頻率是第一階段；第二階段是氣體聲波經過長約17cm的聲道，據說第二階段作為聲音音色音質以及文本無關聲紋識別中的重要部位，也是共振峰產生的場所；第三階段是唇口鼻舌，這些部位影響聲音的發音，比如母音輔音等），為了理解清晰，原始信號的頻率可以默認為是聲帶一張一合的頻率，而音效卡採樣頻率則是1秒鐘採取多少個點，所以這兩者有本質的區別，但都是頻率。

---------------------------------------------------------------------------------------------------------------------------

一，處理.wav格式的音頻文件，得到信號數據和採樣頻率

首先是拿到.wav格式的音頻文件，其它格式比如.mp3需要進行轉換，轉換軟體或方法有很多，不在本文討論範圍。大家很好奇該文件究竟存的什麼東西，我一言以蔽之，其實除了文件頭以外，就是聲音波形圖按照某一刻度刻畫出來的離散點的值，粗暴點可以描述為信號signal = [ 0 0 -1 ..., 627 611 702]，這是我自己demo的數據（接下來都以這一數據作研究）。.wav文件進行處理後就是含有107000個元素的數組，怎麼從 .wav文件得到這個數組，參考python（scipy和numpy庫）的(rate,signal) = scipy.io.wavfile.read("Ansel.wav")，python庫的這一個方法就可以讀取到.wav文件的信號數組signal和該聲音文件的採樣頻率，此處頻率rate=8000HZ。

二，預加重，{A(1*107000)}

預加重我至今沒完全理解，看了一些資料，也專門加了qq群去問人，也沒問到，可能我資質太低，我一直不理解時域上的採樣點為啥有高低頻之分，時域採樣點都是固定的採樣頻率（8000HZ或16000HZ），所以我對下面這段話：

[語音和圖像信號低頻段能量大，高頻段能量小...低頻段信噪比大，高頻信噪比低...增大高頻段信噪比....]

我的理解：

這裡說的高頻低頻是指時域上採樣點分幀後的頻率，針對每一幀，作傅利葉變換得到頻域的N個分量，這N個分量按照頻率為橫坐標，振幅為縱坐標。如果沒有噪音，這N個分量振幅值應該隨著頻率增大逐漸遞減，但因為噪音存在，在某一高頻段處，振幅值出現反常，變得很大，導致信噪比很小，而預加重就是為了把高頻段的信號都放大，從而增大高頻段的信噪比。

我的理解極有可能是錯的，但我知道預加重的做法，還是以我做的demo做數據，signal = [ 0 0 -1 ..., 627 611 702]含有107000個點，預加重做法就是針對這些點，套用公式signal[i]=signal[i+1]-0.97*signal[i]，因此得到一個新的數組A=[ 0. 0. -1. ..., -40.36 2.81 109.33]，數組大小還是107000。係數0.97可以自己選取，據說0.95左右。

三，分幀，{B(1336*200)}

分幀原理就不說了，我說下做法：把上述數組A的元素值，200個為一幀，80個為幀移，總共可以得到1336幀，因為不是整除，1336*80+200=107080，多出的80個點用0填補，所以得到1336*200的二維數組B。

上述200個為一幀，80個為幀移是這樣來的，因為默認語音信號具有短時平穩性，這平穩性不是說短時採樣得到的值相等，而是認為短時間內聲帶、聲道、唇口鼻腔這3個聲音信號源頭具有平穩性，據說是人體肌肉活動短時平穩。短時默認是0.025s，兩個相鄰短時間隔默認是0.01s，因為我分析的.wav文件是8000HZ， 所以0.025s內有0.025*8000=200個採樣點，幀移0.01*8000=80個採樣點。

四，加窗，{C(1336*200)}

需要一個1336*200的窗數組C，這裡默認元素值都初始化為1（為了簡單，用的矩形窗），實際中需用漢明窗，其實漢明窗一樣道理，只不過用漢明窗的的話，這1336*200的二維數組元素值需要根據漢明窗函數計算得到。漢明窗函數是一個餘弦函數，作用是為了使幀和幀之間變得平滑，消除吉布斯效應（傅利葉變換時無法得到邊界值，據說是傅利葉打敗拉格朗日的跨歷史爭執），接下來說下具體怎麼加窗。

分幀後的數組B和窗函數數組C具有相同的維度，把它們對應位置的元素值相乘，即可以得到加窗後的二維數組C，C[i[j]=B[i][j]*C[i][j]]，因為我用的矩形窗，C[i][j]值都為1，所以C=B。

五，離散傅利葉變換（FFT），{D(1336*257)，E(1336*257)，F(1336*1)}

首先還是溫習下[傅里葉變換 掐死教程]，看完這個教程我認為，已知經過加窗後的信號C具有1336*200維，也就是說有1336行（幀），每行200個點，每相鄰兩行之間有120個重合點（幀移80）。那麼對它做傅利葉變換，粗暴的理解為，對有1336幀，每一幀都作N=512的傅利葉變換，這個過程可以認為每一幀分解為257個分量，每個分量其實是一個正/餘弦波的振幅，也就把短時的時域幀（200個點）分解為頻域上的257個分量，有點像力學分解...

因此可以得到1336*257維的頻域信號D，幀數還是1336，對每一幀的257個點的值（個人認為是振幅）取平方，再乘以1/512，便得到能量，這個公式我不理解，我只記得高中時有印象振幅好像代表能量，但是平方再除以512，我就不明白為啥，難道跟面積/邊長計算類似？不管怎樣，得到1336*257的能量普E，然後對每一幀的257個能量值簡單相加，得到該幀的能量總值，一共有1336幀，於是有1336個能量總值，即擁有1336個元素值，記為數組F，數組中每個元素值代表一幀的能量總值。

六，獲得梅爾濾波器{G(26*257)}

如上圖，梅爾值是一個新的量度，據說相比頻率量度，梅爾更接近人耳的聽覺機理，通俗的說，就像納米和米一樣，如果我們用納米衡量我們身邊的事物會是一種什麼感受？所以頻率的某一個值對應著梅爾的某一個值，該印射關係可以用這個公式描述，梅爾值f(f)=2595*lg(1+f/700.0)，如果要反過來，頻率f(m)=700*(10**(m/2595.0)-1)，10**(m/2595.0)意思是10的(m/2595.0)次方。我們看第一個公式，我們的採樣頻率除以2就是真實信號的最大頻率，真實信號的最小頻率為0，依據公式的單調性，我們以這個最大頻率和最小頻率為界限分別得到梅爾刻度的最大最小值，可以把信號的所有頻率值刻畫在這個梅爾區間之內。

梅爾濾波器個數一般默認26個，前期準備工作需要在上述最大最小梅爾區間等間距插入26個值，包括邊界，就是28個值，然後把這28個值的頻率值也算出來，得到28個頻率和梅爾的一一對應關係。接下來對這28個頻率值，依次代入公式y=(512+1)*x/8000)，便可以得到28個y值。我不明白這個公式幹嘛的，只知道，512是傅利葉變換的N，8000是採樣頻率，。

接下來計算26個濾波器的二維數組，先初始化為元素值都為0的26*257二維數組G，然後通過循環填補該二維數組的值，過程大致是針對G的每一行，根據行下標（0~25），再結合上述的28個y值，計算出每行的各個元素值，其計算過程是把相鄰兩個y值相減作分母，分子是每行的元素值下標減去28個元素的下標為行元素值下標的值，再用28個元素的下標為行元素值下標的值減去每行的元素值下標。針對每一行，這樣可以得到三角形形狀的數據分布，三角形頂點為1，除此三角形數據分布之外的點都為0，這個二維數組G經過轉化之後便得到每行只有一部分有值，其餘值為0的二維矩陣，而且有值部分數據呈三角形分布。此處很難解釋，沒圖說個j*，用文字簡直不可描述，還是上圖來說明。如下所示，圖中H1(k)是G的第一行元素值分布，H3(k)是G的第二行元素值分布，H26(k)是G的第26行元素分布，每一行都有257個元素值，比如H1(k)，也就是G的第一行，只有開頭幾個元素值有值，其餘200多個值均為0，且有值部分值大小先線形增大到1，再線性減小到0。

七，得到能量特徵參數的和能量總值{H(1336*26)}

把第五步得到的二維矩陣能量譜E(1336*257)，乘以第六部的二維數組梅爾濾波器G(26*257)的逆，矩陣的逆可得到257*26的矩陣，然後滿足矩陣乘法定律，得到參數H=E*G.T，此處的H其實是1336*26的二維矩陣。還有個參數是第五步計算出來的每幀能量總值F(1336*1)，即擁有1336個元素值的一維數組F。

八，作自然對數運算，離散餘弦變換（DCT）和升倒譜運算{J(1336*13)，K(13*1)，L(1336*13)，feat[1336*13]}

對H的每一個元素值做ln運算，即H[i][j]=ln(H[i][j])，此處我也不明白原理，但很多資料記載說需要做這一部操作，好像是什麼公式。接著對feat的每一行做離散餘弦變換（離散餘弦變換類似傅利葉變換，只不過作用在實數範圍）。據說離散餘弦變換後的數據分布可以把冗餘數據分開，而且大部分信號數據一般集中在變換後的低頻區，所以對每一幀只取前13個數據就好了，於是得到1336*13的二維數組J。

針對1336*13的二維數組J做升倒譜操作，默認升倒譜係數為22，這個過程做法是先產生一個擁有13個元素的一維數組K，這13個元素的值K[i]=1+(22 /2)*sin(pi*i/22)，其中22是升倒譜係數，pi是圓周率3.1415926。得到這個數組K之後，針對1336*13的二維數組J，J[i][j] = J[i][j]*K[j]，得到1336*13的二維數組L，這其實就是mfcc參數的第一組。如果這組參數想要加上能量作為其表示方式，可以把這1336幀，也就是每一行的的第一個元素用一維數組F的每個值替換，即L[i][0] = F[i]。

我們把這經過錯綜複雜得到的L記為feat，它是個二維數組，擁有1336*13個值，這也是mfcc參數的基礎參數，也是第一組，默認是有3組，接下來計算第二和第三組參數。

九，計算第二組和第三組參數{feat[1336*13]，feat"[1336*13]，feat""[1336*13]}

一言以蔽之，第二組參數其實就是在已有的基礎參數下作一階微分操作，第三組參數在第二組參數下作一階微分操作，相當於對基礎參數導數的導數。微分（dx）如果忘了，就把他理解為自變數增長為1的函數值的變化量，準確點描述是要計算離散點之間的變換率，而不是連續函數的導數，但原理類似。

具體操作是這樣的，抽取一個計算一階微分的函數，然後把1336*13的二維數組feat作為參數傳入，返回feat的一階微分feat『，feat"同樣有1336*13個元素值。這個函數有點複雜，我把python函數貼上來。我簡要說兩句，其實這個過程是把feat按照行作循環，每一行擁有13個元素值，如果不考慮邊際效應，feat"[i][j]={feat[i][j+1]-feat[i][j-1] + 2(feat[i][j+2]-feat[i][j-2]) + 3(feat[i][j+3]-feat[i][j-3]) + ... + n(feat[i][j+n]-feat[i][j-n])} / M，M作為分母，是這個函數輸入的big_theta來決定的，big_theta默認值是2，此處的M值可以算出是10，過程參照函數代碼。這個公式可能還是太複雜，如果大家對一階展開式或者泰勒展開式還有印象，對這些應該就不會陌生，上述公式可以再簡化點，大概長成這樣f』(x)={f(x+1)-f(x-1) + 2(f(x+2)-f(x-2)) + 3(f(x+3)-f(x-3)) + ...+ n(f(x+n)-f(x-n)) } / M，這個就平易近人了，所以此處的計算一階微分函數，大致過程就是這樣。最後得到的輸出是1336*13的二維數組feat"，也是mfcc參數的第二組參數。

第三組參數同上，把feat"作為參數運用抽取出來的函數作計算，得到輸出記為feat""，feat""二維數組同樣有1336*13個元素值，這是mfcc參數的第三組參數。

十，得到最後輸出{mfcc(1336*39)}

由前八步和第九步，可以得到feat，feat"和feat""，這3個參數都是擁有1336*13個元素值的二維數組，而且這三個二維數組的每一行第一個元素值可以根據需要，用該行（幀）的能量總值替換。把feat,feat『和feat""拼在一起，即基於feat，每一行橫向追加feat"和feat""每行的元素值，得到擁有1336*39個元素值的一個二維數組，也就是mfcc係數，這就是最後得到的結果。

得到的mfcc係數，可供語音識別或聲紋識別（文本無關，文本相關）等技術，但語音識別需要語料庫，因此還需要建立語音模型來訓練語音，經典的語音模型有HMM，新興的有神經元模型，但不管怎樣，語音模型的建立難度應該比mfcc更大...