梅尔倒频谱系数（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)=\sum _{n=1}^N{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