大學(xué)本科語音信號處理實驗講義8學(xué)時

1、語音信號處理實驗講義時間：2011-12 目錄實驗一 語音信號生成模型分析.3實驗二 語音信號時域特征分析.7實驗三 語音信號頻域特征分析12實驗四 語音信號的同態(tài)處理和倒譜分析16實驗一 語音信號生成模型分析一、 實驗?zāi)康?、了解語音信號的生成機理，了解由聲門產(chǎn)生的激勵函數(shù)、由聲道產(chǎn)生的調(diào)制函數(shù)和由嘴唇產(chǎn)生的輻射函數(shù)。2、編程實現(xiàn)聲門激勵波函數(shù)波形及頻譜，與理論值進(jìn)行比較。3、編程實現(xiàn)已知語音信號的語譜圖，區(qū)分濁音信號和清音信號在語譜圖上的差別。二、實驗原理語音生成系統(tǒng)包含三部分：由聲門產(chǎn)生的激勵函數(shù)、由聲道產(chǎn)生的調(diào)制函數(shù)和由嘴唇產(chǎn)生的輻射函數(shù)。語音生成系統(tǒng)的傳遞函數(shù)由這三個函數(shù)級聯(lián)而成，即

2、 1、激勵模型發(fā)濁音時，由于聲門不斷開啟和關(guān)閉，產(chǎn)生間隙的脈沖。經(jīng)儀器測試它類似于斜三角波的脈沖。也就是說，這時的激勵波是一個以基音周期為周期的斜三角脈沖串。單個斜三角波的頻譜表現(xiàn)出一個低通濾波器的特性?？梢园阉硎境蓏變換的全極點形式這里c是一個常數(shù)，T是脈沖持續(xù)時間。周期的三角波脈沖還得跟單位脈沖串的z變換相乘：這就是整個激勵模型，是一個幅值因子。2、聲道模型當(dāng)聲波通過聲道時，受到聲腔共振的影響，在某些頻率附近形成諧振。反映在信號頻譜圖上，在諧振頻率處其譜線包絡(luò)產(chǎn)生峰值，把它稱為共振峰。一個二階諧振器的傳輸函數(shù)可以寫成實踐表明，用前3個共振峰代表一個元音足夠了。對于較復(fù)雜的輔音或鼻音共振峰

3、要到5個以上。多個疊加可以得到聲道的共振峰模型3、輻射模型從聲道模型輸出的是速度波，而語音信號是聲壓波。二者倒比稱為輻射阻抗，它表征了口唇的輻射效應(yīng)，可用下式表示：三、實驗內(nèi)容1、設(shè)聲門脈沖單個三角波的數(shù)學(xué)表達(dá)式為試畫出三角波波形圖及其頻譜。取。參考程序：%三角波及其頻譜n=linspace(0,25,125);g=zeros(1,length(n);i=0;for i=0:40 if n(i+1)<=5 g(i+1)=0.5*(1-cos(n(i+1)*pi/5); else g(i+1)=cos(n(i+1)-5)*pi/8); endendfigure(1)subplot(1,2,

4、1)plot(n,g)xlabel('時間/ms')ylabel('幅度')axis(0,25,-0.4,1.2)r=fft(g,1024);r1=abs(r);yuanlai=20*log10(r1);signal(1:512)=yuanlai(1:512);pinlv=(0:1:511)*8000/1024;subplot(1,2,2)plot(pinlv,signal);xlabel('頻率/Hz')ylabel('幅度/dB')axis(0,620,0,30) 圖1-1 三角波及其頻譜2、給出語音段“數(shù)字信號處理”（spe

5、ech.wav），畫出它的語譜圖。clear all;x,sr=wavread('speech_dsp.wav');s=length(x);w=round(44*sr/1000);n=w;shift=w/2;h=w-shift;%win=hanning(n)'win=hamming(n)'c=1;ncols=1+fix(s-n)/h);d=zeros(1+n/2),ncols);for b=0:h:(s-n) u=win'.*x(b+1):(b+n); t=fft(u); d(:,c)=t(1:(1+n/2); c=c+1;endtt=0:h:(s-n)

6、/sr;ff=0:(n-2)*sr/n;imagesc(tt,ff/1000,20*log10(abs(d);colormap(gray);axis xyxlabel('時間/s')ylabel('頻率/kHz') 圖1-2 語譜圖四、思考題1、聲門激勵脈沖信號是高頻衰減的還是高頻增強的？2、畫語譜圖時為什么要給語音信號加漢明窗？若加矩形窗會有什么區(qū)別？3、在語譜圖上觀察，濁音信號的和清音信號的頻譜有什么區(qū)別？實驗二 語音信號時域特征分析一、實驗?zāi)康?、了解自相關(guān)函數(shù)及自相關(guān)函數(shù)在語音信號處理中的應(yīng)用。2、編寫程序分析語音信號的短時自相關(guān)特征，計算語音信號的基音

7、周期。3、編寫修正短時自相關(guān)函數(shù)的程序，并與未修正的函數(shù)進(jìn)行比較。二、實驗原理自相關(guān)函數(shù)用于衡量信號自身時間波形的相似性。由前面的討論可知，清音和濁音的發(fā)聲機理不同，因而在波形上也存在著較大的差異。濁音的時間波形呈現(xiàn)出一定的周期性，波形之間相似性較好；清音的時間波形呈現(xiàn)出隨機噪聲的特性，雜亂無章，樣點間的相似性較差，這樣，可以用短時自相關(guān)函數(shù)來測定語音的相似特性。時域離散確定信號的自相關(guān)函數(shù)定義為：對于語音信號來說，采用短時分析方法，可以定義短時自相關(guān)函數(shù)為因為，所以定義，則上式可以寫成 如果長基音周期用窄的窗，將得不到預(yù)期的基音周期；但是如果短的基音周期用長的窗，自相關(guān)函數(shù)將對多個基因周期做

8、平均計算，從而模糊語音的短時特性，這是不希望的。為了解決這個問題，可以采用修正的短時自相關(guān)函數(shù)，選擇的窗長不一定要等于自相關(guān)函數(shù)的最大自變量取值。這種方法可以采用較窄的窗，同時避免了短時自相關(guān)函數(shù)隨k增加而衰減的不足。三、實驗內(nèi)容1、根據(jù)給出的濁音信號，分別畫出濁音信號的時域波形、加矩形窗和加漢明窗后計算短時自相關(guān)歸一化后的結(jié)果。語音的抽樣頻率為8kHz，窗長為320。參考程序：x,fs,nbits=wavread('speech_dsp.wav');s1=x(2500:2819);N=320;A=;for k=1:320 sum=0; for m=1:N-k+1 sum=su

9、m+s1(m)*s1(m+k-1); end A(k)=sumendfor k=1:320 A1(k)=A(k)/A(1);endf=zeros(1,320);n=1,j=1;while j<=320 f(1,j)=s1(n)*0.54-0.46*cos(2*pi*n/319); j=j+1; n=n+1;endB=;for k=1:320 sum=0; for m=1:N-k+1 sum=sum+s1(m)*s1(m+k-1); end B(k)=sumendfor k=1:320 B1(k)=B(k)/B(1);end%畫圖s2=s1/max(s1);figure(1)subplot

10、(3,1,1)plot(s2)title('一幀語音信號');xlabel('樣點數(shù)');ylabel('幅度');axis(0,320,-1,1);subplot(3,1,2)plot(A1)title('加矩形窗的自相關(guān)函數(shù)')xlabel('延時k')ylabel('自相關(guān)函數(shù)R(k)')axis(0,320,-1,1);subplot(3,1,3)plot(B1)title('加漢明窗的自相關(guān)函數(shù)')xlabel('延時k')ylabel('自相關(guān)函數(shù)

11、R(k)')axis(0,320,-1,1); 圖2-1 濁音信號加不同窗時的自相關(guān)函數(shù)2、仍選取上題中的語音信號，改變窗長和截取語音段的長度，計算修正的短時自相關(guān)函數(shù)。取值分別為：(1)N=320,M=640;(2)N=160,M=320;(3)N=70,M=140;x,fs,nbits=wavread('speech_dsp.wav');s1=x(2500:3139);b=s1;%窗長640，自相關(guān)運算取320個點。b1=b(1:640);N=320;A=;for k=1:320 sum=0; for m=1:N sum=sum+b1(m)*b1(m+k-1); e

12、nd A(k)=sum;endfor k=1:320 A1(k)=A(k)/A(1);end%畫圖figure(1)subplot(3,1,1)plot(A1);xlabel('延時k')ylabel('R(k)')legend('N=320')axis(0,320,-0.5,1)%窗長320，自相關(guān)運算取160個點。b2=b(1:320);N=160;A=;for k=1:160 sum=0; for m=1:N sum=sum+b2(m)*b2(m+k-1); end B(k)=sum;endfor k=1:160 B1(k)=B(k)/B(

13、1);end%畫圖figure(1)subplot(3,1,2)plot(B1);xlabel('延時k')ylabel('R(k)')legend('N=160')axis(0,320,-0.5,1)%窗長140，自相關(guān)運算取70個點。b3=b(1:140);N=70;A=;for k=1:70 sum=0; for m=1:N sum=sum+b3(m)*b3(m+k-1); end C(k)=sum;endfor k=1:70 C1(k)=C(k)/C(1);end%畫圖figure(1)subplot(3,1,3)plot(C1);xla

14、bel('延時k')ylabel('R(k)')legend('N=70')axis(0,320,-0.5,1) 圖2-2 修正的自相關(guān)函數(shù)（參加自相關(guān)運算的點數(shù)N取不同值）四、思考題1、自相關(guān)函數(shù)的作用是什么？互相關(guān)函數(shù)的作用是什么？2、濁音信號分別加矩形窗和漢明窗時自相關(guān)函數(shù)有什么不同？3、清音信號的自相關(guān)函數(shù)和濁音信號的有什么區(qū)別？實驗三 語音信號頻域特征分析一、實驗?zāi)康?、了解語音信號進(jìn)行短時傅里葉分析的基礎(chǔ)。2、畫出濁音信號加不同窗時的時域波形和頻譜，并進(jìn)行比較。3、畫出清音信號加不同窗時的時域波形和頻譜，并與濁音信號進(jìn)行比較。二、實驗

15、原理語音信號本質(zhì)上是非平穩(wěn)信號，其非平穩(wěn)特性是由發(fā)聲器官的物理運動產(chǎn)生的。發(fā)聲器官的運動由于存在慣性，所以可以假設(shè)語音信號在1030ms這樣短的時間段內(nèi)是平穩(wěn)的，這是短時傅里葉分析的基礎(chǔ)。某一幀的短時傅里葉變換的定義式如下：式中，是窗函數(shù)，不同的窗函數(shù)，可得到不同的傅里葉變換的結(jié)果。在式中，短時傅里葉變換有兩個變量，即離散時間n及連續(xù)頻率w，若令，則可得到離散的短時傅里葉變換如下： 對于窗來說，它除了具有選出序列中被分析部分的作用外，它的形狀對時變傅里葉變換的特性也有重要作用。采用矩形窗時，基音諧波的各個峰都比較尖銳，而整個頻譜圖顯得比較破碎，這是因為矩形窗的主瓣較窄，具有較高的頻率分辨率，但

16、它也具有較高的旁瓣，因而使基音的相鄰諧波之間的干擾比較嚴(yán)重。相鄰諧波之間的這種“泄露”的現(xiàn)象，抵消了矩形窗主瓣窄的優(yōu)點，因此，在語音短時頻譜分析中極少采用矩形窗。當(dāng)加漢明窗時，得到的短時頻譜要平滑得多，因而在語音分析中漢明窗用得比較普遍。三、實驗內(nèi)容1、根據(jù)給出的濁音語音，畫出它分別加矩形窗和漢明窗的時域波形和短時頻譜。%函數(shù)文件fra.m，對語音信號進(jìn)行幀長為frame，重疊為shift的分幀function output=fra(frame,shift,x)len=length(x);frame_num=floor(len-frame)/shift)+1;y=zeros(frame_num

17、,frame);for i=1:frame_num y(i,:)=x(shift*(i-1)+1):(shift*(i-1)+frame)'endoutput=y;%濁音的波形和短時頻譜圖(窗長256)x,fs,nbits=wavread('speech_dsp.wav');e=fra(256,128,x);ee=e(20,:);subplot(2,2,1)ee1=ee/max(ee);plot(ee1)xlabel('樣點數(shù)')ylabel('幅度')title('原始語音')axis(0,256,-1.5,1.5)%矩

18、形窗傅立葉變換r=fft(ee,1024);r1=abs(r);r1=r1/max(r1);yuanlai=20*log10(r1);signal(1:512)=yuanlai(1:512);pinlv=(0:1:511)*8000/1024;subplot(2,2,2)plot(pinlv,signal);xlabel('頻率/Hz')ylabel('對數(shù)幅度/dB')title('加矩形窗時語音譜')axis(0,4000,-80,15)%加漢明窗f=ee'.*hamming(length(ee);f1=f/max(f);subplo

19、t(2,2,3)plot(f1)axis(0,256,-1.5,1.5)xlabel('樣點數(shù)')ylabel('幅度')title('窗選語音')%加漢明窗傅立葉變換r=fft(f,1024);r1=abs(r);r1=r1/max(r1);yuanlai=20*log10(r1);signal(1:512)=yuanlai(1:512);pinlv=(0:1:511)*8000/1024;subplot(2,2,4)plot(pinlv,signal);xlabel('頻率/Hz')ylabel('對數(shù)幅度/dB

20、9;)title('加漢明窗時語音譜')axis(0,4000,-80,15) 圖3-1 濁音信號加不同窗時的頻譜2、根據(jù)給出的清音語音，畫出它分別加矩形窗和漢明窗的時域波形和短時頻譜。%清音的波形和短時頻譜圖(窗長256)x,fs,nbits=wavread('speech_dsp.wav');e=fra(256,128,x);ee=e(10,:);subplot(2,2,1)ee1=ee/max(ee);plot(ee1)xlabel('樣點數(shù)')ylabel('幅度')title('清音原始語音')axis(

21、0,256,-1.5,1.5)%清音加矩形窗傅立葉變換r=fft(ee,1024);r1=abs(r);r1=r1/max(r1);yuanlai=20*log10(r1);signal(1:512)=yuanlai(1:512);pinlv=(0:1:511)*8000/1024;subplot(2,2,2)plot(pinlv,signal);xlabel('頻率/Hz')ylabel('對數(shù)幅度/dB')title('清音加矩形窗時語音譜')axis(0,4000,-80,15)%清音加漢明窗f=ee'.*hamming(lengt

22、h(ee);f1=f/max(f);subplot(2,2,3)plot(f1)axis(0,256,-1.5,1.5)xlabel('樣點數(shù)')ylabel('幅度')title('清音窗選語音')%清音加漢明窗傅立葉變換r=fft(f,1024);r1=abs(r);r1=r1/max(r1);yuanlai=20*log10(r1);signal(1:512)=yuanlai(1:512);pinlv=(0:1:511)*8000/1024;subplot(2,2,4)plot(pinlv,signal);xlabel('頻率/Hz

23、')ylabel('對數(shù)幅度/dB')title('清音加漢明窗時語音譜')axis(0,4000,-80,15)圖3-1 清音信號加不同窗時的頻譜四、思考題1、濁音信號分別加矩形窗和漢明窗時的頻譜有什么區(qū)別？2、清音信號的頻譜和濁音信號相比有什么不同?3、對語音信號進(jìn)行短時傅里葉分析的基礎(chǔ)是什么？實驗四 語音信號的同態(tài)處理和倒譜分析一、實驗?zāi)康?、掌握同態(tài)處理的原理，倒譜與復(fù)倒譜的含義。2、利用同態(tài)處理對語音信號進(jìn)行基音檢測。3、利用同態(tài)處理對語音信號進(jìn)行共振峰檢測。二、實驗原理1、基音檢測語音的倒譜是將語音的短時譜取對數(shù)后再進(jìn)行IDFT得到的，所以濁

24、音信號的周期性激勵反映在倒譜上是同樣周期的沖激。借此，可從倒譜波形中估計出基音周期。一般把倒譜波形中第二個沖激，認(rèn)為是對應(yīng)激勵源的基頻。下面給出一種倒譜法求基音周期的框圖。先計算倒譜，然后在預(yù)期的基音周期附近尋找峰值。如果倒譜的峰值超出了預(yù)先規(guī)定的門限，則輸入語音段定為濁音，而峰的位置就是基音周期的良好估值。如果沒有超出門限的峰值，則輸入語音段定為清音。如果計算的是一個時變的倒譜，則可估計出激勵源模型及基音周期隨時間的變化。一般每隔1020ms計算一次倒譜，這是因為在一般語音中激勵參數(shù)是緩慢變化的。DFTIDFT峰值檢測 基音檢測 s(n)w(n) 圖4.1 一種倒譜法求基音周期的實現(xiàn)框圖2、

25、共振峰檢測倒譜將基音諧波和聲道的頻譜包絡(luò)分離開來。倒譜的低時部分可以分析聲道、聲門和輻射信息，而高頻部分可用來分析激勵源信息。對倒譜進(jìn)行低時窗選，通過語音倒譜分析系統(tǒng)的最后一級，進(jìn)行DFT后的輸出即為平滑后的對數(shù)模函數(shù)，這個平滑得對數(shù)譜顯示了特定輸入語音段的諧振結(jié)構(gòu)，即譜的峰值基本上對應(yīng)于共振峰頻率，對平滑過的對數(shù)譜中的峰值進(jìn)行定位，即可估計共振峰。原理框圖如圖4.2所示。分幀加窗|FFT|logIFFT加窗|FFT|平滑處理共振峰log圖4.2 共振峰檢測框圖三、實驗內(nèi)容1、利用同態(tài)處理計算給出語音信號的基音周期和基頻。y,fs,nbits=wavread('speech_dsp.w

26、av');time1=1:length(y);time=(1:length(y)/fs;frameSize=floor(30*fs/1000);startIndex=round(2080);endIndex=startIndex+frameSize-1;frame=y(startIndex:endIndex);ylen=frameSize;frameSize=length(frame);frame2=frame.*hamming(length(frame);rwy=rceps(frame2);cepstrum=rwy(1:ylen/2);for i=1:ylen/2 cepstrum1

27、(i)=rwy(ylen/2+1-i);endfor i=(ylen/2+1):ylen cepstrum1(i)=rwy(i+1-ylen/2);end%基音檢測LF=floor(fs/500);HF=floor(fs/70);cn=cepstrum(LF:HF);mx_cep ind=max(cn);if mx_cep>0.08 a=fs/(LF+ind);else a=0;endpitch=a%畫圖figure(1);subplot(3,1,1);plot(time1,y);title('語音波形');axis tightylim=get(gca,'ylim

28、');line(time1(startIndex),time1(startIndex),ylim,'color','r');line(time1(endIndex),time1(endIndex),ylim,'color','r');xlabel('樣點數(shù)');ylabel('幅度');subplot(3,1,2);plot(frame);axis(0,240,-0.5,0.5)title('一幀語音');xlabel('樣點數(shù)');ylabel('幅

29、度');subplot(3,1,3);time2=-119:1:-1,0:1:120;plot(time2,cepstrum1);axis(-120,120,-0.5,0.5)title('一幀語音的倒譜');xlabel('樣點數(shù)');ylabel('幅度'); 圖4-1 對濁音信號進(jìn)行基音檢測2、利用同態(tài)處理計算給出語音信號的共振峰。y,fs,nbits=wavread('speech_dsp.wav');time1=1:length(y);time=(1:length(y)/fs;frameSize=floor(30*fs/1000);startIndex=round(2080);endIndex=startIndex+frameSize-1;frame=y(startIndex:endIndex);%frameSize=length(frame);frame2=frame.*hamming(length(frame);rwy=rceps(frame2);ylen=length(rwy);cepstrum=rwy(1:ylen/2);%基音檢測LF=floor(fs/500);HF=floor(fs/70);cn=cepstrum(LF:HF);mx_cep,ind=ma