1、數(shù)字信號(hào)處理實(shí)驗(yàn)姓名：學(xué)號(hào)：班級(jí)： 專業(yè)：通信工程實(shí)驗(yàn) 1 基 2-FFT 算法實(shí)現(xiàn)一、實(shí)驗(yàn)?zāi)康恼莆栈?2-FFT 的原理及具體實(shí)現(xiàn)方法。編程實(shí)現(xiàn)基 2-FFT 算法。加深理解 FFT 算法的特點(diǎn)。二、實(shí)驗(yàn)設(shè)備與環(huán)境計(jì)算機(jī)、MATLAB 軟件環(huán)境。三、實(shí)驗(yàn)基礎(chǔ)理論FFT 是 DFT 的一種快速算法，能使 DFT 的計(jì)算大大簡(jiǎn)化，運(yùn)算時(shí)間縮短。FFT 利用了的三個(gè)固有特性，即對(duì)稱性、周期性和可約性，將長(zhǎng)序列的 DFT 分解為短序列的 DFT，合并了DFT 運(yùn)算中的某些項(xiàng)，從而減少了 DFT 的運(yùn)算量。FFT 算法基本上可以分為兩大類，即按時(shí)間抽取法和按頻率抽取法。在實(shí)現(xiàn) FFT 算法時(shí)，要重點(diǎn)考

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

3、程實(shí)現(xiàn)序列長(zhǎng)度為N=8 的按時(shí)間抽取的基 2-FFT 算法。給定一個(gè) 8 點(diǎn)序列，采用編寫(xiě)的程序計(jì)算其 DFT，并與 MATLAB 中 fft 函數(shù)計(jì)算的結(jié)果相比較，以驗(yàn)證結(jié)果的正確性。編程實(shí)現(xiàn)序列長(zhǎng)度為 N=8 的按頻率抽取的基 2-FFT 算法。給定一個(gè) 8 點(diǎn)序列，采用編寫(xiě)的程序計(jì)算其 DFT，并與 MATLAB 中 fft 函數(shù)計(jì)算的結(jié)果相比較，以驗(yàn)證結(jié)果的正確性。將上述 FFT 程序推廣到序列長(zhǎng)度為N=2v 的情況，要求利用原位運(yùn)算。五、實(shí)驗(yàn)代碼及分析實(shí)驗(yàn)代碼： 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 兩者結(jié)果完全相同，說(shuō)明結(jié)果正確。實(shí)驗(yàn)代碼：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分析：兩種算法的結(jié)果相同，說(shuō)明結(jié)果相同。實(shí)驗(yàn)代碼：按時(shí)間抽取 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六、 收獲與感悟通過(guò)本次數(shù)字信號(hào)處理實(shí)驗(yàn)，學(xué)習(xí)了基 2-FFT 算法的實(shí)現(xiàn)，通過(guò)對(duì) MATLAB 中 fft 函數(shù)的運(yùn)行結(jié)果與所編寫(xiě)程序結(jié)果的對(duì)比，了解了按時(shí)間抽取算法和按頻率抽取算法的異同，按時(shí)間抽取算法，兩種方式的結(jié)果相同，而按頻率抽取算法，兩種方式的結(jié)果不同。并且知道了碟形圖在 MATLAB 中的實(shí)現(xiàn)方法實(shí)驗(yàn) 2 利用 DFT 分析信號(hào)頻譜一、實(shí)驗(yàn)?zāi)康募由顚?duì) DFT 原理的理解。應(yīng)用 DFT 分析信號(hào)的頻譜。深刻理解利用 DFT 分析信號(hào)頻譜的原理，分析實(shí)現(xiàn)過(guò)程中出現(xiàn)的現(xiàn)象及解決方法。二、實(shí)驗(yàn)設(shè)備與環(huán)境計(jì)算機(jī)、MATLAB 軟件環(huán)境。三、實(shí)驗(yàn)基礎(chǔ)理論DFT 和 DTFT 的關(guān)系有

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

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

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

17、窗函數(shù)均會(huì)引起誤差。因此，要合理選取采樣間隔和截取長(zhǎng)度，必要時(shí)還考慮加適當(dāng)?shù)拇?。?duì)于連續(xù)時(shí)間周期信號(hào)，我們?cè)诓捎糜?jì)算機(jī)進(jìn)行計(jì)算時(shí)，也總是要進(jìn)行截?cái)?，序列總是有限長(zhǎng)的，仍然可以采用上述方法近似計(jì)算。可能用到的 MATLAB 函數(shù)與代碼實(shí)驗(yàn)中 DFT 運(yùn)算可采用 MATLAB 中提供的函數(shù) fft 來(lái)實(shí)現(xiàn)。DTFT 可以利用 MATLAB 矩陣運(yùn)算的方法進(jìn)行計(jì)算X() = = , , , 1 2 四、實(shí)驗(yàn)內(nèi)容=1121.已知 x(n)=2,-1,1,1,完成如下要求：計(jì)算器 DTFT，并畫(huà)出*, +區(qū)間的波形。計(jì)算 4 點(diǎn) DFT，并把結(jié)果顯示在(1)所畫(huà)的圖形中。對(duì) x(n)補(bǔ)零，計(jì)算 64 點(diǎn)

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

19、續(xù)時(shí)間信號(hào) 的頻譜（幅度譜）。分析采用不同的采樣間隔和截取長(zhǎng)度進(jìn)行計(jì)算的結(jié)果，并最終確定適合的參數(shù)。五、實(shí)驗(yàn)結(jié)果及分析1.（1）實(shí)驗(yàn)代碼：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/實(shí)驗(yàn)代碼：fft(x);k=0:3;subplot(211);hold onstem(k,abs(ans);subplot(212);hold onstem(k,angle(ans)/pi);計(jì)算結(jié)果圖示如下：Magnitude321-3-2-10123/ Phase0.50-0.5-3-2-10123/實(shí)驗(yàn)代碼：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 計(jì)算 DTFT。實(shí)現(xiàn)方法：利用對(duì) x(n)補(bǔ)零，提高采樣密度，當(dāng)補(bǔ)零的個(gè)數(shù)足夠多時(shí)，即可用 DFT 的結(jié)果近似 DTFT。2.（1）當(dāng) 0 n 10時(shí)，實(shí)驗(yàn)代碼: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當(dāng) x(n)的長(zhǎng)度變?yōu)?100 時(shí)，實(shí)驗(yàn)代碼：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）當(dāng)0 n 100時(shí)，實(shí)驗(yàn)代碼: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）分析：由頻譜波形可以看出，通過(guò)增加補(bǔ)零的個(gè)數(shù)可以提高頻譜分辨率，但是不能提高分辨力，提高分辨力需要通過(guò)延長(zhǎng)所選信號(hào)的長(zhǎng)度來(lái)實(shí)現(xiàn)。實(shí)驗(yàn)代碼： 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;時(shí)域波形圖為：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 個(gè)頻譜分量。實(shí)驗(yàn)代碼：（1）采樣區(qū)間：0，100，采樣間隔為 1 時(shí)： 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 時(shí)： 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 時(shí)： 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 時(shí)： 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 時(shí)： 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 時(shí)： n=0:2:150; x=exp(-0.1*n); X5=fft(x); k=0:75; stem(k,abs(X5); xlabel(k);波形圖：654321001020304050607080k分析：由

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

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

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

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

32、率為 fc=1kHz。設(shè)采樣頻率為 fs=10kHz，設(shè)計(jì)數(shù)字低通濾波器，滿足如下指標(biāo)通帶截止頻率：fp=1kHz，通帶波動(dòng)：Rp=1dB阻帶截止頻率：fst=1.5kHz，阻帶衰減：As=15dB要求分別采用巴特沃斯、切比雪夫 I 型、切比雪夫 II 型和橢圓模擬原型濾波器及脈沖響應(yīng)不變法進(jìn)行設(shè)計(jì)。結(jié)合實(shí)驗(yàn)結(jié)果，分別討論采用上述設(shè)計(jì)的數(shù)字濾波器是否都能滿足給定指標(biāo)要求，分析脈沖響應(yīng)不變法設(shè)計(jì) IIR 數(shù)字濾波器的優(yōu)缺點(diǎn)及適用范圍。五、實(shí)驗(yàn)結(jié)果及分析1.實(shí)驗(yàn)代碼：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）巴特沃斯濾波器實(shí)驗(yàn)代碼： 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()由圖線可知，結(jié)果滿足設(shè)計(jì)指標(biāo)要求。切比雪夫 I 型實(shí)驗(yàn)代碼：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()由圖線可知，滿足設(shè)計(jì)指標(biāo)要求。切比雪夫 II 型實(shí)驗(yàn)代碼： 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，不滿足設(shè)計(jì)要求。橢圓模擬原型濾波器實(shí)驗(yàn)代碼： 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()由圖線可知，滿足設(shè)計(jì)指標(biāo)要求。脈沖相應(yīng)不變法： 實(shí)驗(yàn)代碼： 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、，滿足設(shè)計(jì)指標(biāo)要求。分析：脈沖響應(yīng)不變法設(shè)計(jì) IIR 數(shù)字濾波器的優(yōu)缺點(diǎn)：優(yōu)點(diǎn)：脈沖響應(yīng)不變法的頻率坐標(biāo)的變換是線性的。因此，如果模擬濾波器的的頻率響應(yīng)是限帶于折疊頻率以內(nèi)的話，則通過(guò)變換后所得的數(shù)字濾波器的頻率響應(yīng)可以不失真的反映原響應(yīng)與頻率的關(guān)系。缺點(diǎn)：脈沖響應(yīng)不變法最大的缺點(diǎn)是有頻譜的周期延拓效應(yīng)。其次，脈沖響應(yīng)不變法所得的數(shù)字濾波器會(huì)出現(xiàn)頻率的混疊失真。第三，數(shù)字濾波器在 T 很?。ɑ蛘哒f(shuō)取樣頻率很高）時(shí)，可能具有不希望的高增益。適用范圍：脈沖響應(yīng)不變法設(shè)計(jì) IIR 數(shù)字濾波器只能用于限帶的頻率響應(yīng)，如衰減特性很好的低通或帶通。六、收獲與體會(huì)通過(guò)本次實(shí)驗(yàn)，學(xué)會(huì)了利用用脈沖響應(yīng)不變法設(shè)計(jì)

53、 IIR 數(shù)字濾波器的原理及具體方法， 通過(guò)比較不同類型的數(shù)字濾波器，了解了它們各自的優(yōu)缺點(diǎn)以及適用范圍，在相同的技術(shù)指標(biāo)下，不同濾波器的效果也有較大差異，不是所有類型都可以滿足要求。了解了如何將所給的技術(shù)指標(biāo)，轉(zhuǎn)化為設(shè)計(jì)過(guò)程中的參數(shù)，包括模擬域和數(shù)字域的轉(zhuǎn)化，加深了對(duì)課本內(nèi)容的理解與認(rèn)識(shí)，進(jìn)一步掌握了 IIR 數(shù)字濾波器的設(shè)計(jì)要點(diǎn)。實(shí)驗(yàn) 7窗函數(shù)法設(shè)計(jì) FIR 數(shù)字濾波器一、實(shí)驗(yàn)?zāi)康恼莆沾昂瘮?shù)法設(shè)計(jì) FIR 數(shù)字濾波器的原理和具體方法二、實(shí)驗(yàn)設(shè)備與環(huán)境計(jì)算機(jī)、Matlab 軟件環(huán)境三、實(shí)驗(yàn)基礎(chǔ)理論1、基本原理窗函數(shù)設(shè)計(jì)法的基本思想為，首先選擇一個(gè)適當(dāng)?shù)睦硐氲臑V波器H()，然后用窗函數(shù)截取它的

54、單位脈沖響應(yīng)h()，得到線性相位和因果的 FIR 濾波器，這種方法的重點(diǎn)是選擇一個(gè)合適的窗函數(shù)和理想濾波器，使設(shè)計(jì)的濾波器的單位脈沖響應(yīng)逼近理想濾波器的單位脈沖響應(yīng)。2、設(shè)計(jì)步驟給定理想濾波器的頻率響應(yīng)H()，在通帶上具有單位增益和線性相位，在阻帶上具有零響應(yīng)。一個(gè)帶寬為( 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，故不滿足設(shè)計(jì)指標(biāo)要求。漢寧窗：實(shí)驗(yàn)代碼： 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分析：由圖可知，滿足設(shè)計(jì)指標(biāo)要求。海明窗：實(shí)驗(yàn)代碼： Wp=0.2*pi;Wst=0.3*pi;tr_width=Wst-Wp;Wn=(Wp+Ws