《語(yǔ)音信號(hào)處理》實(shí)驗(yàn)5-DTW算法實(shí)現(xiàn)及語(yǔ)音模板匹配

1、華南理工大學(xué)語(yǔ)音信號(hào)處理實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)名稱：DTW算法實(shí)現(xiàn)及語(yǔ)音模板匹配姓名：學(xué)號(hào)：班級(jí)：10級(jí)電信5班日期：2013年6月17日一、實(shí)驗(yàn)?zāi)康倪\(yùn)用課堂上所學(xué)知識(shí)以及matlab工具，利用DTW(Dynamic Time Warping，動(dòng)態(tài)時(shí)間規(guī)整)算法，進(jìn)行說(shuō)話者的語(yǔ)音識(shí)別。二、實(shí)驗(yàn)原理1、語(yǔ)音識(shí)別系統(tǒng)概述一個(gè)完整特定人語(yǔ)音識(shí)別系統(tǒng)的方案框圖如圖1所示。輸入的模擬語(yǔ)音信號(hào)首先要進(jìn)行預(yù)處理，包括預(yù)濾波、采樣和量化、加窗、端點(diǎn)檢測(cè)、預(yù)加重等,然后是參數(shù)特征量的提取。提取的特征參數(shù)滿足如下要求：(1)特征參數(shù)能有效地代表語(yǔ)音特征,具有很好的區(qū)分性；(2)參數(shù)間有良好的獨(dú)立性；(3)特征參數(shù)要計(jì)算方便

2、,要考慮到語(yǔ)音識(shí)別的實(shí)時(shí)實(shí)現(xiàn)。圖1語(yǔ)音識(shí)別系統(tǒng)方案框圖語(yǔ)音識(shí)別的過程可以被看作模式匹配的過程，模式匹配是指根據(jù)一定的準(zhǔn)則，使未知模式與模型庫(kù)中的某一個(gè)模型獲得最佳匹配的過程。模式匹配中需要用到的參考模板通過模板訓(xùn)練獲得。在訓(xùn)練階段，將特征參數(shù)進(jìn)行一定的處理后，為每個(gè)詞條建立一個(gè)模型，保存為模板庫(kù)。在識(shí)別階段，語(yǔ)音信號(hào)經(jīng)過相同的通道得到語(yǔ)音特征參數(shù)，生成測(cè)試模板，與參考模板進(jìn)行匹配，將匹配分?jǐn)?shù)最高的參考模板作為識(shí)別結(jié)果。2、語(yǔ)音信號(hào)的處理1、語(yǔ)音識(shí)別的DTW算法本設(shè)計(jì)中，采用DTW算法，該算法基于動(dòng)態(tài)規(guī)劃(DP)的思想解決了發(fā)音長(zhǎng)短不一的模板匹配問題，在訓(xùn)練和建立模板以及識(shí)別階段，都先采用端點(diǎn)檢

3、測(cè)算法確定語(yǔ)音的起點(diǎn)和終點(diǎn)。在本設(shè)計(jì)當(dāng)中，我們建立的參考模板，m為訓(xùn)練語(yǔ)音幀的時(shí)序標(biāo)號(hào)，M為該模板所包含的語(yǔ)音幀總數(shù)，R(m)為第m幀的語(yǔ)音特征矢量。所要識(shí)別的輸入詞條語(yǔ)音稱為測(cè)試模板，n為測(cè)試語(yǔ)音幀的時(shí)序標(biāo)號(hào)，N為該模板所包含的語(yǔ)音幀總數(shù)，T(n)為第n幀的語(yǔ)音特征矢量。參考模板和測(cè)試模板一般都采用相同類型的特征矢量(如LPCC系數(shù))、相同的幀長(zhǎng)、相同的窗函數(shù)和相同的幀移?？紤]到語(yǔ)音中各段在不同的情況下持續(xù)時(shí)間會(huì)產(chǎn)生或長(zhǎng)或短的變化，因而更多地是采用動(dòng)態(tài)規(guī)劃DP的方法。把測(cè)試模板的各個(gè)幀號(hào)n=1N在一個(gè)二維直角坐標(biāo)系中的橫軸上標(biāo)出，把參考模板的各幀號(hào)m=1M在縱軸上標(biāo)出，通過這些形成網(wǎng)格，網(wǎng)格

4、的每一個(gè)交叉點(diǎn)(n,m)即表示測(cè)試模式中某一幀與訓(xùn)練模式中某一幀的交匯點(diǎn)。DP算法即可以歸結(jié)為尋找一條通過此網(wǎng)格中若干個(gè)點(diǎn)的路徑。路徑通過的格點(diǎn)即為此時(shí)與參考模板中進(jìn)行距離計(jì)算的幀號(hào)。應(yīng)當(dāng)注意，路徑不是隨意選擇的，選取的路徑必定是從左下角出發(fā)，在右上角結(jié)束。通常，規(guī)整函數(shù)被限制在一個(gè)平行四邊形的網(wǎng)格內(nèi)，如圖2所示。它的一條邊斜率為2，另一條邊斜率為1/2。規(guī)整函數(shù)的起點(diǎn)是（1, 1），終點(diǎn)為（N，M）。DTW算法的目的是在此平行四邊形內(nèi)由起點(diǎn)到終點(diǎn)尋找一個(gè)規(guī)整函數(shù)，使其具有最小的代價(jià)函數(shù)，保證了測(cè)試模板與參考模板之間具有最大的聲學(xué)相似特性。圖2匹配路徑約束示意圖由于在模板匹配過程中限定了彎折的

5、斜率，因此平行四邊形之外的格點(diǎn)對(duì)應(yīng)的幀匹配距離是不需要計(jì)算的。另外，因?yàn)槊恳涣懈鞲顸c(diǎn)上的匹配計(jì)算只用到了前一列的3個(gè)網(wǎng)格，所以沒有必要保存所有的幀匹配距離矩陣和累積距離矩陣。充分利用這兩個(gè)特點(diǎn)可以減少計(jì)算量和存儲(chǔ)空間的需求，形成一種高效的DTW算法。圖2中，把實(shí)際的動(dòng)態(tài)彎折分為三段，（1，xa），（xa+1，xb），（xb+1，N），其中：xa= (2M-N)/3， xb=2(2N-M)/3xa和xb都取最相近的整數(shù)，由此可得出對(duì)M和N長(zhǎng)度的限制條件：2M-N3， 2N-M2當(dāng)不滿足以上條件時(shí)，認(rèn)為兩者差別太大，則無(wú)法進(jìn)行動(dòng)態(tài)彎折匹配。在x軸上的每一幀不再需要與y軸上的每一幀進(jìn)行比較，而只是與

6、y軸上ymin，ymax間的幀進(jìn)行比較，ymin和ymax的計(jì)算公式為：ymin=x/2，0xxb，2x+(M-2N)，xb< xNymax=2x，0xxa，x/2+(M-N/2)，xa< xN如果出現(xiàn)xa> xb的情況，則彎折匹配的三段為（1，xb），（xb+1，xa），（xa+1，N）。對(duì)于x軸上每前進(jìn)一幀，雖然所要比較的y軸上的幀數(shù)不同，但彎折特性是一樣的，累積距離的更新都是用下式實(shí)現(xiàn)的：D(x，y) = d(x，y)+minD(x-1，y)，D(x-1，y-1)，D(x-1，y-2)號(hào)的短時(shí)能量或短時(shí)平均幅度就能夠把語(yǔ)音段和噪聲背景區(qū)分開。這是僅基于短時(shí)能量的端點(diǎn)檢測(cè)

7、方法。信號(hào)x(n)的短時(shí)能量定義為: 語(yǔ)音信號(hào)的短時(shí)平均幅度定義為: 其中w(n)為窗函數(shù)。2、短時(shí)平均過零率短時(shí)過零表示一幀語(yǔ)音信號(hào)波形穿過橫軸(零電平)的次數(shù)。過零分析是語(yǔ)音時(shí)域分析中最簡(jiǎn)單的一種。對(duì)于連續(xù)語(yǔ)音信號(hào)，過零意味著時(shí)域波形通過時(shí)間軸；而對(duì)于離散信號(hào)，如果相鄰的取樣值的改變符號(hào)稱為過零。過零率就是樣本改變符號(hào)次數(shù)。信號(hào)x(n)的短時(shí)平均過零率定義為: 式中，sgn為符號(hào)函數(shù)，即: 過零率有兩類重要的應(yīng)用:第一，用于粗略地描述信號(hào)的頻譜特性;第二，用于判別清音和濁音、有話和無(wú)話。從上面提到的定義出發(fā)計(jì)算過零率容易受低頻干擾，特別是50Hz交流干擾的影響。解決這個(gè)問題的辦法，一個(gè)是做

8、高通濾波器或帶通濾波，減小隨機(jī)噪聲的影響；另一個(gè)有效方法是對(duì)上述定義做一點(diǎn)修改，設(shè)一個(gè)門限T，將過零率的含義修改為跨過正負(fù)門限。 于是，有定義: 3、檢測(cè)方法利用過零率檢測(cè)清音，用短時(shí)能量檢測(cè)濁音，兩者配合。首先為短時(shí)能量和過零率分別確定兩個(gè)門限，一個(gè)是較低的門限數(shù)值較小，對(duì)信號(hào)的變化比較敏感，很容易超過；另一個(gè)是比較高的門限，數(shù)值較大。低門限被超過未必是語(yǔ)音的開始，有可能是很短的噪聲引起的，高門限被超過并且接下來(lái)的自定義時(shí)間段內(nèi)的語(yǔ)音超過低門限，意味著信號(hào)開始。此時(shí)整個(gè)端點(diǎn)檢測(cè)可分為四段：靜音段、過渡段、語(yǔ)音段、結(jié)束。實(shí)驗(yàn)時(shí)使用一個(gè)變量表示當(dāng)前狀態(tài)。靜音段，如果能量或過零率超過低門限，就開始

9、標(biāo)記起始點(diǎn)，進(jìn)入過渡段。過渡段當(dāng)兩個(gè)參數(shù)值都回落到低門限以下，就將當(dāng)前狀態(tài)恢復(fù)到靜音狀態(tài)。而如果過渡段中兩個(gè)參數(shù)中的任一個(gè)超過高門限，即被認(rèn)為進(jìn)入語(yǔ)音段。處于語(yǔ)音段時(shí)，如果兩參數(shù)降低到門限以下，而且總的計(jì)時(shí)長(zhǎng)度小于最短時(shí)間門限，則認(rèn)為是一段噪音，繼續(xù)掃描以后的語(yǔ)音數(shù)據(jù)，否則標(biāo)一記結(jié)束端點(diǎn)。三、實(shí)驗(yàn)數(shù)據(jù)及平臺(tái)本實(shí)驗(yàn)所采用的數(shù)據(jù)是錄制的語(yǔ)音文件，平臺(tái)是MATLAB。四、實(shí)驗(yàn)過程（步驟）1、語(yǔ)音信號(hào)預(yù)處理語(yǔ)音信號(hào)的預(yù)處理包括預(yù)濾波、采樣和量化、加窗、預(yù)加重、端點(diǎn)檢測(cè)等過程。由于語(yǔ)音信號(hào)在幀長(zhǎng)為10ms30ms之內(nèi)是相對(duì)平穩(wěn)的，同時(shí)為了便于計(jì)算FFT，本系統(tǒng)選取幀長(zhǎng)N為256個(gè)語(yǔ)音點(diǎn)，幀移M為128點(diǎn)

10、。本文采用漢明窗對(duì)語(yǔ)音信號(hào)進(jìn)行分幀處理，如下式：(n) =0.54-0.46cos(2n/(N-1)，0nN-1預(yù)加重用具有6dB/倍頻程的提升高頻特性的一階數(shù)字濾波器實(shí)現(xiàn)：H(z) =1-0.937 5/z端點(diǎn)檢測(cè)采用基于短時(shí)能量和短時(shí)平均過零率法，利用已知為“靜態(tài)”的最初十幀信號(hào)為短時(shí)能量設(shè)置2個(gè)門限ampl和amph，以及過零率閾值z(mì)cr。語(yǔ)音起始點(diǎn)從第11幀開始檢測(cè)，其流程圖如圖4。語(yǔ)音結(jié)束點(diǎn)的檢測(cè)方法與檢測(cè)起點(diǎn)相似，但此時(shí)從后向前搜索。圖3 語(yǔ)音起點(diǎn)檢測(cè)流程圖2、特征參數(shù)提取及語(yǔ)音識(shí)別眾多研究表明，倒譜特征參數(shù)所含的信息量比其他參數(shù)多，能較好地表現(xiàn)語(yǔ)音信號(hào)。本文選取能夠反映人對(duì)語(yǔ)音的

11、感知特性的Mel頻率倒譜系數(shù)(MFCC)作為特征參數(shù)，階數(shù)為12。經(jīng)過MFCC特征參數(shù)提取后，各幀語(yǔ)音信號(hào)就形成了一個(gè)個(gè)特征矢量。識(shí)別時(shí)，將待測(cè)語(yǔ)音與模板庫(kù)中的每一個(gè)模板進(jìn)行模式匹配，找到距離最小的模板作為輸出結(jié)果。五、實(shí)驗(yàn)結(jié)果及討論如圖4為語(yǔ)音信號(hào)“8”的處理結(jié)果，其他語(yǔ)音信號(hào)處理結(jié)果圖不在此一一給出：圖4 語(yǔ)音信號(hào)“8”的處理結(jié)果圖以下為得到的最終處理結(jié)果，10個(gè)數(shù)字識(shí)別正確。經(jīng)測(cè)試，程序等到了較好的語(yǔ)音識(shí)別效果。六、實(shí)驗(yàn)總結(jié)這次實(shí)驗(yàn)利用MATLAB對(duì)信號(hào)進(jìn)行分析和處理，利用DTW算法，對(duì)說(shuō)話人進(jìn)行語(yǔ)音識(shí)別，主要是進(jìn)行1到9這些數(shù)字的識(shí)別，較好地完成了預(yù)期的目標(biāo)。由于這個(gè)實(shí)驗(yàn)難度比較大，加

12、上自身知識(shí)儲(chǔ)備的不足，在實(shí)驗(yàn)過程中還是遇到比較大的困難，代碼是參考網(wǎng)上的，在理解的基礎(chǔ)上進(jìn)行借鑒。重要的是在整個(gè)過程中，自己投入其中，對(duì)這個(gè)過程進(jìn)行了熟悉與體驗(yàn)，收獲是比較大的。在前人的基礎(chǔ)上，自己還是有很多東西可以做的，例如進(jìn)行算法的改進(jìn)，這方面將會(huì)繼續(xù)努力做好。七、 實(shí)驗(yàn)代碼主程序：yuyinshibie.mdisp('正在計(jì)算參考模板的參數(shù).')for i=1:10fname=sprintf('%da.wav',i-1);x=wavread(fname); x1 x2=vad(x); m=mfcc(x); m=m(x1-2:x2-4,:); ref(i).

13、mfcc=m;enddisp('正在分析語(yǔ)音信號(hào).')for i=1:10fname=sprintf('%da.wav',i-1);x,fs,bit=wavread(fname,2000,2512); %采樣%sound(x,fs); %播放語(yǔ)音信號(hào)figure(i);subplot(3,3,1);plot(x(1:256); %原始語(yǔ)音信號(hào)的時(shí)域圖形%title('原始信號(hào)')subplot(3,3,2)h,w=freqz(x,fs) %原始語(yǔ)音信號(hào)的頻率響應(yīng)圖hr=abs(h); %求系統(tǒng)幅頻響應(yīng)plot(w,hr);title('幅

14、頻圖');xlabel('Frequency in rad mple')ylabel('Magnitude in dB')subplot(3,3,3)hphase=angle(h);hphase=unwrap(hphase); %求系統(tǒng)相頻響應(yīng)plot(w,hphase);title('相頻圖');xlabel('Frequency in rad mple')ylabel('Phase in degrees')y=fft(x,512); %傅立葉變換%mag=abs(y);mag1=10*log10(mag

15、);f=fs*(0:255)/512;subplot(3,3,4)plot(f,mag(1:256); %FFT頻譜圖 %title('fft變換后信號(hào)')iff=ifft(y,512); %反傅立葉變換%ifm=abs(iff);subplot(3,3,5)plot(f,ifm(1:256)title('ifft后信號(hào)')% 短時(shí)傅里葉變換Ts=1/fs;%N=T/Ts;N=512;Nw=20; %窗函數(shù)長(zhǎng)L=Nw/2; %窗函數(shù)每次移動(dòng)的樣點(diǎn)數(shù) Tn=(N-Nw)/L+1; %計(jì)算把數(shù)據(jù)x共分成多少段nfft=32; %FFT的長(zhǎng)度TF=zeros(Tn,n

16、fft); %將存放三維譜圖，先清零for i=1:Tn xw=x(i-1)*10+1:i*10+10); %取一段數(shù)據(jù) temp=fft(xw,nfft); %FFT變換 temp=fftshift(temp); %頻譜以0頻為中心 for j=1:nfft; TF(i,j)=temp(j); %把譜圖存放在TF中 endendsubplot(3,3,6)fnew=(1:nfft)-nfft/2)*fs/nfft;tnew=(1:Tn)*L*Ts;F,T=meshgrid(fnew,tnew);mesh(F,T,abs(TF)title('短時(shí)傅立葉變換時(shí)頻圖')subplo

17、t(3,3,7)contour(F,T,abs(TF)title('等高線表示')enddisp('正在計(jì)算測(cè)試模板的參數(shù).')for i=1:10fname=sprintf('%db.wav',i-1);x=wavread(fname); x1 x2=vad(x); m=mfcc(x); m=m(x1-2:x2-4,:); test(i).mfcc=m;enddisp('正在進(jìn)行模板匹配.')dist=zeros(10,10);for i=1:10for j=1:10 dist(i,j)=dtw(test(i).mfcc,ref

18、(j).mfcc);endenddisp('正在計(jì)算匹配結(jié)果.')for i=1:10 d,j=min(dist(i,:); fprintf('測(cè)試信號(hào)%d的識(shí)別結(jié)果為:%dn',i-1,j-1);end各子程序模塊：dtw.mfunction dist=dtw(t,r)n=size(t,1);m=size(r,1);%幀匹配距離矩陣d=zeros(n,m);for i=1:nfor j=1:m d(i,j)=sum(t(i,:)-r(j,:).2);endend%累積距離矩陣D=ones(n,m)*realmax;D(1,1)=d(1,1);%動(dòng)態(tài)規(guī)劃for

19、i=2:nfor j=1:m D1=D(i-1,j); if j>1 D2=D(i-1,j-1); else D2=realmax; end if j>2 D3=D(i-1,j-2); else D3=realmax; end D(i,j)=d(i,j)+min(D1,D2,D3);endenddist=D(n,m);enframe.mfunction f=enframe(x,win,inc)nx=length(x(:);nwin=length(win);if (nwin = 1) len = win;else len = nwin;endif (nargin < 3) in

20、c = len;endnf = fix(nx-len+inc)/inc);f=zeros(nf,len);indf= inc*(0:(nf-1).'inds = (1:len);f(:) = x(indf(:,ones(1,len)+inds(ones(nf,1),:);if (nwin > 1) w = win(:)' f = f .* w(ones(nf,1),:);endmelbankm.mfunction x,mn,mx=melbankm(p,n,fs,fl,fh,w)if nargin < 6 w='tz' if nargin < 5

21、 fh=0.5; if nargin < 4 fl=0;endendend f0=700/fs; fn2=floor(n/2); lr=log(f0+fh)/(f0+fl)/(p+1); % convert to fft bin numbers with 0 for DC term bl=n*(f0+fl)*exp(0 1 p p+1*lr)-f0); b2=ceil(bl(2); b3=floor(bl(3); if any(w='y') pf=log(f0+(b2:b3)/n)/(f0+fl)/lr; fp=floor(pf); r=ones(1,b2) fp fp+

22、1 p*ones(1,fn2-b3); c=1:b3+1 b2+1:fn2+1; v=2*0.5 ones(1,b2-1) 1-pf+fp pf-fp ones(1,fn2-b3-1) 0.5; mn=1; mx=fn2+1;else b1=floor(bl(1)+1; b4=min(fn2,ceil(bl(4)-1; pf=log(f0+(b1:b4)/n)/(f0+fl)/lr; fp=floor(pf); pm=pf-fp; k2=b2-b1+1; k3=b3-b1+1; k4=b4-b1+1; r=fp(k2:k4) 1+fp(1:k3); c=k2:k4 1:k3; v=2*1-pm

23、(k2:k4) pm(1:k3); mn=b1+1; mx=b4+1;end if any(w='n') v=1-cos(v*pi/2);elseif any(w='m') v=1-0.92/1.08*cos(v*pi/2);end if nargout > 1 x=sparse(r,c,v); else x=sparse(r,c+mn-1,v,p,1+fn2);endmfcc.mfunction ccc=mfcc(x)%歸一化mel濾波器組系數(shù)bank=melbankm(24,256,8000,0,0.5,'m');bank=full(b

24、ank);bank=bank/max(bank(:);%DTC系數(shù),12*24for k=1:12 n=0:23; dctcoef(k,:)=cos(2*n+1)*k*pi/(2*24);end%歸一化倒譜提升窗口w=1+6*sin(pi*1:12./12);w=w/max(w);%預(yù)加重濾波器xx=double(x);xx=filter(1 -0.9375,1,xx);%語(yǔ)音信號(hào)分幀xx=enframe(xx,256,80);%計(jì)算每幀的MFCC參數(shù)for i=1:size(xx,1) y=xx(i,:); s=y'.*hamming(256); t=abs(fft(s); t=t.

25、2; c1=dctcoef*log(bank*t(1:129); c2=c1.*w' m(i,:)=c2'end%差分參數(shù)dtm=zeros(size(m);for i=3:size(m,1)-2 dtm(i,:)=-2*m(i-2,:)-m(i-1,:)+m(i+1,:)+2*m(i+2,:);enddtm=dtm/3;%合并mfcc參數(shù)和一階差分mfcc參數(shù)ccc=m dtm;%去除首尾兩幀，因?yàn)檫@兩幀的一階差分參數(shù)為0ccc=ccc(3:size(m,1)-2,:);vad.mfunction x1,x2=vad(x)%幅度歸一化到-1,1x=double(x);x=x/max(abs(x);%常數(shù)設(shè)置FrameLen=240;FrameInc=80;amp1=10;amp2=2;zcr1=10;zcr2=5;maxsilence=3; %3*10ms=30msminlen=15; %15*10ms=150msstatus=0;count=0;silence=0;%計(jì)算過零率tmp1=enframe(x(1:length(x)-1)