1、濾波器設計課程設計題目一、FIR（選取合適的窗函數及窗長度）通帶06000Hz通帶波動士0.01阻帶8500HzFs/2阻帶波動-70dB抽樣頻率44100Hz二、IIR（可選巴特沃斯、切比雪夫I型、切比雪夫II型）1、設計數字低通濾波器，滿足以下指標：通帶邊帶頻率2000Hz通帶波動0.01阻帶邊界頻率2700Hz11025Hz阻帶波動0.05抽樣頻率要求:1、求出h(n)和H(ej)4a)用H(ej)直接畫出幅度響應曲線(dB)及相位響應曲線；(b)做h(n)的FFT,有此畫出幅度響應曲線(dB)及相位響應曲線；(c)并兩種方法的結果進行比較。(d)驗證你的設計結果滿足設計指標。FIR和I

2、IR濾波器設計過程FIR濾波器的設計：1.設計步驟:(1)寫出理想低通濾波器的傳輸函數和單位脈沖響應。?(n=H?/?>?=?2?.-?1*sin?夕?？一、2?-?(2)根據題目已有的模擬濾波器頻率進行歸一化求得數字濾波器的頻率，并對題目所給的通帶波動和阻帶波動進行轉換，可求得通帶截止頻率？=0.8549?阻帶截止頻率？?=1.2110?通帶波動？?=0.174?阻帶衰減？?=70?(3)選擇用布萊克曼窗函數來設計FIR濾波器，通過查閱資料可知道布萊克曼11兀，、,，窗的過渡市范為,最小我減為74?則可以通過matlab中ceil函數向上取整的方式以及阻帶頻率和通帶頻率之差所求的過渡帶

3、頻率來確定濾波器的長度，在經過調試可以取得N=99。當然，在設計中可以之間用blackman()這個函數,但是在這里我們會利用布萊克曼窗函數的公式構造布萊克曼窗：2?4?wn=0.42-0.5cos而%+0.08cos與“?(?？(4)通過求得的矩形窗函數和單位脈沖響應？?？n進行時域上的乘積即可求得FIR濾波器的時域響應，再利用DTFT直接求出H(?X并且可以畫出幅度響應曲線(dB)及相位響應曲線OOX?嚴?=DTFTxn=?=-CO(5)根據時域所求的函數做FFT并由此畫出幅度響應曲線(dB)及相位響應曲線，對這兩種方法的結果進行比較并驗證設計結果是否滿足設計指標。3.設計過程：(1)求出

4、h(n),并做h(n)的DTFT從而直接求出H(?的幅度曲線響應(dB)和相位響應曲線。由設計步驟中可以知道？?？n=嗎箸以及自己通過公式所構造的布萊克曼窗函數即可求得h(n)。然后將h(n)經過DTFT即可求得H(?,進而求得幅度曲線響應和相位響應曲線。波形如下：0.90.BD1D6邕0三D.40.30.2卜廣9/1j1VItRPH*11Vj/L-X尸11-_11020304Q6060703D90100Hi圖一布萊克曼窗函數的波形圖二理想濾波器沖擊響應hdn的波形可|JJJ1L1L11國4口暹0-1001020304050n圖三h(n)的波形?的前100個值0.00010.0003Colum

5、ns1through13-0.0000-0.0000-0.0000-0.00000.00000.0000-0.0001-0.0002-0.00010.00020.0004Columns14through26-0.0003-0.0008-0.00050.00050.00140.0010-0.0008-0.0023-0.00160.00120.00360.0026-0.0016Columns27through39-0.0053-0.00410.00210.00770.0062-0.0026-0.0109-0.00930.00310.01540.0139-0.0035-0.0220Columns40

6、through52-0.02110.00390.03270.0340-0.0043-0.0549-0.06480.00450.13910.27290.32880.27290.1391Columns53through650.0045-0.0648-0.0549-0.00430.03400.03270.0039-0.0211-0.0220-0.00350.01390.01540.0031Columns66through78-0.0093-0.0109-0.00260.00620.00770.0021-0.0041-0.0053-0.00160.00260.00360.0012-0.0016Colu

7、mns79through91-0.0023-0.00080.00100.00140.0005-0.0005-0.0008-0.00030.00030.00040.0002-0.0001-0.0002Columns92through99-0.00010.00000.00010.0000-0.0000-0.0000-0.0000-0.0000山P郛祖垣百順蹙31圖四H(?的幅度響應曲線圖四H(?的相位響應曲線(2)根據時域所求的函數h(n)做FFT并由此畫出幅度響應曲線(dB)及相位響應曲線嗎器竊冏誓聯譬1111rA0701vn舊I111ll1111iii1I1一一!i一J0,3856-乃dBl1

8、ii1j:-:|-i1II!和_11111II14ii!1i1r|Jiffii>1>iiH1i11111111圖五H(?的幅度響應曲線圖六H(?的相位響應曲線由題目可知濾波器的指標：通帶波動為0.174dB,阻帶衰減為-70dB,通帶截止頻率為0.8549Hz,阻帶截止頻率為1.2110Hz。通過DTFT直接求H(?所得到的FIR濾波器的指標是：通帶波動為0.001dB,阻帶衰減為-73.14dB。通過FFT求H(?7所得到的FIR濾波器的指標是：通帶波動為0.0014dB,阻帶衰減為-75dB1.51noci3"1j111II1111mi1h£u1£

9、；111JQJ。名11.2141.6132w?Hz圖七H（?了的幅度響應曲線（|H?/?|）由圖六可以知道這個濾波器兩個正負尖峰頻率差值為：1.2612-0.80381=0.45739Hz設計的的FIR濾波器的過渡帶寬為:1.2110-0.8549=0.3561Hz.I1!1!，4，4一1I*41-1-I-.II-R4riaJ.4Bu4BJ.lld.n.ril->»»»»sKMIW<1*-20r-g亦回加頓g例皿吟+HFT來第應的曲為卜卜，曲珀應的曲線圖八DTFT和FFT所求的幅度響應（dB）曲線由該圖可以看出DTFT所求幅度響應（dB）的曲

10、線和FFT所求的幅度響應(dB)的曲線是基本一致的D1TT求啊應的曲線OFFT未晌應的自裝一，口0050101502025030350404505,JV/p4圖八DTFT和FFT所求的相位響應曲線綜上所述：1.所設計的布萊克曼窗FIR濾波器在阻帶衰減上接近了-70(dB)的指標，但是在通帶波動上0.02(dB)略小于所要求的0.174(dB)的指標。2.FIR濾波器的過渡帶寬小于兩個正負尖峰頻率差值，滿足要求。3.DTFT和FFT兩種方法所設計的濾波器結果及指標基本沒有區別，但是DTFT的計算量會比FFT相對大，而FFT所求的幅度和相位響應的波形是對稱的，需要經過處理才能達到要求。4.DTFT

11、所算出來的指標會相對來說會更接近于題目中要求的指標。IIR濾波器的設計：1.設計步驟：(1)選擇采用沖激響應不變法設計巴特沃斯低通濾波器，已知通帶邊帶頻率為2000Hz,阻帶邊帶頻率為2700Hz,通帶波動為0.01,阻帶波動為0.05,抽樣頻率為11025Hz,根據對頻率的歸一化以及對阻帶波動通帶波動轉換成(dB)的形式可得，Wp=1.1398Hz,Ws=1.5387Hz,Rp=0.087(dB),As=26(dB)。(2)計算樣本模擬濾波器的階次100.1?3?-1?(lg100.1?.1)/2lg壇？可求得N=17o(3)發現Rp與(dB),?月?，由于巴特沃斯濾波器歸一化低通原型的通帶

12、截止頻率？?=1,去歸一化成任意截止頻率時必須是3(dB)衰減處的？為利用通帶或者阻帶滿足衰減指標的公式去求得3(dB)衰減的截止頻率。?)？/?100.1?3?_12?”?/100.1?-1在這里選擇的是阻帶滿足衰減指標的公式，可求得？?=3.8003e+03。(4)根據？?利階數N求模擬低通濾波器的系統函數分子分母系數。(5)沖激響應不變法利用模擬低通濾波器的系統函數和采樣頻率轉換成數字低通濾波器的系統函數。(6)利用數字低通濾波器的系統函數的分子分母系數求得沖激響應即數字低通濾波器的時域響應。(7)根據所求的時域響應作FFT以求得在頻率上的響應。(8)檢測步驟(5)中的求得的模擬濾波器系

13、統函數分子分母系數以及數字濾波器系統函數分子分母系數來檢驗衰減指標。3.設計過程：(1)將要求的各種指標代入巴特沃斯函數中以求得階數和3(dB)衰減的截止頻率，并由此確定模擬濾波器系統函數的分子分母系數，再利用impinvar將模擬濾波器系統函數的分子分母系數轉換成數字濾波器系統函數的分子分母系數。然后利用freqz()函數求出數字濾波器頻域上的幅度響應曲線(dB)和相位響應曲線。i聃字汨轉器帽陰20萼髭祖歸等螞強OoQQQQ口4&凜，K4J1-1144-圖一H(?了的幅度響應曲線II1!,飛.iB-_.V-一I1i11o我字派技界珅吊圖二H(?的相位響應曲線(2)將(1)中所求的數字

14、濾波器系統函數的分子分母的系數利用imz函數即可求得系統的沖激響應即時域函數h(n)。06電也掃u珀應圖三h(n)的波形?的前100個值Columns1through13-0.0000-0.00000.00000.00000.00040.00510.03090.10830.24340.36100.33520.1315-0.1050Columns14through26-0.1764-0.05330.09240.0972-0.0112-0.0796-0.03560.04150.0494-0.0062-0.0414-0.01650.0239Columns27through390.0250-0.006

15、0-0.0226-0.00670.01440.0124-0.0050-0.0123-0.00220.00860.0060-0.0037-0.0066Columns40through52-0.00040.00500.0028-0.0026-0.00350.00030.00290.0012-0.0017-0.00180.00040.00160.0005Columns53through65-0.0010-0.00090.00040.00090.0002-0.0006-0.00040.00030.00050.0000-0.0004-0.00020.0002Columns66through780.000

16、3-0.0000-0.0002-0.00010.00010.0001-0.0000-0.0001-0.00000.00010.0001-0.0000-0.0001Columns79through91-0.00000.00000.0000-0.0000-0.000000000000Columns92through1040000000000000(3)根據時域所求的函數h(n)做FFT并由此畫出幅度響應曲線(dB)及相位響應曲線。圖四H(?的幅度響應曲線圖五H(?的相位響應曲線由題目可知濾波器的指標：通帶波動為0.087dB,阻帶衰減為26dB,通帶截止頻率為1.1398Hz,阻帶截止頻率為1.5

17、387Hz。由于求利用impz函數求h(n)只是數字濾波器的單位沖激響應逼近模擬濾波器的單位沖激響應的過程，再進行FFT變換所得到的數字濾波器的指標誤差相對較大，所以在這里只檢驗直接求頻域響應所得到的數字濾波器的指標。利用書上檢驗衰減指標的代碼并將參數代入可求得：通帶波動為0.0634dB,阻帶衰減為26dB,非常接近題目中要求的指標。20-T-1IIIVI直接啊G|-FFT/由一的舊線，圖六直接求和FFT所求的幅度響應(dB)曲線京田團灣吊_收圖七直接求和FFT所求的相位響應曲線由以上兩幅圖可以知道直接求所得到的幅度響應曲線及相位響應曲線在w/pi=00.7這個范圍內的波形基本保持一致，但是

18、在w/pi>0.7之后FFT所得到的幅度響應的曲線的諧波較多，而且呈現的也不是線性相位。綜上所述，直接求H(?3所設計的低通濾波器基本能夠滿足題目要求的指標，在比較了直接求H(?和FFT設計的數字低通濾波器后，認為在設計低通巴特沃斯濾波器上還是用直接求H(?的方式較為精確。設計過程中遇到的問題和解決方法：1 .FIR濾波器中在求窗長度的時候11pi/(ws-wp)得不到整數。解決方法：用ceil函數向上取整的方式得到整數。2 .在做FFT變換的時候發現結果是呈對稱的形式而且曲線呈現的波形不夠完整。解決方法：首先將波形的范圍延長一倍，然后采用想要的x的范圍去限制matlab畫出的y的波形，

19、就可以得到想要的一般的波形。3 .如何克服矢量維度不一致的問題解決方法：先看兩個矢量的維度是多少，然后可以在那個維度少的矢量后邊用補零函數進行調整使得兩個函數的矢量一致。心得體會：通過這次設計FIR濾波器和IIR濾波器的經歷，我真心覺得自己的matlab水平真的有了一個質的飛躍，我承認網絡上有很多類似的資料，但是要真正達到題目的要求還是需要自己動腦去做的。感覺這個過程就像是建立模型和代碼實現的過程，我也是按照書本上兩個濾波器的模型來設計并且實現的，整個過程發現了不少的問題，有些是關于課程的問題，比如加窗類型的選擇以及調試，在達到過渡帶小于正負尖峰之間那個指標的過程真的是調了很久，最后才確立了布

20、萊克曼窗，參數也是改了好多遍。有些是關于matlab的基礎問題，因為我接觸matlab的時間不多，一開始比較吃力，但是我覺得經過了這個項目之后確實對matlab有了進一步的認識，其實挺希望這個項目能稍微增加難度然后放在課程一開始讓學生研究，我覺得這樣的效果會更好。程序源代碼：FIR濾波器closeallf1=6000;%通帶截止頻率f2=8500;%阻帶截止頻率f=44100;%抽樣頻率wp=2*pi*f1/f%歸一化ws=2*pi*f2/falphap=0.174;%通帶波動dBalphas=70;%阻帶衰減dBc=11*pi;%Blackman窗過渡帶寬wc=(ws+wp)/2;delt_

21、w=ws-wp;M=ceil(c/delt_w);N=1*M+1%窗長度n=0:M;%利用課本公式得到布萊克曼窗win=(0.42-0.5*cos(2*pi*n)/(N-1)+0.08*cos(4*pi*n)/(N-1);win1=win.'figure(1)n=-0.5*M:0.5*M;%低通濾波器的時域響應hd=sin(wc*n)./(pi*n);hd(find(n=0)=wc*cos(wc*0)/pi;ht=hd.*win1'%加窗過程plot(n,ht,'.-')%畫沖激響應系數圖xlabel('n','FontSize'

22、,12);ylabel('ht','fontsize',12);gridonfigure(2)%h,w=freqz(ht,1,1024);%離散時間傅里葉變換得到幅度%W=w/pi;H=20*log10(abs(h);%holdon%plot(W,H);k=1:1024;n=0:M;w=(pi/1024)*k;W=w/pi;h=ht*(exp(-j*pi/1024)A(n'*k);%DTFT變換plot(w/pi,20*log10(abs(h);xlabel('w/pi');ylabel('幅度響應曲線dB');%直接求所

23、得到的頻率幅度響應gridonholdonline(0,1,-70,-70,'color','r');%標出相應的限制范圍line(0,1,-0.174,-0.174,'color','r');dotp=round(mean(find(w>wp-0.1&w<wp+0.1)%找到點wp的坐標dots=round(mean(find(w>ws-0.1&w<ws+0.1)%找到點ws的坐標plot(W(dotp),H(dotp),'.r','MarkerSize',

24、10);%對應處打上顯目的紅點plot(W(dots),H(dots),'.r','MarkerSize',10);text(W(dotp),H(dotp),'','',num2str(round(W(dotp)*1000)/1000),.',',num2str(round(H(dotp)*1000)/1000),'dB','FontSize',15);%標示wp在頻響圖%中的位置，從中可以看來相應增益text(W(dots),H(dots),'',''

25、;,num2str(round(W(dots)*100)/100),.',',num2str(round(H(dots)*100)/100),'dB','fontsize',15);%標示ws在頻響圖%的位置holdofffigure(3)plot(W,angle(h);axis(0,0.5,-4,4);xlabel('w/pi');ylabel('相位響應曲線');%直接求所得到的相位響應曲線gridonfigure(4)y=20*log10(abs(fft(ht,1024);%FFT求得的頻域響應x=(0:10

26、23)/(512);y1=y(x>=0&x<=1);plot(x(x>=0&x<=1),y1);%對所求的對稱圖形進行處理xlabel('w/pi');ylabel('幅度響應曲線dB');%FFT求得的頻域響應曲線gridonholdon;plot(0.3856,-75,'.r','MarkerSize',10);%標出plot(0.2722,-0.0014,'.r','MarkerSize',10);text(0.3856,-75,'0.3856&

27、#39;,'-75','dB','FontSize',15);text(0.2722,-0.0014,'0.2722','-0.0014','dB','FontSize',15);line(0,1,-70,-70,'color','r');%標出相應的限制范圍line(0,1,-0.174,-0.174,'color','r');figure(5)y2=angle(fft(ht,1024);x=(0:1023)/(512

28、);y3=y2(x>=0&x<=1);plot(x(x>=0&x<=1),y3);axis(0,0.5,-4,4);xlabel('w/pi');ylabel('相位響應曲線');%FFT求得的相位響應曲線gridonfigure(10)plot(W,H,'*r',x(x>=0&x<=1),y1,'ob');legend('DTFT求響應的曲線','FFT求響應的曲線);將兩種方法曲線畫在圖上比較xlabel('w/pi');yla

29、bel('幅度響應曲線dB');figure(11)plot(W,angle(h),'.r',x(x>=0&x<=1),y3,'ob');axis(0,0.5,-4,4);legend('DTFT求響應的曲線','FFT求響應的曲線);將兩種方法曲線畫在圖上比較xlabel('w/pi');ylabel('相位響應曲線');figure(12)plot(w,abs(h);%標出幅度曲線在正負尖峰的兩個點，檢測是否滿足過渡帶的要求holdonaxis(0,2,-0.5,1.

30、5);plot(0.8038,1,'.r','MarkerSize',10);plot(1.261,2.88e-05,'.r','MarkerSize',10);text(0.8038,1,'0.8038','1','','FontSize',15);text(1.261,2.88e-05,'1.261','2.88e-05','','FontSize',15);xlabel('w/Hz'

31、);ylabel('H幅度曲線')；figure(13)plot(n,win1);%畫出窗函數的圖像gridon;xlabel('n');ylabel('w(n)');figure(14)n=0:M;plot(n,hd);%畫出理想濾波器沖激響應的圖像gridon;xlabel('n');ylabel('hd(n)');IIR濾波器clc;clearallOP=2*pi*2000;%模擬通帶截止頻率OS=2*pi*2700;%模擬阻帶截止頻率Rp=0.087;%通帶波動dBAs=26;%阻帶衰減dBFs=11025

32、%采樣頻率Wp=OP/Fs;%歸一化Ws=OS/Fs;N,OC=buttord(OP,OS,Rp,As,'s');%求出階數和3dB衰減時的截止頻率%模擬濾波器系統函數分子分母%模擬濾波器系統函數分子分母%數字濾波器系統函數分子分母%直接形式轉化為并聯形式%檢驗衰減指標%數字濾波器的頻域響應b1,a1=u_buttap(N,OC);b,a=butter(N,OC,'s');bz,az=impinvar(b,a,Fs)C,B,A=dir2par(bz,az);w0=Wp,Ws;Hx=freqz(bz,az,w0);H,w=freqz(bz,az);db,mag,p

33、ha,grd,w=freqz_m(bz,az);%求模擬濾波器振幅，相位等ha,x,t=impulse(b1,a1)%模擬濾波器時域響應hn,t=impz(bz,az)%數字濾波器時域響應w1=zeros(1,512-length(w),w;%矢量維度一致dbHx=-20*log10(abs(Hx)/max(abs(H);%檢驗衰減指標figure(1);plot(mag);title('模擬濾波器幅頻');gridonfigure(2);plot(pha);title('模擬濾波器相頻');gridonfigure(3);plot(w1/pi,20*log10

34、(abs(H);title('數字濾波器幅頻,);xlabel('w/pi');ylabel('幅度響應曲線dB');gridonfigure(4);plot(w1/pi,angle(H);title('數字濾波器相頻');xlabel('w/pi');ylabel('相位響應曲線');gridonfigure(5);plot(ha);title('沖擊響應");gridonfigure(6);plot(t,hn);title('單位抽樣響應,);xlabel('n'

35、;,'FontSize',12);ylabel('h(n)','fontsize',12);gridonfigure(7):plot(w/pi,db);hn=hn;zeros(128-length(hn),1);figure(8);y=20*log10(abs(fft(hn,1024);%時域響應的FFT變換求幅度響應x=(0:1023)/(512);%調整FFT生成的波形y1=y(x>=0&x<=1);plot(x(x>=0&x<=1),y1);xlabel('w/pi');ylabel（

36、'幅度響應曲線dB'）;gridonfigure(9) y3=y2(x>=0&x<=1);y2=angle(fft(hn,1024);%時域響應的FFT變換求相位響應x=(0:1023)/(512);%調整FFT生成的波形plot(x(x>=0&x<=1),y3);xlabel('w/pi');ylabel（'相位響應曲線'）;gridonfigure(10)plot(w1/pi,20*log10(abs(H),'*r',x(x>=0&x<=1),y1,'ob&#

37、39;);%將兩種方法曲線畫在圖上比較legend（'直接求響應的曲線','FFT求響應的曲線'）;xlabel('w/pi');ylabel（'幅度響應曲線dB'）;figure(11)plot(w1/pi,angle(H),'.r',x(x>=0&x<=1),y3,'ob');%將兩種方法曲線畫在圖上比較legend('直接求響應的曲線,'FFT求響應的曲線');xlabel('w/pi');ylabel('相位響應曲線');functionC,B,A=dir2par(b,a)M=length(b);N=length(a);r1,p