MFCC理解

MFCC

在语音识别（SpeechRecognition）和话者识别（SpeakerRecognition）方面，最常用到的语音特征就是梅尔倒谱系数（Mel-scaleFrequency Cepstral Coefficients，简称MFCC）。根据人耳听觉机理的研究发现，人耳对不同频率的声波有不同的听觉敏感度。从200Hz到5000Hz的语音信号对语音的清晰度影响对大。两个响度不等的声音作用于人耳时，则响度较高的频率成分的存在会影响到对响度较低的频率成分的感受，使其变得不易察觉，这种现象称为掩蔽效应。由于频率较低的声音在内耳蜗基底膜上行波传递的距离大于频率较高的声音，故一般来说，低音容易掩蔽高音，而高音掩蔽低音较困难。在低频处的声音掩蔽的临界带宽较高频要小。所以，人们从低频到高频这一段频带内按临界带宽的大小由密到疏安排一组带通滤波器，对输入信号进行滤波。将每个带通滤波器输出的信号能量作为信号的基本特征，对此特征经过进一步处理后就可以作为语音的输入特征。由于这种特征不依赖于信号的性质，对输入信号不做任何的假设和限制，又利用了听觉模型的研究成果。因此，这种参数比基于声道模型的LPCC相比具有更好的鲁邦性，更符合人耳的听觉特性，而且当信噪比降低时仍然具有较好的识别性能。
梅尔倒谱系数（Mel-scale Frequency Cepstral Coefficients，简称MFCC）是在Mel标度频率域提取出来的倒谱参数，Mel标度描述了人耳频率的非线性特性，它与频率的关系可用下式近似表示：
在这里插入图片描述
式中f为频率，单位为Hz。下图展示了Mel频率与线性频率的关系：

MFCC，它的13个系数（也许还有13个一阶差分和13个二阶差分）都是通过离散余弦变换（DCT）而来，并取前13个系数。

DCT的作用，为了获得频谱的倒谱，倒谱的低频分量就是频谱的包络，倒谱的高频分量就是频谱的细节，这些都是语音识别上经过科学验证十分有效的语音物理信息。但语音识别是对机器来说的，所以我们必须对物理信息进行某种“编码”，得到的就是MFCC13维特征向量，你是否理解它不重要，只要机器能够识别就行了。
MFCC的物理含义就是——在语音识别领域，将语音物理信息（频谱包络和细节）进行编码运算得到的一组特征向量。简单地说可以理解为语言信号的能量在不同频率范围的分布。
在这里插入图片描述
如果把计算出的系数的低位部分（一般是前12个）作反傅立叶变换（IFFT），就可以得到上图中虚线表示的信号的频谱包络，也就是表示声带特征的那一部分低频信息。
要理解为什么可以这么做，我们先看看倒谱的定义：一种信号的傅里叶变换谱经对数运算后再进行的傅里叶反变换（IFFT）。（IFFT其实就是对频谱再作一次FFT）
人的发声过程可以看作是肺里的气流通过声带这个线性系统。如果用e(t)表示声音的输入激励（音高），h(t)表示声带的响应（也即我们需要获取的特征），那么听到的语音信号为
x(t) = e(t) * h(t)
即两者的卷积。在频域内则可以表示为两者的乘积
X(w) = E(w) * H(w)
通常在频域分析中我们只关注频谱的能量，忽略其相位信息，即
在这里插入图片描述
接下来对频谱作对数运算

最后再作傅里叶反变换便可以得到倒谱系数：

下面两幅图很好地解释了这个过程，语音信号的频谱可以看作是低频的包络和高频的细节相加，在得到倒谱系数后，我们只需要取低位的系数便可以得到包络信息。（这里的x[k]即为倒谱系数）
在这里插入图片描述
注意整个过程中我们并没有明确计算 e(t) 和 h(t) ，而是通过直接对语音信号x(t)作倒谱分析，再提取低位的倒谱系数，就可以获取我们想要的声道特征。

有意思的是，对数频谱作傅里叶反变换后的域被称为quefrency domain（对应频域frequency domain），quefrency domain和时域time domain类似但不完全一样。提取低位系数的操作称为low-time liftering（对应滤波操作filtering）。同样地，我们可以通过high-time liftering来获取激励特征。
在这里插入图片描述

当然在MFCC的实际计算中我们并没有作傅里叶反变换，而是先将频谱通过一组三角滤波器，再作离散余弦变换（DCT）得到MFCC系数。但其物理意义是一样的，即表示信号频谱的能量在不同频率区间的分布。每一个滤波器的作用就是得到对应频率区间的频谱能量，如果我们有26个三角滤波器，就会得到26个MFCC系数，这时候再取低位的系数就可以代表声道的特征。

MFCC特征参数提取过程详解

(1)预处理

预处理包括预加重、分帧、加窗函数。

预加重：预加重的目的是提升高频部分，使信号的频谱变得平坦，保持在低频到高频的整个频带中，能用同样的信噪比求频谱。同时，也是为了消除发生过程中声带和嘴唇的效应，来补偿语音信号受到发音系统所抑制的高频部分，也为了突出高频的共振峰。预加重处理其实是将语音信号通过一个高通滤波器：
在这里插入图片描述
分帧：先将N个采样点集合成一个观测单位，称为帧。通常情况下N的值为256或512，涵盖的时间约为20~30ms左右。为了避免相邻两帧的变化过大，因此会让两相邻帧之间有一段重叠区域，此重叠区域包含了M个取样点，通常M的值约为N的1/2或1/3。通常语音识别所采用语音信号的采样频率为8KHz或16KHz，以8KHz来说，若帧长度为256个采样点，则对应的时间长度是256/8000×1000=32ms。

加窗：将每一帧乘以汉明窗，以增加帧左端和右端的连续性。假设分帧后的信号为S(n), n=0,1…,N-1, N为帧的大小，那么乘上汉明窗
在这里插入图片描述
后，W(n)形式如下：

(2)FFT

由于信号在时域上的变换通常很难看出信号的特性，所以通常将它转换为频域上的能量分布来观察，不同的能量分布，就能代表不同语音的特性。所以在乘上汉明窗后，每帧还必须再经过快速傅里叶变换以得到在频谱上的能量分布。对分帧加窗后的各帧信号进行快速傅里叶变换得到各帧的频谱。

(3)谱线能量

对语音信号的频谱取模平方得到语音信号的谱线能量。

(4)计算通过Mel滤波器的能量

将能量谱通过一组Mel尺度的三角形滤波器组，定义一个有M个滤波器的滤波器组（滤波器的个数和临界带的个数相近），采用的滤波器为三角滤波器，中心频率为f(m) 。M通常取22-26。各f(m)之间的间隔随着m值的减小而缩小，随着m值的增大而增宽，如图所示：
在这里插入图片描述
三角滤波器的频率响应定义为:

对频谱进行平滑化，并消除谐波的作用，突显原先语音的共振峰。（因此一段语音的音调或音高，是不会呈现在MFCC 参数内，换句话说，以MFCC 为特征的语音辨识系统，并不会受到输入语音的音调不同而有所影响）此外，还可以降低运算量。
计算每个滤波器组输出的对数能量为：
在这里插入图片描述
(5)计算DCT倒谱

经离散余弦变换（DCT）得到MFCC系数 :
在这里插入图片描述
将上述的对数能量带入离散余弦变换，求出L阶的Mel参数。L阶指MFCC系数阶数，通常取12-16。这里M是三角滤波器个数。

实现绘图（python）

import numpy as np 
from scipy import signal
from scipy.fftpack import dct
import pylab as pltdef enframe(wave_data, nw, inc, winfunc):'''将音频信号转化为帧。参数含义：wave_data:原始音频型号nw:每一帧的长度(这里指采样点的长度，即采样频率乘以时间间隔)inc:相邻帧的间隔（同上定义）'''wlen=len(wave_data) #信号总长度if wlen<=nw: #若信号长度小于一个帧的长度，则帧数定义为1nf=1else: #否则，计算帧的总长度nf=int(np.ceil((1.0*wlen-nw+inc)/inc))pad_length=int((nf-1)*inc+nw) #所有帧加起来总的铺平后的长度zeros=np.zeros((pad_length-wlen,)) #不够的长度使用0填补，类似于FFT中的扩充数组操作pad_signal=np.concatenate((wave_data,zeros)) #填补后的信号记为pad_signalindices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T  #相当于对所有帧的时间点进行抽取，得到nf*nw长度的矩阵indices=np.array(indices,dtype=np.int32) #将indices转化为矩阵frames=pad_signal[indices] #得到帧信号win=np.tile(winfunc,(nf,1))  #window窗函数，这里默认取1return frames*win   #返回帧信号矩阵Df=5
fs=8000
N=fs/Df
t = np.arange(0,(N-1)/fs,1/fs)      
wave_data=np.sin(2*np.pi*200*t)
#预加重
#b,a = signal.butter(1,1-0.97,'high')
#emphasized_signal = signal.filtfilt(b,a,wave_data)
#归一化倒谱提升窗口
lifts=[]
for n in range(1,13):lift =1 + 6 * np.sin(np.pi * n / 12)lifts.append(lift)
#print(lifts)   #分帧、加窗 
winfunc = signal.hamming(256) 
X=enframe(wave_data, 256, 80, winfunc)    #转置的原因是分帧函数enframe的输出矩阵是帧数*帧长
frameNum =X.shape[0] #返回矩阵行数18，获取帧数
#print(frameNum)
for i in range(frameNum):y=X[i,:]#fftyf = np.abs(np.fft.fft(y)) #print(yf.shape)#谱线能量yf = yf**2#梅尔滤波器系数nfilt = 24low_freq_mel = 0NFFT=256high_freq_mel = (2595 * np.log10(1 + (fs / 2) / 700))  # 把 Hz 变成 Melmel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2)  # 将梅尔刻度等间隔hz_points = (700 * (10**(mel_points / 2595) - 1))  # 把 Mel 变成 Hzbin = np.floor((NFFT + 1) * hz_points / fs)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(yf[0:129], fbank.T)filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks)  # 数值稳定性filter_banks = 10 * np.log10(filter_banks)  # dB filter_banks -= (np.mean(filter_banks, axis=0) + 1e-8)#print(filter_banks)#DCT系数num_ceps = 12c2 = dct(filter_banks, type=2, axis=-1, norm='ortho')[ 1 : (num_ceps + 1)] # Keep 2-13c2 *= lifts
print(c2)
plt.plot(c2)
plt.show()

在这里插入图片描述