數字信號處理實驗

1、數字信號處理實驗姓名：學號：班級： 專業(yè)：通信工程實驗 1 基 2-FFT 算法實現一、實驗目的掌握基 2-FFT 的原理及具體實現方法。編程實現基 2-FFT 算法。加深理解 FFT 算法的特點。二、實驗設備與環(huán)境計算機、MATLAB 軟件環(huán)境。三、實驗基礎理論FFT 是 DFT 的一種快速算法，能使 DFT 的計算大大簡化，運算時間縮短。FFT 利用了的三個固有特性，即對稱性、周期性和可約性，將長序列的 DFT 分解為短序列的 DFT，合并了DFT 運算中的某些項，從而減少了 DFT 的運算量。FFT 算法基本上可以分為兩大類，即按時間抽取法和按頻率抽取法。在實現 FFT 算法時，要重點考

2、慮兩個問題，注意數據的讀取和存儲：（1）輸入輸出的排序；（2）碟形運算的實現。按時間抽取算法中輸入反序輸出順序，按頻率抽取算法中輸入順序輸出反序；運算過程中的每一級都由 N/2 個碟形運算構成，每一個碟形運算單元中，兩個節(jié)點變量運算后得到的結果為下一列相同位置的節(jié)點變量，而和其他節(jié)點變量無關，可以采用原位運算，節(jié)省存儲單元。另外，碟形運算中的復系數可以存儲為能及時查閱的系數表，這樣可以節(jié)約計算量，但是需要 N/2 個復數存儲器。MATLAB 中提供了用于計算 FFT 的函數 fft，可將實驗中所得到的結果與利用 MATLAB 中fft 函數計算的結果相比較，以此驗證結果的正確性。四、實驗內容編

3、程實現序列長度為N=8 的按時間抽取的基 2-FFT 算法。給定一個 8 點序列，采用編寫的程序計算其 DFT，并與 MATLAB 中 fft 函數計算的結果相比較，以驗證結果的正確性。編程實現序列長度為 N=8 的按頻率抽取的基 2-FFT 算法。給定一個 8 點序列，采用編寫的程序計算其 DFT，并與 MATLAB 中 fft 函數計算的結果相比較，以驗證結果的正確性。將上述 FFT 程序推廣到序列長度為N=2v 的情況，要求利用原位運算。五、實驗代碼及分析實驗代碼： x=1+2j 0.5+3j 5+4j -2+3j 6-3j 5+1j 9 -1-2j;n=bin2dec(fliplr(d

4、ec2bin(1:8-1,3)+1;x0=x;for l=1:3for m=1:2(l-1) for k=1:2(3-l)t=x0(k+(m-1)*2(4-l)+2(3-l);x0(k+(m-1)*2(4-l)+2(3-l)=(x0(k+(m-1)*2(4-l)-t)*exp(-2*1j*pi*(k-1)/(8/2(l-1); x0(k+(m-1)*2(4-l)=x0(k+(m-1)*2(4-l)+t;endendendy1=fft(x)for i=1:8y0(i)=x0(n(i);endy0y1 =Columns 1 through 523.5000 + 8.0000i1.4749 +10.

5、7678i-4.0000 -13.5000i-1.5754 + 7.0104i18.5000- 2.0000iColumns 6 through 8-3.4749 + 7.2322i -10.0000 + 3.5000i -16.4246 - 5.0104iy0 =Columns 1 through 523.5000 + 8.0000i1.4749 +10.7678i-4.0000 -13.5000i-1.5754 + 7.0104i18.5000- 2.0000iColumns 6 through 8-3.4749 + 7.2322i -10.0000 + 3.5000i -16.4246

6、- 5.0104i分析：y0 和y1 兩者結果完全相同，說明結果正確。實驗代碼：x=1+2j 0.5+3j 5+4j -2+3j 6-3j 5+1j 9 -1-2j;n=bin2dec(fliplr(dec2bin(1:8-1,3)+1;x0=x;for l=1:3for m=1:2(l-1) for k=1:2(3-l)t=x0(k+(m-1)*2(4-l)+2(3-l);x0(k+(m-1)*2(4-l)+2(3-l)=(x0(k+(m-1)*2(4-l)-t)*exp(-2*1j*pi*(k-1)/(8/2(l-1); x0(k+(m-1)*2(4-l)=x0(k+(m-1)*2(4-l

7、)+t;end endendy1=fft(x)for i=1:8y0(i)=x0(n(i);endy0y1 =Columns 1 through 523.5000 + 8.0000i1.4749 +10.7678i-4.0000 -13.5000i-1.5754 + 7.0104i18.5000- 2.0000iColumns 6 through 8-3.4749 + 7.2322i -10.0000 + 3.5000i -16.4246 - 5.0104iy0 =Columns 1 through 523.5000 + 8.0000i1.4749 +10.7678i-4.0000 -13.5

8、000i-1.5754 + 7.0104i18.5000- 2.0000iColumns 6 through 8-3.4749 + 7.2322i -10.0000 + 3.5000i -16.4246 - 5.0104i分析：兩種算法的結果相同，說明結果相同。實驗代碼：按時間抽取 xn= 1,4,7,5,9,2,3,8,6,5,2,7,6,4,13,2; M=nextpow2(length(xn); N=2M; for m=0:N/2-1WN(m+1)=exp(-j*2*pi/N)m; endDF1=xn,zeros(1,N-length(xn);H=0;for I=0:N-1;if IK=

9、N/2;while H=K;H=H-K;K=K/2;endH=H+K;endfor G=1:M; F=2(G-1); for S=0:F-1;P=2(M-G)*S; for K=S:2G:N-2;T=DF1(K+1)+DF1(K+F+1)*WN(P+1); DF1(K+F+1)=DF1(K+1)-DF1(K+F+1)*WN(P+1); DF1(K+1)=T;end endend DF1 DF1 =Columns 1 through 1084.0000-5.0057-7.8833-5.0141-3.0131-5.2100-7.9867-5.006010.0000-4.5526Columns 11

10、 through 16-8.1429-4.9962-2.9869-5.2093-7.9872-5.0060F1=fft(xn) F1 =Columns 1 through 1084.0000-5.0057-7.8833-5.0141-3.0131-5.2100-7.9867-5.006010.0000-4.5526Columns 11 through 16-8.1429-4.9962-2.9869-5.2093-7.9872-5.0060按頻率抽取 xn= 1,4,7,5,9,2,3,8,6,5,2,7,6,4,13,2;M=nextpow2(length(xn);N=2M; M=log2(N

11、); for k1=0:M-1 D=2k1;E=N/2k1; F=N/2(k1+1); G=N/2(k1+1)-1;Wn=exp(-j*2*pi/E); for g=1:DH1=(g-1)*E;H2=(g-1)*E+F;for r=0:G;k=r+1; xn(k+H1)=xn(k+H1)+xn(k+H2);xn(k+H2)=xn(k+H1)-xn(k+H2)-xn(k+H2)*Wnr;end endendn1=fliplr(dec2bin(0:N-1);n2=bin2dec(n1);for i=1:NXk(i)=xn(n2(i)+1); end XkXk =Columns 1 through

12、1084.0000-5.0057-7.8833-5.0141-3.0131-5.2100-7.9867-5.006010.0000-4.5526Columns 11 through 16-8.1429-4.9962-2.9869-5.2093-7.9872-5.0060F2=fft(xn) F2 =Columns 1 through 1084.0000-5.0057-7.8833-5.0141-3.0131-5.2100-7.9867-5.006010.0000 -4.5526Columns 11 through 16-8.1429-4.9962-2.9869-5.2093-7.9872-5.

13、0060六、 收獲與感悟通過本次數字信號處理實驗，學習了基 2-FFT 算法的實現，通過對 MATLAB 中 fft 函數的運行結果與所編寫程序結果的對比，了解了按時間抽取算法和按頻率抽取算法的異同，按時間抽取算法，兩種方式的結果相同，而按頻率抽取算法，兩種方式的結果不同。并且知道了碟形圖在 MATLAB 中的實現方法實驗 2 利用 DFT 分析信號頻譜一、實驗目的加深對 DFT 原理的理解。應用 DFT 分析信號的頻譜。深刻理解利用 DFT 分析信號頻譜的原理，分析實現過程中出現的現象及解決方法。二、實驗設備與環(huán)境計算機、MATLAB 軟件環(huán)境。三、實驗基礎理論DFT 和 DTFT 的關系有

14、限長序列x(n)(0 n N 1)的離散時間傅里葉變換X()在頻率區(qū)間（0 2） 的 N 個等間隔分布的點k = 2k/N(0 k N 1)上的N 個取樣值可以由下式表示：由上式可知，序列 x(n)的 N 點 DFTX(k),實際上就是 x(n)序列的 DTFT 在 N 個等間隔頻率點 k = 2k/N(0kN 1)上樣本 X(k)。利用 DFT 求 DTFT方法 1：其中(x)為內插函數方法 2：然而在實際 MATLAB 計算中，上述插值運算不見得是最好的辦法。由于 DFT是 DTFT 的取樣值，其相鄰兩個頻率樣本點的間距為 2/N，所以如果我們增加數據的長度N，使得到的 DFT 譜線就更加

15、精細，其包絡就越接近 DTFT 的結果，這樣就可以利用 DFT來近似計算DTFT，如果沒有更多的數據，可以通過補零來增加數據長度。利用 DFT 分析連續(xù)時間信號的頻譜采用計算機分析連續(xù)時間信號的頻譜，第一部就是把連續(xù)時間信號離散化。這里需要兩個操作：一是采樣，二是截斷。對于連續(xù)時間非周期信號()，按采樣間隔T 進行采樣，截取長度為 M，那么對()進行 N 點頻域采樣，得到因此，可以將利用 DFT 分析連續(xù)非周期信號頻譜的步驟歸納如下：（1）確定時域采樣間隔 T，得到離散序列 x(n)；（2）確定截取長度 M，得到 M 點離散序列() = ()(),這里()為窗函數。確定頻域采樣點數 N，要求

16、NM。利用 FFT 計算離散序列的 N 點 DFT，得到()。（5）根據式由()計算()采樣點的近似值。采用上述方法計算()的頻譜，需要注意如下三個問題：頻譜混疊。如果不滿足采樣定理的條件，頻譜會出現混疊誤差。對于頻譜無限寬的信號，應考慮覆蓋大部分主要頻率分量的范圍。柵欄效應和頻譜分辨率。使用 DFT 計算頻譜，得到的結果只是 N 個頻譜樣本值， 樣本值之間的頻譜是未知的，像通過一個柵欄觀察頻譜，稱為“柵欄效應”。頻譜分辨率與記錄長度成反比，要提高頻譜分辨率，就要增加記錄時間。頻譜泄露。對信號階段會把窗函數的頻譜引入信號頻譜，造成頻譜泄露，解決這個問題的主要辦法是采用旁瓣小的窗函數，頻譜泄露和

17、窗函數均會引起誤差。因此，要合理選取采樣間隔和截取長度，必要時還考慮加適當的窗。對于連續(xù)時間周期信號，我們在采用計算機進行計算時，也總是要進行截斷，序列總是有限長的，仍然可以采用上述方法近似計算?？赡苡玫降?MATLAB 函數與代碼實驗中 DFT 運算可采用 MATLAB 中提供的函數 fft 來實現。DTFT 可以利用 MATLAB 矩陣運算的方法進行計算X() = = , , , 1 2 四、實驗內容=1121.已知 x(n)=2,-1,1,1,完成如下要求：計算器 DTFT，并畫出*, +區(qū)間的波形。計算 4 點 DFT，并把結果顯示在(1)所畫的圖形中。對 x(n)補零，計算 64 點

18、 DFT，并顯示結果。根據實驗結果，分析是否可以由 DFT 計算 DTFT，如果可以，如何實現？考察序列0 n10 時，用 DFT 估計 x(n)的頻譜；將 x(n)補零加長到長度為 100 點序列用 DFT估計 x(n)的頻譜，要求畫出相應波形。0 n 100 時，用 DFT 估計 x(n)的頻譜，并畫出波形。根據實驗結果，分析怎樣提高頻譜分辨率。已知信號 ，其中1 =1, 2 = 2, 3 = 3,從 x(n)的表達式可以看出，它包含三個頻率的正弦波，但是，從其時域波形來看，似乎是一個正弦信號，利用 DFT 做頻譜分析，確定適合的參數，使得到的的頻率分辨率復合需要。利用 DFT 近似分析連

19、續(xù)時間信號 的頻譜（幅度譜）。分析采用不同的采樣間隔和截取長度進行計算的結果，并最終確定適合的參數。五、實驗結果及分析1.（1）實驗代碼：x=2 -1 1 1;n=0:3;w=-pi:0.01*pi:pi;X=x*exp(-j*n*w);subplot(211);plot(w,abs(X);xlabel(Omega/pi);title(Magnitude);axis tightsubplot(212);plot(w,angle(X)/pi);xlabel(Omega/pi);title(Phase);axis tight波形圖：Magnitude321-3-2-10123/ Phase0.50

20、-0.5-3-2-10123/實驗代碼：fft(x);k=0:3;subplot(211);hold onstem(k,abs(ans);subplot(212);hold onstem(k,angle(ans)/pi);計算結果圖示如下：Magnitude321-3-2-10123/ Phase0.50-0.5-3-2-10123/實驗代碼：fft(x,64);k=0:63;subplot(211);stem(k,abs(ans);xlabel(k);title(Magnitude);axis tightsubplot(212);stem(k,angle(ans);xlabel(k);tit

21、le(Phase);axis tight波形圖如下：Magnitude32100102030405060k Phase10-10102030405060k分析：可以用 DFT 計算 DTFT。實現方法：利用對 x(n)補零，提高采樣密度，當補零的個數足夠多時，即可用 DFT 的結果近似 DTFT。2.（1）當 0 n 10時，實驗代碼:n=0:10;x=cos(0.48*pi*n)+cos(0.52*pi*n);X=fft(x);subplot(211);k=0:10;stem(k,abs(X);xlabel(k);title(Magnitude);axis tightsubplot(212)

22、;stem(k,angle(X);xlabel(k);title(Phase);axis tightx(n)的頻譜為：Magnitude86420012345678910k Phase10-1012345678910k當 x(n)的長度變?yōu)?100 時，實驗代碼：X1=fft(x,100);subplot(211);k=0:99;stem(k,abs(X1);xlabel(k);title(Magnitude);axis tightsubplot(212);stem(k,angle(X1);xlabel(k);title(Phase);axis tightx(n)的頻譜為：Magnitude1

23、0500102030405060708090k Phase10-10102030405060708090k（2）當0 n 100時，實驗代碼:n=0:100;x=cos(0.48*pi*n)+cos(0.52*pi*n);X=fft(x);subplot(211);k=0:100;stem(k,abs(X);xlabel(k);title(Magnitude);axis tightsubplot(212);stem(k,angle(X);xlabel(k);title(Phase);axis tight頻譜波形：Magnitude402000102030405060708090100k Pha

24、se210-1-20102030405060708090100k（3）分析：由頻譜波形可以看出，通過增加補零的個數可以提高頻譜分辨率，但是不能提高分辨力，提高分辨力需要通過延長所選信號的長度來實現。實驗代碼： t=0:0.01:1f1=1;f2=2;f3=3;x=0.15*sin(2*pi*f1.*t)+sin(2*pi*f2.*t)-0.1*sin(2*pi*f3.*t);plot(t,x);grid;時域波形圖為：1.510.50-0.5-1-1.500.10.20.30.40.50.60.70.80.91 n=0:2000; x=0.15*sin(2*pi*f1.*n)+sin(2*pi

25、*f2.*n)-0.1*sin(2*pi*f3.*n); X=fft(x); figure(2); stem(n,X,filled);-10 x 10420-2-4-6-8-1002004006008001000 1200 1400 1600 1800 2000分析：由圖中頻譜可知，存在 3 個頻譜分量。實驗代碼：（1）采樣區(qū)間：0，100，采樣間隔為 1 時： n=0:100; x=exp(-0.1*n); X1=fft(x); k=0:100; stem(k,abs(X1); xlabel(k);波形圖：1210864200102030405060708090100k（2）采樣區(qū)間：0，1

26、00，采樣間隔為 2 時： n=0:2:100; x=exp(-0.1*n); X2=fft(x); k=0:50; stem(k,abs(X2); xlabel(k);波形圖：654321005101520253035404550k采樣區(qū)間：0，50，采樣間隔為 1 時： n=0:50; x=exp(-0.1*n); X3=fft(x); k=0:50; stem(k,abs(X3); xlabel(k);波形圖：12108642005101520253035404550k采樣區(qū)間：0，50，采樣間隔為 2 時： n=0:2:50; x=exp(-0.1*n); X4=fft(x); k=0

27、:25; stem(k,abs(X4); xlabel(k);波形圖：65432100510152025k（5）采樣區(qū)間：0，150，采樣間隔為 1 時： n=0:150; x=exp(-0.1*n); X5=fft(x); k=0:150; stem(k,abs(X5); xlabel(k);波形圖：121086420050100150k（6）采樣區(qū)間：0，150，采樣間隔為 2 時： n=0:2:150; x=exp(-0.1*n); X5=fft(x); k=0:75; stem(k,abs(X5); xlabel(k);波形圖：654321001020304050607080k分析：由

28、上述對比發(fā)現，采樣區(qū)間取0，100，采樣間隔取 2 時，結果最為接近 x(t)的頻譜。六、收獲及體會通過本次實驗，學習到了如何用 DFT 分析頻譜，知道了如何提高頻譜的分辨力和分辨率，理解力二者的區(qū)別，并且能夠利用 DFT 來近似DTFT 的結果。通過觀察頻譜波形，能夠看到頻譜泄漏現象，使得在時域無法觀察到的波形，能夠在頻域有所體現。不同的采樣間隔和信號截取長度，所得到的 DFT 結果也有所不同，需要選擇合適的參數，才能得到完整清晰的頻譜。實驗中加深了對 Matlab 進行數字信號分析的理解，能夠通過更加直觀的方式，了解到課本上的知識，使得對課本知識的學習有了更深刻的體會。實驗 5 脈沖響應不

29、變法設計 IIR 數字濾波器一、實驗目的掌握利用脈沖響應不變法設計 IIR 數字濾波器的原理及具體方法。加深理解數字濾波器和模擬濾波器之間的技術指標轉化。掌握脈沖響應不變法設計 IIR 數字濾波器的優(yōu)缺點及適用范圍。二、實驗上機環(huán)境計算機、MATLAB 軟件環(huán)境三、實驗原理1、基本原理從時域響應出發(fā)，使數字濾波器的單位脈沖響應 h(n)模仿模擬濾波器的單位沖激響應(),h(n)等于()的取樣值。2、變換方法思路：將H a(s)進行部分分式展開對 H a(s)進行拉式變換對ha (t )時域采樣得到h(n)對h(n)進行 z 變換3.設計步驟確定數字濾波器的性能指標p ,st , Rp , As

30、 。將數字濾波器頻率指標轉換成響應的模擬濾波器頻率指標 ppT ststT根據指標 p ， st ， Rp 和 As 設計模擬濾波器 H a (s) 。將 Ha (s) 展成部分分式形式NH a (s)=k 1Ak。s pkk把模擬極點 p 轉換成數字極點e pkT ，得到數字濾波器NH (z) Akk 1 1 e pkT z1可見 H a (s) 至 H(z)間的變換關系為1s sk11 eskT z1zz eskT在 MATLAB 中有兩種方法可以實現實現上述變換。方法 1：利用 residue 函數和 residuez 函數實現脈沖響應不變變換法，實用方法如下：實現多項式形式和部分分式形

31、式之間的轉換。實現多項式形式和部分分式形式之間的轉換。方法 2：matlab 中提供了 impinvar 函數采用脈沖響應不變法實現模擬濾波器到數字濾波器的變換，其使用如下：bz,az=impinvar(b,a,fs)采用脈沖響應不變法將模擬濾波器系統(tǒng)函數的系數向量 b 和 a 變換成為數字濾波器系統(tǒng)函數的系數向量 bz 和 az，fs 為采樣頻率（默認為 1）。bz,az=impinvar(b,a) 采樣頻率默認為 1 的情況下，采用脈沖響應不變法將模擬濾波器變換為數字濾波器。四、實驗內容設采樣頻率為 fs =4kHz，采用脈沖響應不變法設計一個三階巴特沃斯數字低通濾波器，其 3dB 截止頻

32、率為 fc=1kHz。設采樣頻率為 fs=10kHz，設計數字低通濾波器，滿足如下指標通帶截止頻率：fp=1kHz，通帶波動：Rp=1dB阻帶截止頻率：fst=1.5kHz，阻帶衰減：As=15dB要求分別采用巴特沃斯、切比雪夫 I 型、切比雪夫 II 型和橢圓模擬原型濾波器及脈沖響應不變法進行設計。結合實驗結果，分別討論采用上述設計的數字濾波器是否都能滿足給定指標要求，分析脈沖響應不變法設計 IIR 數字濾波器的優(yōu)缺點及適用范圍。五、實驗結果及分析1.實驗代碼：fs=4000;fc=1000;Wc=fc*2*pi; b=Wc3;a=1 2*Wc 2*Wc2 Wc3;bz,az=impinva

33、r(b,a,fs) bz =0.0000az =0.58130.21141.0000-0.39840.2475-0.0432w=0:500*pi/500; H,w=freqz(bz,az);subplot(221);plot(w/pi,abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(223);plot(w/pi,angle(H)/p

34、i);grid on;xlabel(omega(pi),ylabel(Phase of H(ejomega);grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd);grid on;xlabel(omega(pi),ylabel(Group delay); G=tf(bz,az,0.1,variable,z-1) G =0.5813 z-1 + 0.2114 z-2- 1 - 0.3984 z-1 + 0.2475 z-2 - 0.04321 z-3波形圖如下：1|H(ej)|0.80.60.40.200.51()0|H(ej)|-5-10-1500

35、.51()1Phase of H(ej)0.50-0.5-100.51()2Group delay10-100.51()2.（1）巴特沃斯濾波器實驗代碼： fs=10000; fp=1000; fst=1500; Rp=1; As=15; Ap=Rp; Wp=fp*2*pi; Wst=fst*2*pi; N=ceil(log10(10(Ap/10)-1)/(10(As/10)-1)/(2*log10(Wp/Wst); Wc=Wp/(10(Ap/10)-1)(1/2/N); b,a=butter(N,Wc,s); bz,az=impinvar(b,a,fs) bz =-0.00000.00060

36、.01010.01610.00410.00010az =1.0000-3.36355.0684-4.27592.1066-0.57060.0661 w=0:500*pi/500; H,w=freqz(bz,az); subplot(2,2,1);plot(w/pi,abs(H); grid on;xlabel(omega(pi);ylabel(|H(ejomega)|); subplot(2,2,2);plot(w/pi,20*log10(abs(H)/max(abs(H); grid on;xlabel(omega(pi);ylabel(|H(ejomega)|); subplot(2,2,

37、3);plot(w/pi,angle(H)/pi); grid on;xlabel(omega(pi),ylabel(Phase of H(ejomega); grd=grpdelay(bz,az,w); subplot(2,2,4);plot(w/pi,grd); grid on;xlabel(omega(pi),ylabel(Group delay); G=tf(bz,az,0.1,variable,z-1) G =z-5-5.457e-15 + 0.000631 z-1 + 0.0101 z-2 + 0.01614 z-3 + 0.004101 z-4 + 0.0001033-1 - 3

38、.364 z-1 + 5.068 z-2 - 4.276 z-3 + 2.107 z-4 - 0.5706 z-5 + 0.06607 z-6Sample time: 0.1 seconds Discrete-time transfer function. 波形圖：1|H(ej)|0.5000.51()0|H(ej)|-20-40-60-8000.51()1Phase of H(ej)0.50-0.5-100.51()10Group delay5000.51()由圖線可知，結果滿足設計指標要求。切比雪夫 I 型實驗代碼：fs=10000; fp=1000; fst=1500; Rp=1; As

39、=15; Ap=Rp; epc=sqrt(10(Ap/10)-1); Wp=fp*2*pi; Wc=Wp; Wst=fst*2*pi; N=ceil(acosh(sqrt(10(0.1*As)-1)/epc)/acosh(Wst/Wc); b,a=cheby1(N,Rp,Wp,s); bz,az=impinvar(b,a,fs) bz =0.00000.00540.01810.00400az =1.0000-3.05913.8323-2.29190.5495 w=0:500*pi/500; H,w=freqz(bz,az); subplot(2,2,1);plot(w/pi,abs(H); g

40、rid on;xlabel(omega(pi);ylabel(|H(ejomega)|); subplot(2,2,2);plot(w/pi,20*log10(abs(H)/max(abs(H); grid on;xlabel(omega(pi);ylabel(|H(ejomega)|); subplot(2,2,3);plot(w/pi,angle(H)/pi); grid on;xlabel(omega(pi),ylabel(Phase of H(ejomega); grd=grpdelay(bz,az,w); subplot(2,2,4);plot(w/pi,grd); grid on;

41、xlabel(omega(pi),ylabel(Group delay); G=tf(bz,az,0.1,variable,z-1) G =3.411e-17 + 0.005373 z-1 + 0.0181 z-2 + 0.003985 z-3-1 - 3.059 z-1 + 3.832 z-2 - 2.292 z-3 + 0.5495 z-4Sample time: 0.1 seconds Discrete-time transfer function. 波形圖：1.5|H(ej)|10.5000.51()0|H(ej)|-20-40-60-8000.51()1Phase of H(ej)0

42、.50-0.5-100.51()15Group delay105000.51()由圖線可知，滿足設計指標要求。切比雪夫 II 型實驗代碼： fs=10000; fp=1000; fst=1500; Rp=1; As=15; epc=sqrt(10(Rp/10)-1); Wp=fp*2*pi; Wst=fst*2*pi; Wc=Wst; N=ceil(acosh(sqrt(10(0.1*As)-1)/epc)/acosh(Wc/Wp); b,a=cheby2(N,As,Wst,s); bz,az=impinvar(b,a,fs) bz =-0.42140.9724-0.81190.42050.0

43、000az =1.0000-1.66581.4289-0.51930.0935 w=0:500*pi/500; H,w=freqz(bz,az);subplot(2,2,1);plot(w/pi,abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2,2,2);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2,2,3);plot(w/pi,angle(H)/pi);grid on;xlabe

44、l(omega(pi),ylabel(Phase of H(ejomega);grd=grpdelay(bz,az,w);subplot(2,2,4);plot(w/pi,grd);grid on;xlabel(omega(pi),ylabel(Group delay);G=tf(bz,az,0.1,variable,z-1) G =-0.4214 + 0.9724 z-1 - 0.8119 z-2 + 0.4205 z-3 + 1.663e-06 z-4-1 - 1.666 z-1 + 1.429 z-2 - 0.5193 z-3 + 0.09352 z-4Sample time: 0.1

45、seconds Discrete-time transfer function.波形圖:1.50|H(ej)|H(ej)|1-50.5000.51()-1000.51()1Phase of H(ej)0.50-0.5-100.51()5Group delay0-500.51()由圖線可知，阻帶衰減小于 15dB，不滿足設計要求。橢圓模擬原型濾波器實驗代碼： fs=10000; fp=1000; fst=1500; Rp=1; As=15; epc=sqrt(10(Rp/10)-1); Wp=fp*2*pi; Wst=fst*2*pi; Wc=Wp; A=10(As/20); k1=epc/(s

46、qrt(A2-1); k=Wp/Wst; N=ceil(ellipke(k)*ellipke(sqrt(1-k12)/ellipke(k1)/ellipke(sqrt(1-k2); bz,az=ellip(N,Rp,As,fp/fs*pi) bz =0.18550.04960.04960.1855az =1.0000-1.41071.2357-0.3549 w=0:500*pi/500; H,w=freqz(bz,az);subplot(2,2,1);plot(w/pi,abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2

47、,2,2);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2,2,3);plot(w/pi,angle(H)/pi);grid on;xlabel(omega(pi),ylabel(Phase of H(ejomega);grd=grpdelay(bz,az,w);subplot(2,2,4);plot(w/pi,grd);grid on;xlabel(omega(pi),ylabel(Group delay);G=tf(bz,az,0.1,variable,

48、z-1) G =0.1855 + 0.04956 z-1 + 0.04956 z-2 + 0.1855 z-3-1 - 1.411 z-1 + 1.236 z-2 - 0.3549 z-3Sample time: 0.1 seconds Discrete-time transfer function. 波形圖：1|H(ej)|0.5000.51()0|H(ej)|-20-40-60-8000.51()1Phase of H(ej)0.50-0.5-100.51()15Group delay105000.51()由圖線可知，滿足設計指標要求。脈沖相應不變法： 實驗代碼： fs=10000; fp

49、=1000; fst=1500; Rp=1; As=15; Ap=Rp; Wp=fp*2*pi; Wst=fst*2*pi; N=ceil(log10(10(Ap/10)-1)/(10(As/10)-1)/(2*log10(Wp/Wst); Wc=Wp/(10(Ap/10)-1)(1/2/N); b=Wc3; a=1 2*Wc 2*Wc2 Wc3; bz,az=impinvar(b,a,fs) bz =00.10570.0663az =1.0000-1.64911.0663-0.2450 w=0:500*pi/500; H,w=freqz(bz,az);subplot(2,2,1);plot(

50、w/pi,abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2,2,2);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2,2,3);plot(w/pi,angle(H)/pi);grid on;xlabel(omega(pi),ylabel(Phase of H(ejomega);grd=grpdelay(bz,az,w);subplot(2,2,4);plot(w/pi,grd);gr

51、id on;xlabel(omega(pi),ylabel(Group delay);G=tf(bz,az,0.1,variable,z-1) G =0.1057 z-1 + 0.06633 z-2-1 - 1.649 z-1 + 1.066 z-2 - 0.245 z-3Sample time: 0.1 secondsDiscrete-time transfer function.波形圖：1|H(ej)|0.5000.51()0|H(ej)|-20-40-6000.51()1Phase of H(ej)0.50-0.5-100.51()4Group delay20-200.51()由圖線可知

52、，滿足設計指標要求。分析：脈沖響應不變法設計 IIR 數字濾波器的優(yōu)缺點：優(yōu)點：脈沖響應不變法的頻率坐標的變換是線性的。因此，如果模擬濾波器的的頻率響應是限帶于折疊頻率以內的話，則通過變換后所得的數字濾波器的頻率響應可以不失真的反映原響應與頻率的關系。缺點：脈沖響應不變法最大的缺點是有頻譜的周期延拓效應。其次，脈沖響應不變法所得的數字濾波器會出現頻率的混疊失真。第三，數字濾波器在 T 很?。ɑ蛘哒f取樣頻率很高）時，可能具有不希望的高增益。適用范圍：脈沖響應不變法設計 IIR 數字濾波器只能用于限帶的頻率響應，如衰減特性很好的低通或帶通。六、收獲與體會通過本次實驗，學會了利用用脈沖響應不變法設計

53、 IIR 數字濾波器的原理及具體方法， 通過比較不同類型的數字濾波器，了解了它們各自的優(yōu)缺點以及適用范圍，在相同的技術指標下，不同濾波器的效果也有較大差異，不是所有類型都可以滿足要求。了解了如何將所給的技術指標，轉化為設計過程中的參數，包括模擬域和數字域的轉化，加深了對課本內容的理解與認識，進一步掌握了 IIR 數字濾波器的設計要點。實驗 7窗函數法設計 FIR 數字濾波器一、實驗目的掌握窗函數法設計 FIR 數字濾波器的原理和具體方法二、實驗設備與環(huán)境計算機、Matlab 軟件環(huán)境三、實驗基礎理論1、基本原理窗函數設計法的基本思想為，首先選擇一個適當的理想的濾波器H()，然后用窗函數截取它的

54、單位脈沖響應h()，得到線性相位和因果的 FIR 濾波器，這種方法的重點是選擇一個合適的窗函數和理想濾波器，使設計的濾波器的單位脈沖響應逼近理想濾波器的單位脈沖響應。2、設計步驟給定理想濾波器的頻率響應H()，在通帶上具有單位增益和線性相位，在阻帶上具有零響應。一個帶寬為( Wp=0.2*pi; Wst=0.3*pi; tr_width=Wst-Wp; Wn=(Wp+Wst)/2; Wn1=Wn/pi; Rp=0.25; As=50; N=ceil(1.8*pi/tr_width)+1; h=fir1(N,Wn1,boxcar(N+1); H,w=freqz(h,1); figure(1);

55、plot(w/pi,20*log10(abs(H)/max(abs(H); grid on; xlabel(omegapi);ylabel(dB); figure(2);plot(w/pi,angle(H); xlabel(omegapi);ylabel(rad); grid on;頻率特性曲線：0-10-20-30dB-40-50-60-70-80-9000.10.20.30.40.50.60.70.80.914321rad0-1-2-3-400.10.20.30.40.50.60.70.80.91分析：阻帶衰減小于 50dB，故不滿足設計指標要求。漢寧窗：實驗代碼： Wp=0.2*pi;W

56、st=0.3*pi;tr_width=Wst-Wp;Wn=(Wp+Wst)/2;Wn1=Wn/pi;Rp=0.25;As=50;N=ceil(6.2*pi/tr_width)+1;h=fir1(N,Wn1,hanning(N+1); H,w=freqz(h,1);figure(1);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel(omegapi);ylabel(dB);figure(2);plot(w/pi,angle(H);xlabel(omegapi);ylabel(rad);grid on;頻率特性曲線：0-50dB-100-15000.10.20.30.40.50.60.70.80.914321rad0-1-2-3-400.10.20.30.40.50.60.70.80.91分析：由圖可知，滿足設計指標要求。海明窗：實驗代碼： Wp=0.2*pi;Wst=0.3*pi;tr_width=Wst-Wp;Wn=(Wp+Ws