注意！這個博客里給出的文件是 在matlab2012版本里使用的，而我學習時使用的時matlab2018b，因此發現了很多新版本不兼容的問題，但我沒有給出修改后能直接在新版本中用的代碼。
另外，matlab要安裝voicebox才可以正常進行實驗！
還有一件事，錄音文件要修改好文件名和文件目錄，不然會檢測不到。
這個實驗在沒有同學和老師的指導和幫助下完成，一共花了大概兩周多（也沒有每天做啦），請珍惜勞動成果，不要直接拿去裝逼。
有個關鍵函數提前說一下：
DTW函數：myDTW
功能：利用DTW完成模式匹配
調用格式：
[cost] = myDTW(F,R)
參數說明：輸入參數F為模板的MFCC特征；R為當前待識別語音的MFCC。輸出參數cost為F和R之間的匹配距離。
看看目錄
下面是具體的程序清單：
function [cost] = myDTW(F,R) [r1,c1] = size(F); %模板的維度 [r2,c2] = size(R); %當前待識別語音的維度 distance = zeros(r1,r2); for n = 1:r1 for m=1:r2 FR = F(n,:)-R(m,:)； FR = FR.^2; distance(n,m) = sqrt(sum(FR))/c1; %采用歐式距離度量F和R的相似程度 end end D=zeros(r1+1,r2+1); D(1,:)=inf; D(:,1)=inf; D(1,1)=0; D(2:(r1+1),2:(r2+1)) =distance; %尋找整個過程的最短匹配距離 for i=1:r1; for j= 1:r2; [dmin] = min([D(i,j),D(i,j+1),D(i+1,j)]); D(i+1,j+1) =D(i+1,j+1) +dmin; end end cost =D(r1+1,r2+1); %返回最終距離,沒有反向回溯
這個基於動態時間規整（DTW）的孤立字語音識別的實驗，具體步驟有很多，都寫在下面了：
（1）訓練過程：運行setTemplates.m文件
1）實驗的語音來自三個人，每人說一遍0~9這十個數字，記錄下相應的語音，保存在 目錄“SpeechData”中。以第一個人的語音“1”為例，首先，調用[speechIn1,FS1] =audioread(q)讀入語音“1”，保存在speechIn1中，並獲得采樣率FS1。
2）對speechIn1進行預處理，首先調用my_vad函數進行端點檢測，提取出有效語音段，然后調用mfccf函數對有效語音段提取mfcc特征。
3）以同樣的方式處理其他訓練語音，把第一/二/三個說話人的十個數字0~9語音分別保存s1/s2/s3結構體中，並最終保存在Vectors1.mat/Vectors2.mat/Vectors3.mat文件中。
（2）識別過程：運行matchTemplates.m文件
1） 首先程序提示“Press any key to start 2 seconds of speech recording…”,在MATLAB主窗口按任意鍵開始錄入待識別的語音，程序將錄入的語音保存在speechIn變量中。
2） 對speechIn進行預處理，包括端點檢測和mfcc特征提取，這一過程和訓練中完全相同，得到（NP）矩陣。其中N為幀數，P為維度，接着調用CMN函數，對每一維度分量進行標准化，結果保存在rMatrix中。
3） 調用DTWScores函數，將當前待識別語音的rMatrix與訓練過程中保存的s1、s2、s3中 的0~9各數字匹配，此函數中調用myDTW的匹配過程，返回rMarix和每個說話人說的每一個數字的匹配分值，保存在Sco中（（3*10）的矢量）。
4） 將Sco從大到小進行排序，首先得到每一個模板匹配的兩個最低分值對應的數字序號（有2*3=6個序號），而后返回出現頻率最高的序號作為判決結果，從而完成過程。
但是講這么多，還是很難把實驗完成對不對，現在我把全部用到的都放出來
識別過程matchTemplates.m程序清單：
clear all;
close all;
ncoeff = 12; %MFCC參數階數
N = 10; %10個數字
fs=16000; %采樣頻率
duration2 = 2; %錄音時長
k = 3; %訓練樣本的人數
speech = audiorecorder(fs,16,1);
disp('Press any key to start 2 seconds of speech recording...');
pause
disp('Recording speech...');
recordblocking(speech,duration2) % duration*fs 為采樣點數
speechIn=getaudiodata(speech);
disp('Finished recording.');
disp('System is trying to recognize what you have spoken...');
speechIn = my_vad(speechIn); %端點檢測
rMatrix1 = mfccf(ncoeff,speechIn,fs); %采用MFCC系數作為特征矢量
rMatrix = CMN(rMatrix1); %歸一化處理
Sco = DTWScores(rMatrix,N); %計算DTW值
[SortedScores,EIndex] = sort(Sco,2); %按行遞增排序，並返回對應的原始次序
Nbr = EIndex(:,1:2) %得到每個模板匹配的2個最低值對應的次序
[Modal,Freq] = mode(Nbr(:)); %返回出現頻率最高的數Modal及其出現頻率Freq
Word = char('zero','One','Two','Three','Four','Five','Six','Seven','Eight','Nine');
if mean(abs(speechIn)) < 0.01
fprintf('No microphone connected or you have not said anything.\n');
else if (Freq <2) %頻率太低不確定
fprintf('The word you have said could not be properly recognised.\n');
else
fprintf('You have just said %s.\n',Word(Modal,:));
end
訓練過程setTemplates----程序清單：
%設置模板
ncoeff=12; %mfcc系數的個數
fMatrix1 = cell(1,10);
fMatrix2 = cell(1,10);
fMatrix3 = cell(1,10);
for i = 1:10
q = ['F:\ SpeechData\p1\' num2str(i-1) '.wav'];
[speechIn1,FS1] = audioread(q);
speechIn1 = my_vad(speechIn1);
fMatrix1(1,i) = {mfccf(ncoeff,speechIn1,FS1)};
end
for j = 1:10
q = ['F: \SpeechData\p2\' num2str(j-1) '.wav'];
[speechIn2,FS2] = audioread(q);
speechIn2 = my_vad(speechIn2);
fMatrix2(1,j) = {mfccf(ncoeff,speechIn2,FS2)};
end
for k = 1:10
q = ['F: \SpeechData\p2\' num2str(k-1) '.wav'];
[speechIn3,FS3] = audioread(q);
speechIn3 = my_vad(speechIn3);
fMatrix3(1,k) = {mfccf(ncoeff,speechIn3,FS3)};
end
%將數據保存為mat文件
fields = {'zero','One','Two','Three','Four','Five','Six','Seven','Eight','nine'};
s1 = cell2struct(fMatrix1, fields, 2); %fields項作為行
save Vectors1.mat -struct s1;
s2 = cell2struct(fMatrix2, fields, 2);
save Vectors2.mat -struct s2;
s3 = cell2struct(fMatrix3, fields, 2);
save Vectors3.mat -struct s3;
端點檢測my_vad（）函數----程序清單：
function trimmed_X = my_vad(x) %端點檢測；輸入為錄入語音，輸出為有用信號 Ini = 0.1; %初始靜默時間 Ts = 0.01; %窗的時長 Tsh = 0.005; %幀移時長 Fs = 16000; %采樣頻率 counter1 = 0; %以下四個參數用來尋找起始點和結束點 counter2 = 0; counter3 = 0; counter4 = 0; ZCRCountf = 0; %用於存儲過零率檢測結果 ZCRCountb = 0; ZTh = 40; %過零閾值 w_sam = fix(Ts*Fs); %窗口長度 o_sam = fix(Tsh*Fs); %幀移長度 lengthX = length(x); segs = fix((lengthX-w_sam)/o_sam)+1; %分幀數 sil = fix((Ini-Ts)/Tsh)+1; %靜默時間幀數 win = hamming(w_sam); Limit = o_sam*(segs-1)+1; %最后一幀的起始位置 FrmIndex = 1:o_sam:Limit; %每一幀的起始位置 ZCR_Vector = zeros(1,segs); %記錄每一幀的過零點數 %短時過零點 for t = 1:segs ZCRCounter = 0; nextIndex = (t-1)*o_sam+1; for r = nextIndex+1:(nextIndex+w_sam-1) if (x(r) >= 0) && (x(r-1) >= 0) elseif (x(r) > 0) && (x(r-1) < 0) ZCRCounter = ZCRCounter + 1; elseif (x(r) < 0) && (x(r-1) < 0) elseif (x(r) < 0) && (x(r-1) > 0) ZCRCounter = ZCRCounter + 1; end end ZCR_Vector(t) = ZCRCounter; end %短時平均幅度 Erg_Vector = zeros(1,segs); for u = 1:segs nextIndex = (u-1)*o_sam+1; Energy = x(nextIndex:nextIndex+w_sam-1).*win; Erg_Vector(u) = sum(abs(Energy)); end IMN = mean(Erg_Vector(1:sil)); %靜默能量均值（噪聲均值） IMX = max(Erg_Vector); %短時平均幅度的最大值 I1 = 0.03 * (IMX-IMN) + IMN; %I1，I2為初始能量閾值 I2 = 4 * IMN; ITL = 100*min(I1,I2); %能量閾值下限,前面系數根據實際情況更改得到合適結果 ITU = 10* ITL; %能量閾值上限 IZC = mean(ZCR_Vector(1:sil)); stdev = std(ZCR_Vector(1:sil)); %靜默階段過零率標准差 IZCT = min(ZTh,IZC+2*stdev); %過零率閾值 indexi = zeros(1,lengthX); indexj = indexi; indexk = indexi; indexl = indexi; %搜尋超過能量閾值上限的部分 for i = 1:length(Erg_Vector) if (Erg_Vector(i) > ITU) counter1 = counter1 + 1; indexi(counter1) = i; end end ITUs = indexi(1); %第一個能量超過閾值上限的幀 %搜尋能量超過能量下限的部分 for j = ITUs:-1:1 if (Erg_Vector(j) < ITL) counter2 = counter2 + 1; indexj(counter2) = j; end end start = indexj(1)+1; %第一級判決起始幀 Erg_Vectorf = fliplr(Erg_Vector);%將能量矩陣關於中心左右對稱，若是一行向量相當於逆序 %重復上面過程相當於找結束幀 for k = 1:length(Erg_Vectorf) if (Erg_Vectorf(k) > ITU) counter3 = counter3 + 1; indexk(counter3) = k; end end ITUf = indexk(1); for l = ITUf:-1:1 if (Erg_Vectorf(l) < ITL) counter4 = counter4 + 1; indexl(counter4) = l; end end finish = length(Erg_Vector)-indexl(1)+1;%第一級判決結束幀 %從第一級判決起始幀開始進行第二判決（過零率）端點檢測 BackSearch = min(start,25); for m = start:-1:start-BackSearch+1 rate = ZCR_Vector(m); if rate > IZCT ZCRCountb = ZCRCountb + 1; realstart = m; end end if ZCRCountb > 3 start = realstart; end FwdSearch = min(length(Erg_Vector)-finish,25); for n = finish+1:finish+FwdSearch rate = ZCR_Vector(n); if rate > IZCT ZCRCountf = ZCRCountf + 1; realfinish = n; end end if ZCRCountf > 3 finish = realfinish; end x_start = FrmIndex(start); %最終的起始位置 x_finish = FrmIndex(finish-1); %最終的結束位置 trimmed_X = x(x_start:x_finish);
特征提取mfccf（）函數----程序清單
function FMatrix=mfccf(num,s,Fs) %計算並返回信號s的mfcc參數及其一階和二階差分參數 %參數說明： num 是mfcc系數階數，Fs為采樣頻率 N=512; % FFT 數 Tf=0.02; %窗口的時長 n=Fs*Tf; %每個窗口的長度 M=24; %M為濾波器組數 l=length(s); Ts=0.01; %幀移時長 FrameStep=Fs*Ts; %幀移 a=1; b=[1, -0.97]; %預加重處理的一階FIR高通濾波器系數 noFrames=fix((l-n)/FrameStep)+1; %幀數 FMatrix=zeros(noFrames-2, num); %倒譜系數初始矩陣，一階差分的首尾2幀為0, lifter=1:num; %倒譜加權初始矩陣 lifter=1+floor(num/2)*(sin(lifter*pi/num));%倒譜加權函數 if mean(abs(s)) > 0.01 s=s/max(s); %初始化 end %計算 MFCC 系數 for i=1:noFrames frame=s((i-1)*FrameStep+1:(i-1)*FrameStep+n); %frame用於存儲每一幀數據 framef=filter(b,a,frame); %預加重濾波器 F=framef.*hamming(n); %加漢明窗 FFTo=fft(F,n); %計算FFT melf=melbankm(M,N,Fs); %計算mel濾波器組系數 halfn=1+floor(N/2); spectr=log(melf*(abs(FFTo(1:halfn)).^2)+1e-22);%三角濾波器進行濾波加上1e-22防止進行對數運算后變成無窮 %進行DCT變換 c=zeros(1,num); for p=1:num for m=1:M c(p)=c(p)+spectr(m)*cos(p*(m-0.5)*pi/M); end end ncoeffs=c.*lifter; %對MFCC參數進行倒譜加權 FMatrix(i, :)=ncoeffs; end %調用deltacoeff函數計算MFCC差分系數 d=deltacoeff(FMatrix); %計算一階差分系數 d1=deltacoeff(d); %計算二階差分系數 FMatrix=[FMatrix,d,d1]; %將三組數據作為特征向量
melbankm（）函數---程序清單:
function H=melbankm(M,N,fs,fl,fh) %Mel濾波器組 % 輸入參數: M 濾波器組數量 % N FFT長度 % fs 采樣頻率 % fl -fh 為線性功率譜的有用頻帶(默認為0-0.5*fs) % %輸出參數: H返回濾波器組，每一行為一個濾波器 if nargin<4 fl=0*fs; fh=0.5*fs; end %計算每個濾波器的中心頻率 f=zeros(1,M+1); for m=1:M+2 f(m)=floor((N/fs)*mel2freq(freq2mel(fl)... +(m-1)*(freq2mel(fh)-freq2mel(fl))/(M+1))); end %求濾波器組H c=floor(N/2)+1; y=zeros(1,c); H=zeros(M,c); for m=2:M+1 for k=1:c %由於fh最高為fs/2，那么最多需c位就能存儲H if f(m-1)<=k&&k<=f(m) y(k)=(k-f(m-1))/(f(m)-f(m-1)); elseif f(m)<=k&&k<=f(m+1) y(k)=(f(m+1)-k)/(f(m+1)-f(m)); else y(k)=0; end end H(m,:)=y; end
freq2mel（）函數----程序清單：
function b=freq2mel(f) %%實現線性頻域到mel頻域的轉換 %輸入參數：f為線性頻域頻率 %輸出參數：b為mel頻域頻率 b=2595*log10(1+f/700);
mel2freq（）函數----程序清單：
function f=mel2freq(b) %%實現mel頻域轉換到線性頻域 %輸入參數：b為mel頻域內頻率 %輸出參數：f為線性頻域內頻率 f=700*(10^(b/2595)-1); deltacoeff（）函數----程序清單： function diff = deltacoeff(x) %計算MFCC差分系數 [nr,nc]=size(x); N=2; diff=zeros(nr,nc); for t=3:nr-2 for n=1:N diff(t,:)=diff(t,:)+n*(x(t+n,:)-x(t-n,:)); end diff(t,:)=diff(t,:)/10; %10=2*(1^2+2^2) end
CMN（）函數----程序清單
function NormMatrix = CMN(Matrix) %歸一化處理 [r,c]=size(Matrix); NormMatrix=zeros(r,c); for i=1:c MatMean=mean(Matrix(:,i)); NormMatrix(:,i)=Matrix(:,i)-MatMean; end
DTWScores（）函數----程序清單
function AllScores = DTWScores(rMatrix,N)
%動態時間規整（DTW）尋找最小失真
%輸入參數：rMatrix為當前讀入語音的MFCC參數矩陣,N為每個模板數量詞匯數
%輸出參數：AllScores
%初始化DTW判別矩陣
Scores1 = zeros(1,N);
Scores2 = zeros(1,N);
Scores3 = zeros(1,N);
%加載模板數據
s1 = load('Vectors1.mat');
fMatrixall1 = struct2cell(s1);
s2 = load('Vectors2.mat');
fMatrixall2 = struct2cell(s2);
s3 = load('Vectors3.mat');
fMatrixall3 = struct2cell(s3);
%計算DTW
for i = 1:N
fMatrix1 = fMatrixall1{i,1};
fMatrix1 = CMN(fMatrix1);
Scores1(i) = myDTW(fMatrix1,rMatrix);
end
for j = 1:N
fMatrix2 = fMatrixall2{j,1};
fMatrix2 = CMN(fMatrix2);
Scores2(j) = myDTW(fMatrix2,rMatrix);
end
for k= 1:N
fMatrix3 = fMatrixall3{k,1};
fMatrix3 = CMN(fMatrix3);
Scores3(k) = myDTW(fMatrix3,rMatrix);
end
AllScores = [Scores1;Scores2;Scores3];