數字信號處理上機實驗答案第十章上機實驗數字信號處理是一門理論和實際密切結合的課程，為深入掌握課程內容，最好在學習理論的同時，做習題和上機實驗。上機實驗不僅可以幫助讀者深入的理解和消化基本理論，而且能鍛煉初學者的獨立解決問題的能力。本章在第二版的基礎上編寫了六個實驗，前五個實驗屬基礎理論實驗，第六個屬應用綜合實驗。實驗一系統響應及系統穩定性。實驗二時域采樣與頻域采樣。實驗四IIR數字濾波器設計及軟件實現。實驗五FIR數字濾波器設計與軟件實現實驗六應用實驗——數字信號處理在雙音多頻撥號系統中的應用任課教師根據教學進度，安排學生上機進行實驗。建議自學的讀者在學習完第一章后作實驗一；在學習完第三、四章后作實驗二和實驗三；實驗四IIR數字濾波器設計及軟件實現在。學習完第六章進行；實驗五在學習完第七章后進行。實驗六綜合實驗在學習完第七章或者再后些進行；實驗六為綜合實驗，在學習完本課程后再進行。function?tstem(xn,yn)%時域序列繪圖函數%xn:信號數據序列，yn:繪圖信號的縱坐標名稱(字符串)n=0:length(xn)-1;stem(n,xn,'.');box?onxlabel('n');ylabel(yn);axis([0,n(end),min(xn),*max(xn)])實驗一:系統響應及系統穩定性1.實驗目的(1)掌握求系統響應的方法。(2)掌握時域離散系統的時域特性。(3)分析、觀察及檢驗系統的穩定性。2.實驗原理與方法在時域中，描寫系統特性的方法是差分方程和單位脈沖響應，在頻域可以用系統函數描述系統特性。已知輸入信號可以由差分方程、單位脈沖響應或系統函數求出系統對于該輸入信號的響應，本實驗僅在時域求解。在計算機上適合用遞推法求差分方程的解，最簡單的方法是采用函數conv函數計算輸入信號和系統的單位脈沖響應的線性卷積，求出系系統的時域特性指的是系統的線性時不變性質、因果性和穩定性。重點分析實驗系統的穩定性，包括觀察系統的暫態響應和穩定響應。系統的穩定性是指對任意有界的輸入信號，系統都能得到有界的系統響應。或者系統的單位脈沖響應滿足絕對可和的條件。系統的穩定性由其差分方程的系數決定。實際中檢查系統是否穩定，不可能檢查系統對所有有界的輸入信號，輸出是否都是有界輸出，或者檢查系統的單位脈沖響應滿足絕對可和的條件。可行的方法是在系統的輸入端加入單位階躍序列，如果系統的輸出趨近一個常數(包括零)，就可以斷定系統是穩定的[19]。系統的穩態輸出是指當n時，系統的輸出。如果系統穩定，信號加入系統后，系統輸出的開始一段稱為暫態效應，隨n的加大，幅度趨于穩定，達到穩注意在以下實驗中均假設系統的初始狀態為零。3．實驗內容及步驟(1)編制程序，包括產生輸入信號、單位脈沖響應序列的子程序，制信號波形的功能。(2)給定一個低通濾波器的差分方程為輸入信號x(n)=R(n)18x(n)=u(n)2a)分別求出系統對x(n)=R(n)和x(n)=u(n)的響應序列，并畫182(3)給定系統的單位脈沖響應為1021218用線性卷積法分別求系統h(n)和h(n)對x(n)=1218應，并畫出波形。(4)給定一諧振器的差分方程為000a)用實驗方法檢查系統是否穩定。輸入信號為u(n)時，畫出系b)給定輸入信號為求出系統的輸出響應，并畫出其波形。4．思考題(1)如果輸入信號為無限長序列，系統的單位脈沖響應是有限長序列，可否用線性卷積法求系統的響應如何求(2)如果信號經過低通濾波器，把信號的高頻分量濾掉，時域信號會有何變化，用前面第一個實驗結果進行分析說明。5．實驗報告要求(1)簡述在時域求系統響應的方法。(2)簡述通過實驗判斷系統穩定性的方法。分析上面第三個實驗的穩定輸出的波形。(3)對各實驗所得結果進行簡單分析和解釋。(4)簡要回答思考題。(5)打印程序清單和要求的各信號波形。10.1.2實驗程序清單%實驗1：系統響應及系統穩定性closeall;clearallfilter統對u(n)的響應判斷穩定性======A=[1,];B=[,];%系統差分方程系數向量B和Ax1n=[11111111zeros(1,50)];%產生信號x1(n)=R8(n)x2n=ones(1,128);hn=impz(B,A,58);%產生信號x2(n)=u(n)%求系統單位脈沖響應h(n)subplot(2,2,1);y='h(n)';tstem(hn,y);%調用函數tstem繪圖title('(a)系統單位脈沖響應h(n)');boxony1n=filter(B,A,x1n);%求系統對x1(n)的響應y1(n)subplot(2,2,2);y='y1(n)';tstem(y1n,y);title('(b)系統對R8(n)的響應y1(n)');boxony2n=filter(B,A,x2n);%求系統對x2(n)的響應y2(n)subplot(2,2,4);y='y2(n)';tstem(y2n,y);title('(c)系統對u(n)的響應y2(n)');boxon%===內容2：調用conv函數計算卷積============================x1n=[11111111];%產生信號x1(n)=R8(n)h1n=[ones(1,10)zeros(1,10)];h2n=[11zeros(1,10)];y21n=conv(h1n,x1n);y22n=conv(h2n,x1n);figure(2)subplot(2,2,1);y='h1(n)';tstem(h1n,y);%調用函數tstem繪圖title('(d)系統單位脈沖響應h1(n)');boxonsubplot(2,2,2);y='y21(n)';tstem(y21n,y);title('(e)h1(n)與R8(n)的卷積y21(n)');boxonsubplot(2,2,3);y='h2(n)';tstem(h2n,y);%調用函數tstem繪圖title('(f)系統單位脈沖響應h2(n)');boxonsubplot(2,2,4);y='y22(n)';tstem(y22n,y);title('(g)h2(n)與R8(n)的卷積y22(n)');boxon%=========內容3：諧振器分析========================un=ones(1,256);%產生信號u(n)n=0:255;xsin=sin*n)+sin*n);%產生正弦信號A=[1,,];B=[1/,0,-1/];%系統差分方程系數向量B和Ay31n=filter(B,A,un);%諧振器對u(n)的響應y31(n)y32n=filter(B,A,xsin);%諧振器對u(n)的響應y31(n)figure(3)subplot(2,1,1);y='y31(n)';tstem(y31n,y);title('(h)諧振器對u(n)的響應y31(n)');boxonsubplot(2,1,2);y='y32(n)';tstem(y32n,y);title('(i)諧振器對正弦信號的響應y32(n)');boxon10.1.3實驗程序運行結果及分析討論實驗內容(2)系統的單位沖響應、系統對x(n)=R(n)和x(n)=u(n)182實驗內容(3)系統h(n)和h(n)對x(n)=R(n)的輸出響應分別如圖1218實驗內容(4)系統對u(n)和x(n)=sin(0.014n)+sin(0.4n)的響應序列分別如圖(h)和(i)所示。由圖(h)可見，系統對u(n)的響應逐漸衰減到零，所以系統穩定。由圖(i)可見，系統對x(n)=sin(0.014n)+sin(0.4n)的穩態響rad。0.10.080.060.040.0200.10.080.060.040.020n0.50055n32005n0.60.40.202040n10.80.60.40.2050100n8642001020n8642005n1505n0.0.0200.020.040501001502000.040n10.50505010.1.4簡答思考題n(1)如果輸入信號為無限長序列，系統的單位脈沖響應是有限長序列，可否用線性卷積法求系統的響應。①對輸入信號序列分段；②求單位脈沖響應h(n)與各段的卷積；③將各段卷積結果相加。具體實現方法有第三章介紹的重疊相加法和重疊保留法。(2)如果信號經過低通濾波器，把信號的高頻分量濾掉，時域信號的劇烈變化將被平滑，由實驗內容(1)結果圖10.1.1(a)、(b)和(c)可見，經過系統低通濾波使輸入信號6(n)、x(n)=R(n)和x(n)=u(n)的階182躍變化變得緩慢上升與下降。實驗二時域采樣與頻域采樣10.2.1實驗指導的時域采樣理論與頻域采樣理論是數字信號處理中的重要理論。要求掌握模擬信號采樣前后頻譜的變化，以及如何選擇采樣頻率才能使采樣后的信號不丟失信息；要求掌握頻率域采樣會引起時域周期化的概念，以及頻率域采樣定理及其對頻域采樣點數選擇的指導作用。2.實驗原理與方法a)對模擬信號x(t)以間隔T進行時域等間隔理想采樣，形成的采a樣信號的頻譜X?(j)是原模擬信號頻譜X(j)以采樣角頻率a(=2/T)為周期進行周期延拓。公式為：ssX?(j)=FT[(t)]=1X(jjn)aaTasb)采樣頻率必須大于等于模擬信號最高頻率的兩倍以上，才s能使采樣信號的頻譜不產生頻譜混疊。利用計算機計算上式并不方便，下面我們導出另外一個公式，以便用計機上進行實驗。txt間的關系為：aa(t)=x(t)6(tnT)aa對上式進行傅立葉變換，得到：aan=an在上式的積分號內只有當tnT時，才有非零值，因此：X?(j)x(nT)ejnTaaxnT)＝x(n)，再將T代入，得到：aX?(j)x(n)ejna上式的右邊就是序列的傅立葉變換X(ej)，即X?(j)X(ej)aT上式說明理想采樣信號的傅立葉變換可用相應的采樣序列的傅立葉變換理的要點是：到X(k)X(ej),k0,1,2,,NNN延拓后的主值區序列，公式為：x(n)IDFT[X(k)][x(niN)]R(n)NNNNib)由上式可知，頻域采樣點數N必須大于等于時域離散信號的長度M(即N≥M)，才能使時域不產生混疊，則N點IDFT[X(k)]得到的N序列x(n)就是原序列x(n),即x(n)=x(n)。如果N>M，x(n)比原NNNNM；如果N<M，z則x(n)=IDFT[X(k)]發生NNNN在數字信號處理的應用中，只要涉及時域或者頻域采樣，都必須服從這兩個采樣理論的要點。對比上面敘述的時域采樣原理和頻域采樣原理，得到一個有用的結論，這兩個采樣理論具有對偶性：“時域采樣頻譜周期延拓，頻域采樣時域信號周期延拓”。因此放在一起進行實驗。3.實驗內容及步驟(1)時域采樣理論的驗證。a00圖10.2.1x(t)的幅頻特性曲線a現用DFT(FFT)求該模擬信號的幅頻特性，以驗證時域采樣理安照x(t)的幅頻特性曲線，選取三種采樣頻率，即F=1kHz，as300Hz，200Hz。觀測時間選T=50ms。p為使用DFT，首先用下面公式產生時域離散信號，對三種采樣頻率，采樣序列按順序用x(n)，x(n)，x(n)表示。23a0因為采樣頻率不同，得到的x(n)，x(n)，x(n)的長度不同，長23度(點數)用公式N=T根F計算。選FFT的變換點數為M=64，序列psFFTxnkMkM要求：編寫實驗程序，計算x(n)、x(n)和x(n)的幅度特性，并繪圖23顯示。觀察分析頻譜混疊失真。(2)頻域采樣理論的驗證。(n+1其它頻譜函數X(ejO)=FT[x(n)]在區間[0,2"]上等間隔采樣32322"Ok322"2"2"32163216323232 1616 x(n)和x(n)的波形，進行對比和分析，驗證總結頻域采樣理論。提示：頻域采樣用以下方法容易變程序實現。32也可以按照頻域采樣理論，先將信號x(n)以16為周期進行周期延拓，取其主值區(16點)，再對其進行16點DFT(FFT),得到的就是X(ej)在[0,2]的16點頻率域采樣X(k)。N5.實驗報告及要求a)運行程序打印要求顯示的圖形，。b)分析比較實驗結果，簡述由實驗得到的主要結論c)簡要回答思考題d)附上程序清單和有關曲線。10.2.2實驗程序清單1時域采樣理論的驗證程序清單%時域采樣理論驗證程序2aTp=64/1000;%觀察時間Tp=64微秒%Fs=1000;T=1/Fs;Fs=1000;T=1/Fs;M=Tp*Fs;n=0:M-1;A=;alph=pi*50*2^;omega=pi*50*2^;xnt=A*exp(-alph*n*T).*sin(omega*n*T);yn='xa(nT)';subplot(3,2,1);tstem(xnt,yn);%調用自編繪圖函數tstem繪制序列圖boxon;title('(a)Fs=1000Hz');k=0:M-1;fk=k/Tp;subplot(3,2,2);plot(fk,abs(Xk));title('(a)T*FT[xa(nT)],Fs=1000Hz');xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,*max(abs(Xk))])%=================================================z2頻域采樣理論的驗證程序清單%頻域采樣理論驗證程序M=27;N=32;n=0:M;xa=0:floor(M/2);xb=ceil(M/2)-1:-1:0;xn=[xa,xb];TFx16n=ifft(X16k,N/2);%16點IFFT[X16(k)]得到x16(n)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)]');200])k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),'.');boxontitlec16點頻域采樣n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,'.');boxonIDFT[X_1_6(k)]');xlabel('n');ylabel('x_1_6(n)');axis([0,32,0,20])k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),'.');boxontitlee32點頻域采樣n1=0:N-1;subplot(3,2,6);stem(n1,x32n,'.');boxonIDFT[X_3_2(k)]');xlabel('n');ylabel('x_3_2(n)');axis([0,32,0,20])10.2.3實驗程序運行結果見，采樣序列的頻譜的確是以采樣頻率為周期對模擬信號頻譜的周期延運行結果如圖10.3.3所示。該圖驗證了頻域采樣理論和頻域采樣定理。對信號x(n)的頻譜函數X(ejω)在[0，2π]上等間隔采樣N=16時，N點IDFT[X(k)]得到的序列正是原Nx(n)=IDFT[X(k)]=[x(n+iN)]R(n)NNNNi=由于N<M，所以發生了時域混疊失真，因此。x(n)與x(n)不相同，如圖N圖10.3.3(c)和(d)所示。當N=32時，如圖圖10.3.3(c)和(d)所示，由于N>M，頻域采樣定理，所以不存在時域混疊失真，因此。x(n)與x(n)N10.2.4簡答思考題RnNNi=NNN2NNNN2N10.3.1實驗指導學習用FFT對連續信號和時域離散信號進行譜分析的方法，了解可能出現的分析誤差及其原因，以便正確應用FFT。2.實驗原理用FFT對信號作頻譜分析是學習數字信號處理的重要內容。經常需要進行譜分析的信號是模擬信號和時域離散信號。對信號進行譜分析的重要問題是頻譜分辨率D和分析誤差。頻譜分辨率直接和FFT的變換區間N有關，因為FFT能夠實現的頻率分辨率是2/N，因此要求2/ND。時，得到的是離散譜，而信號(周期信號除外)是連續譜，只有當N較大時離散譜的包絡才能逼近于連續譜，因此N要適當選擇大一些。周期信號的頻譜是離散譜，只有用整數倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜。如果不知道信號周期，可以盡量選擇信號的觀察時間長一些。對模擬信號進行譜分析時，首先要按照采樣定理將其變成時域離散信號。如果是模擬周期信號，也應該選取整數倍周期的長度，經過采樣后形成周期序列，按照周期序列的譜分析進行。3．實驗步驟及內容(1)對以下序列進行譜分析。x(n)=R(n)14(n+1,0n32|0,其它n(4n,0n33|0,其它n(2)對以下周期序列進行譜分析。45分析。分別打印其幅頻特性曲線。并進行對比、分析和討論。(3)對模擬周期信號進行譜分析6選擇采樣頻率F=64Hz，變換區間N=16,32,64三種情況進行譜分析。s分別打印其幅頻特性，并進行分析和討論。4．思考題(1)對于周期序列，如果周期不知道，如何用FFT進行譜分析(2)如何選擇FFT的變換區間(包括非周期信號和周期信號)(3)當N=8時，x(n)和x(n)的幅頻特性會相同嗎為什么N=16呢235．實驗報告要求(1)完成各個實驗任務和要求。附上程序清單和有關曲線。(2)簡要回答思考題。10.3.2實驗程序清單clearall;closeall_%實驗內容(1)===================================================x1n=[ones(1,4)];%產生序列向量x1(n)=R4(n)M=8;xa=1:(M/2);xb=(M/2):-1:1;x2n=[xa,xb];%產生長度為8x3n=[xb,xa];X1k8=fft(x1n,8);X1k16=fft(x1n,16);X2k8=fft(x2n,8);X2k16=fft(x2n,16);X3k8=fft(x3n,8);X3k16=fft(x3n,16);%以下繪制幅頻特性曲線axis([0,2,0,*max(abs(X1k8))])FTaxis([0,2,0,*max(abs(X1k16))])figure(2)subplot(2,2,1);mstem(X2k8);%繪制8點DFT的幅頻特性圖title('(2a)8點DFT[x2(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,*max(abs(X2k8))])FTaxis([0,2,0,*max(abs(X2k16))])axis([0,2,0,*max(abs(X3k8))])FTaxis([0,2,0,*max(abs(X3k16))])%實驗內容(2)周期序列譜分析==================================Nn0:N-1;%FFT的變換區間N=8x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k8=fft(x4n);X5k8=fft(x5n);N=16;n=0:N-1;x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k16=fft(x4n);X5k16=fft(x5n);figure(3)axis([0,2,0,*max(abs(X4k8))])FTaxis([0,2,0,*max(abs(X4k16))])axis([0,2,0,*max(abs(X5k8))])FTaxis([0,2,0,*max(abs(X5k16))])%實驗內容(3)模擬周期信號譜分析===============================figure(4)Fs=64;T=1/Fs;N=16;n=0:N-1;%FFT的變換區間N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);%對X6k16=fftshift(X6k16);%將零頻率移到頻譜中心Tp=N*T;F=1/Tp;%頻率分辨率F(以零頻率為中心)subplotstemfkabsXk');boxon%繪制8點DFT的幅頻特性圖titleaDFTx6(nT)]|');xlabel('f(Hz)');ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0,*max(abs(X6k16))])N=32;n=0:N-1;%FFT的變換區間N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);%對X6k32=fftshift(X6k32);%將零頻率移到頻譜中心Tp=N*T;F=1/Tp;%頻率分辨率F(以零頻率為中心)subplotstemfkabsXk');boxon%繪制8點DFT的幅頻特性圖titlebDFTx6(nT)]|');xlabel('f(Hz)');ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0,*max(abs(X6k32))])N=64;n=0:N-1;%FFT的變換區間N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);%對X6k64=fftshift(X6k64);%將零頻率移到頻譜中心Tp=N*T;F=1/Tp;%頻率分辨率F(以零頻率為中心)subplotstemfkabsXk');boxon%繪制8點DFT的幅頻特性圖titleaDFTx6(nT)]|');xlabel('f(Hz)');ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0,*max(abs(X6k64))])10.3.3實驗程序運行結果程序運行結果分析討論：6請讀者注意，用DFT(或FFT)分析頻譜，繪制頻譜圖時，最好將X(k)的自變量k換算成對應的頻率，作為橫坐標便于觀察頻譜。kN為了便于讀取頻率值，最好關于π歸一化，即以/幾作為橫坐標。1、實驗內容(1)圖(1a)和(1b)說明x(n)=R(n)的8點DFT和16點DFT分別1328832的模相等，如圖(2a)和(3a)。但是，當N=16時，x(n)與x(n)32不滿足循環移位關系，所以圖(2b)和(3b)的模不同。2、實驗內容(2)，對周期序列譜分析4倍，得到正確的單一頻率正弦波的頻譜，僅在π處有1根單一譜線。如圖(4b)和(4b)所示。5期的整數倍，得到的頻譜不正確，如圖(5a)所示。N=16是其一個周期，得到正確的頻譜，僅在π和π處有2根單一譜線,如圖(5b)所3、實驗內容(3)，對模擬周期信號譜分析6123x(t)的周期為。采樣頻率F=64Hz=16f=8f=6.4f。變換區s123間N=16時，觀察時間Tp=16T=，不是x(t)的整數倍周期，所以所得頻譜6不正確，如圖(6a)所示。變換區間N=32,64時，觀察時間Tp=，1s，是x(t)的整數周期，所以所得頻譜正確，如圖(6b)和(6c)所示。圖6(1)用DFT(或FFT)對模擬信號分析頻譜時，最好將X(k)的自變量k換算成對應的模擬頻率fk，作為橫坐標繪圖，便于觀察頻譜。這樣，不管變換區間N取信號周期的幾倍，畫出的頻譜圖中有效離散諧波譜線所在的頻率值不變，如圖(6b)和(6c)所示。F11f=sk=k=k,k=0,1,2,,N1kNNTTp(2)本程序直接畫出采樣序列N點DFT的模值，實際上分析頻譜時最好畫出歸一化幅度譜，這樣就避免了幅度值隨變換區間N變化的缺點。本實驗程序這樣繪圖只要是為了驗證了用DFT對中期序列譜分析的理10.3.4簡答思考題思考題(1)和(2)的答案請讀者在教材3.節找，思考題(3)的答案在程序運行結果分析討論已經詳細回答。R10.4.1實驗指導(1)熟悉用雙線性變換法設計IIR數字濾波器的原理與方法；(2)學會調用MATLAB信號處理工具箱中濾波器設計函數(或濾波器設計分析工具fdatool)設計各種IIR數字濾波器，學會根據濾波需求確定濾波器指標參數。(3)掌握IIR數字濾波器的MATLAB實現方法。(3)通過觀察濾波器輸入輸出信號的時域波形及其頻譜，建立數字2．實驗原理設計IIR數字濾波器一般采用間接法(脈沖響應不變法和雙線性變換法)，應用最廣泛的是雙線性變換法。基本設計過程是：①先將給定的數字濾波器的指標轉換成過渡模擬濾波器的指標；②設計過渡模擬濾波器；③將過渡模擬濾波器系統函數轉換成數字濾波器的系統函數。MATLAB信號處理工具箱中的各種IIR數字濾波器設計函數都是采用雙線性變換法。第六章介紹的濾波器設計函數butter、cheby1、cheby2和橢圓模擬和數字濾波器。本實驗要求讀者調用如上函數直接設計IIR數本實驗的數字濾波器的MATLAB實現是指調用MATLAB信號處理工具箱y(n)。3.實驗內容及步驟(1)調用信號產生函數mstg產生由三路抑制載波調幅信號相加構成信號st，該函數還會自動繪圖顯示st的時域波形和幅頻特性曲但頻域是分離的，所以可以通過濾波的方法在頻域分離，這就是本實驗(2)要求將st中三路調幅信號分離，通過觀察st的幅頻特性曲線，分別確定可以分離st中三路抑制載波單頻調幅信號的三個濾波器 (低通濾波器、帶通濾波器、高通濾波器)的通帶截止頻率和阻帶截止頻率。要求濾波器的通帶最大衰減為,阻帶最小衰減為60dB。提示：抑制載波單頻調幅信號的數學表示式為1s(t)=cos(2ft)cos(2ft)=[cos(2(ff)t)+cos(2(f+f)t)]0c2c0c0其中，cos(2ft)稱為載波，f為載波頻率，cos(2ft)稱為單頻調制信cc0號，f為調制正弦波信號頻率，且滿足f>f。由上式可見，所謂抑制載0c0f+f和差頻ff，這2個頻率成分關于載波頻率f對稱。所以，1路c0c0cc中沒有載頻成分，故取名為抑制載波單頻調幅信號。容易看出，圖10.4.1中三路調幅信號的載波頻率分別為250Hz、500Hz、1000Hz。如果調制信號m(t)具有帶限連續頻譜，無直流成分，則s(t)=m(t)cos(2ft)就cc帶(上下邊帶)，在專業課通信原理中稱為雙邊帶抑制載波(DSB-SC)調幅信號,簡稱雙邊帶(DSB)信號。如果調制信號m(t)有直流成分，則s(t)=m(t)cos(2ft)就是一般的雙邊帶調幅信號。其頻譜圖是關于載波頻c率f對稱的2個邊帶(上下邊帶)，并包含載頻成分。c(3)編程序調用MATLAB濾波器設計函數ellipord和ellip分別設計這三個橢圓濾波器，并繪圖顯示其幅頻響應特性曲線。(4)調用濾波器實現函數filter，用三個濾波器分別對信號產生函123波形，觀察分離效果。functionst=mstg%st=mstg返回三路調幅信號相加形成的混合信號，長度N=1600Fs=10000;T=1/Fs;Tp=N*T;%采樣頻率Fs=10kHz，Tp為采樣時間t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;fc1=Fs/10;%第1路調幅信號的載波頻率fc1=1000Hz,fmfcfm00Hzfc2=Fs/20;fm2=fc2/10;fc3=Fs/40;fm3=fc3/10;xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t);%產生第1路調幅信號xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);%產生第2路調幅信號xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);%產生第3路調幅信號st=xt1+xt2+xt3;fxt=fft(st,N);%三路調幅信號相加%====以下為繪圖部分，繪制st的時域波形和幅頻特性曲線====================subplot(3,1,1)plot(t,st);grid;xlabel('t/s');ylabel('s(t)');axis([0,Tp/8,min(st),max(st)]);title('(a)s(t)的波形')subplot(3,1,2)stem(f,abs(fxt)/max(abs(fxt)),'.');grid;title('(b)s(t)的頻譜')axis([0,Fs/5,0,]);xlabel('f/Hz');ylabel('幅度')ynyn和y3(n)的時域波形和End6．思考題(1)請閱讀信號產生函數mstg，確定三路調幅信號的載波頻率和調mstg樣點數N=800，對st進行N點FFT可以(3)修改信號產生函數mstg，給每路調幅信號加入載波成分，產生調幅(AM)信號，重復本實驗，觀察AM信號與抑制載波調幅信號的時域波形及其頻譜的差別。0c7．實驗報告要求(1)簡述實驗目的及原理。(2)畫出實驗主程序框圖，打印程序清單。(3)繪制三個分離濾波器的損耗函數曲線。(4)繪制經過濾波分理出的三路調幅信號的時域波形。(5)簡要回答思考題。10.4.2濾波器參數及實驗程序清單1、濾波器參數選取觀察圖10.4.1可知，三路調幅信號的載波頻率分別為250Hz、500Hz、1000Hz。帶寬(也可以由信號產生函數mstg清單看出)分別為50Hz、100Hz、200Hz。所以，分離混合信號st中三路抑制載波單頻調幅信號的三個濾波器(低通濾波器、帶通濾波器、高通濾波器)的指標參對載波頻率為250Hz的條幅信號，可以用低通濾波器分離，其指標為帶截止頻率f=280Hz，通帶最大衰減a=0.1dBdB；ppss對載波頻率為500Hz的條幅信號，可以用帶通濾波器分離，其指標為帶截止頻率f=440Hz，f=560Hz，通帶最大衰減plpupzslsus對載波頻率為1000Hz的條幅信號，可以用高通濾波器分離，其指標為ppss說明：(1)為了使濾波器階數盡可能低，每個濾波器的邊界頻率選擇原則是盡量使濾波器過渡帶寬盡可能寬。(2)與信號產生函數mstg相同，采樣頻率Fs=10kHz。(3)為了濾波器階數最低，選用橢圓濾波器。按照圖10.4.2所示的程序框圖編寫的實驗程序為。2、實驗程序清單%IIR數字濾波器設計及軟件實現clearall;closeallFs=10000;T=1/Fs;%采樣頻率%調用信號產生函數mstg產生由三路抑制載波調幅信號相加構成的復st=mstg;%低通濾波器設計與實現=========================================fp=280;fs=450;wp=2*fp/Fs;ws=2*fs/Fs;rp=;rs=60;%DF指標(低通濾波器的通、Fy1t=filter(B,A,st);%濾波器軟件實現%低通濾波器設計與實現繪圖部分figure(2);subplot(3,1,1);myplot(B,A);%調用繪圖函數myplot繪制損耗函數曲線yt='y_1(t)';subplot(3,1,2);tplot(y1t,T,yt);%調用繪圖函數tplot繪制濾波器輸出波形%帶通濾波器設計與實現====================================================fpl=440;fpu=560;fsl=275;fsu=900;wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=;rs=60;wpy2t=filter(B,A,st);%濾波器軟件實現%帶通濾波器設計與實現繪圖部分(省略)1%高通濾波器設計與實現================================================fp=890;fs=600;wp=2*fp/Fs;ws=2*fs/Fs;rp=;rs=60;%DF指標(低通濾波器的通、wpy3t=filter(B,A,st);%濾波器軟件實現%高低通濾波器設計與實現繪圖部分(省略)10.4.3實驗程序運行結果實驗4程序運行結果如圖所示。由圖可見，三個分離濾波器指標參數選取正確，算耗函數曲線達到所給指標。分離出的三路信號y1(n)，y2(n)和y3(n)的波形是抑制載波的單頻調幅波。(a)低通濾波器損耗函數及其分離出的調幅信號y(t)b出的調幅信號y(t)2310.4.4簡要回答思考題思考題(1)已經在10.4.2節解答。思考題(3)很簡單，請讀者按照該題的提示修改程序，運行觀察。思考題(3)因為信號st是周期序列，譜分析時要求觀察時間為整數倍周期。所以，本題的一般解答方法是，先確定信號st的周期，在判煩，故采用下面的方法解答。10.5.1實驗指導(1)掌握用窗函數法設計FIR數字濾波器的原理和方法。(2)掌握用等波紋最佳逼近法設計FIR數字濾波器的原理和方法。(3)掌握FIR濾波器的快速卷積實現原理。(4)學會調用MATLAB函數設計與實現FIR濾波器。(1)認真復習第七章中用窗函數法和等波紋最佳逼近法設計FIR數(2)調用信號產生函數xtg產生具有加性噪聲的信號xt，并自動顯t(3)請設計低通濾波器，從高頻噪聲中提取xt中的單頻調幅信t確定濾波器指標參數。(4)根據濾波器指標選擇合適的窗函數，計算窗函數的長度N，調xt特性曲線、濾波器輸出信號的幅頻特性圖和時域波形圖。(4)重復(3)，濾波器指標不變，但改用等波紋最佳逼近法，調計方法設計的濾波器階數。zTFs根據圖10.6.1(b)和實驗要求，可選擇濾波器指標參數：通帶截止頻率fp=120Hz，阻帶截至頻率fs=150Hz，換算成數字頻率，通帶截止頻率=2fT=0.24，通帶最大衰為，阻帶截至頻率=2fT=0.3，阻帶ppss。FsFs=1000，用窗函數法或等波紋最佳逼近1、計算并繪圖顯示濾波器損耗End4.思考題(1)如果給定通帶截止頻率和阻帶截止頻率以及阻帶最小衰減,如何用窗函數法設計線性相位低通濾波器請寫出設計步驟.(2)如果要求用窗函數法設計帶通濾波器,且給定通帶上、下截止頻率為和，阻帶上、下截止頻率為和，試求理想帶通濾波器plpuslsu的截止頻率和。clcu(3)解釋為什么對同樣的技術指標，用等波紋最佳逼近法設計的濾波器階數低5.實驗報告要求(1)對兩種設計FIR濾波器的方法(窗函數法和等波紋最佳逼近法)進行分析比較，簡述其優缺點。(2)附程序清單、打印實驗內容要求繪圖顯示的曲線圖。(3)分析總結實驗結果。(4)簡要回答思考題。functionxt=xtg(N)%實驗五信號x(t)產生,并顯示信號的幅頻特性曲線%xt=xtg(N)產生一個長度為N,有加性高頻噪聲的單頻調幅信號xt,%載波頻率fc=Fs/10=100Hz,調制正弦波頻率f0=fc/10=10Hz.Fs=1000;T=1/Fs;Tp=N*T;t=0:T:(N-1)*T;fc=Fs/10;f0=fc/10;%載波頻率fc=Fs/10，單頻調制信號頻率為f0=Fc/10;mt=cos(2*pi*f0*t);ct=cos(2*pi*fc*t);xt=mt.*ct;nt=2*rand(1,N)-1;制信號mt，頻率為f0%相乘產生單頻調制信號xt%產生隨機噪聲nt%=======設計高通濾波器hn,用于濾除噪聲nt中的低頻成分,生成高通噪聲=======fp=150;fs=200;Rp=;As=70;fb=[fp,fs];m=[0,1];%濾波器指標f,m,devdev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];[n,fo,mo,W]=remezord(fb,m,dev,Fs);%確定remez函數所需參數hn=remez(n,fo,mo,W);yt=filter(hn,1,10*nt);%濾除隨機噪聲中低頻成分，生成高%================================================================xt=xt+yt;%噪聲加信號fst=fft(xt,N);k=0:N-1;f=k/Tp;subplot(3,1,1);plot(t,xt);grid;xlabel('t/s');ylabel('x(t)');axis([0,Tp/5,min(xt),max(xt)]);title('(a)信號加噪聲波形')subplot(3,1,2);plot(f,abs(fst)/max(abs(fst)));grid;title('(b)信號加噪聲的頻譜')axis([0,Fs/2,0,]);xlabel('f/Hz');ylabel('幅度')10.5.2濾波器參數及實驗程序清單1、濾波器參數選取根據10.5.1節實驗指導的提示③選擇濾波器指標參數：通帶截止頻fpHz頻率fs=150Hz。代入采樣頻率Fs=1000Hz，換算成數字頻率，通帶截止頻率=2fT=0.24，通帶最大衰為，阻帶截至頻pp率=2fT=0.3，阻帶最小衰為60dB。所以選取blackman窗函數。與ss信號產生函數xtg相同，采樣頻率Fs=1000Hz。按照圖10.5.2所示的程序框圖編寫的實驗程序為。2、實驗程序清單%FIR數字濾波器設計及軟件實現clearall;closeall;譜,=========N=1000;xt=xtg(N);fp=120;fs=150;Rp=;As=60;Fs=1000;%輸入給定指標%(1)用窗函數法設計濾波器wc=(fp+fs)/Fs;%理想低通濾波器截止頻率(關于pi歸一化)B=2*pi*(fs-fp)/Fs;%過渡帶寬度指標Nb=ceil(11*pi/B);%blackman窗的長度Nhn=fir1(Nb-1,wc,blackman(Nb));Hw=abs(fft(hn,1024));%求設計的濾波器頻率特性%以下為用窗函數法設計法的繪圖部分(濾波器損耗函數，濾波器輸%省略%(2)用等波紋最佳逼近法設計濾波器fb=[fp,fs];m=[1,0];%確定remezord函數所需參數f,m,devdev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)];[Ne,fo,mo,W]=remezord(fb,m,dev,Fs);%確定remez函數所需參數hn=remez(Ne,fo,mo,W);Hw=abs(fft(hn,1024));yet=fftfilt(hn,xt,N);%求設計的濾波器頻率特性%以下為用等波紋設計法的繪圖部分(濾波器損耗函數，濾波器輸出%省略10.5.3實驗程序運行結果用窗函數法設計濾波器，濾波器長度Nb=184。濾波器損耗函數和濾波器輸出yw(nT)分別如圖10.5.3(a)和(b)所示。用等波紋最佳逼近法設計濾波器，濾波器長度Ne=83。濾波器損耗函數和濾波器輸出ye(nT)分別如圖10.5.3(c)和(d)所示。兩種方法設計的濾波器都能有效地從噪聲中提取信號，但等波紋最佳逼近法設計的濾波器階數低得多，當然濾波實現的運算量以及時延也小得多，從圖10.5.3(b)和(d)可以直觀地看出時延差別。波器階數低10.5.4簡答思考題(1)用窗函數法設計線性相位低通濾波器的設計步驟教材中有詳細的介紹.：clcu=(+)/2,=(+)/2clslplcusupu(3)解釋為什么對同樣的技術指標，用等波紋最佳逼近法設計的濾①用窗函數法設計的濾波器，如果在阻帶截止頻率附近剛好滿足，則離開阻帶截止頻率越遠，阻帶衰減富裕量越大，即存在資源浪費；②幾種常用的典型窗函數的通帶最大衰減和阻帶最小衰減固定，且差別較大，又不能分別控制。所以設計的濾波器的通帶最大衰減和阻帶最小衰減通常都存在較大富裕。如本實驗所選的blackman窗函數，其阻③用等波紋最佳逼近法設計的濾波器，其通帶和阻帶均為等波紋特性，且通帶最大衰減和阻帶最小衰減可以分別控制，所以其指標均勻分布，沒有資源浪費，所以期階數低得多。實驗六數字信號處理在雙音多頻撥號系統中的應用10.6.1實驗指導雙音多頻(DualToneMultiFrequency,DTMF)信號是音頻電就代替了原有的用脈沖計數方式的撥號制式。這種雙音多頻信號制式不僅用在電話網絡中，還可以用于傳輸十進制數據的其它通信系統中，用于語音菜單進行操作。DTMF信號系統是一個典型的小型信號處理系統，它要用數字方法產生模擬信號并進行傳輸，其中還用到了D/A變換器；在接收端用A/D變換器將其轉換成數字信號，并進行數字信號處理與識別。為了系統的檢測速度并降低成本，還開發一種特殊的DFT算法，稱為戈澤爾(Goertzel)算法，這種算法既可以用硬件(專用芯片)實現，也可以用軟件實現。下面首先介紹雙音多頻信號的產生方法和檢測方法，包括戈澤爾算法，行模擬實驗。下面先介紹電話中的DTMF信號的組成。在電話中，數字0~9的中每一個都用兩個不同的單音頻傳輸，所用的8個頻率分成高頻帶和低頻帶兩組，低頻帶有四個頻率：Hz1209Hz,1336Hz,1477Hz和1633Hz.。每一個數字均由高、低頻帶中各一sin(2ft)+sin(2ft)表示，其中f=679Hz，f=1209Hz。這樣8個頻率形212示。表中最后一列在電話中暫時未用。表10.6.1雙頻撥號的頻率分配697Hz770Hz852Hz942Hz147*369#2580633HzABCDDTMF信號在電話中有兩種作用，一個是用撥號信號去控制交換機接通被叫的用戶電話機，另一個作用是控制電話機的各種動作，如播放留言、語音信箱等。2電話中的雙音多頻(DTMF)信號的產生與檢測(1)雙音多頻信號的產生假設時間連續的DTMF信號用x(t)=sin(2ft)+sin(2ft)表示，式中2f和f是按照表10.10.1選擇的兩個頻率，f代表低頻帶中的一個頻率，121f代表高頻帶中的一個頻率。顯然采用數字方法產生DTMF信號，方便而2信號進行采樣，采樣后得到時域離散信號為2形成上面序列的方法有兩種，即計算法和查表法。用計算法求正弦波的序列值容易，但實際中要占用一些計算時間，影響運行速度。查表法是預先將正弦波的各序列值計算出來，寄存在存儲器中，運行時只要按順序和一定的速度取出便可。這種方法要占用一定的存儲空間，但是速序列再送到D/A變換器和平滑濾波器，輸出便是連續時間的DTMF信號。DTMF信號通過電話線路送到交換機。(2)雙音多頻信號的檢測在接收端，要對收到的雙音多頻信號進行檢測，檢測兩個正弦波的頻率是多少，以判斷所對應的十進制數字或者符號。顯然這里仍然要用數字方法進行檢測，因此要將收到的時間連續DTMF信號經過A/D變換，變成數字信號進行檢測。檢測的方法有兩種，一種是用一組濾波器提取所一種是用DFT(FFT)對雙音多頻信號進行頻譜分析，由信號的幅度譜，判斷信號的兩個頻率，最后確定相應的數字或符號。當檢測的音頻數目直接計算DFT的一種線性濾波方法。這里略去Goertzel算法的介紹(請參考文獻[19])，可以直接調用MATLAB信號處理工具箱中戈澤爾算法的T模擬信號進行頻譜分析的問題。根據第三章用DFT對模擬信號進行譜分析的理論，確定三個參數：(1)采樣頻率F，(2)DFT的變換點數N，s(3)需要對信號的觀察時間的長度Tp。這三個參數不能隨意選取，要根據對信號頻譜分析的要求進行確定。這里對信號頻譜分析也有三個要求：(1)頻率分辨率，(2)譜分析的頻譜范圍，(3)檢測頻率的準觀察要檢測的8個頻率，相鄰間隔最小的是第一和第二個頻率，間F=73Hz。DFT的分辨率和對信號的觀察時間T有關，minpT=1/F=1/73=13.7ms。考慮到可靠性，留有富裕量，要求按鍵的時pmin2頻譜分析的頻率范圍要檢測的信號頻率范圍是697~1633Hz，但考慮到存在語音干擾，除了檢測這8個頻率外，還要檢測它們的二次倍頻的幅度大小，波形正常且干擾小的正弦波的二次倍頻是很小的，如果發現二次諧波很大，則不能確定這是DTMF信號。這樣頻譜分析的頻率范圍為697~3266Hz。按照采樣定理，最高頻率不能超過折疊頻率，即0.5F3622Hz，由此要求最s小的采樣頻率應為。因為數字電話總系統已經規定F＝8KHz，因此對頻s譜分析范圍的要求是一定滿足的。按照T=13.7ms，F＝8KHz，算出對pminsminpmins3檢測頻率的準確性DFT率是否準確的問題。序列的N點DFT列，截取周期序列的整數倍周期，進行DFT，其采樣點剛好在周期信號的不可能經過采樣得到周期序列，因此存在檢測頻率的準確性問題。kNkNk模擬域采樣點頻率為f=Fk/N(k=0,1,2,---,N-1)，希望選擇一個合ksNfkkks這樣雖然用幅度最大點檢測的頻率有誤差，但可以準確判斷所對應的DTMF頻率，即可以準確判斷所對應的數字或符號。經過分析研究認為N＝205是最好的。按照F＝8KHz，N＝205，算出8個頻率及其二次諧波對s88個基最近的DFT絕對誤二次諧對應最近的絕對誤頻整數k的差波的整數k差Hz值k值Hzk值值69718139435770201540398522217044394124188247120931241861133634267267147738295474163342326682通過以上分析，確定F＝8KHz，N＝205，T40ms。sp4DTMF信號的產生與識別仿真實驗Xgk=goertzel(xn,K)K是要求計算的DFT[xn]的頻點序號向量，用N表示xn的長度，則要求1K=[18，20，22，24，31，34，38，42]，K=[18，20，22，24，31，34，35，38，39，42，43，47，61，67，74，82]。Xgk是變換結果向量，其中存放的是由K指定的頻率點的DFT[x(n)]的XkDFTxnXgkiXKii=1,2,,length(K)。DTMF信號的產生與識別仿真實驗在MATLAB環境下進行，編寫仿真程序，運行程序，送入6位電話號碼，程序自動產生每一位號碼數字相應的DTMF信號，并送出雙頻聲音，再用DFT進行譜分析，顯示每一位號碼數字的DTMF信號的DFT幅度譜，安照幅度譜的最大值確定對應的頻率，再安照頻率確定每一位對應的號碼數字，最后輸出6位電話號碼。本實驗程序較復雜，所以將仿真程序提供給讀者，只要求讀者讀懂程序，直接運行程序仿真。程序名為exp6。程序分四段：第一段(2—7行)設置參數，并讀入6位電話號碼；第二段(9—20行)根據鍵入的6聲音；第三段(22—25行)對時域離散DTMF信號進行頻率檢測，畫出幅度譜；第四段(26—33行)根據幅度譜的兩個峰值，分別查找并確定輸入6位電話號碼。根據程序中的注釋很容易分析編程思想和處理算法。：%DTMF雙頻撥號信號的生成和檢測程序%clearall;clc;tm=[1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68];%DTMF信號代表的N=205;K=[18,20,22,24,31,34,38,42];f1=[697,770,852,941];f2=[1209,1336,1477,1633];TNr=0;forl=1:6;%行頻率向量%列頻率向量%接收端電話號碼初值為零d=fix(TN/10^(6-l));TN=TN-d*10^(6-l);forp=1:4;forq=1:4;iftm(p,q)==abs(d);break,end%檢測碼相符的iftm(p,q)==abs(d);break,end%檢測碼相符的行號pn=0:1023;%為了發聲，加長序列x=sin(2*pi*n*f1(p)/8000)+sin(2*pi*n*f2(q)/8000);%構成雙頻信號sound(x,8000);pause%接收檢測端的程序X=goertzel(x(1:205),K+1);val=abs(X);%發出聲音subplot(3,2,l);DFT(k)幅度axis([10500120])limit=80;%fors=5:8;ifval(s)>limit,break,endforr=1:4;ifval(r)>limit,break,endTNr=TNr+tm(r,s-4)*10^(6-l);%查找列號%查找行號disp('接收端檢測到的號碼為：')%顯示接收