一、首先讓我們借用並澄清幾個語音學中的概念
1.臨界頻帶與聽覺掩蔽
聽覺臨界頻帶：設純音頻率為，用噪聲（設頻率為）掩蔽純音時，在噪聲湮沒的純音的過程中，起作用的是頻率在以內的噪聲，稱為臨界頻帶。即當噪聲的頻率處於上述區間時，人耳會聽不見該純音，即此頻率的噪聲對該純音的聽覺造成掩蔽。而頻率在區間之外的噪聲，人耳可以正常察覺純音，即不會發生掩蔽。
2.Mel頻率尺度
人耳對音調的感知度，不隨着頻率（Hz）的加倍而加倍，但頻率在Mel尺度內，人對音調的主觀感知度與聲音的頻率則為線性關系。MFCC考慮了人耳的聽覺特性，且沒有任何前提假設[9]。普通頻率轉換為Mel頻率的公式為：
3.濾波器組
將語音信號映射到Mel尺度，並根據人耳所具有的臨界頻帶特性的數學實現，是將每幀語音信號的功率譜，用通過一個如圖所示的濾波器組的方法完成的。
典型的濾波器組是由24個三角形帶通濾波器構成的。每個帶通濾波器具有的中心頻率和頻帶便是人耳的臨界頻帶和聽覺掩蔽特性的反映；且在不同的頻率上，每個帶通濾波器的帶寬是不同的，但在Mel尺度內，則都是等帶寬的。
所以，該濾波器組是通過給每幀的語音信號的功率譜加權而模擬人耳的聽覺特性的。注意：濾波器組的低頻段較密，高頻段較稀疏，這個目的是為了提升低頻段的能量。
二、特征參數提取的目標
特征參數提取的目標，顧名思義，就應該使相同的語音之間的差別盡可能的小，不同的語音之間的差異盡可能的大。在基於語音的線性模型的下，語音的形成可看做為聲門激勵與聲道的耦合——卷積形成的，即：
其中，x(n)為語音信號，h(n)為聲道，e(n)為激勵。
我們的任務是在現有的語音中分離聲道和激勵，即將卷積運算變換成加法運算。
作FFT后，可得到：
在x(n)的頻譜X(k)中，包絡的峰值為共振峰，表示語音的主要頻率成分，共振峰攜帶了的聲道特性，頻譜的細節部分反應了激勵源的信息，對上式取對數得：
若此處用一個低通濾波器，則可以分離聲道與激勵。但事實上，聲道與激勵同是辨識口令的特征。對其再取傅里葉逆變換（IFFT，實際處理時用的是DCT變換，因為DCT的優勢在於它對實數進行變換，得到實數，避免了復數運算）得：
稱為x(n)的倒譜，可作為語音的特征參數。以上的方法為同態處理。
最后，截取所得的
的前12維（一般取12維左右）作為該幀的特征參數，即MFCC，它包含了聲門和聲道的特性。
三、基於MATLAB的實現
function ccc=my_mfcc(x,fs,p)
% x是輸入語音序列，Mel濾波器的個數為p，采樣頻率為fs，frameSize為幀長和FFT點數
%即用統一的數值256，，，，inc為幀移為100，，，；ccc為MFCC參數。
[x1,x2] =my_vad(x); %x為原始語音數據
%frameSize=framelen; %為統一的幀長256
frameSize=256; %frameSize為幀長和FFT點數
inc=100;
% bank=melbankm(p,frameSize,fs);
% % framesize ：length of fft
% % 歸一化Mel濾波器組系數
% bank=full(bank);
% bank=bank/max(bank(:));
%%-------准備工作-------------
%歸一化mel濾波器組系數(24個窗)
bank=melbankm2(p,frameSize,fs,0,0.5,'m');
bank=full(bank);
bank=bank/max(bank(:));
% DCT系數,12*p
for k=1:12
n=0:p-1;
dctcoef(k,:)=cos((2*n+1)*k*pi/(2*p));
end
% 歸一化倒譜提升窗口
w = 1 + 6 * sin(pi * [1:12] ./ 12);
w = w/max(w);
% 語音信號分幀
xx=my_enframe(x,frameSize,inc); %需要再次分幀，使用數據
n2= fix(frameSize/2) +1 ; %計算幀長度的一半
xx=xx(x1:x2,:);
% 計算每幀的MFCC參數
for i=1:(x2-x1+1) %x1,x2只是幀的下標 %幀數，循環
y = xx(i,:);
s = y' .* hamming(frameSize);
t = abs(fft(s));
t = t.^2;
c1=dctcoef * log(bank * t(1:n2));
c2 = c1.*w';
m(i,:)=c2'; %每一行代表每一幀，每一幀里有12個數據
end
% ccc=m;
%
%差分系數
dtm = zeros(size(m)); % 幀數*12的矩陣
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]; %在m的后面再接一個矩陣，共形成了 幀數*24的參數矩陣
%去除首尾兩幀，因為這兩幀的一階差分參數為0
ccc = ccc(3:size(m,1)-2,:);
function [x1,x2] =my_vad(x)
% 預加重濾波器
xx=double(x);
xx=filter([1 -0.9375],1,xx);
%幅度歸一化到[-1,1]
xx = xx / max(abs(xx));
%常數設置
FrameLen = 256;
FrameInc = 100;
NIS=15; %設置前面無聲段的幀數為20，根據自己所采集的樣本，一般有0.5s的無聲段，<30*200/16K=0.4s，可行
amp1 = 6;
% amp2 = 0.3;
amp2 = 0.15;
%zcr1 = 10;
zcr1 = 1;
maxsilence =60; % 30*10ms = 300ms，最大靜默時長，即，未超過maxsilence的長度，仍然為有聲段
minlen = 30; % 15*10ms = 150ms ，用來檢測排查噪聲的
status = 0;
count = 0;
silence = 0;
amp=my_STEn(xx,FrameLen,FrameInc);
zcr=my_STZcr(xx,FrameLen,FrameInc);
ampth=mean(amp(1:NIS)); % 計算初始無話段區間能量和過零率的平均值
zcrth=mean(zcr(1:NIS));
ampth_v=zeros(1,length(zcr));
for i=1:length(zcr)
ampth_v(i)=ampth;
end
zcrth_v=zeros(1,length(zcr));
for j=1:length(zcr)
zcrth_v(j)=zcrth;
end
m_amp=abs(amp-ampth_v);
m_zcr=abs(zcr-zcrth_v); % 每幀的能量和過零率都減去前面都是噪聲段的相應的數值
%計算過零率
% tmp1 = enframe(x(1:end-1), FrameLen, FrameInc);
% tmp2 = enframe(x(2:end) , FrameLen, FrameInc);
% signs = (tmp1.*tmp2)<0;
% diffs = (tmp1 -tmp2)>0.02;
% zcr = sum(signs.*diffs, 2);
%計算短時能量
%amp = sum(abs(enframe(filter([1 -0.9375], 1, x), FrameLen, FrameInc)), 2);
%調整能量門限
amp1 = min(amp1, max(m_amp)/5); %高門限
amp2 = min(amp2, max(m_amp)/15); %低門限
zcr2 = min(zcr1, max(m_zcr)/5);
%開始端點檢測
x1 = 0;
x2 = 0;
for n=1:length(zcr) %幀數循環
% goto = 0;
switch status
case {0,1} % 0 = 靜音, 1 = 可能開始
if amp(n) > amp1 % 確信進入語音段
% if amp(n) > amp1 | ... % 確信進入語音段
% zcr(n) > zcr2
x1 = max(n-count-1,1);
status = 2;
silence = 0;
count = count + 1;
elseif amp(n) > amp2 | ... % 可能處於語音段
zcr(n) > zcr2
status = 1;
count = count + 1;
else % 靜音狀態
status = 0;
count = 0;
end
case 2, % 2 = 語音段
if amp(n) > amp2 | ... % 保持在語音段
zcr(n) > zcr2
count = count + 1;
else % 語音將結束
silence = silence+1;
if silence < maxsilence % 靜音還不夠長，尚未結束
count = count + 1;
elseif count < minlen % 語音長度太短，認為是噪聲
status = 0;
silence = 0;
count = 0;
else % 語音結束
status = 3;
end
end
case 3,
break;
end
end
count =count-silence/1.2; %根據經驗設定的長度，silence/3也可以
x2 = x1+count-1; %x1,x2為語音的起止幀號
% inc=FrameInc;
function frameout=my_enframe(x,win,inc)
nx=length(x(:)); % 取數據長度
nwin=length(win); % 取窗長
if (nwin == 1) % 判斷窗長是否為1，若為1，即表示沒有設窗函數
len = win; % 是，幀長=win
else
len = nwin; % 否，幀長=窗長
end
if (nargin < 3) % 如果只有兩個參數，設幀inc=幀長
inc = len;
end
nf = fix((nx-len+inc)/inc); % 計算幀數
frameout=zeros(nf,len); % 初始化
indf= inc*(0:(nf-1)).'; % 設置每幀在x中的位移量位置
inds = (1:len); % 每幀數據對應1:len
frameout(:) = x(indf(:,ones(1,len))+inds(ones(nf,1),:)); % 對數據分幀
if (nwin > 1) % 若參數中包括窗函數，把每幀乘以窗函數
w = win(:)'; % 把win轉成行數據
frameout = frameout .* w(ones(nf,1),:); % 乘窗函數
end
%短時能量計算函數
function para=my_STEn(x,win,inc)
X=my_enframe(x,win,inc)'; % 分幀
fn=size(X,2); % 求出幀數
for i=1 : fn
u=X(:,i); % 取出一幀
u2=u.*u; % 求出能量
para(i)=sum(u2); % 對一幀累加求和
end
end
%短時能量計算函數
function para=my_STEn(x,win,inc)
X=my_enframe(x,win,inc)'; % 分幀
fn=size(X,2); % 求出幀數
for i=1 : fn
u=X(:,i); % 取出一幀
u2=u.*u; % 求出能量
para(i)=sum(u2); % 對一幀累加求和
end
end
%短時過零率計算函數
function zcr=my_STZcr(x,win,inc)
X=my_enframe(x,win,inc)'; % 分幀
fn=size(X,2); % 求出幀數
if length(win)==1
wlen=win; % 求出幀長
else
wlen=length(win);
end
zcr=zeros(1,fn); % 初始化
delta=0.01; % 設置一個很小的閾值
for i=1:fn
z=X(:,i); % 取得一幀數據
for k=1 : wlen % 中心截幅處理
if z(k)>=delta
ym(k)=z(k)-delta;
elseif z(k)<-delta
ym(k)=z(k)+delta;
else
ym(k)=0;
end
end
zcr(i)=sum(ym(1:end-1).*ym(2:end)<0); % 取得處理后的一幀數據尋找過零率
end
