一、簡介

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

三、執行結果

四、備註

版本：2014a
完整程式碼或代寫加1564658423