梅爾倒頻譜系數（Mel-frequency cepstrum Coefficents）

簡介

梅爾倒頻譜（MFC）
在聲音處理中，梅爾倒頻譜（MFC）表示了聲音短時功率譜。它基於非線性梅爾刻度頻率的對數功率譜的一個線性余弦變換。
梅爾頻率倒譜系數（MFCC）
梅爾頻率倒譜系數（MFCC）是所有構成MFC的系數。

倒譜和梅爾頻率倒譜的區別
在梅爾倒頻譜中，頻帶是等距地分布在Mel尺度上的，相比於在正常倒譜中線性間隔的頻帶，這種等距分布的頻帶其更接近於人類的聽覺系統。這種頻帶彎曲能更好地表現聲音。

MFCC 提取步驟

（1）把信號切分成幀
（2）計算每幀的功率譜周期圖估計（the periodogram estimate of the power spectrum）
（3）使用梅爾濾波器對每一幀的功率譜濾波， 再把每一幀濾波后的能量相加得到能量總和
（4）求所有濾波器組能量總和的對數
（5）對（4）得到的結果做離散余弦變換
（6）保留2-13的離散余弦對數，舍棄其他的

對每一個步驟進一步解釋

（1）把信號切分成幀
語音信號是時變的。在很短的一個時間間隔上，假定信號沒有發生變化或者變化很小。這個時間間隔設為20-30ms。過短不能保證有足夠的采樣點獲取可靠的能譜估計，過長是信號改變太多。

（2）計算每幀的功率譜周期圖估計
這是一個類似人耳蝸的工作的步驟。不同的頻率的聲音在耳蝸不同的地方振動，這些不同地方的振動通過神經告訴大腦聲音的不同振動。周期譜估計就是為了區分每一幀聲音中包含的頻率。

（3）使用梅爾濾波器對每一幀的功率譜濾波， 再把每一幀濾波后的能量相加
周期譜能量估計包含了很多語音識別不需要的信息。由於人的耳蝸很難辨別相近的兩個頻率，尤其是在高頻。通過集合周期圖譜並把他們相加得到各個頻段的能量。梅爾濾波器第一個得到 0Hz 附近的能量，越往后的濾波器頻帶越寬。梅爾濾波器是一組包含20-40個（標准26個）三角濾波器的濾波器。根據梅爾尺度（Mel scale）設置每一個濾波器，下面會詳細介紹如何設置。
頻率到梅爾尺度（Mel scale）轉換公式：

M(f)=1125ln(1+f/700) $$M(f)=1125ln(1+f/700)$$

梅爾尺度轉為頻率：

M−1=700(exp(m/1125)−1) $$M^{-1}=700(exp(m/1125)-1)$$

（4）求濾波后每一幀的能量的對數
這也是類似人的聽覺，我們不能在線性范圍聽到響度。使用對數而不是使用立方根是因為對數能夠讓我們使用倒數均值減法。

（5）對（4）得到的結果做離散余弦變換，保留2-13變換系數
因為濾波器組是有重疊的，濾波器的能量彼此相關；離散余弦變換能夠去相關，比如在HMM分類器，對角協方差矩陣就能構造。但是只有12-26 離散余弦變換系數被保留，由於高階系數表示快速變化的濾波器能量，其對信號的識別不好，所以將它們舍棄。

具體實現

(1)信號切分幀. 假設信號的采樣頻率為 16kHz，每一幀切分為25ms， 那么每一幀有400個采樣點。每一幀開始的步長通常設為160個采樣點，就是說第一幀從0采樣點開始，第二幀從第160個采樣點開始，一次類推，每一幀包含400個采樣點，使用 si(n) $s_i(n)$表示， i $i$表示每一幀， n $n$表示每一幀的采樣點。

Nframe = floor( (length(x) - wlen) / nstep) + 1; 
for k = 1:Nframe 
  idx = (1:wlen) + (k-1) * nstep; 
  x_frame(1,k) = x(idx); 
end 

有些情況下這里還有一個預加重的步驟，目的是為了高頻部分，凸顯高頻共振峰。實現方法如下公式：

si(n)∗′=si(n)−β∗si(n−1) $$s_i(n){\ast }^{\prime }=s_i(n)-\beta \ast s_i(n-1)$$

(2)離散傅里葉變換求能量譜

Si(n)=∑n=1Nsi(n)h(n)e−j2πkn/N,1≤k≤K $$S_i(n)=\underset{n=1}{\overset{N}{\sum }}{s}_i(n)h(n)e^{-j2\pi kn/N},1\le k\le K$$

K $K$ 離散傅里葉變換長度, h(n) $h(n)$表示N采樣長度的分析窗，比如 hamming窗。加窗的目的是為了平滑信號，在這里還能減弱傅里葉變換后旁瓣大小和頻譜泄漏。hamming窗函數如下， α $\alpha$通常取0.46，

w(n,α)=(1−α)−αcos(2∗π∗n/(N−1))，0≤n≤N−1 $$w(n,\alpha )=(1-\alpha )-\alpha cos(2\ast \pi \ast n/(N-1))，0\le n\le N-1$$

直接相乘給信號加窗，具體實現如下：

si(n)∗′=w(n,α)×si(n)=[(1−α)−αcos(2∗π∗n/(N−1))]×sn，0≤n≤N−1 $$s_i(n){\ast }^{\prime }=w(n,\alpha )\times s_i(n)=[(1-\alpha )-\alpha cos(2\ast \pi \ast n/(N-1))]\times s_n，0\le n\le N-1$$

每一幀的周期能量譜估計，

Pi(k)=1N|Si(k)|2 $$P_i(k)=\frac{1}{N}|S_i(k){|}^2$$

(3)使用梅爾濾波器對每一幀的功率譜濾波
當選擇濾波器的個數為26時，濾波器表現為以26個長度為 K/2 $K/2$ 向量的三角濾波器。傅里葉的窗為512個采樣點，那么這一組濾波器如下圖（a）。單獨一個濾波器如第8個和第20個濾波器如圖（c）（e），將每一個梅爾濾波器乘以能量譜就得到了每一個濾波器濾波后的能量（how much energy was in each filterbank）。
下面會詳細介紹如何計算梅爾濾波器能量。

(4)對3中得到的每個能量取對數
(5) 對4得到的26個對數做離散余弦變換得到26個倒譜系數，2-13的保留

設計梅爾濾波器組

（1）先選定一個上限頻率和下限頻率。將上、下限頻率在梅爾尺度下表示，如300Hz對應401.25Mels，8000Hz對應2834.99Mels。

M(f)=1125ln(1+f/700) $$M(f)=1125ln(1+f/700)$$

（2）假如我們使用10個濾波器，那么需要12個點。上面我們給出了上下限兩個點，中間還需要10個等間距的點，如下：

m(i) = 401.25, 622.50, 843.75, 1065.00, 1286.25, 1507.50, 1728.74, 1949.99, 2171.24, 2392.49, 2613.74, 2834.99

%code
midpoint(1,:) = lower_frequency;
for c=1:(num_filter+2)
    midpoint(1,c+1) = midpoint(1,c) + (upper_frequency-lower_frequency)/(num_filter+2);
end

（3）將梅爾尺度表示的值轉換會頻域

M−1=700(exp(m/1125)−1) $$M^{-1}=700(exp(m/1125)-1)$$

h(i) = 300, 517.33, 781.90, 1103.97, 1496.04, 1973.32, 2554.33, 3261.62, 4122.63, 5170.76, 6446.70, 8000

（4）沒有頻率精度要求去把濾波器放在精確到點，所以我們需要把這些頻率取在 the nearest FFT bin。 To convert the frequncies to fft bin numbers，我們需要知道FFT的大小和采樣率：

f(i) = floor((nfft+1)*h(i)/samplerate)

把h(i)代入得到f(i)

f(i) = 9, 16, 25, 35, 47, 63, 81, 104, 132, 165, 206, 256

最后的濾波器在256， 與8kHz， 512個點的FFT大小對應。

（5）建立濾波器
濾波器是三角濾波器，第一個濾波器從第一點開始，第二個時取得最大值，第三個點又歸零。第二個濾波器從第二個點開始，第三點時達到最大值，第四點歸零，依次類推。由下面的公式表達：

m 表示濾波器的數量， f() $f()$表示m+2梅爾間隔頻率（Mel-spaced frequencies）列表 。

Deltas and Delta-Deltas

這兩個系數也稱為微分系數和加速系數。MFCC特征向量只描述了一幀的功率譜包絡。另外隨着時間的推移，MFCC系數的軌跡如何？ delta系數可以通過下面的公式計算，多少個MFCC系數將得到多少個delta系數。將delta系數補充MFCC向量上得到一個長度是原來MFCC向量兩倍的一個新的特征向量。

dt $d_t$表示delta系數，N的典型取值為2。

delta-delta系數（加速度系數）由delta系數計算得到，把上面的計算公式的 ct $c_t$換成 dt $d_t$就可以。

極力推薦閱讀參考2.
在python上的應用實例：https://github.com/jameslyons/python_speech_features
在Matlab上的應用實例：https://labrosa.ee.columbia.edu/matlab/rastamat/

1. https://en.wikipedia.org/wiki/Mel-frequency_cepstrum
2. http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/
3. http://www.speech.cs.cmu.edu/15-492/slides/03_mfcc.pdf
4. http://haythamfayek.com/2016/04/21/speech-processing-for-machine-learning.html
5. http://blog.csdn.net/xmdxcsj/article/details/51228791
6. https://tspace.library.utoronto.ca/bitstream/1807/44123/1/Mohamed_Abdel-rahman_201406_PhD_thesis.pdf
7. http://mirlab.org/jang/books/audioSignalProcessing/speechFeatureMfcc_chinese.asp?title=12-2%20MFCC
8. http://ibillxia.github.io/blog/2012/07/18/MFCC-feature-extraction/#disqus_thread