Mel频率倒谱系数-MFCC

MFCC：Mel频率倒谱系数的缩写。

目的：模拟人耳对不同频率语音的感知

Mel频率和Hz频率的关系

人类对不同频率语音有不同的感知能力:

1kHz以下，与频率成线性关系。
1kHz以上，与频率成对数关系。

Mel频率倒谱系数(MFCC)则是利用它们之间的这种关系，计算得到的Hz频谱特征。

由于Mel频率与Hz频率之间非线性的对应关系，使得MFCC随着频率的提高，其计算精度随之下降。因此，在应用中常常只使用低频MFCC，而丢弃中高频MFCC。

Fmel=2595lg(1+f/700)F_{mel}=2595lg(1+f/700) Fmel=2595lg(1+f/700)
FmelF_{mel}Fmel是以美尔（Mel）为单位的感知频率；f是以Hz为单位的实际频率。

Mel滤波器组

具体频率组的公式后面再说

在matlab中有个voicebox工具箱中有一个melbankm函数可用于计算Mel滤波器组。

h=melbankm(p,n,fs,fl,fh,w);
%{
fs是采样频率
fh是设计滤波器最高频率(用fs进行归一，一般为0)
fl是设计滤波器最低频率(用fs进行归一，一般为0.5)
p是在fl和fh之间设计mel滤波器的个数
n是一帧FFT后数据的长度
w是窗函数，'t'是三角窗，'n'和'm'是海明窗和汉宁窗
输出h是滤波器的频域响应，是一个p*(n/2+1)的数组，p个滤波器就有p个滤波器的频域响应，每个滤波器的响应曲线长为n/2+1，相当于取正频率的部分
}%

我们可以用代码的形式画出mel滤波器组的样子：

% pr3_3_1
clear all; clc; close all;% 调用melbankm函数,在0-0.5区间设计24个Mel滤波器,用三角形窗函数
bank=melbankm(24,256,8000,0,0.5,'t');
bank=full(bank);
bank=bank/max(bank(:));              % 幅值归一化df=8000/256;                         % 计算分辨率
ff=(0:128)*df;                       % 频率坐标刻度
for k=1 : 24                         % 绘制24个Mel滤波器响应曲线plot(ff,bank(k,:),'k','linewidth',2); hold on;
end
hold off; grid;
xlabel('频率/Hz'); ylabel('相对幅值')
title('Mel滤波器组的频率响应曲线')

得出下图

mfcc特征参数的提取

步骤：

1、预处理（预加重、分帧、加窗）

2、快速傅里叶变换

3、计算谱线能量

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

5、计算DCT倒谱

1、预处理

之前讲过分帧分窗了

预加重：目的是为了补偿高频分量的损失，提升高频分量。预加重的滤波器常设为
H(z)=1−az−1H(z)=1-az^{-1} H(z)=1−az−1
式中，a为一个常数。

分帧：由于语音信号是一个准稳态信号，把它分成较短的帧，在每一帧中可将其看做稳态信号，可用处理稳态信号的方法来处理。同时，为了使一帧到另一帧之间参数能较平稳的过渡，在两帧之间互相有部分重叠。

加窗：目的是减少频域中的泄漏，将对每一帧语音乘以汉明窗或海宁窗。语音信号处理x(n)经处理后为xi(m)x_i(m)xi(m),其中下标为第i帧。

2、快速傅里叶变换

将时域数据变成频域数据
X(i,k)=FFT(xi(m))X(i,k)=FFT(x_i(m)) X(i,k)=FFT(xi(m))

3、计算谱线能量

对每一帧FFT后的数据计算谱线的能量：
E(i,k)=[X(i,k)]2E(i,k)=[X(i,k)]^2 E(i,k)=[X(i,k)]2

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

把求出的每帧谱线能量谱通过Mel滤波器，并计算在该Mel滤波器中的能量。在频域中相当于把每帧的能量谱E(i,k)(其中i表示第i帧，k表示频谱中第k条谱线)与滤波器的频域响应Hm(k)H_m(k)Hm(k)相乘并相加：
S(i,m)=∑k=0N−1E(i,k)HmK,0≤m<MS(i,m)=\sum_{k=0}^{N-1}E(i,k)H_mK, 0\leq m<M S(i,m)=k=0∑N−1E(i,k)HmK,0≤m<M

5、计算DCT倒谱

FFT倒谱x^(n)\hat x(n)x^(n)为:
x^(n)=FT−1[X^(k)]\hat x(n)=FT^{-1}[\hat X(k)] x^(n)=FT−1[X^(k)]
式中，
X^(k)=lnFT[x(n)]=lnX(k)\hat X(k)=ln{FT[x(n)]}=ln{X(k)} X^(k)=lnFT[x(n)]=lnX(k)
而序列x(n)的DCT为：
X(k)=2N∑n=0N−1C(k)x(n)cos[π(2n+1)k2N],k=0,1,⋯,N−1X(k)=\sqrt{\frac{2}{N}}\sum^{N-1}_{n=0}C(k)x(n)cos[\frac{\pi(2n+1)k}{2N}],k=0,1,\cdots,N-1 X(k)=N2n=0∑N−1C(k)x(n)cos[2Nπ(2n+1)k],k=0,1,⋯,N−1
式中，参数N是序列x(n)的长度，C(k)是正交因子，可表示为：
C(k)={22，k=01，k=1,2,⋯,N−1C(k)=\begin{cases} \frac{\sqrt2}{2}，k=0\\ 1， k=1,2,\cdots,N-1\end{cases} C(k)={22，k=01，k=1,2,⋯,N−1
DCT倒谱和FFT倒谱很相似。把mel滤波器的能量取对数后，再计算DCT：
mfcc(i,n)=2M∑m=0M−1log[S(i,m)]cos[πn(2m−1)2M],k=0,1,⋯,N−1mfcc(i,n)=\sqrt{\frac{2}{M}}\sum^{M-1}_{m=0}log[S(i,m)]cos[\frac{\pi n(2m-1)}{2M}],k=0,1,\cdots,N-1 mfcc(i,n)=M2m=0∑M−1log[S(i,m)]cos[2Mπn(2m−1)],k=0,1,⋯,N−1
式中S(i,m)是第4步求出的Mel滤波器能量。m是指第m个Mel滤波器，i是指第i帧，n是DCT后的谱线。

这样就计算出了MFCC参数。

6、倒谱提升窗口

在为每帧数据计算出K阶MFCC参数后，通常还要为这K个系数分别乘以不同的权系数ω，实际上是一个短的窗口。
ωm=1+K2sin(πmK),1≤m≤K\omega_m=1+\frac{K}{2}sin(\frac{\pi m}{K}),1\leq m\leq K ωm=1+2Ksin(Kπm),1≤m≤K

c^m=ωmCm\hat c_m=\omega_mC_m c^m=ωmCm

7、差分倒谱参数

标准的MFCC参数只反映了语音参数的静态特性，而人耳对语音的动态特征更为敏感，如何获得反映语音动态变化的参数，通常是用差分倒谱参数来描述这种特性。

差分参数的计算公式如下：
d(n)=1∑i=−kki2∑i=−kki⋅c(n+i)d(n)=\frac{1}{\sqrt{\sum^k_{i=-k}i^2}}\sum^k_{i=-k}i·c(n+i) d(n)=∑i=−kki21i=−k∑ki⋅c(n+i)
这里c和d都表示一帧语音参数，k为常数，通常取2，这时差分参数就称为当前帧的前两帧和后两帧参数的线性组合。

由此公式可求出差分参数的一阶MFCC差分参数，用同样的公式对一阶差分参数进行计算，可以得到二阶差分MFCC参数。实际应用中，将MFCC参数和各阶差分参数合并为一个完整的矢量，作为一帧语音信号的参数。

mfcc代码

代码如下：

function ccc=mfcc_m(x,fs,p,frameSize,inc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 function ccc=mfcc_m(x);
%对输入的语音序列x进行MFCC参数的提取，返回MFCC参数和一阶
%差分MFCC参数，Mel滤波器的个数为p，采样频率为fs
%对x每frameSize点分为一帧，相邻两帧之间的帧移为inc
% FFT的长度为帧长
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 按帧长为frameSize，Mel滤波器的个数为p，采样频率为fs
% 提取Mel滤波器参数，用汉明窗函数
bank=melbankm(p,frameSize,fs,0,0.5,'m');
% 归一化Mel滤波器组系数
bank=full(bank);
bank=bank/max(bank(:));% DCT系数,12*p
for k=1:12n=0:p-1;dctcoef(k,:)=cos((2*n+1)*k*pi/(2*p));
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,frameSize,inc);
n2=fix(frameSize/2)+1;
% 计算每帧的MFCC参数
for i=1:size(xx,1)y = xx(i,:);s = y' .* hamming(frameSize);t = abs(fft(s));t = t.^2;c1=dctcoef * log(bank * t(1:n2));c2 = c1.*w'; % 倒谱提升m(i,:)=c2';
end%差分系数
dtm = zeros(size(m));
for i=3:size(m,1)-2dtm(i,:) = -2*m(i-2,:) - m(i-1,:) + m(i+1,:) + 2*m(i+2,:);
end
dtm = dtm / 3;
%合并MFCC参数和一阶差分MFCC参数
ccc = [m dtm];
%去除首尾两帧，因为这两帧的一阶差分参数为0
ccc = ccc(3:size(m,1)-2,:);