北京理工大學數(shù)字信號處理實驗報告

1、 本科實驗報告實驗名稱： 數(shù)字信號處理實驗 課程名稱：數(shù)字信號處理實驗實驗時間：任課教師：實驗地點：4423實驗教師：實驗類型： 原理驗證 綜合設(shè)計 自主創(chuàng)新學生姓名：學號/班級：組 號：學 院：信息與電子學院同組搭檔：專 業(yè)：信息工程成 績：實驗1 利用DFT分析信號頻譜一、實驗?zāi)康?.加深對DFT原理的理解。2.應(yīng)用DFT分析信號頻譜。3.深刻理解利用DFT分析信號頻譜的原理，分析現(xiàn)實過程現(xiàn)象及解決辦法。二、實驗原理1、DFT和DTFT的關(guān)系 有限長序列的離散時間傅里葉變換在頻率區(qū)間的N個等分點上的N個取樣值可以由下式表示： 由上式可知，序列的N點DFT，實際上就是序列的DTFT在N個等間

2、隔頻率點上樣本。2、利用DFT求DTFT 方法1：由恢復出的方法如圖2.1所示： 圖 2.1.由 N點DFT恢復頻譜DTFT的流程由圖2.1所示流程圖可知：由式2-2可以得到其中為內(nèi)插函數(shù) 方法2：然而在實際MATLAB計算中，上訴插值公式不見得是最好的方法。由于DFT是DTFT的取樣值，其相鄰的兩個頻率樣本點的間距為,所以如果我們增加數(shù)據(jù)的長度N，使得得到的DFT譜線就更加精細，其包絡(luò)就越接近DTFT的結(jié)果，這樣可以利用DFT來近似計算DTFT。如果沒有更多的數(shù)據(jù)，可以通過補零來增加數(shù)據(jù)長度。3、利用DFT分析連續(xù)時間信號的頻譜 采用計算機分析連續(xù)時間信號的頻譜，第一步就是把連續(xù)時間信號離散

3、化，這里需要進行連個操作：一是采樣，二是截斷。對于連續(xù)非周期信號，按采樣間隔T進行采樣，截取長度為M，那么對進行N點的頻率采樣，得到因此，可以將利用DFT分析連續(xù)非周期信號頻譜的步驟歸納如下： （1）確定時域采樣間隔T，得到離散序列； （2）確定截取長度M，得到M點離散序列，這里的為窗函數(shù)。 （3）確定頻域采樣點數(shù)N，要求。 （4）利用FFT計算離散序列的N點DFT，得到。 （5）根據(jù)式（2-6）由計算采樣點的近似值。 采用上訴方法計算的頻譜，需要注意如下三點問題：（1）頻譜混疊。如果不滿足采樣定理的條件，頻譜會很出現(xiàn)混疊誤差。對于頻譜無限寬的信號，應(yīng)考慮覆蓋大部分主要頻率的范圍。（2）柵欄效

4、應(yīng)和頻譜分辨率。使用DFT計算頻譜，得到的結(jié)果只是N個頻譜樣本值，樣本值之間的頻譜是未知的，就像通過一個柵欄觀察頻譜，稱為“柵欄效應(yīng)”。頻譜分辨率與記錄長度成正比，提高頻譜分辨率，就要增加記錄時間。（3）頻譜泄露。對于信號截斷會把窗函數(shù)的頻譜會引入到信號頻譜中，造成頻譜泄露。解決這問題的主要辦法是采用旁瓣小的窗函數(shù)，頻譜泄露和窗函數(shù)均會引起誤差。 因此，要合理選取采樣間隔和截取長度，必要時還需考慮適當?shù)拇啊?對于連續(xù)周期信號，我們在采用計算機進行計算時，也總是要進行截斷，序列總是有限長的，仍然可以采用上訴方法近似計算。4、可能用到MATLAB函數(shù)與代碼實驗中的DFT運算可以采用MATLAB中提

5、供的FFT來實現(xiàn)。DTFT可以利用MATLAB矩陣運算的方法進行計算。三、實驗內(nèi)容1. ，完成如下要求：（1）計算其DTFT，并畫出區(qū)間的波形。（2）計算4 點DFT，并把結(jié)果顯示在（1）所畫的圖形中。（3）對補零，計算64 點DFT，并顯示結(jié)果。（4）是否可以由DFT 計算DTFT，如果可以，請編程實現(xiàn)。2. 考察序列（1）時，用DFT 估計的頻譜；將補零加長到長度為100點序列用DFT估計的頻譜。要求畫出相應(yīng)波形。（2）時，用DFT 估計x(n)的頻譜，并畫出波形。3. 已知信號 ，其中，。從的表達式可以看出，它包含三個頻率的正弦波，但是，從其時域波形來看，似乎是一個正弦信號，利用DFT做

6、頻譜分析，確定適合的參數(shù)，使得到的頻譜的頻率分辨率符合需要。4.利用DFT近似分析連續(xù)時間信號xt=e-0.1tu(t)的頻譜（幅度譜）。分析采用不同的采樣間隔和截取長度進行計算的結(jié)果，并最終確定合適的參數(shù)。四、實驗代碼及實驗結(jié)果實驗1.1.1實驗結(jié)果：實驗代碼：>> n=0:3;>> x=2 -1 1 1;>> w=-pi:0.01*pi:pi;>> X=x*exp(-j*n'*w);>> subplot(211);>> plot(w,abs(X);>> title('幅度');xla

7、bel('w');ylabel('|X|');>> axis tight;>> subplot(212);>> plot(w,angle(X);>> title('相位');xlabel('w');ylabel('Angle（X）');>> axis tight;實驗1.1.2實驗結(jié)果：實驗代碼： n=0:3;x=2 -1 1 1;w=-pi:0.01*pi:pi;X=x*exp(-j*n'*w);subplot(211);plot(w,abs(X

8、);title('幅度');xlabel('w');ylabel('|X|');axis tight;hold on;subplot(212);plot(w,angle(X);title('相位');xlabel('w');ylabel('Angle（X）');axis tight;hold on;H=fft(x);subplot(211);>> stem(n,abs(H),'filled');>> subplot(212);>> stem(n,a

9、ngle(H),'filled');實驗1.1.3實驗結(jié)果：實驗代碼：>> x=2 -1 1 1 zeros(1,60);>> X=fft(x);>> subplot(211);>> n=0:63;>> stem(n,abs(X),'filled');>> title('幅度');xlabel('n');ylabel('|X|');>> subplot(212);>> stem(n,angle(X),'fille

10、d');title('相位');xlabel('n');ylabel('angle(X)');分析：可以由DFT計算DTFT。通過補零加長序列，提高采樣密度，可以由DFT近似計算DTFT。實驗1.2.1實驗結(jié)果：實驗代碼：n=0:10;x=cos(0.48*n*pi)+cos(0.52*n*pi);X=fft(x);subplot(211);stem(n,abs(X),'filled');title('幅度');xlabel('n');ylabel('|X|');subplo

11、t(212);stem(n,angle(X),'filled');title('相位');xlabel('n');ylabel('angle(X)');補零加長：實驗結(jié)果：實驗代碼：h=x zeros(1,89);H=fft(h);n=0:99;subplot(211);stem(n,abs(H),'filled');title('幅度');xlabel('n');ylabel('|H|');subplot(212);stem(n,angle(H),'fill

12、ed');title('相位');xlabel('n');ylabel('angle(H)');實驗1.2.2實驗結(jié)果：實驗代碼：n=0:100;>> x=cos(0.48*n*pi)+cos(0.52*n*pi);>> X=fft(x);>> subplot(211);stem(n,abs(X),'filled');title('幅度');xlabel('n');ylabel('|X|');subplot(212);stem(n,angl

13、e(X),'filled');title('相位');xlabel('n');ylabel('angle(X)');分析：可以通過增大截取長度和增加補零的個數(shù)來提高頻譜分辨率，但是補零不能夠增加分辨力。實驗1.3實驗結(jié)果：實驗代碼：f1=1;f2=2;f3=3;Fs=500;Tp=1;t=0:1/Fs:Tp;x=0.15*sin(2*pi*f1*t)+sin(2*pi*f2*t)-0.1*sin(2*pi*f3*t);N=501;F=Fs/N;f=0:Fs/(N-1):Fs;X=fft(x,N)；stem(f,abs(X)/250

14、,'filled');axis(0,4,0,1);xlabel('f/Hz');title('Magnitude');分析：通過選取合適的采樣周期，可以完整的恢復出原信號的頻譜波形。實驗1.4Tp=1s，0n100實驗結(jié)果：實驗代碼：n=0:1:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=5s,0n100實驗結(jié)果：實驗代碼：n=0:5:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=25s，0

15、n100實驗結(jié)果：實驗代碼：n=0:25:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=1s，0n50實驗結(jié)果：實驗代碼：n=0:1:50;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=1,0n25實驗結(jié)果：實驗代碼：n=0:1:25;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');分析： 最終確定參數(shù)：tp=5s，0n100。五、心得與體會通過本次實驗，我們掌握并加深對DFT原理的理解

16、并且學會應(yīng)用DFT分析信號頻譜，在此基礎(chǔ)上利用DFT分析信號頻譜的原理，掌握了利用matlab分析現(xiàn)實問題的步驟及辦法。并且，通過這次的實驗，對信號序列有了更加深刻的認識，單純的一個信號序列是沒有意義的，只有配合他本身的時間序列才是一個完整的信號序列，才可以對其進行分析。實驗2 利用FFT計算線性卷積（選作）一、實驗?zāi)康?.掌握利用FFT計算線性卷積的原理及具體實現(xiàn)方法。2.加深理解重疊相加法和重疊保留法。3.考察利用FFT計算線性卷積各種方法的適用范圍。二、實驗原理1.線性卷積與圓周卷積設(shè)x(n)為L點序列，h(n)為M點序列，x(n)和h(n)的線性卷積為 （3-1）的長度為L+M-1x(

17、n)和h(n)的圓周卷積為 （3-2）圓周卷積與線性卷積相等而不產(chǎn)生交疊的必要條件為NL+M+1 (3-3)圓周卷積定理：根據(jù)DFT性質(zhì)，x(n)和h(n)的N點圓周卷積的DFT等于它們的DFT的乘積： (3-4)2.快速卷積快速卷積發(fā)運用圓周卷積實現(xiàn)線性卷積，根據(jù)圓周卷積定理利用FFT算法實現(xiàn)圓周卷積?？蓪⒖焖倬矸e運算的步驟歸納如下：(1)必須選擇；為了能使用基-2算法，要求。采用補零的辦法使得x(n)和h(n)的長度均為N。(2)計算x(n)和h(n)的N點FFT。(3)組成乘積(4)利用IFFT計算Y(k)的IDFT，得到線性卷積y(n)3.分段卷積我們考察單位取樣響應(yīng)為h(n)的線性系

18、統(tǒng)，輸入為x(n)，輸出為y(n)，則yn=xn*h(n)當輸入序列x(n)極長時，如果要等x(n)全部集齊時再開始進行卷積，會使輸出有較大延時；如果序列太長，需要大量存儲單元。為此，我們把x(n)分段，為別求出每段的卷積，合在一起得到最后的總輸出。這稱為分段卷積。分段卷積可以細分為重疊保留法和重疊相加法。重疊保留法：設(shè)x(n)的長度為，h(n)的長度為M。把序列x(n)分成多段N點序列，每段雨前一段重寫M-1個樣本。并在第一個輸入段前面補M-1個零。計算每一段與h(n)的圓周卷積，其結(jié)果中前M-1個不等與線性卷積，應(yīng)當舍去，只保留后面N-M+1個正確的輸出樣本，把它們合起來得到總的輸出。利用

19、FFT實現(xiàn)重疊保留法的步驟如下：(1)在x(n)前面填充M-1個零，擴大以后的序列為(2)將x(n)分為若干段N點子段，設(shè)L=N-M+1為每一段的有效長度，則第i段的數(shù)據(jù)為：(3)計算每一段與h(n)的N點圓周卷積，利用FFT計算圓周卷積(4)舍去每一段卷積結(jié)果的前M-1個樣本，連接剩下的樣本得到卷積結(jié)果y(n)。重疊相加法：設(shè)h(n)長度為M，將信號x(n)分解成長為L的子段。以表示沒斷信號，則：每一段卷積的長度為L+M-1，所以在做求和時，相鄰兩段序列由M-1個樣本重疊，即前一段的最后M-1個樣本和下一段前M-1個樣本序列重疊，這個重疊部分相加，再與不重疊的部分共同組成y(n)。利用FFT

20、實現(xiàn)重疊保留法的步驟如下：(1)將x(n)分為若干L點子段。(2)計算每一段與h(n)的卷積，根據(jù)快速卷積法利用FFT計算卷積。(3)將各段相加，得到輸出y(n)。4、可能得到的MATLAB函數(shù) 實驗中FFT運算可采用MATLAB中提供的函數(shù)fft來實現(xiàn)。三、實驗內(nèi)容假設(shè)要計算序列x(n)=u(n)-u(n-L),0nL和h(n)=cos(0.2n),0nM的線性卷積完成以下實驗內(nèi)容。1.設(shè)L=M,根據(jù)線性卷積的表達式和快速卷積的原理分別編程實現(xiàn)計算兩個序列線性卷積的方法,比較當序列長度分別為8,16,32,64,256,512,1024時兩種方法計算線性卷積所需時間。2當L=2048且M=2

21、56時比較直接計算線性卷積和快速卷積所需的時間,進一步考察當L=4096且M=256時兩種算法所需的時間。3. 編程實現(xiàn)利用重疊相加法計算兩個序列的線性卷積,考察L=2048且M=256時計算線性卷積的時間,與2題的結(jié)果進行比較。4. 編程實現(xiàn)利用重疊保留法計算兩個序列的線性卷積,考察L=2048且M=256時計算線性卷積的時間,與2題的結(jié)果進行比較。四、實驗代碼及實驗結(jié)果實驗2.1實驗代碼：for M=8 16 32 64 256 512 1024;L=M;n=0:1:L-1;x=ones(1,L);h=cos(0.2.*pi.*n);disp('L=M=');disp(M)

22、disp('線性卷積用時：')ticy1=conv(x,h);tocdisp('快速卷積用時：');ticX1=fft(x);H=fft(h);Y=X1.*H;y=ifft(Y);tocend實驗結(jié)果：L=M= 8線性卷積用時：時間已過 0.000122 秒。快速卷積用時：時間已過 0.000042 秒。L=M= 16線性卷積用時：時間已過 0.000045 秒。快速卷積用時：時間已過 0.000026 秒。L=M= 32線性卷積用時：時間已過 0.000051 秒。快速卷積用時：時間已過 0.000030 秒。L=M= 64線性卷積用時：時間已過 0.0000

23、50 秒?？焖倬矸e用時：時間已過 0.000030 秒。L=M= 256線性卷積用時：時間已過 0.000089 秒?？焖倬矸e用時：時間已過 0.000066 秒。L=M= 512線性卷積用時：時間已過 0.000540 秒?？焖倬矸e用時：時間已過 0.001388 秒。L=M= 1024線性卷積用時：時間已過 0.000617 秒。快速卷積用時：時間已過 0.000115 秒。分析：可見在在相同長度下，快速卷積比線性卷積快差不多一倍的時間。實驗2.2實驗代碼：L=2048;M=256;n1=0:1:L;n2=0:1:M;x=ones(1,L);h=cos(0.2.*pi.*n2);disp(

24、'線性卷積：')ticy1=conv(x,h);tocdisp('快速卷積：')N=M+L-1;ticX1=fft(x,N);H=fft(h,N);Y=X1.*H;y=ifft(Y,N);toc實驗結(jié)果： L=2048，M=256線性卷積：時間已過 0.000437 秒?？焖倬矸e：時間已過 0.000407 秒。L=4096，M=256線性卷積：時間已過 0.000445 秒?？焖倬矸e：時間已過 0.001395 秒。分析：當序列長度較短時，快速卷積比線性卷積快；當序列長度過長時，快速卷積比線性卷積慢了許多。實驗2.3實驗代碼：L=2048;M=256;disp

25、('L=');disp(L);disp('M=');disp(M);x=ones(1,L);n=0:1:M-1;h=cos(0.2.*pi.*n);Li=M;Num=L/Li;Ni=Li+M-1;N_yu=L+M-1;y = zeros(1,N_yu);disp('重疊相加法卷積：');ticH = fft(h,Ni);for i=1:Num;i_low=(i-1)*M+1;xi=x(i_low:i_low+M-1);Xi=fft(xi,Ni);Yi=Xi.*H;yi=ifft(Yi,Ni);y(i_low:i_low+Ni-1)=y(i_low

26、:i_low+Ni-1)+yi;endtoc實驗結(jié)果：L= 2048M= 256重疊相加法卷積時間時間已過 0.000626 秒。分析：較第2題結(jié)果快。實驗2.4實驗代碼：L=2048;M=256;disp('L=');disp(L);disp('M=');disp(M);x=ones(1,L);n=0:1:M-1;h=cos(0.2.*pi.*n);Li=M;Num=L/Li;Ni=Li+M-1;N_yu=L+M-1;y=zeros(1,N_yu);x_add=zeros(1,M-1),x;disp('Öصþ&

27、#177;£Áô·¨¾í»ýʱ¼ä');ticH = fft(h,Ni);for i=1:Num;i_low=(i-1)*M+1;xi=x_add(i_low:i_low+M-1);Xi=fft(xi,Ni);Yi=Xi.*H;yi=ifft(Yi,Ni);y(i_low:i_low+M-1)=yi(M:Ni);endxi=x_add(i*M+1):N_yu);Xi=fft(xi,Ni);Y=Xi.*H;yi=ifft(Y,Ni);y(i*M+1):N_y

28、u)=yi(M:Ni-1);toc實驗結(jié)果：L= 2048M= 256重疊保留法卷積時間時間已過 0.000508 秒。分析：較第2題結(jié)果快。五、心得與體會本次實驗要求我們掌握利用FFT計算線性卷積的原理及具體實現(xiàn)方法，以及通過對快速卷積和線性卷積運算速度的比較，更加直觀的看到利用FFT快速卷積的優(yōu)點。本次實驗讓我切實看到了FFT算法的高效性及對于龐大的數(shù)據(jù)的處理實時性，同時對重疊保留法、重疊相加法也有了更深的認識，對以后的實驗有了很好的理論基礎(chǔ)，受益頗多。實驗3 IIR數(shù)字濾波器設(shè)計一、實驗?zāi)康?.掌握利用脈沖響應(yīng)不變法和雙線性變換法設(shè)計IIR數(shù)字濾波器的原理及具體方法。2.加深理解數(shù)字濾波

29、器和模擬濾波器之間的技術(shù)指標轉(zhuǎn)化。3.掌握脈沖響應(yīng)不變法和雙線性變換法設(shè)計IIR數(shù)字濾波器的優(yōu)缺點及適用范圍。二、實驗內(nèi)容1.設(shè)采樣頻率為fs =4kHz，采用脈沖響應(yīng)不變法設(shè)計一個三階巴特沃斯數(shù)字低通濾波器，其3dB截止頻率為fc=1kHz。2.設(shè)采樣頻率為fs=10kHz，設(shè)計數(shù)字低通濾波器，滿足如下指標 通帶截止頻率：fp=1kHz，通帶波動：Rp=1dB 阻帶截止頻率：fst=1.5kHz，阻帶衰減：As=15dB要求分別采用巴特沃斯、切比雪夫I型、切比雪夫II型和橢圓模擬原型濾波器及脈沖響應(yīng)不變法進行設(shè)計。結(jié)合實驗結(jié)果，分別討論采用上述設(shè)計的數(shù)字濾波器是否都能滿足給定指標要求，分析脈

30、沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的優(yōu)缺點及適用范圍。四、實驗代碼及實驗結(jié)果巴特沃斯：實驗代碼：fs=10*1000;fp=1*1000;fst=1.5*1000;Rp=1;As=15;Wp = 2*pi*fp/fs;Ws = 2*pi*fst/fs;N = ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws);Wc = Wp /(10(Rp/10)-1)(1/(2*N);b,a = butter(N,Wc/pi,'s');bz,az = impinvar(b,a,1)w=0:500*pi/500;H,w=freqz(bz,a

31、z); %頻率響應(yīng)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)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);grid on;xlabel('omega(pi)');y

32、label('Phase of H(ejomega)(pi)');grd=grpdelay(bz,az,w); %群延時響應(yīng)subplot(224);plot(w/pi,grd);grid on;xlabel('omega(pi)');ylabel('Group delay');實驗結(jié)果：bz = 1.0e-04 * -0.0000 0.0091 0.2033 0.4463 0.1524 0.0051 0az =1.0000 -5.1378 11.0528 -12.7383 8.2921 -2.8898 0.4211切比雪夫型：實驗代碼：fs=1

33、0*1000; fp=1*1000;fst=1.5*1000;Rp=1;As=15;Wp = 2*pi*fp/fs; Ws = 2*pi*fst/fs; e = sqrt(10(Rp/10)-1);A = 10(As/20);N = ceil(acosh( sqrt(A2-1)/e) / ( acosh(Ws/Wp);b,a = cheby1(N,Rp,Wp/pi,'s');bz,az = impinvar(b,a,1)w=0:500*pi/500;H,w=freqz(bz,az); subplot(221);plot(w/pi,abs(H);grid on;xlabel(&#

34、39;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)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);grid on;xlabel('omega(pi)');ylabel('Phase of H(ejomega)(pi)');grd=grpdelay(bz

35、,az,w); subplot(224);plot(w/pi,grd);grid on;xlabel('omega(pi)');ylabel('Group delay');實驗結(jié)果：bz = 1.0e-03 * 0.0000 0.0623 0.2372 0.0567 0az =1.0000 -3.7710 5.3742 -3.4293 0.8265切比雪夫型：實驗代碼：fs=10*1000;fp=1*1000;fst=1.5*1000;Rp=1;As=15;Wp = 2*pi*fp/fs; Ws = 2*pi*fst/fs; e = sqrt(10(Rp/10)

36、-1);A = 10(As/20);N = ceil( acosh( sqrt(A2-1)/e) / ( acosh(Ws/Wp) );b,a = cheby2(N,As,Wp/pi,'s');bz,az = impinvar(b,a,1)w=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

37、(abs(H);grid on;xlabel('omega(pi)');ylabel('|H(ejomega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);grid on;xlabel('omega(pi)');ylabel('Phase of H(ejomega)(pi)');grd=grpdelay(bz,az,w); subplot(224);plot(w/pi,grd);grid on;xlabel('omega(pi)');ylabel('Group de

38、lay');實驗結(jié)果：bz = 0.0884 -0.3342 0.5117 -0.3718 0.1076az = 1.0000 -3.4904 4.6022 -2.7148 0.6048橢圓型：實驗代碼：fs=10*1000;fp=1*1000;fst=1.5*1000;Rp=1;As=15;Wp = 2*pi*fp/fs; Ws = 2*pi*fst/fs; e = sqrt(10(Rp/10)-1);A = 10(As/20);N = ceil(acosh( sqrt(A2-1)/e) / ( acosh(Ws/Wp);b,a = ellip(N,Rp,As,Wp/pi,'

39、s');bz,az = impinvar(b,a,1)w=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)|(dB)');subplot(223);plot(w

40、/pi,angle(H)/pi);grid on;xlabel('omega(pi)');ylabel('Phase of H(ejomega)(pi)');grd=grpdelay(bz,az,w); subplot(224);plot(w/pi,grd);grid on;xlabel('omega(pi)');ylabel('Group delay');實驗結(jié)果：bz = 0.1470 -0.5636 0.8366 -0.56az =1.0000 -3.7726 5.3933 -3.4604 0.8408分析：脈沖響應(yīng)不變法的

41、優(yōu)點是頻率坐標的變換是線性的，因此如果模擬濾波器的頻響是限帶于折疊頻率以內(nèi)的話，通過變換后的數(shù)字濾波器的頻響可以不失真地反映原響應(yīng)與頻率之間的關(guān)系。但是，其最大的缺點是有頻譜的周期延拓效應(yīng)，會產(chǎn)生混疊失真。故，脈沖響應(yīng)不變法只能用于限帶的頻響特性，如衰減特性較好的低通，或帶通，而且高頻衰減越大，頻響的混疊效應(yīng)就越小。雙線性變換法的優(yōu)點是消除了脈沖響應(yīng)不變法的固有的混疊失真，但缺點是頻率變換之間的非常嚴重的非線性，會使變換后時域上的圖像產(chǎn)生非常嚴重的失真。五、心得與體會 通過實驗，學習了利用脈沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的基本原理為：從時域響應(yīng)出發(fā)，使數(shù)字濾波器的單位脈沖響應(yīng)h(n)模仿模擬

42、濾波器的單位沖激響應(yīng)ha(t)，h(n)等于ha(t)的取樣值。所以設(shè)計中，應(yīng)先根據(jù)數(shù)字指標轉(zhuǎn)換來的模擬指標，設(shè)計模擬濾波器，再利用impinvar函數(shù)實現(xiàn)脈沖響應(yīng)不變法模擬濾波器到數(shù)字濾波器的變換。脈沖響應(yīng)不變法和雙線性法設(shè)計IIR數(shù)字濾波器設(shè)計的基礎(chǔ)，掌握了這倆種方法能夠使我們對數(shù)字濾波器設(shè)計有更加深刻的認識。實驗4 頻率取樣法設(shè)計FIR數(shù)字濾波器一、實驗?zāi)康恼莆疹l率取樣法設(shè)計FIR數(shù)字濾波器的原理及具體方法二、實驗原理1、基本原理 頻率取樣法從頻域出發(fā)，把理想的濾波器等間隔采樣得到,將作為實際設(shè)計濾波器的。 （8-1）得到以后可以由來確定唯一確定濾波器的單位脈沖響應(yīng),可以由求得： （8-

43、2、3）其中為內(nèi)插函數(shù) （8-4）由求得的頻率響應(yīng)將逼近 如果我們設(shè)計的是線性相位FIR濾波器，則的幅度和相位滿足線性相位濾波器的約束條件。我們將表示為如下形式： （8-5）當為實數(shù)，則由此得到 （8-6） 即為中心偶對稱。在利用線性相位條件可知，對于1型和2型線性相位濾波器 （8-7）對于3型和4型線性相位濾波器 （8-8）其中，x表示取小于該數(shù)的最大的整數(shù)。設(shè)計步驟（1）由給定的理想濾波器給出和。（2）由求得（3) 根據(jù)求得或三、實驗內(nèi)容1. 設(shè)計一個數(shù)字低通FIR濾波器，其技術(shù)指標如下：_p=0.2，R_p=0.25dB_st=0.3，A_s=50dB分別采用矩形窗、漢寧窗、海明窗、布萊

44、克曼窗、凱瑟窗設(shè)計該濾波器。結(jié)合實驗結(jié)果分別討論上述方法設(shè)計的數(shù)字濾波器是否符合指標。2. 設(shè)計一個數(shù)字帶通FIR濾波器，其技術(shù)指標如下：3. 采用頻率取樣設(shè)計法設(shè)計FIR數(shù)字低通濾波器，滿足以下指標_p=0.2，R_p=0.25dB_st=0.3，A_s=50dB（1）取N=20，過渡帶沒有樣本。（2）取N=40，過渡帶有一個樣本，T=0.39。（3）取N=60，過渡帶有兩個樣本，T1=0.5925, T2=0.1099。（4）分別采用上述方法設(shè)計的數(shù)字濾波器是否都能滿足給定的指標要求。4. 采用頻率取樣技術(shù)設(shè)計下面的高通濾波器_st=0.6，A_s=50dB_p=0.8，R_p=1dB對于

45、高通濾波器，N必須為奇數(shù)（或1型濾波器）。選擇N=33，過渡帶有兩個樣本，過渡帶樣本最優(yōu)值為T1=0.1095，T2=0.598。四、實驗代碼及實驗結(jié)果實驗4.1矩形窗：實驗代碼：wp = 0.2 * pi; wst = 0.3 * pi;tr_width = wst - wp;N = ceil(1.8 * pi / tr_width) + 1;n = 0 : N - 1;wc = (wp + wst) / 2;alpha = (N - 1) / 2;hd = (wc / pi) * sinc(wc / pi) * (n - alpha);w_boxcar = boxcar(N)'h

46、= hd .* w_boxcar;subplot(221);stem(n, hd, 'filled');axis tight; xlabel('n'); ylabel('hd(n)');Hr, w1 = zerophase(h);subplot(222); plot(w1/pi, Hr);axis tight; xlabel('omega/pi'); ylabel('H(omega)');subplot(223); stem(n, h, 'filled');axis tight; xlabel(&#

47、39;n'); ylabel('h(n)');H, w = freqz(h, 1);subplot(224); plot(w / pi, 20 * log10(abs(H) / max(abs(H);xlabel('omega/pi'); ylabel('dB');grid on;實驗結(jié)果：分析：阻帶的衰減小于50db，所以，不符合指標。漢寧窗：實驗代碼：wp = 0.2 * pi; wst = 0.3 * pi;tr_width = wst - wp;N = ceil(6.2 * pi / tr_width) + 1;n = 0 : N

48、 - 1;wc = (wp + wst) / 2;alpha = (N - 1) / 2;hd = (wc / pi) * sinc(wc / pi) * (n - alpha);w_hanning = hanning(N)'h = hd .* w_hanning;subplot(221);stem(n, hd, 'filled');axis tight; xlabel('n'); ylabel('hd(n)');Hr, w1 = zerophase(h);subplot(222); plot(w1/pi, Hr);axis tight;

49、 xlabel('omega/pi'); ylabel('H(omega)');subplot(223); stem(n, h, 'filled');axis tight; xlabel('n'); ylabel('h(n)');H, w = freqz(h, 1);subplot(224); plot(w / pi, 20 * log10(abs(H) / max(abs(H);xlabel('omega/pi'); ylabel('dB');grid on;實驗結(jié)果：分析：阻帶衰

50、減小于50db，故，不符合指標。海明窗：實驗代碼：wp = 0.2 * pi; wst = 0.3 * pi;tr_width = wst - wp;N = ceil(6.6 * pi / tr_width) + 1;n = 0 : N - 1;wc = (wp + wst) / 2;alpha = (N - 1) / 2;hd = (wc / pi) * sinc(wc / pi) * (n - alpha);w_hamming = hamming(N)'h = hd .* w_hamming;subplot(221);stem(n, hd, 'filled');ax

51、is tight; xlabel('n'); ylabel('hd(n)');Hr, w1 = zerophase(h);subplot(222); plot(w1/pi, Hr);axis tight; xlabel('omega/pi'); ylabel('H(omega)');subplot(223); stem(n, h, 'filled');axis tight; xlabel('n'); ylabel('h(n)');H, w = freqz(h, 1);subplot(

52、224); plot(w / pi, 20 * log10(abs(H) / max(abs(H);xlabel('omega/pi'); ylabel('dB');grid on;實驗結(jié)果：分析：阻帶衰減大于50db，故，該數(shù)字濾波器符合指標。布萊克曼窗：實驗代碼：wp = 0.2 * pi; wst = 0.3 * pi;tr_width = wst - wp;N = ceil(11 * pi / tr_width) + 1;n = 0 : N - 1;wc = (wp + wst) / 2;alpha = (N - 1) / 2;hd = (wc / pi

53、) * sinc(wc / pi) * (n - alpha);w_blackman = blackman(N)'h = hd .* w_blackman;subplot(221);stem(n, hd, 'filled');axis tight; xlabel('n'); ylabel('hd(n)');Hr, w1 = zerophase(h);subplot(222); plot(w1/pi, Hr);axis tight; xlabel('omega/pi'); ylabel('H(omega)')

54、;subplot(223); stem(n, h, 'filled');axis tight; xlabel('n'); ylabel('h(n)');H, w = freqz(h, 1);subplot(224); plot(w / pi, 20 * log10(abs(H) / max(abs(H);xlabel('omega/pi'); ylabel('dB');grid on;實驗結(jié)果：分析：阻帶衰減大于50db，故，該數(shù)字濾波器符合要求。凱瑟窗：實驗代碼：wp = 0.2 * pi; wst = 0.3

55、* pi;tr_width = wst - wp;As = 50;N = ceil(As - 7.95) / (2.285 * tr_width) + 1;beta = 0.1102 * (As - 8.7);n = 0 : N - 1;wc = (wp + wst) / 2;alpha = (N - 1) / 2;hd = (wc / pi) * sinc(wc / pi) * (n - alpha);w_kaiser = kaiser(N, beta)'h = hd .* w_kaiser;subplot(221);stem(n, hd, 'filled');axi

56、s tight; xlabel('n'); ylabel('hd(n)');Hr, w1 = zerophase(h);subplot(222); plot(w1/pi, Hr);axis tight; xlabel('omega/pi'); ylabel('H(omega)');subplot(223); stem(n, h, 'filled');axis tight; xlabel('n'); ylabel('h(n)');H, w = freqz(h, 1);subplot(2

57、24); plot(w / pi, 20 * log10(abs(H) / max(abs(H);xlabel('omega/pi'); ylabel('dB');grid on;實驗結(jié)果：分析：阻帶衰減大于50db，故，該濾波器符合指標。實驗4.2實驗代碼：ws1 = 0.2 * pi; ws2 = 0.8 * pi;wp1 = 0.35 * pi; wp2 = 0.65 * pi;tr_width = wp1 - ws1;wc1 = (ws1 + wp1) / 2;wc2 = (ws2 + wp2) / 2;N = ceil(11 * pi / tr_width) + 1;n = 0 : N - 1;alpha = (N - 1) / 2