北京理工大學數字信號處理實驗報告

本科實驗報告實驗名稱：數字信號處理實驗課程名稱：數字信號處理實驗實驗時間：任課教師：實驗地點：4—423實驗教師：實驗類型：□原理驗證□綜合設計□自主創(chuàng)新學生姓名：組號：學院：信息與電子學院同組伙伴：專業(yè)：信息工程成績：實驗1利用DFT分析信號頻譜一、實驗目的1.加深對DFT原理的理解。2.應用DFT分析信號頻譜。3.深刻理解利用DFT分析信號頻譜的原理，分析現實過程現象及解決方法。二、實驗原理1、DFT和DTFT的關系有限長序列的離散時間傅里葉變換在頻率區(qū)間的N個等分點上的N個取樣值可以由下式表示：由上式可知，序列的N點DFT，實際上就是序列的DTFT在N個等間隔頻率點上樣本。2、利用DFT求DTFT方法1：由恢復出的方法如圖2.1所示：圖2.1.由N點DFT恢復頻譜DTFT的流程由圖2.1所示流程圖可知：由式2-2可以得到其中為內插函數方法2：然而在實際MATLAB計算中，上訴插值公式不見得是最好的方法。由于DFT是DTFT的取樣值，其相鄰的兩個頻率樣本點的間距為,所以如果我們增加數據的長度N，使得得到的DFT譜線就更加精細，其包絡就越接近DTFT的結果，這樣可以利用DFT來近似計算DTFT。如果沒有更多的數據，可以通過補零來增加數據長度。3、利用DFT分析連續(xù)時間信號的頻譜采用計算機分析連續(xù)時間信號的頻譜，第一步就是把連續(xù)時間信號離散化，這里需要進行連個操作：一是采樣，二是截斷。對于連續(xù)非周期信號，按采樣間隔T進行采樣，截取長度為M，那么對進行N點的頻率采樣，得到因此，可以將利用DFT分析連續(xù)非周期信號頻譜的步驟歸納如下：〔1〕確定時域采樣間隔T，得到離散序列；〔2〕確定截取長度M，得到M點離散序列，這里的為窗函數。〔3〕確定頻域采樣點數N，要求?！?〕利用FFT計算離散序列的N點DFT，得到?！?〕根據式〔2-6〕由計算采樣點的近似值。采用上訴方法計算的頻譜，需要注意如下三點問題：〔1〕頻譜混疊。如果不滿足采樣定理的條件，頻譜會很出現混疊誤差。對于頻譜無限寬的信號，應考慮覆蓋大局部主要頻率的范圍?！?〕柵欄效應和頻譜分辨率。使用DFT計算頻譜，得到的結果只是N個頻譜樣本值，樣本值之間的頻譜是未知的，就像通過一個柵欄觀察頻譜，稱為“柵欄效應〞。頻譜分辨率與記錄長度成正比，提高頻譜分辨率，就要增加記錄時間?！?〕頻譜泄露。對于信號截斷會把窗函數的頻譜會引入到信號頻譜中，造成頻譜泄露。解決這問題的主要方法是采用旁瓣小的窗函數，頻譜泄露和窗函數均會引起誤差。因此，要合理選取采樣間隔和截取長度，必要時還需考慮適當的窗。對于連續(xù)周期信號，我們在采用計算機進行計算時，也總是要進行截斷，序列總是有限長的，仍然可以采用上訴方法近似計算。4、可能用到MATLAB函數與代碼實驗中的DFT運算可以采用MATLAB中提供的FFT來實現。DTFT可以利用MATLAB矩陣運算的方法進行計算。三、實驗內容1.，完成如下要求：〔1〕計算其DTFT，并畫出區(qū)間的波形。〔2〕計算4點DFT，并把結果顯示在〔1〕所畫的圖形中。〔3〕對補零，計算64點DFT，并顯示結果?！?〕是否可以由DFT計算DTFT，如果可以，請編程實現。2.考察序列〔1〕時，用DFT估計的頻譜；將補零加長到長度為100點序列用DFT估計的頻譜。要求畫出相應波形?！?〕時，用DFT估計x(n)的頻譜，并畫出波形。3.信號，其中，，。從的表達式可以看出，它包含三個頻率的正弦波，但是，從其時域波形來看，似乎是一個正弦信號，利用DFT做頻譜分析，確定適合的參數，使得到的頻譜的頻率分辨率符合需要。4.利用DFT近似分析連續(xù)時間信號xt=e-0.1tu(t)的頻譜〔四、實驗代碼及實驗結果實驗實驗結果：實驗代碼：>>n=0:3;>>x=[2-111];>>w=-pi:0.01*pi:pi;>>X=x*exp(-j*n'*w);>>subplot(211);>>plot(w,abs(X));>>title('幅度');xlabel('w');ylabel('|X|');>>axistight;>>subplot(212);>>plot(w,angle(X));>>title('相位');xlabel('w');ylabel('Angle〔X〕');>>axistight;實驗實驗結果：實驗代碼：n=0:3;x=[2-111];w=-pi:0.01*pi:pi;X=x*exp(-j*n'*w);subplot(211);plot(w,abs(X));title('幅度');xlabel('w');ylabel('|X|');axistight;holdon;subplot(212);plot(w,angle(X));title('相位');xlabel('w');ylabel('Angle〔X〕');axistight;holdon;H=fft(x);subplot(211);>>stem(n,abs(H),'filled');>>subplot(212);>>stem(n,angle(H),'filled');實驗實驗結果：實驗代碼：>>x=[2-111zeros(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),'filled');title('相位');xlabel('n');ylabel('angle(X)');分析：可以由DFT計算DTFT。通過補零加長序列，提高采樣密度，可以由DFT近似計算DTFT。實驗實驗結果：實驗代碼：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|');subplot(212);stem(n,angle(X),'filled');title('相位');xlabel('n');ylabel('angle(X)');補零加長：實驗結果：實驗代碼：h=[xzeros(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),'filled');title('相位');xlabel('n');ylabel('angle(H)');實驗實驗結果：實驗代碼：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,angle(X),'filled');title('相位');xlabel('n');ylabel('angle(X)');分析：可以通過增大截取長度和增加補零的個數來提高頻譜分辨率，但是補零不能夠增加分辨力。實驗1.3實驗結果：實驗代碼：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,'filled');axis([0,4,0,1]);xlabel('f/Hz');title('Magnitude');分析：通過選取適宜的采樣周期，可以完整的恢復出原信號的頻譜波形。實驗1.4Tp=1s，0《n《100實驗結果：實驗代碼：n=0:1:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=5s,0《n《100實驗結果：實驗代碼：n=0:5:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=25s，0《n《100實驗結果：實驗代碼：n=0:25:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=1s，0《n《50實驗結果：實驗代碼：n=0:1:50;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=1,0《n《25實驗結果：實驗代碼：n=0:1:25;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');分析：最終確定參數：tp=5s，0《n《100。五、心得與體會通過本次實驗，我們掌握并加深對DFT原理的理解并且學會應用DFT分析信號頻譜，在此根底上利用DFT分析信號頻譜的原理，掌握了利用matlab分析現實問題的步驟及方法。并且，通過這次的實驗，對信號序列有了更加深刻的認識，單純的一個信號序列是沒有意義的，只有配合他本身的時間序列才是一個完整的信號序列，才可以對其進行分析。實驗2利用FFT計算線性卷積〔選作〕一、實驗目的1.掌握利用FFT計算線性卷積的原理及具體實現方法。2.加深理解重疊相加法和重疊保存法。3.考察利用FFT計算線性卷積各種方法的適用范圍。二、實驗原理1.線性卷積與圓周卷積設x(n)為L點序列，h(n)為M點序列，x(n)和h(n)的線性卷積為〔3-1〕的長度為L+M-1x(n)和h(n)的圓周卷積為〔3-2〕圓周卷積與線性卷積相等而不產生交疊的必要條件為N≥L+M+1(3-3)圓周卷積定理：根據DFT性質，x(n)和h(n)的N點圓周卷積的DFT等于它們的DFT的乘積：(3-4)2.快速卷積快速卷積發(fā)運用圓周卷積實現線性卷積，根據圓周卷積定理利用FFT算法實現圓周卷積?？蓪⒖焖倬矸e運算的步驟歸納如下：(1)必須選擇；為了能使用基-2算法，要求。采用補零的方法使得x(n)和h(n)的長度均為N。(2)計算x(n)和h(n)的N點FFT。(3)組成乘積(4)利用IFFT計算Y(k)的IDFT，得到線性卷積y(n)3.分段卷積我們考察單位取樣響應為h(n)的線性系統(tǒng)，輸入為x(n)，輸出為y(n)，則y當輸入序列x(n)極長時，如果要等x(n)全部集齊時再開始進行卷積，會使輸出有較大延時；如果序列太長，需要大量存儲單元。為此，我們把x(n)分段，為別求出每段的卷積，合在一起得到最后的總輸出。這稱為分段卷積。分段卷積可以細分為重疊保存法和重疊相加法。重疊保存法：設x(n)的長度為，h(n)的長度為M。把序列x(n)分成多段N點序列，每段雨前一段重寫M-1個樣本。并在第一個輸入段前面補M-1個零。計算每一段與h(n)的圓周卷積，其結果中前M-1個不等與線性卷積，應當舍去，只保存后面N-M+1個正確的輸出樣本，把它們合起來得到總的輸出。利用FFT實現重疊保存法的步驟如下：(1)在x(n)前面填充M-1個零，擴大以后的序列為(2)將x(n)分為假設干段N點子段，設L=N-M+1為每一段的有效長度，則第i段的數據為：(3)計算每一段與h(n)的N點圓周卷積，利用FFT計算圓周卷積(4)舍去每一段卷積結果的前M-1個樣本，連接剩下的樣本得到卷積結果y(n)。重疊相加法：設h(n)長度為M，將信號x(n)分解成長為L的子段。以表示沒斷信號，則：每一段卷積的長度為L+M-1，所以在做求和時，相鄰兩段序列由M-1個樣本重疊，即前一段的最后M-1個樣本和下一段前M-1個樣本序列重疊，這個重疊局部相加，再與不重疊的局部共同組成y(n)。利用FFT實現重疊保存法的步驟如下：(1)將x(n)分為假設干L點子段。(2)計算每一段與h(n)的卷積，根據快速卷積法利用FFT計算卷積。(3)將各段相加，得到輸出y(n)。4、可能得到的MATLAB函數實驗中FFT運算可采用MATLAB中提供的函數fft來實現。三、實驗內容假設要計算序列x(n)=u(n)-u(n-L),0≤n≤L和h(n)=cos(0.2πn),0≤n≤M的線性卷積完成以下實驗內容。1.設L=M,根據線性卷積的表達式和快速卷積的原理分別編程實現計算兩個序列線性卷積的方法,比擬當序列長度分別為8,16,32,64,256,512,1024時兩種方法計算線性卷積所需時間。2．當L=2048且M=256時比擬直接計算線性卷積和快速卷積所需的時間,進一步考察當L=4096且M=256時兩種算法所需的時間。3.編程實現利用重疊相加法計算兩個序列的線性卷積,考察L=2048且M=256時計算線性卷積的時間,與2題的結果進行比擬。4.編程實現利用重疊保存法計算兩個序列的線性卷積,考察L=2048且M=256時計算線性卷積的時間,與2題的結果進行比擬。四、實驗代碼及實驗結果實驗2.1實驗代碼：forM=[81632642565121024];L=M;n=0:1:L-1;x=ones(1,L);h=cos(0.2.*pi.*n);disp('L=M=');disp(M)disp('線性卷積用時：')ticy1=conv(x,h);tocdisp('快速卷積用時：');ticX1=fft(x);H=fft(h);Y=X1.*H;y=ifft(Y);tocend實驗結果：L=M=8線性卷積用時：時間已過0.000122秒?？焖倬矸e用時：時間已過0.000042秒。L=M=16線性卷積用時：時間已過0.000045秒?？焖倬矸e用時：時間已過0.000026秒。L=M=32線性卷積用時：時間已過0.000051秒。快速卷積用時：時間已過0.000030秒。L=M=64線性卷積用時：時間已過0.000050秒?？焖倬矸e用時：時間已過0.000030秒。L=M=256線性卷積用時：時間已過0.000089秒。快速卷積用時：時間已過0.000066秒。L=M=512線性卷積用時：時間已過0.000540秒。快速卷積用時：時間已過0.001388秒。L=M=1024線性卷積用時：時間已過0.000617秒?？焖倬矸e用時：時間已過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('線性卷積：')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實驗結果：L=2048，M=256線性卷積：時間已過0.000437秒。快速卷積：時間已過0.000407秒。L=4096，M=256線性卷積：時間已過0.000445秒?？焖倬矸e：時間已過0.001395秒。分析：當序列長度較短時，快速卷積比線性卷積快；當序列長度過長時，快速卷積比線性卷積慢了許多。實驗2.3實驗代碼：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);disp('重疊相加法卷積：');ticH=fft(h,Ni);fori=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:i_low+Ni-1)+yi;endtoc實驗結果：L=2048M=256重疊相加法卷積時間時間已過0.000626秒。分析：較第2題結果快。實驗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('??μt±￡á?·¨?í?yê±??');ticH=fft(h,Ni);fori=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_yu)=yi(M:Ni-1);toc實驗結果：L=2048M=256重疊保存法卷積時間時間已過0.000508秒。分析：較第2題結果快。五、心得與體會本次實驗要求我們掌握利用FFT計算線性卷積的原理及具體實現方法，以及通過對快速卷積和線性卷積運算速度的比擬，更加直觀的看到利用FFT快速卷積的優(yōu)點。本次實驗讓我切實看到了FFT算法的高效性及對于龐大的數據的處理實時性，同時對重疊保存法、重疊相加法也有了更深的認識，對以后的實驗有了很好的理論根底，受益頗多。實驗3IIR數字濾波器設計一、實驗目的1.掌握利用脈沖響應不變法和雙線性變換法設計IIR數字濾波器的原理及具體方法。2.加深理解數字濾波器和模擬濾波器之間的技術指標轉化。3.掌握脈沖響應不變法和雙線性變換法設計IIR數字濾波器的優(yōu)缺點及適用范圍。二、實驗內容1.設采樣頻率為fs=4kHz，采用脈沖響應不變法設計一個三階巴特沃斯數字低通濾波器，其3dB截止頻率為fc=1kHz。2.設采樣頻率為fs=10kHz，設計數字低通濾波器，滿足如下指標通帶截止頻率：fp=1kHz，通帶波動：Rp=1dB阻帶截止頻率：fst=1.5kHz，阻帶衰減：As=15dB要求分別采用巴特沃斯、切比雪夫I型、切比雪夫II型和橢圓模擬原型濾波器及脈沖響應不變法進行設計。結合實驗結果，分別討論采用上述設計的數字濾波器是否都能滿足給定指標要求，分析脈沖響應不變法設計IIR數字濾波器的優(yōu)缺點及適用范圍。四、實驗代碼及實驗結果巴特沃斯：實驗代碼：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,az);%頻率響應subplot(221);plot(w/pi,abs(H));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');grd=grpdelay(bz,az,w);%群延時響應subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');實驗結果：bz=1.0e-04*-0.00000.00910.20330.44630.15240.00510az=1.0000-5.137811.0528-12.73838.2921-2.88980.4211切比雪夫Ⅰ型：實驗代碼：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(A^2-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));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');實驗結果：bz=1.0e-03*0.00000.06230.23720.05670az=1.0000-3.77105.3742-3.42930.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)-1);

A=10^(As/20);

N=ceil(acosh(sqrt(A^2-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));

gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');

subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));

gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');

subplot(223);plot(w/pi,angle(H)/pi);

gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');

grd=grpdelay(bz,az,w);

subplot(224);plot(w/pi,grd);

gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');實驗結果：bz=0.0884-0.33420.5117-0.37180.1076az=1.0000-3.49044.6022-2.71480.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(A^2-1)/e)/(acosh(Ws/Wp)));[b,a]=ellip(N,Rp,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));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');實驗結果：bz=0.1470-0.56360.8366-0.56az=1.0000-3.77265.3933-3.46040.8408分析：脈沖響應不變法的優(yōu)點是頻率坐標的變換是線性的，因此如果模擬濾波器的頻響是限帶于折疊頻率以內的話，通過變換后的數字濾波器的頻響可以不失真地反映原響應與頻率之間的關系。但是，其最大的缺點是有頻譜的周期延拓效應，會產生混疊失真。故，脈沖響應不變法只能用于限帶的頻響特性，如衰減特性較好的低通，或帶通，而且高頻衰減越大，頻響的混疊效應就越小。雙線性變換法的優(yōu)點是消除了脈沖響應不變法的固有的混疊失真，但缺點是頻率變換之間的非常嚴重的非線性，會使變換后時域上的圖像產生非常嚴重的失真。五、心得與體會通過實驗，學習了利用脈沖響應不變法設計IIR數字濾波器的根本原理為：從時域響應出發(fā)，使數字濾波器的單位脈沖響應h(n)模仿模擬濾波器的單位沖激響應ha(t)，h(n)等于脈沖響應不變法和雙線性法設計IIR數字濾波器設計的根底，掌握了這倆種方法能夠使我們對數字濾波器設計有更加深刻的認識。實驗4頻率取樣法設計FIR數字濾波器一、實驗目的掌握頻率取樣法設計FIR數字濾波器的原理及具體方法二、實驗原理1、根本原理頻率取樣法從頻域出發(fā)，把理想的濾波器等間隔采樣得到,將作為實際設計濾波器的。〔8-1〕得到以后可以由來確定唯一確定濾波器的單位脈沖響應,可以由求得：〔8-2、3〕其中為內插函數〔8-4〕由求得的頻率響應將逼近如果我們設計的是線性相位FIR濾波器，則的幅度和相位滿足線性相位濾波器的約束條件。我們將表示為如下形式：〔8-5〕當為實數，則由此得到〔8-6〕即為中心偶對稱。在利用線性相位條件可知，對于1型和2型線性相位濾波器〔8-7〕對于3型和4型線性相位濾波器〔8-8〕其中，x設計步驟〔1〕由給定的理想濾波器給出和?！?〕由求得〔3)根據求得或三、實驗內容1.設計一個數字低通FIR濾波器，其技術指標如下：ω_p=0.2π，R_p=0.25dBω_st=0.3π，A_s=50dB分別采用矩形窗、漢寧窗、海明窗、布萊克曼窗、凱瑟窗設計該濾波器。結合實驗結果分別討論上述方法設計的數字濾波器是否符合指標。2.設計一個數字帶通FIR濾波器，其技術指標如下：3.采用頻率取樣設計法設計FIR數字低通濾波器，滿足以下指標ω_p=0.2π，R_p=0.25dBω_st=0.3π，A_s=50dB〔1〕取N=20，過渡帶沒有樣本。〔2〕取N=40，過渡帶有一個樣本，T=0.39?！?〕取N=60，過渡帶有兩個樣本，T1=0.5925,T2=0.1099?！?〕分別采用上述方法設計的數字濾波器是否都能滿足給定的指標要求。4.采用頻率取樣技術設計下面的高通濾波器ω_st=0.6π，A_s=50dBω_p=0.8π，R_p=1dB對于高通濾波器，N必須為奇數〔或1型濾波器〕。選擇N=33，過渡帶有兩個樣本，過渡帶樣本最優(yōu)值為T1=0.1095，T2=0.598。四、實驗代碼及實驗結果實驗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=hd.*w_boxcar;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/pi,Hr);axistight;xlabel('\omega/\pi');ylabel('H(\omega)');subplot(223);stem(n,h,'filled');axistight;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');gridon;實驗結果：分析：阻帶的衰減小于50db，所以，不符合指標。漢寧窗：實驗代碼：wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp;N=ceil(6.2*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_hanning=hanning(N)';h=hd.*w_hanning;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/pi,Hr);axistight;xlabel('\omega/\pi');ylabel('H(\omega)');subplot(223);stem(n,h,'filled');axistight;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');gridon;實驗結果：分析：阻帶衰減小于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');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/pi,Hr);axistight;xlabel('\omega/\pi');ylabel('H(\omega)');subplot(223);stem(n,h,'filled');axistight;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');gridon;實驗結果：分析：阻帶衰減大于50db，故，該數字濾波器符合指標。布萊克曼窗：實驗代碼：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)*sinc((wc/pi)*(n-alpha));w_blackman=blackman(N)';h=hd.*w_blackman;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/pi,Hr);axistight;xlabel('\omega/\pi');ylabel('H(\omega)');subplot(223);stem(n,h,'filled');axistight;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');gridon;實驗結果：分析：阻帶衰減大于50db，故，該數字濾波器符合要求。凱瑟窗：實驗代碼：wp=0.2*pi;wst=0.3*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');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/pi,Hr);axistight;xlabel('\omega/\pi');ylabel('H(\omega)');subplot(223);stem(n,h,'filled');axistight;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');gridon;實驗結果：分析：阻帶衰減大于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;hd=(wc2/pi)*sinc((wc2/pi)*(n-alpha))-(wc1/pi)*sinc((wc1/pi)*(n-alpha));w_blackman=blackman(N)';h=hd.*w_blackman;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/pi,Hr);axistight;xlabel('\omega/\pi');ylabel('H(\omega)');subplot(223);stem(n,h,'filled');axistight;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');gridon;實驗結果：實驗代碼：N=20;alpha=(N-1)/2;l=0:N-1;wl=(2*pi/N)*l;Hrs=[1,1,1,zeros(1,15),1,1];Hdr=[1,1,0,0];wdl=[0,0.25,0.25,1];k1=0:floor((N-1)/2);k2=floor((N-1)/2)+1:N-1;angH=[-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)];H=Hrs.*exp(j*angH);h=ifft(H,N);w=[0:500]*pi/500;H=freqz(h,1,w);[Hr,wr]=zerophase(h);subplot(221);plot(wdl,Hdr,wl(1:11)/pi,Hrs(1:11),'o');axis([0,1,-0.1,1.1]);xlabel('\omega(\pi)');ylabel('Hr(k)');subplot(222);stem(l,h,'filled');axis([0,N-1,-0.1,0.3]);xlabel('n');ylabel('h(n)');subplot(223);plot(wr/pi,Hr,wl(1:11)/pi,Hrs(1:11),'o');axis([0,1,-0.2,1.2]);xlabel('\omega(\pi)');ylabel('Hr(\omega)');subplot(224);plot(w/pi,20*log10((abs(H)/max(abs(H)))));axis([0,1,-50,5]);grid;xlabel('\omega(\pi)');ylabel('dB');實驗結果：實驗實驗代碼：N=40;alpha=(N-1)/2;l=0:N-1;wl=(2*pi/N)*l;Hrs=[1,1,1,1,1,0.39,zeros(1,29),0.39,1,1,1,1];Hdr=[1,1,0,0];wdl=[0,0.25,0.25,1];k1=0:floor((N-1)/2);k2=floor((N-1)/2)+1:N-1;angH=[-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)];H=Hrs.*exp(j*angH);h=ifft(H,N);w=[0:500]*pi/500;[H,w]=freqz(h,1,w);Hr=abs(H);subplot(221);plot(wdl,Hdr,wl(1:21)/pi,Hrs(1:21),'o');axis([0,1,-0.1,1.1]);xlabel('\omega(\pi)');ylabel('Hr(k)');subplot(222);stem(l,h,'filled');axis([0,N-1,-0.1,0.3]);xlabel('n');ylabel('h(n)');subplot(223);plot(w/pi,Hr,wl(1:21)/pi,Hrs(1:21),'o');axis([0,1,-0.2,1.2]);xlabel('\omega(\pi)');ylabel('Hr(\omega)');subplot(224);plot(w/pi,20*log10((abs(H)/max(abs(H)))));axis([0,1,-50,5]);grid;xlabel('\omega(\pi)');ylabel('dB');實驗結果：實