计算MFCC系数
由于信号在时域上的变换很难看出特征，因此，我们通常将它转换为频域上的能量分布以便于观察。不同的能量分布，代表不同语音的特征。语音原信号在与窗函数（如汉明窗）相乘后，每帧还必须再经过快速傅里叶变换以得到频谱上的能量分布。对语音信号分帧加窗后的各帧的频谱，然后对频谱进行取模平方运算后即为语音信号的功率谱。
对信号幅度谱、功率谱以及对数功谱的计算实例代码如下：
import numpy import logging def msgspec(frames,NFFT): """ 计算帧中每个帧的幅度谱。如果帧为N*D，则输出N*（NFFT/2+1） """ if numpy.shape(frames)[1]>NFFT: logging.warn('frame length (%d)is greater than FFT size(%d),frame will be truncated .Increase NFFT to avoid.',numpy.shape(frames)[1],NFFT) complex_spec=numpy.fft.rfft(frames,NFFT) return numpy.absolute(complex_spec) def power_spectrum(frames,NFFT): return 1.0/NFFT*numpy.square(spectrum_magnitude(frames,NFFT)) def log_power_spectrum(frames,NFFT,norm=1): spec_power=power_spectrum(frames,NFFT) spec_power[spec_power<1e-30] log_spec_power=10*numpy.log10(spec_power) if norm: return log_spec_power-numpy.max(log_spec_power) else: return log_spec_power
此外，信号的每一帧的音量（即能量），也是语音的特征，而且非常容易计算。因此，通常会再加上一帧的能量，使得每一帧基本的语音特征增加一个维度，包括一个对数能量和倒谱参数。标准的倒谱参数MFCC，只反映了语音参数的静态特征，语音参数的动态特征可以用这些静态特征的差分普来描述。
MFCC的全部组成如下：N维MFCC系数（N/3 MFCC系数+N/3 一阶差分系数+N/3二阶差分系数）+帧能量。以语音识别中常用的39维MFCC为例，即为：13个静态系数+13个一阶差分系数（Delta系数）+13个二阶差分系数（Delta-Delta系数）。其中，差分系数用来描述动态特征，即声学特征在相邻帧间的变化情况。
在MFCC计算中还涉及频率与梅尔刻度之间的转换，其转换方式如下：
m=2595lg(1+f700)m=2595lg(1+f700) def hz2mel(hz): return 2595*numpy.log10(1+hz/700.0)
同样，我们也可以推出下列公式：
f=700(10m/2595−1)f=700(10m/2595−1)
Delta系数的计算公式为：
dt=∑Nn=1n(ct+n−ct−n)2∑Nn=1n2dt=∑n=1Nn(ct+n−ct−n)2∑n=1Nn2
其中，dtdt为Delta系数，从帧t根据静态系数ct−N到ct+Nct−N到ct+N 计算而得。N一般取值为2。Delta-Delta（加速度）系数的计算方法相同，但他们是根据Delta而不是静态系数来进行计算得到的。计算Deltaa系数的示例代码如下：
def delta(feat,N): if N<1: raise ValueError('N must be an integer>=1') NUMFRAMES=len(feat) denominator=2*sum([i**2 for i in range(1,N+1)]) delta_feat=numpy.pad(feat,((N,N),(0,0)),mode='edge') for t in range(NUMFRAMES): delta_feat[t]=numpy.dot(numpy.arange(-N,N+1),padded[t:t+2*N+1])/denominator return delta_feat
当然除了自己定义函数，也可以直接使用工具包中的API。
对语音信号进行预加重
import numpy as np import matplotlib.pyplot as plt from python_speech_features.sigproc import * from python_speech_features import * from scipy.fftpack import dct import scipy.io.wavfile as wav sample_rate,signal=wav.read('./test.wav') #保留语音的前3.5秒 signal=signal[0:int(3.5*sample_rate)] #信号预加重 emphasized_signal=preemphasis(signal,coeff=0.95) #显示信号 plt.plot(signal) plt.title("Original Signal") plt.plot(emphasized_signal) plt.title("Preemphasis Signal") plt.show()
上述示例代码，对信号进行预加重处理的是preemphasis（signal，coeff）函数，除了这个函数，也可以使用以下代码实现：
pre_emphasis=0.95
emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1])
源代码：
import numpy as np import matplotlib.pyplot as plt from python_speech_features.sigproc import * from python_speech_features import * from scipy.fftpack import dct import scipy.io.wavfile as wav sample_rate,signal=wav.read('./test.wav') pre_emphasis=0.95 emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1]) #保留语音的前3.5秒 #signal=signal[0:int(3.5*sample_rate)] #信号预加重 #emphasized_signal=preemphasis(signal,coeff=0.95) #显示信号 plt.plot(signal) plt.title("Original Signal") plt.plot(emphasized_signal) plt.title("Preemphasis Signal") plt.show()
通过上面的程序可知，两种函数都可以进行预加重处理，可以自行选择合适的方法。
对语音信号进行短时傅里叶变换
在对语音信号进行处理之前，我们需要对不稳定的语音信号进行短时分帧以获取傅里叶变换必需的稳定信号。语音处理范围内的典型帧大小范围为20ms~40ms，连续帧之间重叠50%左右。因此一般将帧长度设置为25ms。短时傅里叶变换（Short-Time Fourier Transform，SIFT）在MFCC计算过程中主要用于短时分帧处理后，通过对信号进行时域到频域的转换来获取语音信号的频谱。
#对信号进行短时分帧处理 frame_size=0.025 #设置帧长 #计算帧对应采样数（frame_length)以及步长对应采样数（frame_step) frame_length,frame_step=frame_size*sample_rate,frame_stride*sample_rate signal_length=len(emphasized_signal) #信号总采样数 frame_length=int(round(frame_length)) #帧采样数 frame_step=int(round(frame_step)) #num_frames为总帧数，确保我们至少有一个帧 num_frames=int(np.ceil(float(np.abs(signal_length-frame_length))/frame_step)) pad_signal_length=num_frames*frame_step+frame_length z=np.zeros((pad_signal_length-signal_length)) #填充信号以后确保所有的帧的采样数相等 pad_signal=np.append(emphasized_signal,z) indices=np.tile(np.arange(0,frame_length),(num_frames,1))+np.tile(np.arange(0,num_frames*frame_step,frame_step),(frame_length,1)).T frames=pad_signal[indices.astype(np.int32,copy=False)]
信号经过短时分帧之后，可通过短时傅里叶变换得到各种频谱
NFFT=512 mag_frames=np.absolute(np.fft.rfft(frames,NFFT)) pow_frames=((1.0/NFFT)*((mag_frames)**2)) log_pow_frames=logpowspec(pow_frames,NFFT,norm=1) #保留语音的前3.5秒 #signal=signal[0:int(3.5*sample_rate)] #信号预加重 #emphasized_signal=preemphasis(signal,coeff=0.95) #显示信号 plt.plot(mag_frames) plt.title("Mag_Spectrum") plt.plot(emphasized_signal) plt.show() plt.plot(pow_frames) plt.title("Power_Spectrum") plt.show() plt.plot(pow_frames) plt.title("Log_Power_Spectrum") plt.show()
运行上面的程序，就可以得到处理结果，下面展示原有的所有代码：
import numpy as np import matplotlib.pyplot as plt from python_speech_features.sigproc import * from python_speech_features import * from scipy.fftpack import dct import scipy.io.wavfile as wav sample_rate,signal=wav.read('./test.wav') pre_emphasis=0.95 emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1]) #对信号进行短时分帧处理 frame_size=0.025 #设置帧长 frame_stride=0.1 #计算帧对应采样数（frame_length)以及步长对应采样数（frame_step) frame_length,frame_step=frame_size*sample_rate,frame_stride*sample_rate signal_length=len(emphasized_signal) #信号总采样数 frame_length=int(round(frame_length)) #帧采样数 frame_step=int(round(frame_step)) #num_frames为总帧数，确保我们至少有一个帧 num_frames=int(np.ceil(float(np.abs(signal_length-frame_length))/frame_step)) pad_signal_length=num_frames*frame_step+frame_length z=np.zeros((pad_signal_length-signal_length)) #填充信号以后确保所有的帧的采样数相等 pad_signal=np.append(emphasized_signal,z) indices=np.tile(np.arange(0,frame_length),(num_frames,1))+np.tile(np.arange(0,num_frames*frame_step,frame_step),(frame_length,1)).T frames=pad_signal[indices.astype(np.int32,copy=False)] NFFT=512 mag_frames=np.absolute(np.fft.rfft(frames,NFFT)) pow_frames=((1.0/NFFT)*((mag_frames)**2)) log_pow_frames=logpowspec(pow_frames,NFFT,norm=1) #保留语音的前3.5秒 #signal=signal[0:int(3.5*sample_rate)] #信号预加重 #emphasized_signal=preemphasis(signal,coeff=0.95) #显示信号 plt.plot(mag_frames) plt.title("Mag_Spectrum") plt.plot(emphasized_signal) plt.show() plt.plot(pow_frames) plt.title("Power_Spectrum") plt.show() plt.plot(log_pow_frames) plt.title("Log_Power_Spectrum") plt.show()
(a)幅度谱
(b)功率谱
（c)功率对数谱
音频文件使用不同，最终结果也会不同，大家自己使用自己的音频，注意音频格式为“.wav”
定义滤波器组
将信号通过一组梅尔刻度的三角形滤波器组，采用的滤波器为三角形滤波器，中心频率为f(m),m=1,2,3,```````,M,M通常取22~26. 各f(m)之间的间隔随着m值的减少而减少。随着m值的增大而增大。如图：
三角形滤波器的频率响应定义公式：4
Hm(k)=⎧⎩⎨⎪⎪2(k−f(m−1))(f(m+1)−f(m−1)(f(m)−f(m−1)))2(f(m+1)−k)(f(m+1)−f(m−1)(f(m+1)−f(m−1)))f(m−1)≤k≤f(m+1)f(m)≤k≤f(m+1)Hm(k)={2(k−f(m−1))(f(m+1)−f(m−1)(f(m)−f(m−1)))f(m−1)≤k≤f(m+1)2(f(m+1)−k)(f(m+1)−f(m−1)(f(m+1)−f(m−1)))f(m)≤k≤f(m+1)
对于其他的情况，例如，k<f(m-1)和k>=f(m+1)则为0，当k=f(m)时为1.
定义梅尔刻度的三角形滤波器组的示例代码为：
low_freq_MEL=0 #将频率转换为梅尔刻度 nfilt=40 #窗的数目 #计算m=2595*log10(1+f/700) high_freq_mel=(2595*np.log10(1+(sample_rate/2)/700)) mel_points=np.linspace(low_freq_MEL,high_freq_mel,nfilt+2) #梅尔刻度的均匀分布 #计算f=700(10**(m/2595)-1) hz_points=(700*(10**(mel_points/2595)-1)) bin=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]) #三角形滤波器左边频率f(m-1) f_m=int(bin[m]) #三角形滤波器中间频率fm f_m_plus=int(bin[m+1]) #三角形滤波器右边频率f(m-1) for k in range(f_m_minus,f_m): fbank[m-1,k]=(k-bin[m-1])/(bin[m+1]-bin[m]) plt.plot(fbank.T) plt.show()
三角形滤波器有两个主要功能，其一，对频谱进行平滑并消除谐波的作用，突显原先语音的共振峰；其二，用以降低运算量。如图所示的滤波器组中的每个滤波器在中心频率处响应为1，并朝着0线性减少，直至达到响应为0的两个相邻滤波器的中心频率。
计算MFCC系数
如果计算出的滤波器组系数高度相关，则在某些机器学习算法中可能会存在问题。我们可用离散余弦变换对滤波器组系数进行去相关，并产生滤波器组的压缩表示。滤波器组输出的对数能量经离散余弦变换后，即可得到MFCC系数。示例代码如下：
import numpy as np import matplotlib.pyplot as plt from python_speech_features.sigproc import * from python_speech_features import * from scipy.fftpack import dct import scipy.io.wavfile as wav sample_rate,signal=wav.read('./test.wav') pre_emphasis=0.95 emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1]) #对信号进行短时分帧处理 frame_size=0.025 #设置帧长 frame_stride=0.1 #计算帧对应采样数（frame_length)以及步长对应采样数（frame_step) frame_length,frame_step=frame_size*sample_rate,frame_stride*sample_rate signal_length=len(emphasized_signal) #信号总采样数 frame_length=int(round(frame_length)) #帧采样数 frame_step=int(round(frame_step)) #num_frames为总帧数，确保我们至少有一个帧 num_frames=int(np.ceil(float(np.abs(signal_length-frame_length))/frame_step)) pad_signal_length=num_frames*frame_step+frame_length z=np.zeros((pad_signal_length-signal_length)) #填充信号以后确保所有的帧的采样数相等 pad_signal=np.append(emphasized_signal,z) indices=np.tile(np.arange(0,frame_length),(num_frames,1))+np.tile(np.arange(0,num_frames*frame_step,frame_step),(frame_length,1)).T frames=pad_signal[indices.astype(np.int32,copy=False)] NFFT=512 mag_frames=np.absolute(np.fft.rfft(frames,NFFT)) pow_frames=((1.0/NFFT)*((mag_frames)**2)) log_pow_frames=logpowspec(pow_frames,NFFT,norm=1) #保留语音的前3.5秒 #signal=signal[0:int(3.5*sample_rate)] #信号预加重 #emphasized_signal=preemphasis(signal,coeff=0.95) #显示信号 ''' plt.plot(mag_frames) plt.title("Mag_Spectrum") plt.plot(emphasized_signal) plt.show() plt.plot(pow_frames) plt.title("Power_Spectrum") plt.show() plt.plot(log_pow_frames) plt.title("Log_Power_Spectrum") plt.show() ''' low_freq_MEL=0 #将频率转换为梅尔刻度 nfilt=40 #窗的数目 #计算m=2595*log10(1+f/700) high_freq_mel=(2595*np.log10(1+(sample_rate/2)/700)) mel_points=np.linspace(low_freq_MEL,high_freq_mel,nfilt+2) #梅尔刻度的均匀分布 #计算f=700(10**(m/2595)-1) hz_points=(700*(10**(mel_points/2595)-1)) bin=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]) #三角形滤波器左边频率f(m-1) f_m=int(bin[m]) #三角形滤波器中间频率fm f_m_plus=int(bin[m+1]) #三角形滤波器右边频率f(m-1) for k in range(f_m_minus,f_m): fbank[m-1,k]=(k-bin[m-1])/(bin[m+1]-bin[m]) plt.plot(fbank.T) plt.show() filter_banks=np.dot(pow_frames,fbank.T) filter_banks=np.where(filter_banks==0,np.finfo(float).eps,filter_banks) filter_banks=20*np.log10(filter_banks) num_ceps=12 #取12个系数 #通过DCT计算MFCC系数 mfcc=dct(filter_banks,type=2,axis=1,norm='ortho')[:,1:(num_ceps+1)] #对MFCC进行倒谱提升可以改善噪声信号中的语音识别 (nframes,ncoeff)=mfcc.shape n=np.arange(ncoeff) cep_lifter=22 #倒谱滤波系数，定义倒谱所用到的滤波器组内滤波器个数 lift=1+(cep_lifter/2)*np.sin(np.pi*n/cep_lifter) mfcc*=lift mfcc-=(np.mean(mfcc,axis=0)+1e-8) plt.imshow(np.flipud(mfcc.T),cmap=plt.cm.jet,aspect=0.2,extent=[0,mfcc.shape[0],0,mfcc.shape[1]]) #绘制MFCC热力图 plt.show()
对MFCC进行如下的归一化操作，运行操作，其相应的热力图如下：
filter_banks-=(np.mean(filter_banks,axis=0)+1e-8) 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.show()
归一化的MFCC热力图
