背景引入
特征值提取，在模式识别领域是很常见的一种算法和手段。特征值看起来好像很陌生，其实在我们日常生活中也很常见。我们使用的身份认证，ID，都可以视为不同系统下的特征值。
MFCC在语音识别领域就是一组特征向量，它通过对语音信号（频谱包络与细节）进行编码运算来得到。MFCC有39个系数，其中包括13个静态系数，13个一阶差分系数，以及13个二阶差分系数。差分系数用来描述动态特征，也就是声学特征在相邻帧间的变化情况。这些系数都是通过离散余弦变换（Discrete Cosine Transform，DCT）计算而来。
MFCC语音特征值提取算法简介
MFCC意为梅尔频率倒谱系数，顾名思义，MFCC语音特征提取包含两个关键步骤；将语音信号转化为梅尔频率，然后进行倒谱分析。梅尔频谱是一个可用来代表短期音频的频谱，梅尔刻度（Mel Scale）则是一种基于人耳对等距的音高变化的感官判断而确定的非线性频率刻度。梅尔频率和正常的频率f之间的关系：
mel(f)=25951g(1+f/700)mel(f)=25951g(1+f/700)
当梅尔刻度均匀分布，则对应的频率之间的距离会越来越大。梅尔刻度的滤波器组在低频部分的分辨率高，跟人耳的听觉特性比较相符，这也是梅尔刻度的物理意义。在梅尔频域内，人对音调的感知度为线性关系，如果两段语音的梅尔频率相差两倍，则人耳听起来两者的音调也相差两倍。
转化为梅尔频率时，首先对时域信号进行离散傅里叶变换，将信号转换到频域，然后再利用梅尔刻度的滤波器组对频域信号进行切分，使每个频率段对应一个数值。倒谱（Cepstrum）通过对一个时域信号进行傅里叶变换后取对数，并再次进行反傅里叶变换（Inverse Fast Fourier Transform，IFFT）得到。倒谱可分为复倒谱（Complex Cepstrum），实倒谱（Real Cepstrum）和功率倒谱（Power Cepstrum）。倒谱分析可用于信号分解，也就是将两个信号的卷积转化为两个信号的相加。
MFCC的物理含义，简而言之，可理解为语音信号的能量在不同频率范围的分布。
人的发声过程可以看作是肺里的气流通过声带这个线性系统。如果用e（t）表示输入声音的音高，h(t)表示声带的响应（也即我们需要获取的语音特征），那么听到的语音信号x(t)即为二者的卷积：
x(t)=e(t)∗h(t)x(t)=e(t)∗h(t)
x（t)为时域信号，对其进行离散傅里叶变换后可得到频域信号X（K），亦即频谱：
X(K)=DFT(x(t))X(K)=DFT(x(t))
时域信号的卷积在频域内则可表示为二者的乘积：
X(K)=E(k)∗H(k)X(K)=E(k)∗H(k)
通常，在频域分析中我们只关注频谱的能量而忽略其相位信息，即：
||X(k)||=||E(k)||||H(k)||||X(k)||=||E(k)||||H(k)||
对频谱进行对数运算：
log||X(k)||=log||E(k)||+log||H(k)||log||X(k)||=log||E(k)||+log||H(k)||
然后进行反傅里叶计算：
c(n)=IFFT(log||X(k)||)=IFFT(log||E(k)||)+IFFT(log||H(k)||))c(n)=IFFT(log||X(k)||)=IFFT(log||E(k)||)+IFFT(log||H(k)||))
c(n)即为倒谱系数，已经和原始的时域信号x(t)不一样。并且时域信号的卷积关系已经转化为频域信号的线性相加关系。
语音信号分帧
语音信号属于准稳态信号，这也意味着，在一定的时间内，信号会保持稳定。这个时常对于我们人类来说，很短，一般只有10ms~30ms。在这一区间（即帧）内，可将语音信号看成稳态信号，只有稳态信号才能进行信号处理。
信号分帧一般会涉及到一个加窗的操作，即将原始信号与一个窗函数相乘。我们用计算机处理信号的时候，一般不会取无限长的信号，而是会取其中间的一段信号，这将会减少工作量，也会加快程序法分析的时间。
无限长的信号被截断后，其频谱会发生畸变，从而导致频谱能量泄露。 为了减少这种能量泄露，我们可采用不同的截取函数对信号进行截断。执行截断操作的函数称为窗函数，简称为窗。常用的窗函数有矩形窗，三角窗，汉明（Hamming）窗及汉宁窗等。
汉宁窗也叫升余弦窗，是很有用的窗函数。如果测试测试信号有多个频率分量，频谱表现非常复杂，测试目的更多在于关注频率点而非能量大小，则用汉宁窗。汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减少泄漏的观点出发，汉宁窗明显优于矩形窗。但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨率下降，他与矩形窗相比，泄露以及波动较小，选择性则相应较高。
汉明窗是用来加权余弦形成的锥形窗，也称之为改进的升余弦窗，只是加权系数不同，其旁瓣更小，但其旁瓣衰减速度比汉宁窗要慢。汉明窗是以著名的美国数学家理查德·卫斯理·汉明（Richard Wesley Hamming）的名字来命名：
w(n)=0.54−0.46cos(2πnM−1)(0≤n⩽M−1)w(n)=0.54−0.46cos(2πnM−1)(0≤n⩽M−1)
下面的代码就是用python来生成汉明窗和汉宁窗：
import matplotlib.pyplot as plt import scipy #信号处理工具包 plt.figure(figsize=(6,2)) plt.plot(scipy.hanning(512),"b--",label="Hanning") #绘制汉宁窗 plt.plot(scipy.hamming(512),"r--",label="Hamming") #绘制汉明窗 plt.title("Demo Hanning & Hamming Window") plt.legend() plt.show()
除了scipy模块可以实现汉宁窗和汉明窗，我们也可以用NumPy来实现汉宁窗和汉明窗。示例代码如下：
import numpy as np import matplotlib.pyplot as plt hanWing=np.hanning(512)#定义汉宁窗 hamWin=np.hamming(512) #定义汉明窗 plt.plot(hanWing,'y--',label="Hanning") plt.plot(hamWin,'b--',label="Hamming") plt.title("Hamming & Hanning window") plt.ylabel("Amplitude") plt.xlabel("Sample") plt.legend() plt.show()
信号加窗，从本质上而言，就是将原始信号与一个窗函数相乘。进行加窗操作之后，我们就可以对信号进行傅里叶展开。加窗的代价就是，一帧信号的两端部分将会被消弱。所以在进行信号分帧处理时，帧与帧之间需要有部分重叠。相邻两帧重叠后，其起始位置的时间差称之为帧移，即步长（Stride） 。
以下为简单的信号加窗操作示意图：
import numpy as np import matplotlib.pyplot as plt import scipy x=np.linspace(0,10,1000) originWav=np.sin(x**2) #示例原信号 win=scipy.hamming(1000) #定义一个窗函数，这里使用的汉明窗 winFrame=originWav*win #结果可视化 plt.title("Signal Chunk with Hamming Windows") plt.plot(originWav) plt.plot(win) plt.plot(winFrame) plt.legend() plt.show()
运行程序，其中蓝色波形为原信号，橙色波形为窗函数，绿色为加窗操作之后的信号。
假设x为语音信号，w为窗函数，则分帧信号为：
y(n)=∑N/2n=−(N/2)+1x(m)w(n−m)y(n)=∑n=−(N/2)+1N/2x(m)w(n−m)
其中，w(n-m)为窗口序列，当n去不同的值时，窗口w(n-m)沿x(m)。因此，w(n-m)是一个“滑动的”窗口。y(n)为短时傅里叶变换（SIFT）。由于窗口是有限长度的，满足绝对可和条件，所以这个变幻的前提条件是存在的，这也是信号分帧的理论依据。
以下示例代码从指定文件夹读取一个音频文件，然后将该音频文件分帧并显示其中一个分帧信号的波形：
#读取指定音频文件 import matplotlib.pyplot as plt import numpy as np import wave #导入波形处理工具包 import os import soundfile def audioSignalFrame(signal,nw,inc): ''' signal:原始音频信号 nw:每一帧的长度 inc:相邻帧的间隔 ''' #信号总长度 signal_length=len(signal) #若信号长度小于一个帧的长度，则帧数定义为1 if signal_length<=nw: nf=1 else: nf=int(np.ceil((1.0*signal_length-nw+inc)/inc)) #所有帧加起来的总的铺平的长度 pad_length=int((nf-1)*inc+nw) #长度不够时，使用0填补，类似于FFT中的扩充数组长度 zeros=np.zeros((pad_length-signal_length,)) #填补后的信号 pad_signal=np.concatenate((signal,zeros)) #相当于对所有帧的时间点进行抽取，得到nf*nw的长度的矩阵 indices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T #将indices转化为矩阵 indices=np.array(indices,dtype=np.int32) #得到帧信号 frames=pad_signal[indices] #窗函数，这里默认取1 return frames def readSignalWave(filename): f=wave.open(filename,'rb') params=f.getparams() nchannels,sampwidth,frammerate,nframes=params[:4] #读取音频，字符串格式 strData=f.readframes(nframes) #将字符串转化为int waveData=np.fromstring(strData,dtype=np.int16) f.close() #信号幅值归一化 waveData=waveData*1.0/(max(abs(waveData))) waveData=np.reshape(waveData,[nframes,nchannels]).T return waveData if __name__=='__main__': filepath="./test.wav" #dirname=os.listdir(filepath) #filename=filepath+dirname[3] data=readSignalWave(filepath) #初始化每帧长度及帧间隔 nw=512 inc=128 Frame=audioSignalFrame(data[0],nw,inc) #显示原始信号 plt.plot(data[0]) plt.title("Original Signal") plt.show() #显示第一帧信号 plt.plot(Frame[0]) plt.title("First Frame") plt.show()
上面的代码中，没有对信号进行加窗处理，若要执行信号加床操作，只需将分帧函数稍作修改，
def audioSignalFrame(signal,nw,inc,winfunc): ''' signal:原始音频信号 nw:每一帧的长度 inc:相邻帧的间隔 ''' #信号总长度 signal_length=len(signal) #若信号长度小于一个帧的长度，则帧数定义为1 if signal_length<=nw: nf=1 else: nf=int(np.ceil((1.0*signal_length-nw+inc)/inc)) #所有帧加起来的总的铺平的长度 pad_length=int((nf-1)*inc+nw) #长度不够时，使用0填补，类似于FFT中的扩充数组长度 zeros=np.zeros((pad_length-signal_length,)) #填补后的信号 pad_signal=np.concatenate((signal,zeros)) #相当于对所有帧的时间点进行抽取，得到nf*nw的长度的矩阵 indices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T #将indices转化为矩阵 indices=np.array(indices,dtype=np.int32) #得到帧信号 frames=pad_signal[indices] #窗函数，这里默认取1 win=np.tile(winfunc,(nf,1)) return frames*win
当然，随着函数的改变，主函数中对函数的调用也需要改变，只需要改变参数即可。除了调用工具包中的汉明窗函数，也可以使用公式来定义。
def hamming(n): return 0.54-0.46*cos(2*pi/n*(arange(n)+0.5))
语音信号在进行分帧之前，一般需要进行一个与加重操作。语音信号的预加重，是为了对语音的高频部分进行加重，使信号变得平坦，保持在地频到高频的整个频带中能用同样的信噪比求频谱。同时也为了消除发声过程中声带和口唇辐射效应，补偿语音信号受到发音系统所抑制的高频部分，增加语音的高频分辨率。
我们一般通过一阶有限长单位冲激响应（Finite Impulse Response，FIR）高通数字滤波器来实现预加重。FIR滤波器H(z)=1−az−1H(z)=1−az−1作为传递函数，其中a为预加重系数，0.9<a<1.0。假设t时刻的语音采样值为x(t)，经过预加重处理后的结果为y(t)=x(t)−ax(t−1)y(t)=x(t)−ax(t−1).a一般默认取0.95.
信号的预加重处理示例代码：
def preemphasis(signal,coeff=0.95): ''' signal:要滤波的输入信号 coeff:预加重系数。0表示无过滤，默认为0.95 返回值：滤波信号 ''' return numpy.append(signal[0],signal[1:]-coeff*signal[:-1])
