一、简介
MFCC(Mel-frequency cepstral coefficients):梅尔频率倒谱系数。梅尔频率是基于人耳听觉特性提出来的， 它与Hz频率成非线性对应关系。梅尔频率倒谱系数(MFCC)则是利用它们之间的这种关系，计算得到的Hz频谱特征。主要用于语音数据特征提取和降低运算维度。例如：对于一帧有512维(采样点)数据，经过MFCC后可以提取出最重要的40维(一般而言)数据同时也达到了将维的目的。
MFCC一般会经过这么几个步骤：预加重，分帧，加窗，快速傅里叶变换(FFT)，梅尔滤波器组，离散余弦变换(DCT).其中最重要的就是FFT和梅尔滤波器组，这两个进行了主要的将维操作。
1.预加重
将经采样后的数字语音信号s(n)通过一个高通滤波器(high pass filter)：其中a一般取0.95左右。经过预加重后的信号为：
预加重的目的是提升高频部分，使信号的频谱变得平坦，保持在低频到高频的整个频带中，能用同样的信噪比求频谱。同时，也是为了消除发生过程中声带和嘴唇的效应，来补偿语音信号受到发音系统所抑制的高频部分，也为了突出高频的共振峰。
2.分帧
为了方便对语音分析，可以将语音分成一个个小段，称之为：帧。先将N个采样点集合成一个观测单位，称为帧。通常情况下N的值为256或512，涵盖的时间约为20~30ms左右。为了避免相邻两帧的变化过大，因此会让两相邻帧之间有一段重叠区域，此重叠区域包含了M个取样点，通常M的值约为N的1/2或1/3。通常语音识别所采用语音信号的采样频率为8KHz或16KHz，以8KHz来说，若帧长度为256个采样点，则对应的时间长度是256/8000×1000=32ms。
3.加窗
语音在长范围内是不停变动的，没有固定的特性无法做处理，所以将每一帧代入窗函数，窗外的值设定为0，其目的是消除各个帧两端可能会造成的信号不连续性。常用的窗函数有方窗、汉明窗和汉宁窗等，根据窗函数的频域特性，常采用汉明窗。
将每一帧乘以汉明窗，以增加帧左端和右端的连续性。假设分帧后的信号为S(n), n=0,1…,N-1, N为帧的大小，那么乘上汉明窗后 ，W(n)形式如下：
不同的a值会产生不同的汉明窗，一般情况下a取0.46.
4.快速傅里叶变换
由于信号在时域上的变换通常很难看出信号的特性，所以通常将它转换为频域上的能量分布来观察，不同的能量分布，就能代表不同语音的特性。所以在乘上汉明窗后，每帧还必须再经过快速傅里叶变换以得到在频谱上的能量分布。对分帧加窗后的各帧信号进行快速傅里叶变换得到各帧的频谱。并对语音信号的频谱取模平方得到语音信号的功率谱。设语音信号的DFT为：
式中x(n)为输入的语音信号，N表示傅里叶变换的点数。
这里需要先介绍下Nyquist频率，奈奎斯特频率（Nyquist频率）是离散信号系统采样频率的一半，因哈里·奈奎斯特（Harry Nyquist）或奈奎斯特－香农采样定理得名。采样定理指出，只要离散系统的奈奎斯特频率高于被采样信号的最高频率或带宽，就可以避免混叠现象。在语音系统中我通常采样率取16khz，而人发生的频率在300hz~3400hz之间，按照Nyquist频率的定义就有Nyquist频率等于8khz高于人发生的最高频率，满足Nyquist频率的限制条件。FFT就是根据Nyquist频率截取采样率的一半来计算，具体来说就是，假设一帧有512个采样点，傅里叶变换的点数也是512，经过FFT计算后输出的点数是257(N/2+1)，其含义表示的是从0(Hz)到采样率/2(Hz)的N/2+1点频率的成分。也就是说在经过FFT计算时不仅把信号从时域转到了频域并且去除了高于被采样信号的最高频率的点的影响，同时也降低了维度。
5.梅尔滤波器组
由于人耳对不同频率的敏感程度不同，且成非线性关系，因此我们将频谱按人耳敏感程度分为多个Mel滤波器组，在Mel刻度范围内，各个滤波器的中心频率是相等间隔的线性分布，但在频率范围不是相等间隔的，这个是由于频率与Mel频率转换的公式形成的，公式如下：
式中的log是以log10为底，也就是lg。
将能量谱通过一组Mel尺度的三角形滤波器组，定义一个有M个滤波器的滤波器组（滤波器的个数和临界带的个数相近），采用的滤波器为三角滤波器，中心频率为f(m),m=1,2,…,M。M通常取22-26。各f(m)之间的间隔随着m值的减小而缩小，随着m值的增大而增宽，如图所示：
式中的k指经过FFT计算后的点的下标，也就是前面例子中的0~257，f(m)也对应点的下标，具体求法如下：
1.确定语音信号最低(一般是0hz)最高(一般是采样率的二分之一)频率以及Mel滤波器个数
2.计算对应最低最高频率的mel频率
3.计算相邻两个mel滤波器中心频率的距离：(最高mel频率-最低mel频率)/(滤波器个数+1)
4.将各个中心Mel频率转成频率
5.计算频率对应FFT中点的下标
例如：假设采样率为16khz，最低频率为0hz，滤波器个数为26，帧大小为512，则傅里叶变换点数也为512，那么带入Mel频率与实际频率的转换公式中得到最低Mel频率为0，最高Mel频率为2840.02.中心频率距离为：(2840.02-0)/(26+1)=105.19，这样我们就可以得到Mel滤波器组的中心频率：[0，105.19，210.38，…，2840.02]，然后再将这组中心频率转成实际频率组(按公式操作即可，这里不列出来了)，最后计算实际频率组对应FFT点的下标，计算公式为：实际频率组中的每个频率/采样率*(傅里叶变换点数 + 1)。这样就得到FFT点下标组：[0,2,4,7,10,13,16，…，256]，也就是f(0),f(1),…,f(27)。
有了这些，我们在计算每个滤波器的输出，计算公式如下：
式中的M指滤波器的个数，N指FFT中的点数(上述的例子中是257)。经过上面的计算后每帧数据我们得到一个与滤波器个数相等的维数，降低了维数(本例中是26维)。
6.离散余弦变换
离散余弦变换经常用于信号处理和图像处理，用来对信号和图像进行有损数据压缩，这是由于离散余弦变换具有很强的"能量集中"特性：大多数的自然信号（包括声音和图像）的能量都集中在离散余弦变换后的低频部分，实际就是对每帧数据在进行一次将维。其公式如下：
将上述每个滤波器的对数能量带入离散余弦变换，求出L阶的Mel-scale Cepstrum参数。L阶指MFCC系数阶数，通常取12-16。这里M是三角滤波器个数。
7.动态差分参数的提取
标准的倒谱参数MFCC只反映了语音参数的静态特性，语音的动态特性可以用这些静态特征的差分谱来描述。实验证明：把动、静态特征结合起来才能有效提高系统的识别性能。差分参数的计算可以采用下面的公式：
式中,dt表示第t个一阶差分，Ct表示第t个倒谱系数，Q表示倒谱系数的阶数，K表示一阶导数的时间差，可取1或2。将上式的结果再代入就可以得到二阶差分的参数。
因此，MFCC的全部组成其实是由： N维MFCC参数（N/3 MFCC系数+ N/3 一阶差分参数+ N/3 二阶差分参数）+帧能量（此项可根据需求替换）。
这里的帧能量是指一帧的音量（即能量），也是语音的重要特征，而且非常容易计算。因此，通常再加上一帧的对数能量（定义：一帧内信号的平方和，再取以10为底的对数值，再乘以10）使得每一帧基本的语音特征就多了一维，包括一个对数能量和剩下的倒频谱参数。另外，解释下最开始说的40维是怎么回事，假设离散余弦变换的阶数取13，那么经过一阶二阶差分后就是39维了再加上帧能量总共就是40维，当然这个可以根据实际需要动态调整。
二、源代码
clc;
clear;
load traindata Myfeature
A1=zeros(1,30);
A2=ones(1,30);
Group=[A1,A2];
TrainData=Myfeature;
SVMStruct = svmtrain(TrainData,Group);
N=5.3;
Tw = 25; % analysis frame duration (ms)
Ts = 10; % analysis frame shift (ms)
alpha = 0.97; % preemphasis coefficient
R = [ 300 3700 ]; % frequency range to consider
M = 20; % number of filterbank channels
C = 13; % number of cepstral coefficients
L = 22; % cepstral sine lifter parameter
fs = 16000;
hamming = @(N)(0.54-0.46*cos(2*pi*[0:N-1].'/(N-1)));
[filename, pathname] = uigetfile({'*.*';'*.flac'; '*.wav'; '*.mp3'; }, '选择语音');
% %没有图像
if filename == 0
return;
end
[speech,fs] = audioread([pathname, filename]);
[voice,fs]=extractvoice_simple(speech,-30, -20,0.2);
voicex=voice(1:N*16000);
[ mfccs, FBEs, frames ] = ...
mfcc( voicex, fs, Tw, Ts, alpha, hamming, R, M, C, L );
ceps_mfccx=mfccs(:);
[cep,ER]=lpces(voicex,17,256,256); ceps_lpc=cep(2:17,:);%LPC
%[lpc,ER]=lpces(voice,12,256,256);
%ceps_lpcc=lpc2lpcc(cep);%LPCC
ceps_lpcx=ceps_lpc(:);
ceps=[ceps_mfccx(1000:2000);ceps_lpcx(1:2000)];
TestData = ceps';
languagex=svmclassify(SVMStruct,TestData);
if languagex == 1
language='Chinese'
else
language='English'
end
% t=[1:2000];
% figure
% scatter(t,ceps_lpcx(1:2000),50,'r');
% xlabel('sample point');
% ylabel('LPC');
% title('LPC features');
% hold on
% [filename, pathname] = uigetfile({'*.*';'*.flac'; '*.wav'; '*.mp3'; }, '选择语音');
% % %没有图像
% if filename == 0
% return;
% end
% [speech,fs] = audioread([pathname, filename]);
% [voice,fs]=extractvoice_simple(speech,-30, -20,0.2);
% voicex=voice(1:N*16000);
% [ mfccs, FBEs, frames ] = ...
% mfcc( voicex, fs, Tw, Ts, alpha, hamming, R, M, C, L );
% ceps_mfccx=mfccs(:);
% [cep,ER]=lpces(voicex,17,256,256); ceps_lpc=cep(2:17,:);%LPC
%
function [ H, f, c ] = trifbank( M, K, R, fs, h2w, w2h )
% TRIFBANK Triangular filterbank.
%
% [H,F,C]=TRIFBANK(M,K,R,FS,H2W,W2H) returns matrix of M triangular filters
% (one per row), each K coefficients long along with a K coefficient long
% frequency vector F and M+2 coefficient long cutoff frequency vector C.
% The triangular filters are between limits given in R (Hz) and are
% uniformly spaced on a warped scale defined by forward (H2W) and backward
% (W2H) warping functions.
%
% Inputs
% M is the number of filters, i.e., number of rows of H
%
% K is the length of frequency response of each filter
% i.e., number of columns of H
%
% R is a two element vector that specifies frequency limits (Hz),
% i.e., R = [ low_frequency high_frequency ];
%
% FS is the sampling frequency (Hz)
%
% H2W is a Hertz scale to warped scale function handle
%
% W2H is a wared scale to Hertz scale function handle
%
% Outputs
% H is a M by K triangular filterbank matrix (one filter per row)
%
% F is a frequency vector (Hz) of 1xK dimension
%
% C is a vector of filter cutoff frequencies (Hz),
% note that C(2:end) also represents filter center frequencies,
% and the dimension of C is 1x(M+2)
%
% Example
% fs = 16000; % sampling frequency (Hz)
% nfft = 2^12; % fft size (number of frequency bins)
% K = nfft/2+1; % length of each filter
% M = 23; % number of filters
%
% hz2mel = @(hz)(1127*log(1+hz/700)); % Hertz to mel warping function
% mel2hz = @(mel)(700*exp(mel/1127)-700); % mel to Hertz warping function
%
% % Design mel filterbank of M filters each K coefficients long,
% % filters are uniformly spaced on the mel scale between 0 and Fs/2 Hz
% [ H1, freq ] = trifbank( M, K, [0 fs/2], fs, hz2mel, mel2hz );
%
% % Design mel filterbank of M filters each K coefficients long,
% % filters are uniformly spaced on the mel scale between 300 and 3750 Hz
% [ H2, freq ] = trifbank( M, K, [300 3750], fs, hz2mel, mel2hz );
%
% % Design mel filterbank of 18 filters each K coefficients long,
% % filters are uniformly spaced on the Hertz scale between 4 and 6 kHz
% [ H3, freq ] = trifbank( 18, K, [4 6]*1E3, fs, @(h)(h), @(h)(h) );
%
% hfig = figure('Position', [25 100 800 600], 'PaperPositionMode', ...
% 'auto', 'Visible', 'on', 'color', 'w'); hold on;
% subplot( 3,1,1 );
% plot( freq, H1 );
% xlabel( 'Frequency (Hz)' ); ylabel( 'Weight' ); set( gca, 'box', 'off' );
%
% subplot( 3,1,2 );
% plot( freq, H2 );
% xlabel( 'Frequency (Hz)' ); ylabel( 'Weight' ); set( gca, 'box', 'off' );
%
% subplot( 3,1,3 );
% plot( freq, H3 );
% xlabel( 'Frequency (Hz)' ); ylabel( 'Weight' ); set( gca, 'box', 'off' );
%
% Reference
% [1] Huang, X., Acero, A., Hon, H., 2001. Spoken Language Processing:
% A guide to theory, algorithm, and system development.
% Prentice Hall, Upper Saddle River, NJ, USA (pp. 314-315).
% Author Kamil Wojcicki, UTD, June 2011
if( nargin~= 6 ), help trifbank; return; end; % very lite input validation
f_min = 0; % filter coefficients start at this frequency (Hz)
f_low = R(1); % lower cutoff frequency (Hz) for the filterbank
f_high = R(2); % upper cutoff frequency (Hz) for the filterbank
f_max = 0.5*fs; % filter coefficients end at this frequency (Hz)
f = linspace( f_min, f_max, K ); % frequency range (Hz), size 1xK
fw = h2w( f );
% filter cutoff frequencies (Hz) for all filters, size 1x(M+2)
c = w2h( h2w(f_low)+[0:M+1]*((h2w(f_high)-h2w(f_low))/(M+1)) );
cw = h2w( c );
H = zeros( M, K ); % zero otherwise
for m = 1:M
% implements Eq. (6.140) on page 314 of [1]
% k = f>=c(m)&f<=c(m+1); % up-slope
% H(m,k) = 2*(f(k)-c(m)) / ((c(m+2)-c(m))*(c(m+1)-c(m)));
% k = f>=c(m+1)&f<=c(m+2); % down-slope
% H(m,k) = 2*(c(m+2)-f(k)) / ((c(m+2)-c(m))*(c(m+2)-c(m+1)));
% implements Eq. (6.141) on page 315 of [1]
k = f>=c(m)&f<=c(m+1); % up-slope
H(m,k) = (f(k)-c(m))/(c(m+1)-c(m));
k = f>=c(m+1)&f<=c(m+2); % down-slope
H(m,k) = (c(m+2)-f(k))/(c(m+2)-c(m+1));
end
三、运行结果
四、备注
完整代码添加QQ1575304183