语音信号处理—MFCC特征提取(matlab代码)

        MFCC 参数考虑了人耳的听觉特性，将频谱转化为基于梅尔频标的非线性频谱，然后转换到倒谱域上。由于充分考虑了人的听

觉特性，而且没有任何前提假设，MFCC 参数具有良好的识别性能和抗噪能力。

       由于人类对于声音高低的的感知强度与该声音的频率的对数近似成正比，梅尔频率正是体现出了这种声音频率与人类感知声音高

低的对应关系。在梅尔频率域内，人对声音高低的感知程度是线性的，梅尔频率与语音频率的关系式如下： 

                                                                                                  f(mel)=2595*lg ( 1 + f / 700 )

主要过程有：

                    输入语音——预加重——分帧、加窗——FFT、|FFT|^2——mel滤波器组——ln( )、DCT——输出系数

(1).预加重：人说话时唇端辐射使高频能量丢失，预加重可以提升高频部分，使信号频谱变得平坦，利于分析。一般通过传递函数为的一阶高通滤波器来实现预加重，对应的差分方程为。其中为 预加重系数，一般0.9<<1.0。这里取0.98。

(2).分帧、加窗：由于发音器官的惯性运动，可认为语音在一小段时间（10-30ms）内频谱特征和一些物理特征参数基本保持不变，即语音的短时平稳性，因此可以将平稳过程的处理方法和理论引入到语音信号的短时处理中，将语音信号划分为很多短时的语音段，每个短时的语音段称为一帧。为了减小吉布斯效应，需要加窗处理，一般采用hamming窗。

(3).FFT、|FFT|^2:对每一帧进行FFT变换，再取绝对值，求平方，得幅度谱。

(4).mel滤波器组：将这些幅度谱通过mel滤波器组(mel滤波器组系数需先计算好)。

(5).然后取对数，再做DCT变换。

详细原理过程见：http://write.blog.csdn.net/postedit/44623627

matlab代码一：

x=wavread('E:\\测试语音\\kugou.wav');
%%%%%%%%%MEL三角滤波参数%%%%%%%%%%%%%%%%%%%%%%%%%%%
fh=4000; %fs=8000Hz,fh=4000Hz    语音信号的频率一般在300-3400Hz，所以一般情况下采样频率设为8000Hz即可。
max_melf=2595*log10(1+fh/700);
M=24;%三角滤波器的个数。
N=256;%256是语音片段的长度，也是fft的点数，也是24个滤波器(0-4000Hz)总点数
i=0:25;
f=700*(10.^(max_melf/2595*i/(M+1))-1);%将mei频域中的 各滤波器的中心频率 转到实际频率
F=zeros(24,256);
for m=1:24
    for k=1:256
        i=fh*k/N;
        if (f(m)<=i)&&(i<=f(m+1))
            F(m,k)=(i-f(m))/(f(m+1)-f(m));
        else if (f(m+1)<=i)&&(i<=f(m+2))
                F(m,k)=(f(m+2)-i)/(f(m+2)-f(m+1));
            else
                F(m,k)=0;
            end
        end
    end
end
plot((1:256)*4000/256,F);
%%%%%%%%%%%%%%%DCT系数%%%%%%%%%%%
dctcoef=zeros(12,24);
for k=1:12
  n=1:24;
  dctcoef(k,:)=cos((2*n-1)*k*pi/(2*24));
end

%%%%%%%%%%%对语音信号进行预加重处理%%%%%%%%%%
len=length(x);
alpha=0.98;
y=zeros(len,1);
for i=2:len
    y(i)=x(i)-alpha*x(i-1);
end

%%%%%%%%%%%%%%%MFCC特征参数的求取%%%%%%%%%%%%
h=hamming(256);%256*1
num=N/2;
count=floor(len/num)-1;
c1=zeros(count,12);
for i=1:count
    x_frame=y(num*(i-1)+1:num*(i-1)+N);%256*1
    w = x_frame.* h;%
    Fx=abs(fft(x_frame));
    s=log(F*Fx.^2);%ln
    c1(i,:)=(dctcoef*s)';  
end

%%%%%%%%%%%%差分系数%%%%%%%%%%%
dtm = zeros(size(c1));
for i=3:size(c1,1)-2
  dtm(i,:) = -2*c1(i-2,:) - c1(i-1,:) + c1(i+1,:) + 2*c1(i+2,:);
end
dtm = dtm / 3;

%%%%合并mfcc参数和一阶差分mfcc参数%%%%%
ccc = [c1 dtm];
%去除首尾两帧，因为这两帧的一阶差分参数为0
ccc = ccc(3:size(c1,1)-2,:);

matlab代码二:

其中有个地方不太懂

function ccc = mfcc(x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 function ccc=mfcc(x);
%对输入的语音序列x进行MFCC参数的提取，返回MFCC参数和一阶
%差分MFCC参数，Mel滤波器的阶数为24
%fft变换的长度为256，采样频率为8000Hz，对x 256点分为一帧
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bank=melbankm(24,256,8000,0,0.5,'m');%24:滤波器个数  256: length of fft  8000:采样频率
%size(bank);%24*129   为什么不是24*256？？
% 归一化mel滤波器组系数
bank=full(bank);
bank=bank/max(bank(:));
tem=(0:floor(256/2))*8000/256;
plot(tem,bank)   %fs=sample frequency

% DCT系数,12*24
for k=1:12
  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);%y(n)=x(n)-0.9375x(n-1);

% 语音信号分帧
xx=enframe(xx,256,80);

% 计算每帧的MFCC参数
frame_num=size(xx,1);%size(xx,1)返回xx的行数  249*256
h=hamming(256);
for i=1:frame_num
  y = xx(i,:); 
  s = y'.* h;%乘窗  x.*y是对应位置的元素相乘  例如：A=[1 2; 3 4];B=A;A.*B=[1*1 2*2;3*3 4*4]=[1 4;9 16]而A*B=[1*1+2*3,1*2+2*4;3*1+4*3,3*2+4*4]=[7 10;15 22]
  t = abs(fft(s));
  t = t.^2;      %计算能量%功率谱
  mel=log(bank * t(1:129));
  c1=dctcoef * mel;%这里选择的帧长为256点，然后FFT的点数也为256，由于是对称的，所以只取前面一半的点计算频谱。然后加入到三角滤波器中。
  c2 = c1.*w'; % 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;

%合并mfcc参数和一阶差分mfcc参数
ccc = [m dtm];
%去除首尾两帧，因为这两帧的一阶差分参数为0
ccc = ccc(3:size(m,1)-2,:);