声事件识别之mfcc特征提取的研究

一、    基本介绍

MFCC（Mel-Frequency Cepstral Coefficients）包含两步：第一将频率转化为梅尔频率；第二进行倒谱分析。它是基于人耳听觉特性提出来的，它与Hz频率成非线性对应关系，充分考虑了人耳对声音的听觉感知特性，目前主流的声事件分类方法之一，论文中常用其做对比算法。

二、    主要应用

该特征主要应用于语音特征提取和降低运算维度（例如将1帧512维的数据经过mfcc算法后降低到12维）。常常作为语音特征用于语音识别，对于基于生活环境中的声事件分类，应用也比较广泛。

三、    算法步骤介绍

1、         梅尔频率

梅尔刻度是一种基于人耳的听觉特性而定的非线性频率刻度。和频率（赫兹）的关系如下: 

当在梅尔刻度上面上是均匀分度的话，对于的赫兹之间的距离将会越来越大，所以梅尔刻度的滤波器组的尺度变化如下：

梅尔刻度的滤波器组在低频部分的分辨率高，跟人耳的听觉特性是相符合的，这也是梅尔刻度的物理意义所在。 

这一步的意义在于，首先对时域信号进行傅里叶变换转换到频域，然后再利用梅尔频率滤波器组将频域信号进行切分，最后在每个频率段得到一个梅尔频率。

2、         倒谱分析

倒谱的含义是：对时域信号做傅里叶变换，然后取log，然后再进行反傅里叶变换即得到倒谱，倒谱可以分为复倒谱、实倒谱和功率倒谱，我们所用的是功率倒谱。倒谱分析用于将信号分解，即将两个信号的卷积转化为两个信号的相加。

 

上面第一幅图为原始信号的频谱X(k)，第二幅图为该信号的频谱包络H（k）,第三幅图为该信号频谱的细节E(K)。那么这个信号的频谱可以看成频谱包络与频谱细节的乘积：

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

然后假设频谱包络和频谱细节的时域信号是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)，于此只需要进行一次离散余弦反变换即可达到以上两步效果之和。

四、    基本流程

1、         预加重

预加重在语音信号是为了消除发声过程中，声带和嘴唇造成的效应，以此来补偿语音信号受到发音系统所压抑的高频部分。并且能突显高频的共振峰。 我们的识别对象虽然是环境音，该步骤也课用来抵消麦克风对高频信号压制。

2、         分帧

目的：大多数声音信号属于非平稳信号，为了处理此类信号需要使用短时傅里叶变换工具，将信号看成短时平稳的，因此将信号分成很多帧便于处理。帧长一般取10-30ms，最后的帧若不够一帧所定的点数要补零，因为做fft需要信号长度为2^n。

关键问题：分帧在能得到局部时频特征的同时还会导致频域分辨率降低，甚至可能会导致频谱扩展。假如采样频率fs=44.1K，每帧包含的点数N=22050，那么频谱分辨率为fs/N=2hz（正常来说每帧不会有如此多的点，这势必会影响频谱分辨率）。若输入为纯音999hz，此时由于频谱分辨率的不足将导致在998hz  996hz  1000hz 1002hz处出现幅值，造成频谱扩展。

3、         加窗

加窗的目的是平滑时域信号更好地满足FFT处理的周期性要求，减少泄漏。因而窗函数的选择也非常关键，使用汉明窗加以平滑的话，相比于矩形窗函数，会减弱FFT以后旁瓣大小以及频谱泄露。频谱泄露：做短时傅里叶变换时需要截取其中一部分（截断）信号，然后将该信号做周期延拓，此时加窗函数后相当与在时域乘上该窗函数，而在频域就相当于频域卷积，于是频谱中除了本来该有的主瓣之外，还会出现本不该有的旁瓣，这就是频谱泄露。为了减弱频谱泄露，可以采用加权的窗函数如平顶窗、汉宁窗、高斯窗等等。而未加权的矩形窗导致的泄露最为严重。在这一般使用汉明窗（a=0.46）。

4、         时频变换

目的：快速傅里叶变换将时域信号映射到频域，然后转换为功率谱。便于后面梅尔尺度的变换。

5、梅尔滤波器组过滤信号

   因为频域信号有很多冗余，三角滤波器组不仅可以对频谱进行平滑还可以对频域的幅值进行精简，每一个频段用一个值来表示。

对于FFT得到的功率谱，分别跟每一个滤波器进行频率相乘累加，得到的值即为该帧数据在该滤波器对应频段的能量值。如果滤波器的个数为22，那么此时应该得到22个能量值，注意：使用22个三角滤波器得到的梅尔滤波器组并不意味着得到的mfcc参数是22个值或者11个值，其mfcc的特征数量取决于余弦离散变换，得到最终的mfcc数量。

6、能量值取log

由于人耳对声音的感知并不是线性的，用log这种非线性关系更好描述。取完log才可以进行倒谱分析。

7、离散余弦变换

按照倒谱的定义，该步需要进行反傅里叶变换然后通过低通滤波器获得最后的低频信号。然而在这里直接使用DCT就可以获取频率谱的低频信息。由于滤波器之间是有重叠的，所以前面的获得的能量值之间是具有相关性的，DCT还可以对数据进行降维压缩和抽象，获得最后的特征参数。相比于傅里叶变换，离散余弦变换的结果没有虚部，更好计算。

    将每个滤波器的对数能量带入离散余弦变换后，求出L阶的Mel-scale Cepstrum参数。L阶指MFCC系数阶数，通常取12-16。。

8、差分

由于声音信号是时域连续的，mfcc提取的特征信息只反应了本帧声音的静态特性，为了使特征更能体现声音的动态特性，可以在特征维度增加一阶差分和二阶差分。

五、    总结

MFCC参数作为一种经典的声音特征不但能做到对声音信号的压缩还能体现声音信号中重要的特征，所以，我对此作了详细研究，并在网上查找了相关程序，已经在MATLAB中实现对该特征的提取。

本文对于快速傅里叶变换、梅尔三角滤波器组、离散余弦变换并未做进行详细的原理研究。

六、    所遇到的问题

1、   环境音的声音库太少，虽然声音种类多，但是每一种的样本数少。远远达不到深度学习的数据数量要求。

七、MATLAB程序

 

% MFCC特征提取
[x fs]=audioread('C:\Users\lifangzu\Desktop\car crash\play(4).mp3');  %读取音频信号 x为信号样本，fs信号采样频率
%% 判断并提取单声道
if size(x,2)>1
    x=x(:,1);
else
    x=x;
end
bank=melbankm(24,1024,fs,0,0.4,'t'); %Mel滤波器的阶数为24，FFT变换的长度为256，采样频率为fs  
%归一化Mel滤波器组系数  
bank=full(bank); %将稀疏矩阵转化为完全矩阵
bank=bank/max(bank(:));  %滤波器组归一化
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);%预加重处理
%语音信号分帧  
xx=enframe(xx,1024,512);%对xx1024点分为一帧  
%计算每帧的MFCC参数  
for i=1:size(xx,1)  
    y=xx(i,:);  
    s=y'.*hamming(1024);  
    t=abs(fft(s));%FFT快速傅里叶变换  
    t=t.^2;  
    c1=dctcoef*log(bank*t(1:513));  
    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,:);  

（文章经多方总结得到，若有错误，望批评指正）