一、MFCC概述[1]
在语音识别（SpeechRecognition）和话者识别（SpeakerRecognition）方面，最常用到的语音特征就是梅尔倒谱系数（Mel-scaleFrequency Cepstral Coefficients，简称MFCC）。根据人耳听觉机理的研究发现，人耳对不同频率的声波有不同的听觉敏感度。从200Hz到5000Hz的语音信号对语音的清晰度影响对大。两个响度不等的声音作用于人耳时，则响度较高的频率成分的存在会影响到对响度较低的频率成分的感受，使其变得不易察觉，这种现象称为掩蔽效应。由于频率较低的声音在内耳蜗基底膜上行波传递的距离（速度）大于频率较高的声音，故一般来说，低音容易掩蔽高音，而高音掩蔽低音较困难。在低频处的声音掩蔽的临界带宽较高频要小。所以，人们从低频到高频这一段频带内按临界带宽的大小由密到疏安排一组带通滤波器，对输入信号进行滤波。将每个带通滤波器输出的信号能量作为信号的基本特征，对此特征经过进一步处理后就可以作为语音的输入特征。由于这种特征不依赖于信号的性质，对输入信号不做任何的假设和限制，又利用了听觉模型的研究成果。因此，这种参数比基于声道模型的LPCC相比具有更好的鲁邦性，更符合人耳的听觉特性，而且当信噪比降低时仍然具有较好的识别性能。
MFCC是在Mel标度频率域提取出来的倒谱参数，Mel标度描述了人耳频率的非线性特性，它与频率的关系可用下式近似表示：
其中
为频率，单位Hz。下图展示了Mel频率与线性频率之间的关系：
附：
import
二、MFCC提取过程[2]
基本流程：
连续语音 预加重 分帧 加窗 FFT Mel滤波器组 对数运算 DCT
0. 连续语音
我们这里截取了3.5s的原始声音数据
import numpy as npimport scipy.io.wavfilefrom matplotlib import pyplot as pltfrom scipy.fftpack import dct# 原始数据,读取前3.5s 的数据sample_rate, signal = scipy.io.wavfile.read('OSR_us_000_0010_8k.wav')original_signal = signal[0:int(3.5*sample_rate)]signal_num = np.arange(len(signal))sample_num = np.arange(len(original_signal))# 绘图 01plt.figure(figsize=(11,7), dpi=500)plt.subplot(211)plt.plot(signal_num/sample_rate, signal, color='black')plt.plot(sample_num/sample_rate, original_signal, color='blue')plt.ylabel("Amplitude")plt.title("signal of Voice")plt.subplot(212)plt.plot(sample_num/sample_rate, original_signal, color='blue')plt.xlabel("Time (sec)")plt.ylabel("Amplitude") plt.title("3.5s signal of Voice ")plt.savefig('mfcc_01.png')
预加重
预加重处理其实是将语音信号通过一个高通滤波器：其中 ，我们通常取为0.97。预加重的目的是提升高频部分，使信号的频谱变得平坦，保持在低频到高频的整个频带中，能用同样的信噪比求频谱。同时，也是为了消除发生过程中声带和嘴唇的效应，来补偿语音信号受到发音系统所抑制的高频部分，也为了突出高频的共振峰。
# 预加重
2. 分帧
先将N个采样点集合成一个观测单位，称为帧。通常情况下N的值为256或512，涵盖的时间约为20~30ms左右。为了避免相邻两帧的变化过大，因此会让两相邻帧之间有一段重叠区域，此重叠区域包含了M个取样点，通常M的值约为N的1/2或1/3。通常语音识别所采用语音信号的采样频率为8KHz或16KHz，以8KHz来说，若帧长度为256个采样点，则对应的时间长度是256/8000 1000=32ms。
预加重后，我们需要将信号分成短时帧。这一步背后的基本原理是信号中的频率随时间而变化，所以在大多数情况下，对整个信号进行傅立叶变换是没有意义的，因为我们会随着时间的推移丢失信号的频率轮廓。为了避免这种情况，我们可以安全地假设信号中的频率在很短的时间内是平稳的。因此，通过在这个短时间帧内进行傅里叶变换，我们可以通过连接相邻帧来获得信号的频率轮廓的良好近似。
语音处理范围内的典型帧大小范围为20毫秒到40毫秒，连续帧之间重叠50％（+/- 10％）。流行设置25毫秒的帧大小，frame_size = 0.025和-10毫秒的步幅（15毫秒重叠）， frame_stride = 0.01。
# 分帧
3. 加窗
将每一帧乘以汉明窗，以增加帧左端和右端的连续性。假设分帧后的信号为S(n), n=0,1,…,N-1, N为帧的大小，那么乘上汉明窗后 形式如下：不同的 值会产生不同的汉明窗，一般情况下取为0.46。
# 汉明窗
在将信号切分为帧后，我们将一个窗口函数（如Hamming窗口）应用于每个帧。有几个原因需要对窗框应用窗函数，特别是为了抵消FFT所假设的数据是无限的，并减少频谱泄漏。
# 加汉明窗
4. 快速傅里叶变换 FFT
由于信号在时域上的变换通常很难看出信号的特性，所以通常将它转换为频域上的能量分布来观察，不同的能量分布，就能代表不同语音的特性。所以在乘上汉明窗后，每帧还必须再经过快速傅里叶变换以得到在频谱上的能量分布。对分帧加窗后的各帧信号进行快速傅里叶变换得到各帧的频谱。并对语音信号的频谱取模平方得到语音信号的功率谱。设语音信号的DFT为： ，式中 为输入的语音信号， 表示傅里叶变换的点数。
# 傅里叶变换和功率谱NFFT = 512mag_frames = np.absolute(np.fft.rfft(frames, NFFT)) # Magnitude of the FFTpow_frames = (1.0 / NFFT) * (mag_frames ** 2)
5. 三角带通滤波器 （Mel滤波）
将能量谱通过一组Mel尺度的三角形滤波器组，定义一个有M个滤波器的滤波器组（滤波器的个数和临界带的个数相近），采用的滤波器为三角滤波器，中心频率为 。M通常取22-26。各f(m)之间的间隔随着m值的减小而缩小，随着m值的增大而增宽，如图所示
三角滤波器的频率响应定义为:式中，
三角带通滤波器有两个主要目的：对频谱进行平滑化，并消除谐波的作用，突显原先语音的共振峰。（因此一段语音的音调或音高，是不会呈现在MFCC 参数内，换句话说，以MFCC 为特征的语音辨识系统，并不会受到输入语音的音调不同而有所影响）此外，还可以降低运算量。
6. 计算每个滤波器组输出的对数能量为
，
# 将频率转换为Mel频率low_freq_mel = 0nfilt = 40high_freq_mel = (2595 * np.log10(1 + (sample_rate / 2) / 700))mel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2) # Equally spaced in Mel scalehz_points = (700 * (10**(mel_points / 2595) - 1)) # Convert Mel to Hzbin = np.floor((NFFT + 1) * hz_points / sample_rate)fbank = np.zeros((nfilt, int(np.floor(NFFT / 2 + 1))))for m in range(1, nfilt + 1):f_m_minus = int(bin[m - 1]) # leftf_m = int(bin[m]) # centerf_m_plus = int(bin[m + 1]) # rightfor k in range(f_m_minus, f_m):fbank[m - 1, k] = (k - bin[m - 1]) / (bin[m] - bin[m - 1])for k in range(f_m, f_m_plus):fbank[m - 1, k] = (bin[m + 1] - k) / (bin[m + 1] - bin[m])filter_banks = np.dot(pow_frames, fbank.T)filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks) # Numerical Stabilityfilter_banks = 20 * np.log10(filter_banks) # dB
7. 经离散余弦变换（DCT）得到MFCC系数
，将上述的对数能量带入离散余弦变换，求出L阶的Mel-scale Cepstrum参数。L阶指MFCC系数阶数，通常取12-16。这里M是三角滤波器个数。
事实证明，前一步计算出的滤波器组系数高度相关，这在某些机器学习算法中可能存在问题。因此，我们可以应用离散余弦变换（DCT）去相关滤波器组系数并产生滤波器组的压缩表示。通常，对于自动语音识别（ASR），所得到的倒谱系数2-13被保留，其余的被丢弃; num_ceps = 12。丢弃其他系数的原因是它们表示滤波器组系数的快速变化，并且这些细节不会有助于自动语音识别（ASR）。
num_ceps
可以将正弦升降1应用于MFCC以降低已被声称在噪声信号中改善语音识别的较高MFCC.
n = np.arange(ncoeff)cep_lifter =22lift = 1 + (cep_lifter / 2) * np.sin(np.pi * n / cep_lifter)mfcc *= lift
如前所述，为了平衡频谱并改善信噪比（SNR），我们可以简单地从所有帧中减去每个系数的平均值。平均归一化滤波器组：
filter_banks
对于MFCC也是如此，平均归一化MFCC：
mfcc
最后结果为
# 绘图 04plt.figure(figsize=(11,7), dpi=500)plt.subplot(211)plt.imshow(np.flipud(filter_banks.T), cmap=plt.cm.jet, aspect=0.2, extent=[0,filter_banks.shape[1],0,filter_banks.shape[0]]) #画热力图plt.title("MFCC")plt.subplot(212)plt.imshow(np.flipud(mfcc.T), cmap=plt.cm.jet, aspect=0.2, extent=[0,mfcc.shape[0],0,mfcc.shape[1]])#热力图plt.title("MFCC")plt.savefig('mfcc_04.png')
8. 对数能量
此外，一帧的音量（即能量），也是语音的重要特征，而且非常容易计算。因此，通常再加上一帧的对数能量（定义：一帧内信号的平方和，再取以10为底的对数值，再乘以10）使得每一帧基本的语音特征就多了一维，包括一个对数能量和剩下的倒频谱参数。注：若要加入其它语音特征以测试识别率，也可以在此阶段加入，这些常用的其它语音特征包含音高、过零率以及共振峰等。
9.动态查分参数的提取（包括一阶差分和二阶差分）
标准的倒谱参数MFCC只反映了语音参数的静态特性，语音的动态特性可以用这些静态特征的差分谱来描述。实验证明：把动、静态特征结合起来才能有效提高系统的识别性能。差分参数的计算可以采用下面的公式：式中,dt表示第t个一阶差分；Ct表示第t个倒谱系数；Q表示倒谱系数的阶数；K表示一阶导数的时间差，可取1或2。将上式中结果再代入就可以得到二阶差分的参数。
总结：因此，MFCC的全部组成其实是由：N维MFCC参数（N/3MFCC系数+ N/3一阶差分参数+ N/3二阶差分参数）+帧能量（此项可根据需求替换）。
