语音识别中的声学特征提取：梅尔频率倒谱系数MFCC | 老炮儿改名PPLOVELL | 5th

MFCC是Mel-Frequency Cepstral Coefficients的缩写，全称是梅尔频率倒谱系数。它是在1980年由Davis和Mermelstein提出来的，是一种在自动语音和说话人识别中广泛使用的特征。顾名思义，MFCC特征提取包含两个关键步骤：梅尔频率分析和倒谱分析，下面分别进行介绍。

梅尔（Mel）频率分析

对于人类听觉感知的实验表明，人类听觉的感知只聚焦在某些特定的区域，而不是整个频谱包络，而Mel频率分析就是基于人类听觉感知实验的。实验观测发现人耳就像一个滤波器组一样，它只关注某些特定的频率分量（人的听觉对频率是有选择性的）。也就说，它只让某些频率的信号通过，而压根就直接无视它不想感知的某些频率信号。但是这些滤波器在频率坐标轴上却不是统一分布的，在低频区域有很多的滤波器，他们分布比较密集，但在高频区域，滤波器的数目就变得比较少，分布很稀疏。

梅尔刻度是一种基于人耳对等距的音高(pitch)变化的感官判断而定的非线性频率刻度，和频率的赫兹的关系如下:

m=2595*log10(1+f/700)

所以当在梅尔刻度上面上是均匀分度的话，对于的赫兹之间的距离将会越来越大。梅尔刻度的滤波器组在低频部分的分辨率高，跟人耳的听觉特性是相符的，这也是梅尔刻度的物理意义所在。

倒谱分析

对于一个语音的频谱图，峰值表示语音的主要频率成分，也称为共振峰，而共振峰携带了声音的辨识属性，在语音识别中，我们需要把共振峰的位置和它们转变的过程提取出来，这个变化的过程是一条连接这些共振峰点的平滑曲线，我们称为频谱的包络（Spectral Envelope）。

原始频谱可以看成由两部分组成：包络和频谱细节，如果我们把这两部分分离开，我们就可以得到包络。那怎么把他们分离开呢？也就是，怎么在给定log X[k]的基础上，求得log H[k] 和 log E[k]以满足log X[k] = log H[k] + log E[k]呢？

这就用到了倒谱分析。倒谱（cepstrum）是一种信号的傅里叶变换经对数运算后再进行傅里叶反变换得到的谱。倒谱分析可用于将信号分解，两个信号的卷积转化为两个信号的相加。

假设上面的频率谱X(k)，时域信号为x(n)，那么满足

X(k)=DFT(x(n))

考虑将频域X(k)拆分为两部分的乘积：

X(k)=H(k)E(k)

假设两部分对应的时域信号分别是h(n)和e(n)，那么满足：

x(n)=h(n)∗e(n)

此时我们是无法区分开h(n)和e(n)，对频域两边取log：

log(X(k))=log(H(k))+log(E(k))

然后进行反傅里叶变换：

IDFT(log(X(k)))=IDFT(log(H(k)))+IDFT(log(E(k)))

得到的时域信号如下：

X’(n)=h’(n)+e’(n)

此时获得时域信号x’(n)即为倒谱，已经和原始的时域信号x(n)不一样，但是可以把时域信号的卷积关系转化为了线性加关系。倒谱分析已经将两部分对应的时域信号的卷积关系转化为了线性加关系，所以只需要将倒谱通过一个低通滤波器即可获得包络部分对应的时域信号h’(t)。

梅尔频率倒谱系数MFCC

梅尔频率倒谱系数MFCC考虑到了人类的听觉特征，先将线性频谱映射到基于听觉感知的Mel非线性频谱中，然后转换到倒谱上。即将频谱通过一组Mel滤波器就得到Mel频谱，公式表述就是log X[k] = log (Mel-Spectrum)，这时候我们再在log X[k]上进行倒谱分析，在Mel频谱上面获得的倒谱系数h[k]就称为Mel频率倒谱系数，简称MFCC。通常，计算MFCC之前，还会通过预加重、分帧和加窗、短时FFT等手段将原始原始声音信号的spectrogram声谱图，MFCC对声谱信号进行分析。

这里总结一下提取MFCC特征的过程：

1）先对语音进行预加重、分帧和加窗；

2）对每一个短时分析窗，通过FFT得到对应的频谱；

3）将上面的频谱通过Mel滤波器组得到Mel频谱；

4）在Mel频谱上面进行倒谱分析（取对数，做逆变换，实际逆变换一般是通过DCT离散余弦变换来代替上文的IDFT，取DCT后的第2个到第13个系数作为MFCC系数），获得Mel频率倒谱系数MFCC。

举个栗子：MFCC的matlab实现

[x fs]= wavread ('test.wav');
bank=melbankm(24,256,fs,0,0.4,'t');%Mel滤波器的阶数为24，fft变换的长度为256，采样频率为16000Hz
%归一化mel滤波器组系数
bank=full(bank);
bank=bank/max(bank(:));
for k=1:12                            %归一化mel滤波器组系数
n=0:23;
dctcoef(k,:)=cos((2*n+1)*k*pi/(2*24));
end
w=1+6*sin(pi*[1:12]./12);%归一化倒谱提升窗口
w=w/max(w);%预加重滤波器
xx=double(x);
xx=filter([1-0.9375],1,xx);%语音信号分帧
xx=enframe(xx,256,80);%对x 256点分为一帧
%计算每帧的MFCC参数
for i=1:size(xx,1)
y=xx(i,:);
s=y'.*hamming(256);
         t=abs(fft(s));%fft快速傅立叶变换
t=t.^2;
         c1=dctcoef*log(bank*t(1:129));
c2=c1.*w';
m(i,:)=c2';
end
%求取一阶差分系数
dtm=zeros(size(m));
for i=3:size(m,1)-2
dtm(i,:)=-2*m(i-2,:)-m(i-1,:)+m(i+1,:)+2*m(i+2,:);
end
dtm=dtm/3;
%求取二阶差分系数
dtmm=zeros(size(dtm));
for i=3:size(dtm,1)-2
dtmm(i,:)=-2*dtm(i-2,:)-dtm(i-1,:)+dtm(i+1,:)+2*dtm(i+2,:);
end
dtmm=dtmm/3;
%合并mfcc参数和一阶差分mfcc参数
ccc=[m dtm dtmm];
%去除首尾两帧，因为这两帧的一阶差分参数为0
ccc=ccc(3:size(m,1)-2,:);
ccc
% subplot(2,1,1)
% ccc_1=ccc(:,1);
% plot(ccc_1);title('MFCC');ylabel('幅值');
% [h,w]=size(ccc);
% A=size(ccc);
% subplot(212)   
% plot([1,w],A);
% xlabel('维数');
% ylabel('幅值');
% title('维数与幅值的关系')

（matlab2014之后版本，取消了wavread函数，改用audioread函数。）

【参考文章】

1、https://blog.csdn.net/zouxy09/article/details/9156785

2、https://blog.csdn.net/xmdxcsj/article/details/51228791

3、https://blog.csdn.net/cxroom/article/details/26989275