1、實驗 4基于 MATLAB的 FIR數字濾波器設計實驗目的： 加深對數字濾波器的常用指標和設計過程的理解。通帶邊緣頻率P ，阻帶邊緣頻率S，通帶起伏P，SPPSPassbandStopband TransitionbandFig 1 Typical magnitude通帶峰值起伏p20 log 10 (1p ) dB ，specification for a digital LPF阻帶起伏s ，最小阻帶衰減S20 log10 (s ) dB 。實驗原理： 低通濾波器的常用指標：H ()1P1PH ()1P , forP1PH ()S , forS?s數字濾波器有 IIR 和 FIR 兩種類型，

2、它們的特點和設計方法不同。在 MATLAB中，可以用b=fir1(N,Wn, ftyp,etaper)等函數輔助 設計 FIR 數字濾波器。 N 代表濾波器階數； Wn 代表濾波器的截止頻率(歸一化頻率 )，當設計帶通和帶阻濾波器時，Wn為雙元素相量； ftype代表濾波器類型，如 high高通， stop帶阻等； taper 為窗函數類 型 ， 默 認為 海 明 窗 ， 窗系 數 需 要 實現 用 窗 函 數blackman, hamming,hanning chebwin, kaiser產生。例 1用凱塞窗設計一FIR 低通濾波器，通帶邊界頻率p0.3，阻帶邊界頻率s0.5，阻帶衰減不小于

3、50dB。解 首先由過渡帶寬和阻帶衰減來決定凱塞窗的 N 和sp0.2，上圖給出了以上設計的頻率特性，(a) 為 N=30 直接截取的頻率特性 (b)為凱塞窗設計的頻率特性。凱塞窗設計對應的MATLAB程序為：wn=kaiser(30,4.55); nn=0:1:29;alfa=(30-1)/2;hd=sin(0.4*pi*(nn-alfa)./(pi*(nn-alfa); h=hd.*wn;h1,w1=freqz(h,1);或者： b = fir1(29,0.4,kaiser(30,4.55); h1,w1=freqz(b,1);plot(w1/pi,20*log10(abs(h1); ax

4、is(0,1,-80,10);grid;xlabel(歸一化頻率 /) ; ylabel(幅度/dB) ;還可以使用 n,Wn,beta,ftype = kaiserord(f,a,dev)函數來估計濾波器階數等，得到凱塞窗濾波器：fcuts = 0.30.5; %歸一化頻率 omega/pimags = 1 0;devs = 0.05 10(-2.5);n,Wn,beta,ftype = kaiserord(fcuts,mags,devs);%計算出凱塞窗 N，beta 的值hh = fir1(n,Wn,ftype,kaiser(n+1,beta),noscale);freqz(hh);實際

5、中，一般調用 MATLAB 信號處理工具箱函數 remezord 來計算等波紋濾波器階數 N 和加權函數 W()，調用函數 remez 可進行等波紋濾波器的設計，直接求出濾波器系數。 函數 remezord中的數組 fedge 為通帶和阻帶邊界頻率， 數組 mval 是兩個邊界處的幅值， 而數組 dev 是通帶和阻帶的波動， fs 是采樣頻率單位為 Hz。例 2 利用雷米茲交替算法設計等波紋濾波器，設計一個線性相位低通 FIR 數字濾波器，其指標為：通帶邊界頻率 f c=800Hz，阻帶邊界 f r=1000Hz，通帶波動阻帶最小衰減 At=40dB ，采樣頻率fs=4000Hz。解在 MAT

6、LAB中可以用 remezord 和 remez 兩個函數設計，其結果如圖2， MATLAB程序如下：fedge=800 1000;mval=1 0;dev=0.0559 0.01; fs=4000;N,fpts,mag,wt=remezord(fedge,mval,dev,fs); b=remez(N,fpts,mag,wt); h,w=freqz(b,1,256); plot(w*2000/pi,20*log10(abs(h);grid;xlabel(頻率/Hz) ;ylabel(幅度/dB);一、實驗內容：利用 MATLAB編程設計一個數字帶通濾波器，指標要求如下：通帶邊緣頻率：P10.

7、45，P 20.65，通帶峰值起伏：p1 dB 。阻帶邊緣頻率：S10.3，S20.75，最小阻帶衰減：S40 dB 。分別用窗函數法和等波紋濾波器法設計兩種FIR 數字濾波器。實驗要求：給出 FIR 數字濾波器的沖激響應，繪出它們的幅度和相位頻響曲線， 討論它們各自的實現形式和特點。1-1）用窗函數法實現：調用函數 n,wn,bta,ftype=kaiserord（ f，a， dev，fs）參數： f=0.3 0.45 0.65 0.8為對應數字頻率S10.3，P10.45，P 20.65，S20.75a=0 1 0為由 f 指定的各個頻帶上的幅值向量，一般只有0 和 1表示；和 f 長度關

8、系為（ 2*a 的長度） 2=（f 的長度）devs=0.01 0.1087 0.01用于指定各個頻帶輸出濾波器的頻率響應與其期望幅值之間的最大輸出誤差或偏差，長度與 a 相等， 計算公式：阻帶衰減誤差 =，通帶波動衰減誤差=fs 缺省值為 2HZn,wn,bta,ftype=kaiserord(0.30.45 0.650.8,010,0.010.1087 0.01);%用 kaiserord 函數估計出濾波器階數n 和 beta 參數 h1=fir1(n,wn,ftype,kaiser(n+1,bta),noscale); hh1,w1=freqz(h1,1,256); figure(1)

9、subplot(2,1,1) plot(w1/pi,20*log10(abs(hh1) grid xlabel(歸一化頻率 w);ylabel( 幅度/db); subplot(2,1,2) plot(w1/pi,angle(hh1) grid xlabel(歸一化頻率 w);ylabel( 相位/rad);h1 =Columns 1 through 80.00410.0055-0.0091-0.0018-0.0056-0.00000.0391-0.0152Columns 9 through 16-0.03810.0077-0.02930.09400.0907-0.2630-0.05170.3

10、500Columns 17 through 24-0.0517-0.26300.09070.0940-0.02930.0077-0.0381-0.0152Columns 25 through 310.0391-0.0000-0.0056-0.0018-0.00910.00550.0041圖 4-1如果直接用freqz(h1,1,256);畫圖得：1-2）用等波紋法設計：調用函數 n,fpts,mag,wt=remezord(f， a，dev) f=0.3 0.45 0.65 0.8a=0 1 0dev=0.01 0.1087 0.01其含義同 函數n,wn,bta,ftype=kaiseror

11、d（f ，a， dev，fs）中的參數相同。n,fpts,mag,wt=remezord(0.30.45 0.65 0.8,010,0.010.10870.01);%用 remezord 函數估算出 remez 函數要用到的階n、歸一化頻帶邊緣矢量 fpts、頻帶內幅值響應矢量mag 及加權矢量 w，使 remez 函數設計出的濾波器滿足f、a 及 dev 指定的性能要求。 h2=remez(n,fpts,mag,wt);%設計出等波紋濾波器 hh2,w2=freqz(h2,1,256); figure(2) subplot(2,1,1) plot(w2/pi,20*log10(abs(hh2

12、) grid xlabel(歸一化頻率 w);ylabel( 幅度/db); subplot(2,1,2) plot(w2/pi,angle(hh2) grid xlabel(歸一化頻率 w);ylabel( 相位/rad); h2h2 =Columns 1 through 9-0.00130.0092-0.0255-0.06420.11770.0922-0.2466-0.04660.3116Columns 10 through 17-0.0466-0.24660.09220.1177-0.0642-0.02550.0092-0.0013圖 4-2用 freqz(h2,1,256);直接得圖：

13、二、對課本作業 9.23 畫圖。2-1）用漢寧窗實現：1）手動計算 n=0.001:58.001; hd=sin(0.18125*pi*(n-29)./(pi*(n-29); win=0.5+0.5*cos(2*pi*(n-29)/58); h1=2*cos(pi*(n-29)/2).*hd.*win; hh1,w1=freqz(h1,1,256); figure(1) subplot(2,1,1) plot(w1,20*log10(abs(hh1) grid xlabel( 數字頻率w/rad);ylabel( 幅度 /db); subplot(2,1,2) plot(w1,angle(hh

14、1) grid xlabel( 數字頻率w/rad);ylabel( 相位 /rad);圖 4-32-2）用自帶的 fir1 函數： n=59; wn=3/8 5/8; h2=fir1(n,wn,bandpass,hann(n+1); hh2,w2=freqz(h2,1,256); figure(2) subplot(2,1,1) plot(w2,20*log10(abs(hh2) xlabel( 數字頻率w/rad);ylabel( 幅度 /db); grid subplot(2,1,2) plot(w2,angle(hh2) xlabel( 數字頻率w/rad);ylabel( 相角 /r

15、ad); grid圖 4-42-3）用等波紋法設計：調用函數 n,fpts,mag,wt=remezord(f ， a， dev) f=0.2625 0.375 0.625 0.7375a=0 1 0dev=0.01 0.1087 0.01n,fpts,mag,wt=remezord(0.26250.3750.6250.7375,010,0.010.10870.01);h2=remez(n,fpts,mag,wt); hh2,w2=freqz(h2,1,256); figure(2) subplot(2,1,1) plot(w2/pi,20*log10(abs(hh2) grid xlabel

16、( 歸一化頻率w);ylabel( 幅度 /db); subplot(2,1,2) plot(w2/pi,angle(hh2) grid xlabel( 歸一化頻率w);ylabel( 相位 /rad); h2畫圖為：另外帶通濾波器還可以用低通和高通級聯的方法實現，因為步驟過多，這里不作討論。三、實驗分析：總結 ： FIR濾波器實現一般采用窗函數法和等紋波設計法。窗函數法還包含兩個分支，一種是用公式先手動算出N 值和其他對應得窗函數參數值，再代入窗函數和fir1 實現，一種是用函數 *rord 估算出 N 和相應參數再用fir1 實現。不過要注意 *rord 會低估或高估階次n，可能會使濾波器

17、達不到指定的性能，這時應稍微增加或降低階次。如果截止頻率在0 或Nyquist 頻率附近，或者設定的dev 值較大，則得不到正確結果。濾波器實現形式及特點：由于一般的濾波器在利用窗函數是其通帶波紋和阻帶波紋不同（一般為第一個阻帶波紋最大）因此， 在滿足第一個阻帶衰減旁瓣時，比其頻率高的旁瓣，它們的衰減都大大超出要求。 而根據阻帶衰減與項數的近似關系越多。NP( At)fs可得當阻帶衰減越大，所需項數TW等波紋設計和窗函數設計不同之處在于實現形式及特點：窗函數設計是通過最小平方積分辦法來設計的即該濾波器的誤差為：E窗| H 理想 ()H 實際 () | 2 d即要求E窗 最小方法來設計濾波器，這

18、樣的濾波器更忠實于理想濾波器（即濾波系數更接近于理想濾波器。證明如下：*E窗| H 理想 ()H 實際 () |2 d=(H 理想 ()H 實際 ()( H *()H 實際 () d=(| H 理想 () |2H 理想 ()*理想Hhh()H()實際理想實際理想H 實際 ()| H 實際 () |2 )d= 2(| h理想n |2h理想 n*n* nh實際n| h實際 n |)= 2(| h理想nh實際n |) 2因此，幅度頻譜差值越小，實際濾波器就越接近理想濾波器。而等波紋濾波器是通過最大加權誤差最小化來實現，其誤差為：E等波紋W ()( H 理想 ()H 實際 ()要求該誤差最小來實現濾波

19、器，得出來的濾波系數較窗函數設計相差較遠。以下通過對題目中的h1 及 h2 作比較。%sigsum 是用來對數組各元素進行求和function y=sigsum(n1,n2,n,x); y=0;for i=n1+1-min(n):n2+1-min(n)y=y+x(i);end n=0.001:30.001; h=2*cos(0.55*pi*(n-15).*sin(0.175*pi*(n-15)./(pi*(n-15); delta1=h-h1; n=0.001:16.001; h=2*cos(0.55*pi*(n-15).*sin(0.175*pi*(n-15)./(pi*(n-15); delta2=h-h2; y1=sigsum(0,30,0:30,(abs(delta1).2)/31; y2=sigsum(0,16,0:16,(abs(delta2).2)/17; y1y1 =1.9099e-004 y2y2 =0.0278由此得到用窗函數實現的濾波系數比用等波紋濾波器系數的每一項更接近于理想濾波器（ y1 為用窗函數實現的與理想濾波器的差值，y2 為用等波紋濾波器實現的與理想濾波器的差值）；對比圖 4-1 及圖 4-2 可見在幅度頻譜上等波紋濾波器阻帶邊緣比用窗函數實現的更平滑（理想濾波器為