1、南京郵電大學通達學院實驗報告實驗名稱語音時頻參數的提取和應用實驗語音編碼實驗課程名稱信息處理技術專業綜合實驗班級學號姓名開課時間2016/2017學年，第二學期實驗一語音時頻參數的提取和應用一、實驗目的熟悉根本的語音時頻參數提取方法，觀察并比擬各類參數采用不同提取方法的效果，了解這些時頻參數在語音處理中的應用。二、實驗內容1、編寫并調試語音時域參數提取程序，包括短時功率、短時過零率等函數。2、編寫并調試語音頻域參數提取程序，包括頻譜、語譜圖、基音頻率、共振峰等。3、利用提取出的時頻參數對語音信號進展靜/清/濁音的劃分。三、實驗原理語音信號具有短時平穩性，需對其進展分段幀）處理幀長一般取10至3

2、0ms之間）。根據鼓勵方式的不同，語音信號可以分成濁音和清音兩大類：濁音：鼓勵源可模擬成準周期性的脈沖串。該周期稱為基音周期，其倒數就是基音頻率，基音頻率一般在幾百赫茲以內。清音：鼓勵源可模擬成隨機白噪聲。簡化的語音生成數學模型如下圖：對語音時頻波形進展繪制與觀察，有助于了解語音的清/濁音等特性。對基音周期等聲門參數和共振峰等聲道參數進展提取與分析，是語音處理與應用的重要內容。對于語音信號，確定了各線性預測系數后，根據H(z)可得其頻率響應的估值即LPC譜。5. LPC頻譜與基于DFT求得的語音對數幅度譜相比有何不同？LPC譜優點在于可以很好地表示共振峰構造而不出現額外的蜂起和起伏，根本不含有

3、噪聲，其頻率靈敏度與人耳不相匹配；DET得到的對數幅度譜受基頻諧波的影響，最大值只出現在諧波頻率上，共振峰測定誤差較大。6. LPC合成語音的質量如何？有何改良措施？LPC合成技術的優點是簡單直觀。其合成過程實質上只是一種簡單的解碼和拼接過程。另外，由于波形拼接技術的合成基元是語音的波形數據，保存了語音的全部信息，因而對于單個合成基元來說能夠獲得很高的自然度。但是，由于自然語流中的語音和孤立狀況下的語音有著極大的區別，如果只是簡單地把各個孤立的語音生硬地拼接在一起，其整個語流的質量勢必是不太理想的。而LPC技術從本質上來說只是一種錄音+重放，對于合成整個連續語流LPC合成技術的效果是不理想的。

4、因此，LPC合成技術必須和其他技術相結合，才能明顯改善LPC合成的質量。實驗過程記錄與結果分析1、計算短時自相關函數1.1分別計算正弦信號和白噪聲的短時自相關函數，估計正弦信號的基音周期。其Matlab代碼如下:x=sin(2*pi*0.01*(0:499);r,eta=xcorr(x,l00,unbiased);stem(eta,r);w=randn(500,1);r,eta=xcorr(w,100,unbiased);stcm(cta,r);1.2計算正弦信號疊加白噪聲的短時自相關函數，試估計正弦信號的基音周期。其Matlab代碼如下:x=sin(2*pi*0.01*(0:499);w=r

5、andn(500,l);x1=x+w;r,eta=xcorr(x1,1()0,unbiased);stem(eta,r);1.3分別畫出一幀濁音和一幀清音的語音自相關波形(采樣頻率為10kHz,幀長為25.6ms,每幀有256個樣點)，估計濁音的基音周期。loaddigits;x=digits.three1;m=2756;N=256;n=m-N+1:m;r,eta=xcorr(x(n),250,unbiased);plot(eta,r);m=500;N=256;n=m-N+l:ni;r,eta=xcorr(x(n),250,unbiased);plot(eta,r);2、求取LPC系數2.1加

6、載“digits.threel”語音數據。基于自相關法，求出其中一幀數據采樣頻率為10kHz,幀長為25.6ms,每幀有256個樣點)的14階LPC系數。畫出相應的LPC譜。loaddigits;x=digits.three1;m=2756:N=256;n=m-N+1:m;M=14;|r,cta=xcorr(x(n),250,unbiased);Rx=tocplitz(r(M+1:2*M);rx=r(M+2:2*M+l);a=Rxrx;NFFT=IO24;k=l:NFFT/2;X=fft(x(n).*hann(N),NFFT);Thcta=l./fft(l;-a,NFFT);plot(k,20

7、*logl0(abs(353*Theta(k)X(k)axis(0NFFT/2-infinf)2.2編寫一個用以實現杜賓遞推算法的函數“durbin.m”,其Matlab代碼如下：functiona,xi,kappa=durbin(r,M)kappa=zeros(M,l);a=zeros(M,1);xi=r(l);zeros(MJ);for(j=l:M)kappa(j)=r(j+l)-a(l:j-l)*r(j:-l:2)/xi(j);a(j)=kappa；a(1:j-1)=a(1:j-l)-kappa(j)*a(j-l:-l:1);xi(j+l)=a(l:j-l)-kappa(j)*a(j-l

8、xi(j+l)=xi(j)*(l-kappa(j)A2);end2.3利用函數“durbin.m”,計算2.1中語音數據幀的14階LPC系數，并與2.1中的結果進展比擬：x=digils.three1;m=2756;N=256;n=m-N+l:m;M=14;Ir,eta=xcorr(x(n),250,unbiased);aLD,xi,kappa=durbin(r(M+l:2*M+1),M);a,aLD,norm(a-aLD)a=0.35120.20720.38080.31630.0522-0.1639-0.2903-0.2898-0.12220.19600.34920.30240.1713-0

9、.8258aLD=0.35120.20720.38080.31630.0522-0.1639-0.2903-0.2898-0.12220.19600.34920.30240.1713-0.8258ans=l. (X)47e-142.4編寫一個用以實現反射系數轉換為LPC系數的函數“rf21pc.m”,其Matlab代碼如下：functiona=rf21pc(kappa)M=length(kappa);a=zcros(MJ);for(j=l:M)a(j)=kappa(j);a(1:j-1)=a(1:j-1)-kappa(j)*a(j-11:1);end2.5編寫一個用以實現LPC系數轉換為反射系

10、數的函數“lpc2rf.m”，其Matlab代碼如下：functionkappa=lpc2rf(a)M=lengih(a);kappa=zeros(M,1);for(j=M:-l:l)kappa=a(j);a(l：j-D=(a(1:j-1)+a(j)*a(j-l:-l:l)/(1-kappa(j)A2);end2.6使用函數“rf21pc.m”和Tpc2rf.m”,檢驗反射系數和LPC系數相互轉換結果。norm(kappa-lpc2rf(aLD)norm(aLD-rf2Ipc(kappa)ans=1.4760e-14ans=3、語音信號的逆濾波3.1利用2.3中求出的一幀語音數據的LPC系數，

11、構造逆濾波器，并畫出該幀語音信號的殘差波形。ehat=filter(11,1,x(n);plot(x(n)ehal)3.2編寫一個用以實現語音信號LPC分析的函數“Ipcaut。，其Matlab代碼如下：functionar,xi,e,m=lpcauto(x,M,win,Olap)Nx=lcngth(x);N=length(win);if(N=1)N=win;win=ones(NJ);endF=fix(Nx-Olap)/(N-OIap);ar=zeros(M+l,F);xi=zeros(M+l,F);e=zcros(Nx,l);m=zeros(F,I);n=1:N;nl=1:Olap;n2=N

12、-Olap+1:N;n3=Olap+1:N;win1=win(nl)./(win(n1)+win(n2)+eps);win2=win(n2)./(win(n1)+win(n2)+eps);for(f=l:F)r,eta=xcorr(x(n).*win.M,biased);a,xi(:,f),kappa=durbin(r(M+1:2*M+1),M);ar(:,f)=l;-a;ehat=filtcr(ar(:,f).l,x(n);e(n)=e(n(n1).*win2+ehat(n1).*win1;ehat(n3);%Overlap-add.m(f)=n(N);n=n+(N-Olap);end3.3

13、加載“timitl”語音數據，利用“Ipcauio”函數對該句語音進展LPC分析。畫出語音信號及其殘差波形。loadlimitI;x=timit1;M=14;N=256;ar,xi,e,m=lpcauto(x,M,hann(N),N/2);subplot(121),plot(x)subplot(122),plot(e)axis(0,60000,-1,1)plot(x,e)soundsc(x)soundsc(c)soundsc(x-e)4、LPC譜估計4.1編寫一個用以計算和顯示LPC譜的函數Ipcplot”,其Matlab代碼如下：functionlpcplot(A,Nfft,Fs,m)M,N

14、=size(A);if(N=l)Theta,F=freqz(l,A,Nfft,Fs);plot(F,20*logl0(abs(Thcta);xlabel(Frequency,itFHz);ylabel(Magnitude,|theta(oniega)|dB);elseif(length(m)=N)error(ThecolumndimensionofAmustbeequaltothelengthofm.)endTheta=zeros(Nfft,N);for(n=l:N)Theta(:,n),F=freqz(1,A(:,n),Nfft,Fs);endMeshHndl=meshz(m,F,20*Io

15、g10(abs(Theta);axisij;view(-45,45);set(MeshHndl,MeshStyle,Column);axistight;axisautoy;axisautoz;xlabcl(SamplcNumber,(itn);ylabelCFrequency,itF)Hz);zlabcl(Magnitudc,|thcta(onicga)|dB);end4.2利用“Ipcplot”函數，畫出3.3中語音數據的幾幀采樣頻率為16kHz,幀長為32ms,每幀有512個樣點)LPC譜。plot(x(m(71):m(147)lpcplot(ar(:,71:147),512,16000,

16、m(71:147)lpcplot(ar(:,82),512,16000,0)lpcplot(ar(:,ll0),512,16000,0)lpcplot(ar(:,l34),512,16000,0)對于時變語音信號雙)，可以計算它的短時功率乙(7)和短時過零率ZA.(/72)o第，”幀幀長為N個樣點)的功率計算公式為：過零率計算公式為：對于濁音語音，可以利用其頻譜X(C)具有豐富的諧波分量的特點，求出其諧波乘積譜，計算公式為：式中，R般取為5。在諧波乘積譜中，基頻分量變得很大，更易于估計基音周期。四、實驗方法及程序1.調用mallab中的wavread和wavplay讀入并播放語音文件。2.調用

17、fft命令計算語音的頻譜。3.調用specgram命令得到語音的語譜圖。4.調用plot命令繪制并觀察語音的時頻域波形。5.根據短時功率、短時過零率、諧波乘積譜的計算公式，編寫分別用以計算短時功率、短時過零率、基音頻率的函數：functionPx=stpowerfunctionZx=stzerocross(x.N)和functionHPSx=hpspectrum(x,N,R)。6.利用短時功率和短時過零率這兩個參數可以對語音信號進展濁音/清音分類。編寫實現此分類功能的函數voiunvoi.mo五、實驗結果與分析1. 清音和濁音的短時功率、短時過零率各有何特點？清音段的能量一般比濁音段的小得多。

18、濁音時能量集中于較低頻率段內，具有較低的過零率，而清音時能量集中于較高頻率段內，具有較高的過零率。2. 如何利用短時功率和短時過零率這兩個參數對語音信號進展濁音/清音的分類？發濁音時，由于聲門波引起了譜的高頻跌落。所以其語音能量集中在3kHz以下；發清音時，多數能量出現在高頻。高頻率意味著高的平均過零率，低頻率意味著低的平均過零率，-般可認為濁音具有較低的平均過零率，清音具有較高的平均過零率。3. 如何根據時域波形估計各幀元音的基音周期？看圖可知四幀語音的周期均為0.01s,第一幀()l()()Hz有10個峰值點，基頻約為10Hz,第二幀0l()()Hz有5個峰值點，基頻約為2()Hz,第三幀

19、()l()()Hz有5個峰值點，基頻約為20Hz,第四幀()l()0Hz有4個峰值點，基頻約為25Hzo4. 如何根據基于DFT的對數幅度譜估計某帔濁音的基音周期？01172為一個基因周期。5. 如何根據基于DFT的對數幅度譜，估計出共振峰頻率？有7個峰點，因此基頻約為167。6. 時域對語音信號進展加窗，反映在頻域，其窗譜對基于DFT的對數幅度譜有何影響？如何估計出窗譜的主瓣寬度？因為時間窗幅度的傅立葉變換為脈沖狀態，那么頻譜偏差很小，當定常過程x(n)的平均值6為零且N很大時ETN(o)=f(a)(l2)所以加窗的周期圖TN(3)是頻譜密度函數的無偏估計。固有TN(o)=l+頃如果原始信號

20、的頻譜成份與FFT中的譜線完全一致，這種情況下采樣數據的長度為信號周期的整數倍，頻譜中只有主瓣。7. 如何根據窄帶、寬帶語譜圖，提取出語音的基音和共振峰軌跡？在基音提取中，廣泛采用語音波形或誤差信號波形的低通濾波。提取共振峰特性最簡便的手段是使用語譜儀。提取共振峰還有倒譜法、LPC分析法等更為有效、準確的方法。做完試驗后，記住要答復這些問題)實驗過程記錄與結果分析1、進入matlab1.1在Windows環境中，創立一個名為speech的文件夾，存放語音數據和與實驗相關的Matlab文件。如：E:speecho用于實驗的語音數據(.mat)及相應說明(.txt)包括：digits.matdig

21、its.txt%gliss.mat%兩個包含/i/的滑音gliss.txtletters.matletlers.txt%英文字母表中26個字母的發音mal.txtmal-1.matma1-2.matma1-3.matmal-4.mat%一句連續語音%一句連續語音%一句連續語音%一句連續語音timit.txtlimit1.mattimit2.mattimit3.mattimit4.ma(vowels.matvowels.txtwords,matwords.txt%一句連續語音%一句連續語音%一句連續語音%一句連續語音%元音/a/,/i/,/o/,/u/的發音%十個孤立字的發音1.2在Window

22、s桌面上，雙擊Matlab的圖標，進入Matlab的工作環境。在命令窗口提示符后，鍵入cdE:spcech指令，進入實驗文件夾路徑。1.3在命令窗口鍵入clearall指令，去除工作空間的所有變量。2、加載語音數據，熟悉matlab命令2.1加載timitl.mat語音數據，命令為：loadlimitI2.2用命令“who”和“whos列出當前工作空間中的變量。變量“timitl”會以包含55911個元素的矩陣列向量)形式出現Yourvariablesare:timitlNameSizeBytesClassAttributestimitl55911x1447288double2.3利用命令“l

23、ength”和“size”查看語音信號timitl的長度和維數：m=length(timitl)lm,nj=size(timitl)(粘貼實驗結果)in=55911m=55911n=I2.4通過以下命令，分別查看語音信號timitl中間的1個、2個和20個數據值:limit1(5001)timitl(5001,5003)timitl(5001:5020)(粘貼實驗結果)ans=-0.0646ans=-0.06460.0111ans=-0.06460.03070.0111-0.02000.20210.40120.1869-0.1756-0.19020.08230.23630.17120.0637

24、-0.0209-0.0683-0.03940.00810.01050.03110.06962.5Matlab中，符號/可以實現轉置運算，運行以下命令進展觀察：size(timill)粘貼實驗結果)ans=1559112.6使用函數max和min可以分別得到信號的最大值和最小值：max(timit1)(粘貼實驗結果ans=Ians=-0.80683、語音信號的時域分析3.1使用以下命令，畫出語音信號timitl(采樣頻率為16kHz,長度為3.5秒)的時域波形。plot(timitl)title(4imitr)xlabel(TimeIndex,iln)ylabel(Amplitude)axis(

25、0length(tiniitI)-1l|)grid3.2分別畫出一幀濁音和一幀清音的語音時域波形(采樣頻率為16kHz,幀長為25ms,每幀有400個樣點)。subplot(211),plot(timitl(14501:14900)%/a/indark.subplot(212),plot(timitl(3550l:35900)%/s/inwash.3.3編寫個一用以計算信號短時功率的函數“stpower.m”,其Matlab代碼如下：functionPx=stpower(x,N)M=length(x);Px=zeros(M,1);Px(N)=x(1:N)*x(1:N)/N;For(m=(N+l

26、):M)Px(m)=Px(m-1)4-(x(m)A2-x(m-N)A2)/N;3.4編寫一個用以計算信號短時過零率的函數“stzerocross.ni,其Matlab代碼如下：functionZx=stzerocross(x,N)M=length(x);Zx=zeros(M,l);Zx(N+1)=suni(abs(sign(x(2:N+1)-sign(x(1:N)/(2*N);for(m=(N+2):M)Zx(m)=Zx(m-1)+(abs(sign(x(m)-sign(x(m-l).abs(sign(x(m-N)-sign(x(m-N-1)/(2*N);end3.5加載“digits.mat

27、”語音數據，該數據包含英文數字“0”至U“9”的發音。計算其中的單詞“four”的短時功率和過零率(采樣頻率為10kHz,幀長為30ms,每幀有300個樣點)。loaddigits;N=300;x=digits.fourl;Px=stpower(x.N);Zx=stzerocross(x,N);plot(Px*le-5Zxx/20001)3.6加載vowels.mat”語音數據,該數據包含元音/a/,/i/o/,/u/的發音。分別畫出一幀/a/、一幀/i/、一幀/o/和一幀/u/的時域波形采樣頻率為10kHz,幀長為30ms,每幀有300個樣點)。loadvowelssubplot(221)p

28、lot(vowels.a_1(2001:2300)subplot(222)plot(vowels.i_1(2001:2300)subplot(223)plot(vowels.o_1(2001:2300)subplot(224)plot(vowels.u_1(2001:23OO)3.7編寫一個利用短時功率和短時過零率這兩個參數對語音信號進展濁音/清音分類的函數voiunvoi.m”,其Matlab代碼如下：functionvoi=voiunvoi(x,N,Pth,Zth)%Short-timepowerPxandzerocrossingZxmeasures.Px=stpower(x,N);Zx=

29、stzerocross(x,N);%Compareestimateswiththresholdvalues.voi=(PxPth*max(Px)&(ZxZth);%Shiftthevoi-flagN/2samplestotheleft(middlesampleinwindow).voi=voi(fix(N/2)+1:length(voi);voi(length(voi)*ones(fix(N/2),1);4、語音信號的頻域分析4.1加載“mal_l”語音數據。基于DFT變換，畫出其中一幀數據采樣頻率為8kHz,幀長為37.5ms,每幀有30()個樣點)的頻域波形(對數幅度譜)。loadmal_

30、l;x=mal_l(4161:4460);subplot(121),plot(x)N=1024;k=0:N/2-l;X=fftshift(fft(x.*hann(length(x).N);subplot。22),plot(k,20*logl0(abs(X(N/2:-1:1),axis(0N/2-1-infinf)4.2編寫一個用以計算信號諧波乘積譜的函數“hpspectrum.m,其Matlab代碼如下：functionHPSx=hpspcctrum(x,N,R)k=l:R:N/2;K=length(k);X=fft(x.*hann(length(x),N);HPSx=X(k);for(r=H

31、PSx=HPSx.*X(l:r:r*K);end4.3利用函數“hpspectrum.m”,來估計一幀濁音的基音周期。代碼如下：X=ma1_1(4161:4460);N=1024;R=5;HPSx=hpspectrum(x,N,R);plot(20*logl0(abs(HPSx)4.4加載“vowels.mat”語音數據，分別畫出一*幀/i/和一幀/u/采樣頻率為10kHz,幀長為30ms,每幀有300個樣點)的基于DFT的對數幅度譜。其Matlab代碼如下：loadvowelsx=vowels.i_l(2001:2300);N=1024;k=-N/2:N/2-l;X=fftshift(fft

32、(x.*hann(length(x),N);plot(k,20*log!0(abs(X),axis(0fix(N/2)0100J)x=vowels.u_l(2001:2300);N=1024;k=-N/2:N/2-l;X=fftshift(fft(x.*hann(length(x),N);plot(k,20*logl0(abs(X),axis(0fix(N/2)0100)4.5畫出一幀清音語音的基于DFT的對數幅度譜。語音數據為mal.l中的第15701-15860個樣點采樣頻率為8kHz,幀長為20ms,每幀有160個樣點)。loadma1_1;x=mal-l(4161:4460);plot

33、(x)N=1024;k=-N/2;N/21;X=fftshift(fft(x.*hann(length(x),N);plot(k,20*log!0(abs(X),axis(0fix(N/2)-infinf)4.6加載“timitl”語音數據。利用函數specgram,畫出該句語音的語譜圖。其Matlab代碼如下:loadtimitl;NFFT=256;Fs=16000;Win=256;Noverlap=128;specgram(timitl,NFFT,Fs,Win,Noverlap);4.7加載“gliss語音數據。畫出該句語音的語譜圖。loadgliss;NFFT=256;Fs=10000;

34、Win=256;Noverlap=128;specgram(gliss.i_2,NFFT,Fs,Win,Noverlap);4.8基于濁音/清音分類函數voiunvoi.m和諧波乘積譜函數hpspectrum.m,編寫一個用以計算濁音短時基音周期的函數“slpilch.m”。其Matlab代碼如下：functionFp=stpitch(x,N,Pth,Zth,NFFT,R,Fs)M=length(x);N=2*fix(N/2);win=hann(N);Fp=zeros(M,l);voi=voiunvoi(x,N,Pth,Zth);invoi=0;for(m=N:N/2:M)n1-m-N+l:m

35、-N/2;n2=m-N/2+1:m;n=nln2;if(any(*voi(n)if(invoi)Fp(n1)=Fmax;invoi=0;endelseHPSmax,Fmax=max(abs(hpspectrum(x(n),NFFT,R);if(invoi)Fp(n)=Fp(n)+win*Fmax;elseFp(n)=Fmax*ones(N/2,1)win(N/2+l:N);invoi=1;endendendwsave=warning;warning(off);Fp=(Fp*Fs/NFFT).*(Fp./Fp);warning(wsave);1. 調用matlab中的wavread和wavplay讀入并播放語音文件。2. 調用fft命令計算語音的頻譜。3. 調用speegram命令得到語音的語譜圖。4. 調用plot命令繪制并觀察語音的時頻域波形。5. 根據短時功率、短時過零率、諧波乘枳諾的計算公式，編寫分別用以計算短時功率、短時過零率、基音頻率的函數：functionPx=stpowcr、functionZx=stzerocross(x,N)和functionHPSx=hpspectr