時(shí)域采樣與頻域采樣定理的驗(yàn)證實(shí)驗(yàn)

實(shí)驗(yàn)一時(shí)域采樣與頻域采樣定理的驗(yàn)證實(shí)驗(yàn)實(shí)驗(yàn)?zāi)康臅r(shí)域采樣理論與頻域采樣理論是數(shù)字信號(hào)處理中的重要理論。要求掌握模擬信號(hào)采樣前后頻譜的變化，以及如何選擇采樣頻率才能使采樣后的信號(hào)不丟失信息；要求掌握頻率域采樣會(huì)引起時(shí)域周期化的概念，以及頻率域采樣定理及其對(duì)頻域采樣點(diǎn)數(shù)選擇的指導(dǎo)作用。實(shí)驗(yàn)原理與方法時(shí)域采樣定理的要點(diǎn)是：①對(duì)模擬信號(hào)xa()以T進(jìn)行時(shí)域等間隔理想采樣，形成的采樣信號(hào)的頻譜父。(]Q)會(huì)以采樣角頻率Qs(Qs=2n/T)為周期進(jìn)行周期延拓。公式為X(jQ)=FT[X(t)]=1工X(jQ—jkQ)

aaTas②采樣頻率Qs必須大于等于模擬信號(hào)最高頻率的兩倍以上，才能使采樣信號(hào)的頻譜不產(chǎn)生頻譜混疊。利用計(jì)算機(jī)計(jì)算上式并不方便，下面我們導(dǎo)出另外一個(gè)公式，以便在計(jì)算機(jī)上進(jìn)行實(shí)驗(yàn)。理想采樣信號(hào)X(t)和模擬信號(hào)xa()之間的關(guān)系為：x(t)=x(t)藝6(t—nT)n=一8對(duì)上式進(jìn)行傅里葉變換，得到：X(jQ)=「[x(t)工6(t—nT)]e-jQtdt

a—3an=—3=￡「x(t)6(t—nT)e-jQdt

a—3n=—3上式中，在數(shù)值上x(chóng)(nT)=x(n)，再將?=QT代入，得到：aX(jQ)=Sx(nT)e—jQnTn=—3上式的右邊就是序列的傅里葉變換，即人一X(jQ)=X即)agQT上式說(shuō)明采樣信號(hào)的傅里葉變換可用相應(yīng)序列的傅里葉變換得到，只要將自變量3用QT代替即可。頻域采樣定理的要點(diǎn)是：①對(duì)信號(hào)x(n)的頻譜函數(shù)在[0,2丸]上等間隔采樣N點(diǎn)，得到：X(k)=X(ej?)Nk=0,1,2,…，N—1則N點(diǎn)IDFT[XN(k)]得到的序列就是原序列x(n)以N為周期進(jìn)行周期延拓后的主值區(qū)序列，公式為%(n)=IDFT[XJk)七=[X(k)=X(ej?)Nk=0,1,2,…，N—1i=-s②由上式可知，頻域采樣點(diǎn)數(shù)N必須大于等于時(shí)域離散信號(hào)的長(zhǎng)度M(即NNM)，才能使時(shí)域不產(chǎn)生混疊，則N點(diǎn)IDFT[XN(k)]得到的序列XN(n)就是原序列x(n),即XN(n)=x(n)。如果N>M，XN(n)比原序列尾部多N-M個(gè)零點(diǎn)；如果N<M，則XN(n)=IDFT[XN(k)]發(fā)生了時(shí)域混疊失真，而且XN(n)的長(zhǎng)度N也比x(n)的長(zhǎng)度M短，因此，％(n)與x(n)不相同。對(duì)比上面敘述的時(shí)域采樣原理和頻域采樣原理，得到一個(gè)有用的結(jié)論：這兩個(gè)采樣理論具有對(duì)偶性，即“時(shí)域采樣頻譜周期延拓，頻域采樣時(shí)域信號(hào)周期延拓”。因此，將它們放在一起進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)內(nèi)容及步驟(1)時(shí)域采樣理論的驗(yàn)證。給定模擬信號(hào)x(t)=Ae-atsin(Qt)u(t)式中，A=444.128，a=50t&Q0=50%'2nrad/s，它的幅頻特性曲線(xiàn)如圖1.1所示。//Hz圖1.1xa(t)的幅頻特性曲線(xiàn)現(xiàn)用DFT(FFT)求該模擬信號(hào)的幅頻特性，以驗(yàn)證時(shí)域采樣理論。按照Xa(t)的幅頻特性曲線(xiàn)，選取三種采樣頻率，即Fs=1kHz，300Hz，200Hz。觀(guān)測(cè)時(shí)間選T=64ms。p為使用DFT，首先用下面的公式產(chǎn)生時(shí)域離散信號(hào)，對(duì)三種采樣頻率，采樣序列按順序用X](n)、X2(n)、X3(n)表示。

x(n)=x(nT)=Ae-楓tsin(QnT)u(nT)因?yàn)椴蓸宇l率不同，得到的x1(n)、x2(n)、x3(n)的長(zhǎng)度不同，長(zhǎng)度(點(diǎn)數(shù))用公式N=TpXFs計(jì)算。選FFT的變換點(diǎn)數(shù)為M=64，序列長(zhǎng)度不夠64的尾部加零。X(k)=FFT[x(n)]，k=0,1,2,3，…，M—1式中，k代表的頻率為°2n3k=—k要求：編寫(xiě)實(shí)驗(yàn)程序，計(jì)算xi(n)>x2(n)和x3(n)的幅度特性，并繪圖顯示。觀(guān)察分析頻譜混疊失真。(2)頻域采樣理論的驗(yàn)證。給定信號(hào)如下：'n+10<n<13x(n)=<27-n14<n<260其它編寫(xiě)程序分別對(duì)頻譜函數(shù)X(ejw)=FT[x(n)]在區(qū)間[0,2n]上等間隔采樣32點(diǎn)和16點(diǎn)，得到X32(k)和xi6(k)：X(k)=X(ej3)323=習(xí)k32k=0,1,2,???31X(k)=X(ej3),k=0,1,2,...15163=k16再分別對(duì)X32(k)和X16(k)進(jìn)行32點(diǎn)和16X(k)=X(ej3)323=習(xí)k32k=0,1,2,???31X32(n)=IFFT[X32(k)]32,n=0,1,2,???,31x16(n)=IFFT[X16(k)]16,n=0,1,2,???,15分別畫(huà)出、X32(k)和X/k)的幅度譜，并繪圖顯示x(n)、x32(n)和xjn)的波形，進(jìn)行對(duì)比和分析，驗(yàn)證總結(jié)頻域采樣理論。時(shí)域：A=444.128;a=64*sqrt(2)*pi;w0=64*sqrt(2)*pi;Tp=64/1000;F1=1000;F2=300;F3=200;T1=1/F1;T2=1/F2;T3=1/F3;n1=0:Tp*F1-1;n2=0:Tp*F2-1;n3=0:Tp*F3-1;x1=A*exp(-a*n1*T1).*sin(w0*n1*T1);x2=A*exp(-a*n2*T2).*sin(w0*n2叮2);x3=A*exp(-a*n3*T3).*sin(w0*n3*T3);f1=fft(x1,length(n1));f2=fft(x2,length(n2));f3=fft(x3,length(n3));k1=0:length(f1)-1;fk1=k1/Tp;k2=0:length(f2)-1;fk2=k2/Tp;k3=0:length(f3)-1;fk3=k3/Tp;subplot(3,2,1)stem(n1,x1,'.')title('(a)Fs=1000HZ');xlabel('n');ylabel('x1(n)');subplot(3,2,3)stem(n2,x2,'.')title('(b)Fs=300HZ');xlabel('n');ylabel('x2(n)');subplot(3,2,5)stem(n3,x3,'.')title('(c)Fs=200HZ');xlabel('n');ylabel('x3(n)');subplot(3,2,2)>>plot(fk1,abs(f1));>>title('(a)FT[xa(nT)],Fs=1000Hz')>>xlabel('f(Hz)');>>ylabel('幅度')；>>subplot(3,2,4);>>plot(fk2,abs(f2));>>title('(a)FT[xa(nT)],Fs=300Hz');>>xlabel('f(Hz)');ylabel('幅度')>>plot(fk3,abs(f3));>>title('(a)FT[xa(nT)],Fs=200Hz')>>xlabel('f(Hz)');ylabel('幅度')

頻域：M=27;N=32;n=0:M;xa=0:floor(M/2);xb=ceil(M/2)-1:-1:0;xn=[xa,xb];Xk=fft(xn,1024);X32k=fft(xn,32);X32n=ifft(X32k);X16k=X32k(1:2:N);X16n=ifft(X16k,N/2);subplot(3,2,2);stem(n,xn,'.');boxontitle('(b)三角形序列x(n)');xlabel('n');ylabel('x(n)');axis([0,32,0,20])k=0:1023;wk=2*k/1024;subplot(3,2,1);plot(wk,abs(Xk));title('(a)FT[x(n)]');xlabel('\omegaApi');ylabel('IX(eAjA\omega)l');axis([0,1,0,200])k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),'.');boxon;title('(c)16點(diǎn)頻率采樣')；xlabel('k');ylabel('IX_1_6(k)I');axis([0,8,0,200]);n1=0:N/2-1;subplot(3,2,4);stem(n1,abs(X16n),'.');boxon;title('(d)16IDFTIX_1_6(k)I');xlabel('n');ylabel('IX_1_6(n)I');axis([0,32,0,20]);k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),'.');boxon;title('(e)32點(diǎn)頻率采樣')；xlabel('k');ylabel('IX_3_2(k)I');axis([0,16,0,200]);n1=0:N-1;subplot(3,2,6);stem(n1,abs(X32n),'.');boxon;title('(f)32點(diǎn)IDFTIX_3_2(k)I');xlabel('n');ylabel('IX_3_2(n)I');axis([0,16,0,20])思考題如果序列x(n)的長(zhǎng)度為M,希望得到其頻譜在[0,2丸]上的N點(diǎn)等間隔采樣，當(dāng)N<M時(shí)，如何用一次最少點(diǎn)數(shù)的DFT得到該頻譜采樣？答：從實(shí)驗(yàn)的結(jié)果中，可以看出對(duì)于求頻域采樣點(diǎn)數(shù)N小于原時(shí)域序列長(zhǎng)度M的N點(diǎn)離散頻譜是，可以先對(duì)元序列x(n)以N為周期進(jìn)行周期延拓后去其主值區(qū)序列x(n)=IDFT[X(k)]=[黨x(n+iN)]R(n)NNNNi=-s再計(jì)算N點(diǎn)DFT則得到N點(diǎn)頻域采樣：X(k)=X(ej@),k=0,1,2，…，N-1N1N其所求的N點(diǎn)離散頻譜對(duì)應(yīng)的時(shí)域離散序列是元序列x(n)以N為周期進(jìn)行周期延拓后取主值區(qū)序列，而不是原序列x(n).5、實(shí)驗(yàn)小結(jié)通過(guò)這次實(shí)驗(yàn)，我對(duì)時(shí)域采樣和頻域采樣的理論、定理的理解更加透徹，以前只是課堂上老師說(shuō)的，現(xiàn)在通過(guò)自己親手做了，也知道了其內(nèi)在的原理是怎么樣的了。在這次實(shí)驗(yàn)中，無(wú)論是在時(shí)域還是頻域，對(duì)信號(hào)采樣必須仔細(xì)考慮采樣的參數(shù)：采樣頻譜、采樣周期、采樣點(diǎn)數(shù)，因?yàn)槲以趯?shí)驗(yàn)過(guò)程中，把stem()函數(shù)的參數(shù)設(shè)置錯(cuò)了，其直接導(dǎo)致結(jié)果出錯(cuò)，我對(duì)代碼檢查了好幾遍才改過(guò)來(lái)的，所以，我們要做，就是選取好采樣的參數(shù)，避免另一個(gè)域周期延拓時(shí)發(fā)生混疊，否則，我們采樣所得的數(shù)據(jù)肯定丟失一部分原信號(hào)的信息，我們便無(wú)法對(duì)原信號(hào)對(duì)原信號(hào)進(jìn)行恢復(fù)和正確分析。此次實(shí)驗(yàn)所遇到的問(wèn)題：主要是時(shí)域和頻域的理論通過(guò)自己親自設(shè)計(jì)的時(shí)候，對(duì)一些小知識(shí)點(diǎn)理解的不夠徹底，使得在用MATLAB命令時(shí)，有點(diǎn)難下手。所以，理論應(yīng)用于實(shí)踐，實(shí)踐源于理論，理論很重要的。設(shè)計(jì)二正余弦信號(hào)的譜分析實(shí)驗(yàn)?zāi)康膶W(xué)習(xí)用FFT對(duì)連續(xù)信號(hào)和時(shí)域離散信號(hào)進(jìn)行譜分析的方法，了解可能出現(xiàn)的分析誤差及其原因，以便正確應(yīng)用FFT。實(shí)驗(yàn)原理用FFT對(duì)信號(hào)作頻譜分析是學(xué)習(xí)數(shù)字信號(hào)處理的重要內(nèi)容。經(jīng)常需要進(jìn)行譜分析的信號(hào)是模擬信號(hào)和時(shí)域離散信號(hào)。對(duì)信號(hào)進(jìn)行譜分析的重要問(wèn)題是頻譜分辨率尸和分析誤差。頻譜分辨率直接和FFT的變換區(qū)間N有關(guān)，因?yàn)镕FT能夠?qū)崿F(xiàn)的頻率分辨率是2n/N，因此要求2n/NWF。可以根據(jù)此式選擇FFT的變換區(qū)間N。誤差主要來(lái)自于用FFT作頻譜分析時(shí)，得到的是離散譜，而信號(hào)(周期信號(hào)除外)是連續(xù)譜，只有當(dāng)N較大時(shí)，離散譜的包絡(luò)才能逼近于連續(xù)譜，因此N要適當(dāng)選擇大一些。周期信號(hào)的頻譜是離散譜，只有用整數(shù)倍周期的長(zhǎng)度作FFT，得到的離散譜才能代表周期信號(hào)的頻譜。如果不知道信號(hào)周期，可以盡量選擇信號(hào)的觀(guān)察時(shí)間長(zhǎng)一些。對(duì)模擬信號(hào)進(jìn)行譜分析時(shí)，首先要按照采樣定理將其變成時(shí)域離散信號(hào)。如果是模擬周期信號(hào)，也應(yīng)該選取整數(shù)倍周期的長(zhǎng)度，經(jīng)過(guò)采樣后形成周期序列，按照周期序列的譜分析進(jìn)行。3?實(shí)驗(yàn)步驟及內(nèi)容對(duì)一個(gè)頻率為10Hz，采樣頻率為64Hz的32點(diǎn)余弦序列進(jìn)行譜分析，畫(huà)出其頻譜圖；若將頻率改為11Hz，其他參數(shù)不變，重新畫(huà)出該序列的頻譜圖，觀(guān)察頻率泄漏現(xiàn)象，分析原因。

1DHZ余弦序列11HZ余弦序列n10田時(shí)32點(diǎn)FFT幅頻曲線(xiàn)n11HZ時(shí)32點(diǎn)FFT幅頻曲線(xiàn)M文件：n=0:31;x1=cos(pi*20*n/64);%x2=cos(pi*22*n/64);subplot(2,2,1);%stem(n,x1)%繪制序列的幅頻特性曲線(xiàn)xlabel('n');ylabel('x1(n)');title('10HZ余弦序列')；subplot(2,2,2);stem(n,x2)xlabel('n');ylabel('x2(n)');title('11HZ余弦序列');X1=abs(fft(x1,32));%求x1余弦序歹0subplot(2,2,3)k=0:31;stem(k,X1)%繪制序列的幅頻特性xlabel('k');ylabel('X1(k)');title('10HZ時(shí)32點(diǎn)FFT幅頻曲線(xiàn)');X2=abs(fft(x2,32));%求x2余弦序歹0subplot(2,2,4)k=0:31;stem(k,X2)xlabel('k');ylabel('X2(k)');title('11HZ時(shí)32點(diǎn)FFT幅頻曲線(xiàn)')

泄漏的原因：通過(guò)圖可看出：頻率為10Hz的余弦曲線(xiàn)DFT只有兩個(gè)點(diǎn)不等于零，位于k=5和k=27處，k=5對(duì)應(yīng)于頻率10Hz，k=27對(duì)應(yīng)于頻率54Hz。這樣DFT確實(shí)正確的分辨了余弦信號(hào)的頻率。將頻率改為11Hz，采樣頻率和窗長(zhǎng)度依然為32點(diǎn)，計(jì)算圖像可看出：頻譜圖上k=5和k=27處都有較大的峰值，而其它的點(diǎn)上幅度不再為零。這兩個(gè)峰值對(duì)應(yīng)的頻率為1DHZ余弦序列11HZ余弦序列n10田時(shí)32點(diǎn)FFT幅頻曲線(xiàn)n11HZ時(shí)32點(diǎn)FFT幅頻曲線(xiàn)（2）考察DFT的長(zhǎng)度對(duì)雙頻率信號(hào)頻譜分析的影響。設(shè)待分析的信號(hào)為x（n）-0.5sin（2兀fn）+sin（2兀fn）0<n<1512令兩個(gè)長(zhǎng)度為16的正余弦序列的數(shù)字頻率為f、=0.22及f2=0.34。取N為四個(gè)不同值16,32，64，128。畫(huà)出四個(gè)DFT幅頻圖。分析DFT長(zhǎng)度對(duì)頻譜分辨率的影響。32點(diǎn)FFT幅頻曲線(xiàn)程序：n=0:16;x=0.5*sin(pi*2*0.22*n)+sin(pi*2*0.34*n);X1=abs(fft(x,16));%求余弦序列的16點(diǎn)FFTsubplot(2,2,1)k=0:15;stem(k,X1)%繪制序列的幅頻特性曲線(xiàn)xlabel('k');ylabel('X1(k)');title('16點(diǎn)FFT幅頻曲線(xiàn)')X2=abs(fft(x,32));%求余弦序列的32點(diǎn)FFTsubplot(2,2,2)k=0:31;stem(k,X2)%繪制序列的幅頻特性曲線(xiàn)xlabel('k');ylabel('X2(k)');title('32點(diǎn)FFT幅頻曲線(xiàn)')X3=abs(fft(x,64));%求余弦序歹0的64點(diǎn)FFTsubplot(2,2,3)k=0:63;stem(k,X3)%繪制序列的幅頻特性曲線(xiàn)xlabel('k');ylabel('X3(k)');title('64點(diǎn)FFT幅頻曲線(xiàn)')X4=abs(fft(x,128));%求余弦序歹0的128點(diǎn)FFTsubplot(2,2,4)k=0:127;stem(k,X4)%繪制序列的幅頻特性曲線(xiàn)xlabel('k');ylabel('X4(k)');title('128點(diǎn)FFT幅頻曲線(xiàn)')DFT長(zhǎng)度對(duì)頻譜分辨率的影響：DFT樣本值就是其DTFT在相應(yīng)位置的采樣。在圖中很難看出兩個(gè)峰值，因此要提高它的分辨率，故把R增大，逐漸可以看出它有兩個(gè)峰值，將k換算成數(shù)字頻率f=w/2*pi=k/R.這樣可確定峰值的位置大體在f=0.21和0.35之附近與信號(hào)的給定頻率有一定的誤差，這也是截?cái)嗪托孤稁?lái)的問(wèn)題，在這圖上還可以看到一些較小的峰，這是很難判斷是輸入信號(hào)固有的還是由泄露引起的。這說(shuō)明了增加DFT長(zhǎng)度R減小了相鄰樣本間的頻率間距，提高頻譜的視在分辨率，因而可以提高樣本未知的測(cè)定精度。32點(diǎn)FFT幅頻曲線(xiàn)(3)在上題中若把兩個(gè)正弦波的頻率取得較近，令f=0.22，f2=0.25，試問(wèn)怎樣選擇FFT參數(shù)才能在頻譜分析中分辨出這兩個(gè)分量？要能分清兩個(gè)頻率，分辨率至少應(yīng)達(dá)到f=0.03。因?yàn)榇颂幍臄?shù)字頻率是對(duì)采樣頻率Fs進(jìn)行歸一化，因此總的樣本數(shù)至少要達(dá)到1/0.03=33。加窗以后可以使頻譜函數(shù)更加光滑便于分辨峰值位置和準(zhǔn)確的數(shù)值，為了提高實(shí)際分辨率，應(yīng)該盡量增加信號(hào)的長(zhǎng)度N及DF長(zhǎng)度R,當(dāng)受到條件限制不能提高N，則單獨(dú)提高R可以提高視在分辨率。4.思考題對(duì)于周期序列，如果周期不知道，如何用FFT進(jìn)行譜分析？答：如果X(n)的周期預(yù)先不知道，可先截取M點(diǎn)進(jìn)行DFT，即x(n)=x(n)R(n)Xm(k)=DFT[Xm(n)]0<=K<=M-1再將截取長(zhǎng)度擴(kuò)大一倍，截取x(n)=x(n)R(n)X2m(k)=DFT[X2m(n)]0<=K<=2M-1比較XM(k)和X2M(k)，如果二者的主譜差別滿(mǎn)足分析誤差要求，則以XM(k)或X2m(k)近似表示x(n)的頻譜，否則繼續(xù)將截取長(zhǎng)度加倍，直至前后再次分析所得主譜頻率差別滿(mǎn)足誤差要求，設(shè)最后截取長(zhǎng)度為1M，則X^(K0)表示w=[2K/(iM)]k0點(diǎn)的譜線(xiàn)強(qiáng)度如何選擇FFT的變換區(qū)間？(包括非周期信號(hào)和周期信號(hào))答：對(duì)于非周期信號(hào)：有頻譜分辨率F，而頻譜分辨率直接和FFT的變換區(qū)間有關(guān)，因?yàn)镕FT能夠?qū)崿F(xiàn)的頻率分辨率是2n/N...因此有最小的N>2n/F。就可以根據(jù)此式選擇FFT的變換區(qū)間。對(duì)于周期信號(hào)，周期信號(hào)的頻譜是離散譜，只有用整數(shù)倍周期的長(zhǎng)度作FFT，得到的離散譜才能代表周期信號(hào)的頻譜。設(shè)計(jì)三語(yǔ)音信號(hào)濾波處理實(shí)驗(yàn)?zāi)康牧私庹Z(yǔ)音信號(hào)的產(chǎn)生、采集，能繪制語(yǔ)音信號(hào)的頻率響應(yīng)曲線(xiàn)及頻譜圖；學(xué)會(huì)用MATLAB對(duì)語(yǔ)音信號(hào)進(jìn)行分析和處理；掌握用濾波器去除語(yǔ)音信號(hào)噪聲的方法，觀(guān)察去噪前后的語(yǔ)音信號(hào)。實(shí)驗(yàn)原理語(yǔ)音信號(hào)的采集在MATLAB軟件平臺(tái)下，利用函數(shù)wavread()對(duì)語(yǔ)音信號(hào)采集，并記錄采樣頻率和采樣點(diǎn)數(shù)。將語(yǔ)音信號(hào)轉(zhuǎn)換成計(jì)算機(jī)能夠運(yùn)算的有限長(zhǎng)序列。wavread函數(shù)的調(diào)用格式如下：y=wavread(file)讀取file所規(guī)定的wav文件，返回采樣值放在向量y中。[y,fs]=wavread(file)采樣值放在向量y中，fs表示采樣頻率(Hz)。用FFT作譜分析FFT即快速傅立葉變換，它是從DFT運(yùn)算中發(fā)展起來(lái)的，利用系數(shù)Wf的對(duì)稱(chēng)性和周期性減少運(yùn)算量。長(zhǎng)度為N的序列直接計(jì)算DFT需要N2次復(fù)乘和N(N1)次復(fù)加，N、…而用FFT進(jìn)行運(yùn)算一般需要ylog2N次復(fù)乘和Nlog2N次復(fù)加，從而使DFT的運(yùn)算大大簡(jiǎn)化。用FFT對(duì)連續(xù)信號(hào)進(jìn)行譜分析的步驟如下:圖3.1連續(xù)信號(hào)譜分析過(guò)程引入前置低通濾波器LPF是為了消除或減少時(shí)域連續(xù)信號(hào)轉(zhuǎn)換成序列時(shí)可能出現(xiàn)的頻譜混疊現(xiàn)象。w(n)表示時(shí)域有限的窗函數(shù)。設(shè)計(jì)濾波器去除語(yǔ)音信號(hào)的噪聲通過(guò)wavread()函數(shù)將語(yǔ)音信號(hào)讀入，通過(guò)頻率采樣及fft()產(chǎn)生信號(hào)，并對(duì)之加噪，通過(guò)窗函數(shù)法設(shè)計(jì)濾波器濾掉該語(yǔ)音信號(hào)的噪聲，并對(duì)比濾波前后的語(yǔ)音波形和頻譜。實(shí)驗(yàn)內(nèi)容⑴利用Windows下的錄音機(jī)錄制一段自己的話(huà)音，時(shí)間在1s內(nèi)。然后在Matlab軟件平臺(tái)下，利用wavread函數(shù)對(duì)語(yǔ)音信號(hào)進(jìn)行采樣，記住采樣頻率和采樣點(diǎn)數(shù)；答：程序：x=wavread('D:\test.wav');[x,fs]=wavread('D:\test.wav');%把語(yǔ)音信號(hào)進(jìn)行加載入Matlab仿真軟件平臺(tái)中sound(x,fs);%對(duì)加載的語(yǔ)音信號(hào)進(jìn)行回放stem(x);title('語(yǔ)音信號(hào)的時(shí)域波形')(2)畫(huà)出語(yǔ)音信號(hào)的時(shí)域波形，對(duì)采樣后的語(yǔ)音信號(hào)進(jìn)行快速傅立葉變換，得到信號(hào)的頻譜特性，畫(huà)出采樣信號(hào)的時(shí)域波形和頻譜圖；答：程序：Fs=8000;%給出抽樣頻率x=wavread('D:\test.wav');[x,fs,nbits]=wavread('D:\test.wav');%把語(yǔ)音信號(hào)進(jìn)行加載入Matlab仿真軟件平臺(tái)中sound(x,fs,nbits);n=length(x);%求出語(yǔ)音信號(hào)的長(zhǎng)度X=fft(x,n);%傅里葉變換subplot(3,2,1);plot(x);title('原始信號(hào)時(shí)域波形');subplot(3,2,2);plot(abs(X));title('原始域波形')[x,fs,nbits]=wavread('D:\test.wav');n=length(x);%求出語(yǔ)音信號(hào)的長(zhǎng)度noise=0.04火sin(10000*pi*x);%sin函數(shù)產(chǎn)生正弦噪聲s=x+noise;%語(yǔ)音信號(hào)加入正弦噪聲sound(s);subplot(3,2,3);plot(s);title('加正弦噪語(yǔ)音信號(hào)的時(shí)域波形');S=fft(s);%傅里葉變換subplot(3,2,4);plot(abs(S));title('加正弦噪語(yǔ)音信號(hào)的頻域波形')[x,fs,nbits]=wavread('D:\test.wav');n=length(x);%求出語(yǔ)音信號(hào)的長(zhǎng)度noise=rand(size(x/5));%產(chǎn)生白噪聲s=x+noise;%語(yǔ)音信號(hào)加入白噪聲sound(s);subplot(3,2,5);plot(s);title('加白噪聲語(yǔ)音信號(hào)的時(shí)域波形');S=fft(s);%傅里葉變換subplot(3,2,6);plot(abs(S));title('加白噪聲語(yǔ)音信號(hào)的頻域波形')根據(jù)對(duì)語(yǔ)音信號(hào)譜分析的結(jié)果，確定濾除噪聲所需濾波器的技術(shù)指標(biāo)，設(shè)計(jì)合適的數(shù)字濾波器，并畫(huà)出濾波器的頻域響應(yīng)；答：根據(jù)對(duì)加噪語(yǔ)音信號(hào)譜的分析得濾除噪聲所需的濾波器應(yīng)為低通濾波器。技術(shù)指標(biāo)如下：Fp=1000HZ;Fs=1200HZ;rp=0.1;rs=60程序:Fp=1000;Fs=1200;wp=2火pi*Fp;ws=2火pi*Fs;rp=0.1;rs=60[n,wn]=ellipord(wp,ws,rp,rs,'s')[b,a]=ellip(n,rp,rs,wn,'low','s')fk=linspace(1,3000,100)w=2火pi*fkh=freqs(b,a,w)magh=abs(h)subplot(2,1,1);plot(fk,magh);title('低通幅頻響應(yīng)')phah=unwrap(angle(h))subplot(2,1,2);plot(fk,phah);title('低通相頻響應(yīng)')(4)用所設(shè)計(jì)的濾波器對(duì)采集的信號(hào)進(jìn)行濾波，在同一個(gè)窗口畫(huà)出濾波前后信號(hào)的時(shí)域波形和頻譜，并對(duì)濾波前后的信號(hào)進(jìn)行對(duì)比，分析信號(hào)的變化；答：程序：Ft=8000;Fp=1000;Fs=1200;wp=2火pi*Fp/Ft;ws=2火pi*Fs/Ft;fp=2火Ft*tan(wp/2);fs=2火Fs*tan(wp/2);[n11,wn11]=buttord(wp,ws,1,50,'s');%求低通濾波器的階數(shù)和截止頻率[b11,a11]=butter(n11,wn11,'s');%求S域的頻率響應(yīng)的參數(shù)[num11,den11]=bilinear(b11,a11,0.5);%利用雙線(xiàn)性變換實(shí)現(xiàn)頻率響應(yīng)S域到Z域的變換[x,fs,nbits]=wavread('D:\test.wav');n=length(x);%求出語(yǔ)音信號(hào)的長(zhǎng)度noise=0.04火sin(10000*pi*x);%sin函數(shù)產(chǎn)生正弦噪聲s=x+noise;%語(yǔ)音信號(hào)加入正弦噪聲z11=filter(num11,den11,s);sound(z11);m11=fft(z11);%求濾波后的信號(hào)subplot(2,2,1);plot(z11);title('加入正弦噪聲濾波后的時(shí)域波形');subplot(2,2,2);plot(abs(m11),'r');title('加入正弦噪聲濾波后的頻域波形');[x,fs,nbits]=wavread('D:\test.wav');n=length(x);%求出語(yǔ)音信號(hào)的長(zhǎng)度noise=rand(size(x/5));%產(chǎn)生白噪聲s=x+noise;%語(yǔ)音信號(hào)加入白噪聲z11=filter(num11,den11,s);sound(z11);m11=fft(z11);%求濾波后的信號(hào)subplot(2,2,3);plot(z11);title('加入白噪聲濾波后的時(shí)域波形');subplot(2,2,4);plot(abs(m11),'r');title('加入白噪聲濾波后的頻域波形');

(5)回放語(yǔ)音信號(hào)；sound(x)(6)為使編制的程序操作方便，設(shè)計(jì)一個(gè)信號(hào)處理用戶(hù)界面。在所設(shè)計(jì)的系統(tǒng)界面上可以選擇濾波器的類(lèi)型，輸入濾波器的參數(shù)，顯示濾波器的頻率響應(yīng)，選擇信號(hào)等。數(shù)字濾波器的設(shè)計(jì)及實(shí)現(xiàn)1.設(shè)計(jì)目的熟悉IIR數(shù)字濾波器的設(shè)計(jì)原理和方法;學(xué)會(huì)調(diào)用MATLAB信號(hào)處理工具箱中的濾波器設(shè)計(jì)函數(shù)設(shè)計(jì)各種IIR數(shù)字濾波器，學(xué)會(huì)根據(jù)濾波要求確定濾波器指標(biāo)參數(shù);掌握用IIR數(shù)字濾波器的MATLAB實(shí)現(xiàn)方法，并能繪制濾波器的幅頻特性、相頻特性;(4)2.(4)2.設(shè)計(jì)原理抑制載波單頻調(diào)幅信號(hào)的數(shù)學(xué)表達(dá)式為T(mén)OC\o"1-5"\h\zS0)=cos(2兀f01)cos(2兀ft)=1lcos(2兀(f-f0)t)+cos(2兀(f+f0)t)](4.1)uccucu其中，cos(2兀ft)稱(chēng)為載波，f為載波頻率，cos(2兀ft)稱(chēng)為單頻調(diào)制信號(hào)，f為調(diào)制正弦ccUU波信號(hào)頻率，且滿(mǎn)足fc>f0。由(4.1)式可見(jiàn)，所謂抑制載波單頻調(diào)制信號(hào)，就是兩個(gè)正弦信號(hào)相乘，它有2個(gè)頻率成分：和頻f+f°，差頻fc-f°，這兩個(gè)頻率成分關(guān)于載波頻率f對(duì)稱(chēng)。所以，1路抑制載波單頻調(diào)幅信號(hào)的頻譜圖是關(guān)于載波頻率fc對(duì)稱(chēng)的兩根譜線(xiàn)。設(shè)計(jì)內(nèi)容調(diào)用信號(hào)產(chǎn)生函數(shù)mstg產(chǎn)生由三路抑制載波調(diào)幅信號(hào)相加構(gòu)成的復(fù)合信號(hào)st，觀(guān)察st的時(shí)域波形和幅頻特性曲線(xiàn)；functionst=mstg%產(chǎn)生信號(hào)序列st，并顯示st的時(shí)域波形和頻譜%st=mstg返回三路調(diào)幅信號(hào)相加形成的混合信號(hào)，長(zhǎng)度N=800N=800;%信號(hào)長(zhǎng)度N為800Fs=10000;T=1/Fs;Tp=N*T;%采樣頻率Fs=10kHz，Tp為采樣時(shí)間t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;fc1=Fs/10;%第1路調(diào)幅信號(hào)載波頻率fc1=1000Hzfm1=fc1/10;%第1路調(diào)幅信號(hào)的調(diào)制信號(hào)頻率fm1=100Hzfc2=Fs/20;%第2路調(diào)幅信號(hào)載波頻率fc2=500Hzfm2=fc2/10;%第2路調(diào)幅信號(hào)的調(diào)制信號(hào)頻率fm2=50Hzfc3=Fs/40;%第3路調(diào)幅信號(hào)載波頻率fc3=250Hzfm3=fc3/10;%第3路調(diào)幅信號(hào)的調(diào)制信號(hào)頻率fm3=25Hzxt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t);%產(chǎn)生第1路調(diào)幅信號(hào)xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);%產(chǎn)生第2路調(diào)幅信號(hào)xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);%產(chǎn)生第3路調(diào)幅信號(hào)st=xt1+xt2+xt3;%三路信號(hào)相加，得到復(fù)合信號(hào)fxt=fft(st,N);%計(jì)算信號(hào)st的頻譜%以下為繪圖命令subplot(2,1,1);plot(t,st);grid;xlabel('t/s');ylabel('s(t)');axis([0,Tp,min(st),max(st)]);title('(a)s(t)的波形')subplot(2,1,2);stem(f,abs(fxt)/max(abs(fxt)),'.');grid;title('(b)s(t)的頻譜')axis([0,Fs/8,0,1.2]);xlabel('f/Hz');ylabel('幅度')；該函數(shù)產(chǎn)生由三路已知載波調(diào)幅信號(hào)相加構(gòu)成的復(fù)合信號(hào)st，并繪圖顯示st的時(shí)域波形和幅頻特性曲線(xiàn)如圖4.1所示。t/s

1O(a)s(t)的波形f/Hz(b)s(t)的頻譜圖4.1三路調(diào)幅信號(hào)st的時(shí)域波形和幅頻特性曲線(xiàn)由圖4.1可見(jiàn)，三路信號(hào)時(shí)域混疊無(wú)法在時(shí)域進(jìn)行分離，但頻域是分離的。容易看出，這三路調(diào)幅信號(hào)的載波頻率分別為250Hz、500Hz和1000Hz。1O(a)s(t)的波形f/Hz要求將st中的三路調(diào)幅信號(hào)分離，通過(guò)觀(guān)察st的幅頻特性曲線(xiàn)，分別確定可以分離st中三路抑制載波單頻調(diào)幅信號(hào)的三個(gè)濾波器(低通濾波器、帶通濾波器、高通濾波器)的通帶截止頻率和阻帶截止頻率，要求濾波器的通帶最大衰減為0.1dB，阻帶最小衰減為60dB；，低通濾波器的通帶截止頻率應(yīng)大于275HZ所以選300HZ，阻帶截止頻率應(yīng)小于450HZ，選400HZ，低通濾波器程序與幅頻特性圖要求通帶截止300HZ，最大衰減0.1dB；通帶截止400HZ，最小衰減60dB。,帶通濾波器的通帶截止頻率應(yīng)包含450HZ~550HZ故選400HZ~600HZ阻帶截止頻率應(yīng)包含275HZ和900HZ,故選300HZ~800HZ，帶通濾波器程序與幅頻特性圖，要求通帶截止400HZ和600HZ,最大衰減0.1dB;阻帶截止300HZ和800HZ,最小衰減60dB。，高通濾波器的通帶截止頻率應(yīng)小于900HZ，故選800HZ，阻帶截止頻率應(yīng)大于550HZ，故選600HZ高通濾波器程序與幅頻特性圖要求通帶截止800HZ，最大衰減0.1dB;阻帶截止600H最小衰減60dB編程調(diào)用MATLAB濾波器設(shè)計(jì)函數(shù)分別設(shè)計(jì)這三個(gè)數(shù)字濾波器，并繪圖顯示其幅頻特性曲線(xiàn)；調(diào)用濾波函數(shù)filter，用所設(shè)計(jì)的三個(gè)濾波器分別對(duì)復(fù)合信號(hào)st進(jìn)行濾波，分離出st中的三路不同載波頻率的調(diào)幅信號(hào)，并繪圖顯示濾波后信號(hào)的時(shí)域波形和頻譜，觀(guān)察分離效果。⑶(4)答案:，低通濾波器的通帶截止頻率選300HZ,阻帶截止頻率選400HZ,低通濾波器程序與幅頻特性圖要求通帶截止300HZ，最大衰減0.1dB；通帶截止400HZ，最小衰減60dB。程序：functionst=mstg%st=mstg返回三路調(diào)幅信號(hào)相加的形成的混合信號(hào)％長(zhǎng)度N=800N=800;Fs=10000;T=1/Fs;Tp=N*T;%采樣頻率Fs=10kHz,Tp為采樣時(shí)間t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;fc1=Fs/10;fm1=fc1/10;%第1路調(diào)幅信號(hào)載波頻率fc1=1000Hz,調(diào)制頻fm1=100Hzfc2=Fs/20;fm2=fc2/10;%第2路調(diào)幅信號(hào)載波頻率fc2=500Hz,調(diào)制頻率fm2=50Hzfc3=Fs/40;fm3=fc3/10;%第3路調(diào)幅信號(hào)載波頻率fc3=250Hz,調(diào)制頻率fm3=25Hzxt1=cos(2火pi*fm1*t).*cos(2火pi*fc1*t);%產(chǎn)生第1路調(diào)幅信號(hào)xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);%產(chǎn)生第2路調(diào)幅信號(hào)xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);%產(chǎn)生第3路調(diào)幅信號(hào)st=xt1+xt2+xt3;%3路信號(hào)相加得到復(fù)合信號(hào)fp=300;fs=400;Fs=10000;Wp=2火fp/Fs;Ws=2火fs/Fs;Rp=0.1;Rs=60;%低通濾波器指標(biāo)[N,Wpo]=ellipord(Wp,Ws,Rp,Rs);[B,A]=ellip(N,Rp,Rs,Wpo);%計(jì)算濾波器系統(tǒng)函數(shù)系數(shù)subplot(3,1,1);W=0:0.01:pi;[H,W]=freqz(B,A,W);%計(jì)算頻率向量W上的濾波器頻率響應(yīng)H=20*log(abs(H));%縱坐標(biāo)單位取dBplot(W/pi,H)%繪制幅頻響應(yīng)曲線(xiàn)gridon%添加網(wǎng)格axis([0,1,-200,10])xlabel('\omega/\pi(rad/s)');ylabel('|H(e八j八\omega)|/dB');subplot(3,1,2);y3=filter(B,A,st);plot(t,y3);axis([0,Tp,min(st),max(st)]);title('(a)低通濾波器的波形');subplot(3,1,3);Y3=abs(fft(y3,800));stem(f,Y3/max(Y3),'.');grid;title('(b)低通濾波器的頻譜')axis([0,Fs/8,0,1.2]);xlabel('f/Hz');ylabel('幅度');圖4-2:低通濾波器幅頻特性曲線(xiàn),帶通濾波器的通帶截止頻率應(yīng)包含450HZ~550HZ故選400HZ~600HZ阻帶截止頻率應(yīng)包含275HZ和900HZ,故選300HZ~800HZ，帶通濾波器程序與幅頻特性圖，要求通帶截止400HZ和600HZ,最大衰減0.1dB;阻帶截止300HZ和800HZ,最小衰減60dB。程序：functionst=mstg%st=mstg返回三路調(diào)幅信號(hào)相加的形成的混合信號(hào)％長(zhǎng)度N=800N=800;Fs=10000;T=1/Fs;Tp=N*T;%采樣頻率Fs=10kHz,Tp為采樣時(shí)間t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;fc1=Fs/10;fm1=fc1/10;%第1路調(diào)幅信號(hào)載波頻率fc1=1000Hz,調(diào)制頻fm1=100Hzfc2=Fs/20;fm2=fc2/10;%第2路調(diào)幅信號(hào)載波頻率fc2=500Hz,調(diào)制頻率fm2=50Hzfc3=Fs/40;fm3=fc3/10;%第3路調(diào)幅信號(hào)載波頻率fc3=250Hz,調(diào)制頻率fm3=25Hzxt1=cos(2火pi*fm1*t).*cos(2火pi*fc1*t);%產(chǎn)生第1路調(diào)幅信號(hào)xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);%產(chǎn)生第2路調(diào)幅信號(hào)xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);%產(chǎn)生第3路調(diào)幅信號(hào)st=xt1+xt2+xt3;%3路信號(hào)相加得到復(fù)合信號(hào)Fs=10000;Wp=[2*400/Fs,2*600/Fs];Ws=[2*300/Fs,2*800/Fs];Rp=0.1;Rs=60;%帶通濾波器指標(biāo)[N,Wpo]=ellipord(Wp,Ws,Rp,Rs);[B,A]=ellip(N,Rp,Rs,Wpo);%計(jì)算濾波器系統(tǒng)函數(shù)系數(shù)subplot(3,1,1);W=0:0.01:pi;[H,W]=freqz(B,A,W);%計(jì)算頻率向量W上的濾波器頻率響應(yīng)H=20*log(abs(H));%縱坐標(biāo)單位取dB

plot(W/pi,H)%繪制幅頻響應(yīng)曲線(xiàn)axis([0,1,-200,10])xlabel('\omega/\pi(rad/s)');ylabel('|H(e八j八\omega)|/dB');subplot(3,1,2);y3=filter(B,A,st);plot(t,y3);axis([0,Tp,min(st),max(st)]);title('(a)plot(W/pi,H)%繪制幅頻響應(yīng)曲線(xiàn)axis([0,1,-200,10])xlabel('\omega/\pi(rad/s)');ylabel('|H(e八j八\omega)|/dB');subplot(3,1,2);y3=filter(B,A,st);plot(t,y3);axis([0,Tp,min(st),max(st)]);title('(a)帶通濾波器的波形');subplot(3,1,3);Y3=abs(fft(y3,800));stem(f,Y3/max(Y3),'.');grid;title('(b)帶通濾波器的頻譜')axis([0,Fs/8,0,1.2]);xlabel('f/Hz');ylabel('幅度');FileEditView[nsertT&oIeDesktopWindowHelp您％0?豎屆、且□E□03210-10.10.20.3040.5O.E0.70.80.91w/p[rad閭

同帶通濾'波器的波形0.010.020.030.040一050.0&0.070.0S㈣帶通濾波器的頻譜IIiiiiIIIIIIiiii■■■IiiiiIIIIiiiiii暨0.520040060080010001200圖4-3:帶通濾波器幅頻特性曲線(xiàn)，高通濾波器的通帶截止頻率應(yīng)小于900HZ,故選800HZ，阻帶截止頻率應(yīng)大于550HZ,故選600HZ高通濾波器程序與幅頻特性圖要求通帶截止800HZ，最大衰減0.1dB;阻帶截止600H最小衰減60dB程序：functionst=mstg%st=mstg返回三路調(diào)幅信號(hào)相加的形成的混合信號(hào)％長(zhǎng)度N=800N=800;Fs=10000;T=1/Fs;Tp=N*T;%采樣頻率Fs=10kHz,Tp為采樣時(shí)間t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;fc1=Fs/10;fm1=fc1/10;%第1路調(diào)幅信號(hào)載波頻率fc1=1000Hz,調(diào)制頻fm1=100Hzfc2=Fs/20;fm2=fc2/10;%第2路調(diào)幅信號(hào)載波頻率fc2=500Hz,調(diào)制頻率fm2=50Hzfc3=Fs/40;fm3=fc3/10;%第3路調(diào)幅信號(hào)載波頻率fc3=250Hz,調(diào)制頻率fm3=25Hzxt1=cos(2火pi*fm1*t).*cos(2火pi*fc1*t);%產(chǎn)生第1路調(diào)幅信號(hào)xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);%產(chǎn)生第2路調(diào)幅信號(hào)xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);%產(chǎn)生第3路調(diào)幅信號(hào)st=xt1+xt2+xt3;%3路信號(hào)相