1.概述
Mel倒谱系数（Mel-frequency cepstral coefficients，MFCC）是受人的听觉系统研究成果推动而导出的声学特征。 研究发现，当两个音调的频率差小于临界带宽时，人就会把两个音调听成一个(屏蔽效应)。Mel刻度是对这一临界带宽的 度量方法之一, MFCC在语音识别领域应用广泛。本文详细介绍了Mel频率倒谱系数参数的6大提取步骤。
2.什么是Mel频率倒谱系数？
Mel频率倒谱系数（Mel Frequency Cepstrum Coefficient）的缩写是MFCC，Mel频率是基于人耳听觉特性提出来的， 它与Hz频率成非线性对应关系。Mel频率倒谱系数(MFCC)则是利用它们之间的这种关系，计算得到的Hz频谱特征。
用录音设备录制一段模拟语音信号后，经由自定的取样频率(如8000 Hz、16000 Hz等)采样后转换(A/D)为数字语音信号。 由于在时域(time domain)上语音信号的波形变化相当快速、不易观察，因此一般都会在频域(frequency domain)上来观察， 其频谱是随着时间而缓慢变化的，因此通常可以假设在一较短时间中，其语音信号的特性是稳定的，通常我们定义这个较短 时间为一帧(frame)，根据人的语音的音调周期值的变化，一般取10~20ms。
3.Mel频率倒谱系数(MFCC)参数的提取步骤
(1) 预加重(pre-emphasis)
将经采样后的数字语音信号s(n)通过一个高通滤波器(high pass filter)：
(2) 分帧(frame blocking)
一般取10-30ms为一帧，为了避免窗边界对信号的遗漏，因此对帧做偏移时候，要有帧迭(帧与帧之间需要重叠一部分)。 一般取帧长的一半作为帧移，也就是每次位移一帧的二分之一后再取下一帧，这样可以避免帧与帧之间的特性变化太大。
(3) 计算短时能量(energy)
短时能量代表着音量的高低，亦即声音振幅的大小，可以根据此能量的值来过滤掉语音信号中的一些细微噪声。当一帧的能量 值低于我们定的门槛值(threshold)时，则将此帧作为静音段(silence)。
(4) 加窗(hamming window)
语音在长范围内是不停变动的，没有固定的特性无法做处理，所以将每一帧代入窗函数，窗外的值设定为0，其目的是消除各个 帧两端可能会造成的信号不连续性。常用的窗函数有方窗、汉明窗和汉宁窗等，根据窗函数的频域特性，常采用汉明窗。公式是在 加窗范围内，$w(i)=0.54-0.46*cos(2* \pi * \frac{i}{n-1}), i \in [0,n-1]$。
(5) 快速傅立叶变换(FFT transform)
由于语音信号在时域上的变化快速而不稳定，所以通常都将它转换到频域上来观察，此时它的频谱会随着时间作缓慢的变化。 所以通常将加窗后的帧经过FFT (Fast Fourier Transform)求出每帧的频谱参数。
(6) 三角形带通滤波器(triangular band-pass filter)
将每帧的频谱参数通过一组N个三角形带通滤波器(N一般为20~30个)所组成的梅尔刻度滤波器，将每个频带的输出取对数，求出每一个 输出的对数能量(log energy)，k = 1,2,… N。 再将此N个参数进行余弦变换(cosine transform)求出L阶的Mel-scale cepstrum参数。
4.Matlab程序实现
function r = mfcc(s, fs)
% MFCC参数提取
% Reference: 论文《MFCC和LPCC特征参数在说话人识别中的研究》
%
% Inputs: s contains the signal to analize
% fs is the sampling rate of the signal
%
% Output: r contains the transformed signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = 256; % 帧长
m = 100; % 帧移
l = length(s); % 信号总长度
nbFrame = floor((l - n) / m) + 1; % 信号总帧数
for i = 1:n
for j = 1:nbFrame
M(i, j) = s(((j - 1) * m) + i); % 分帧
end
end
h = hamming(n); % Hamming窗w = 0.54 - 0.46*cos(2*pi*x);
M2 = diag(h) * M; % 对M加窗，形成对角矩阵M2
for i = 1:nbFrame
frame(:,i) = fft(M2(:, i)); % 进行快速傅里叶变换，将其转换到频域上
end
t = n / 2;
% tmax = l / fs;
m = melfb(20, n, fs); % 调用20阶MEL滤波器组进行滤波
n2 = 1 + floor(t);
z = m * abs(frame(1:n2, : )).^2; % 取前n2行帧列向量的平方
r = dct(log(z)); % 取对数后进行反余弦变换
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function m = melfb(p, n, fs)
% MELFB Determine matrix for a mel-spaced filterbank
%
% Inputs: p number of filters in filterbank
% n length of fft
% fs sample rate in Hz
%
% Outputs: x a (sparse) matrix containing the filterbank amplitudes
% size(x) = [p, 1+floor(n/2)]
%
% Usage: For example, to compute the mel-scale spectrum of a
% colum-vector signal s, with length n and sample rate fs:
% f = fft(s);
% m = melfb(p, n, fs);
% n2 = 1 + floor(n/2);
% z = m * abs(f(1:n2)).^2;
%
% z would contain p samples of the desired mel-scale spectrum
%
% To plot filterbanks e.g.:
% plot(linspace(0, (12500/2), 129), melfb(20, 256, 12500)'),
% title('Mel-spaced filterbank'), xlabel('Frequency (Hz)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f0 = 700 / fs;
fn2 = floor(n/2);
lr = log(1 + 0.5/f0) / (p+1);
% convert to fft bin numbers with 0 for DC term
bl = n * (f0 * (exp([0 1 p p+1] * lr) - 1));
b1 = floor(bl(1)) + 1;
b2 = ceil(bl(2));
b3 = floor(bl(3));
b4 = min(fn2, ceil(bl(4))) - 1;
pf = log(1 + (b1:b4)/n/f0) / lr;
fp = floor(pf);
pm = pf - fp;
r = [fp(b2:b4) 1+fp(1:b3)];
c = [b2:b4 1:b3] + 1;
v = 2 * [1-pm(b2:b4) pm(1:b3)];
m = sparse(r, c, v, p, 1+fn2);
参考文献
[1] http://www.semxi.com/TechnologyDetail.aspx?nID=27 [2] MFCC和LPCC特征参数在说话人识别中的研究