基于Matlab的FIR濾波器設(shè)計(jì)與實(shí)現(xiàn) - sunev - 博客園

1、sunev· 博客園· 閃存· 首頁(yè)· 新隨筆· 聯(lián)系· 管理· 訂閱 隨筆- 43  文章- 0  評(píng)論- 162  基于Matlab的FIR濾波器設(shè)計(jì)與實(shí)現(xiàn) 一、摘要前面一篇文章介紹了通過(guò)FDATool工具箱實(shí)現(xiàn)濾波器的設(shè)計(jì)，見(jiàn)“基于Matlab中FDATool工具箱的濾波器設(shè)計(jì)及相關(guān)文件的生成”，這里通過(guò)幾個(gè)例子說(shuō)明采用Matlab語(yǔ)言設(shè)計(jì)FIR濾波器的過(guò)程。二、實(shí)驗(yàn)平臺(tái)Matlab7.1三、實(shí)驗(yàn)原理以低通濾波器為例，其常用的設(shè)計(jì)指標(biāo)有：1. 通帶邊緣頻率fp(數(shù)字頻率為p)2. 阻帶邊緣

2、頻率fst (數(shù)字頻率為st)3. 通帶內(nèi)最大紋波衰減p=-20log10(1-p)，單位為 dB4. 阻帶最小衰減s=-20log10(s)，單位為 dB5. 阻帶起伏s6. 通帶峰值起伏p其中，以1、2、3、4條最為常用。5、6條在程序中估算濾波器階數(shù)等參數(shù)時(shí)會(huì)用到。數(shù)字頻率 = 模擬頻率/采樣頻率四、實(shí)例分析例1  用凱塞窗設(shè)計(jì)一FIR低通濾波器，通帶邊界頻率p=0.3pi，阻帶邊界頻率 s=0.5pi，阻帶衰減s不小于50dB。方法一：手動(dòng)計(jì)算濾波器階數(shù)N和值，之后在通過(guò)程序設(shè)計(jì)出濾波器。第一步：通過(guò)過(guò)渡帶寬度和阻帶衰減，計(jì)算濾波器的階數(shù)B和值。

3、第二步：通過(guò)程序設(shè)計(jì)濾波器。程序如下：b = fir1(29,0.4,kaiser(30,4.55);h1,w1=freqz(b,1);  plot(w1/pi,20*log10(abs(h1); axis(0,1,-80,10); grid;xlabel('歸一化頻率/p') ylabel('幅度/dB')  波形如下： 方法二：采用n,Wn,beta,ftype = kaiserord(f,a,dev)函數(shù)來(lái)估計(jì)濾波器階數(shù)等，得到凱塞窗濾波

4、器。這里的函數(shù)kaiserord(f,a,dev)或者kaiserord(f,a,dev,fs)：f為對(duì)應(yīng)的頻率，fs為采樣頻率；當(dāng)f用數(shù)字頻率表示時(shí)，fs則不需要寫(xiě)。a=1 0為由f指定的各個(gè)頻帶上的幅值向量，一般只有0和1表示；a和f長(zhǎng)度關(guān)系為（2*a的長(zhǎng)度）- 2=（f的長(zhǎng)度）devs=0.05 10(-2.5)用于指定各個(gè)頻帶輸出濾波器的頻率響應(yīng)與其期望幅值之間的最大輸出誤差或偏差，長(zhǎng)度與a相等，計(jì)算公式：阻帶衰減誤差=s，通帶衰減誤差=p，可有濾波器指標(biāo)中的3、4條得到。fs缺省為2Hz。 程序如下：fcuts = 0.3 

5、 0.5; %歸一化頻率omega/pi，這里指通帶截止頻率、阻帶起始頻率mags = 1 0;devs = 0.05 10(-2.5);n,Wn,beta,ftype = kaiserord(fcuts,mags,devs);  %計(jì)算出凱塞窗N，beta的值hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale'); freqz(hh); 波形如下： 實(shí)際中，一般調(diào)用M

6、ATLAB信號(hào)處理工具箱函數(shù)remezord來(lái)計(jì)算等波紋濾波器階數(shù)N和加權(quán)函數(shù)W()，調(diào)用函數(shù)remez可進(jìn)行等波紋濾波器的設(shè)計(jì)，直接求出濾波器系數(shù)。函數(shù)remezord中的數(shù)組fedge為通帶和阻帶邊界頻率，數(shù)組mval是兩個(gè)邊界處的幅值，而數(shù)組dev是通帶和阻帶的波動(dòng)，fs是采樣頻率單位為Hz。 例2  利用雷米茲交替算法設(shè)計(jì)等波紋濾波器，設(shè)計(jì)一個(gè)線性相位低通FIR數(shù)字濾波器，其指標(biāo)為：通帶邊界頻率fc=800Hz，阻帶邊界fr=1000Hz，通帶波動(dòng) 阻帶最小衰減At=40dB，采樣頻率fs=4000Hz。  解：在MATLAB中可以用r

7、emezord 和remez兩個(gè)函數(shù)設(shè)計(jì) 程序如下：  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

8、(h);  grid;  xlabel('頻率/Hz')  ylabel('幅度/dB'); 波形如下：例3 利用MATLAB編程設(shè)計(jì)一個(gè)數(shù)字帶通濾波器，指標(biāo)要求如下：通帶邊緣頻率：p1=0.45pi，p2=0.65pi，通帶峰值起伏：1<=1dB。阻帶邊緣頻率：s1=0.3pi，s2=0.8pi，最小阻帶衰減：2>=40dB 。方法一：窗函數(shù)法 程序如下：n,wn,bta,ftype=kaiserord(0.3 0.45 0.65

9、 0.8,0 1 0,0.01 0.1087 0.01);%用kaiserord函數(shù)估計(jì)出濾波器階數(shù)n和beta參數(shù)h1=fir1(n,wn,ftype,kaiser(n+1,bta),'noscale');hh1,w1=freqz(h1,1,256);figure(1)subplot(2,1,1)plot(w1/pi,20*log10(abs(hh1)gridxlabel('歸一化頻率w');ylabel('幅度/db');subplot(2,1,2)plot(w1/pi,angle(hh1)g

10、ridxlabel('歸一化頻率w');ylabel('相位/rad'); 波形如下： 濾波器系數(shù)為：h1 =   Columns 1 through 8      0.0041    0.0055   -0.0091   -0.0018   -0.0056  

11、60;-0.0000    0.0391   -0.0152   Columns 9 through 16     -0.0381    0.0077   -0.0293    0.0940    0.0907   -0.2630 

12、;  -0.0517    0.3500   Columns 17 through 24     -0.0517   -0.2630    0.0907    0.0940   -0.0293    0.0077   -

13、0.0381   -0.0152   Columns 25 through 31      0.0391   -0.0000   -0.0056   -0.0018   -0.0091    0.0055    0.0041 如果直接用fr

14、eqz(h1,1,256)，得幅頻特性和相頻特性曲線： 方法二：等波紋法設(shè)計(jì) 程序如下：n,fpts,mag,wt=remezord(0.3 0.45 0.65 0.8,0 1 0,0.01 0.1087 0.01);%用remezord函數(shù)估算出remez函數(shù)要用到的階n、歸一化頻帶邊緣矢量fpts、頻帶內(nèi)幅值響應(yīng)矢量mag及加權(quán)矢量w，使remez函數(shù)設(shè)計(jì)出的濾波器滿足f、a及dev指定的性能要求。h2=remez(n,fpts,mag,wt);%設(shè)計(jì)出等波紋濾波器hh2,w2=freqz(h2,1,256);figure(2)subplot(2,1,1)plot(w2

15、/pi,20*log10(abs(hh2)gridxlabel('歸一化頻率w');ylabel('幅度/db');subplot(2,1,2)plot(w2/pi,angle(hh2)gridxlabel('歸一化頻率w');ylabel('相位/rad');h2 波形如下：濾波器系數(shù)如下：h2 =   Columns 1 through 9     -0.0013   

16、 0.0092   -0.0255   -0.0642    0.1177    0.0922   -0.2466   -0.0466    0.3116   Columns 10 through 17     -0.0466 &

17、#160; -0.2466    0.0922    0.1177   -0.0642   -0.0255    0.0092   -0.0013 如果直接用freqz(h2,1,256);得幅頻特性和相頻特性曲線： 方法三：采用FDATool工具 這種方法需要事先計(jì)算出濾波器的階數(shù)，bate值，然后設(shè)置相應(yīng)參數(shù)，最后生成濾波器。 設(shè)置界

18、面如下圖所示：  將上述圈圈的區(qū)域設(shè)置好之后，生成濾波器，最后通過(guò)analysis菜單可以觀察生成的濾波器的各種特性曲線和濾波器系數(shù)。這里的濾波器系數(shù)跟方法一的一樣。 波形如下： 五、結(jié)果分析5.1  濾波器設(shè)計(jì)總結(jié)FIR濾波器實(shí)現(xiàn)一般采用窗函數(shù)法和等紋波設(shè)計(jì)法。窗函數(shù)法還包含兩個(gè)分支，一種是用公式先手動(dòng)算出N值和其他對(duì)應(yīng)得窗函數(shù)參數(shù)值，再代入窗函數(shù)和fir1實(shí)現(xiàn)，一種是用函數(shù)*rord估算出N和相應(yīng)參數(shù)再用fir1實(shí)現(xiàn)。不過(guò)要注意*rord會(huì)低估或高估階次n，可能會(huì)使濾波器達(dá)不到指定的性能，這時(shí)應(yīng)稍微增加或降低階次。如果截止頻率在0或Nyqu

19、ist頻率附近，或者設(shè)定的dev值較大，則得不到正確結(jié)果。濾波器實(shí)現(xiàn)形式及特點(diǎn)：由于一般的濾波器在利用窗函數(shù)是其通帶波紋和阻帶波紋不同（一般為第一個(gè)阻帶波紋最大）因此，在滿足第一個(gè)阻帶衰減旁瓣時(shí)，比其頻率高的旁瓣，它們的衰減都大大超出要求。而根據(jù)阻帶衰減與項(xiàng)數(shù)的近似關(guān)系N = P(2)*fs/TW，可得當(dāng)阻帶衰減越大，所需項(xiàng)數(shù)越多。5.2  窗函數(shù)法和等波紋設(shè)計(jì)的不同之處窗函數(shù)設(shè)計(jì)是通過(guò)最小平方積分辦法來(lái)設(shè)計(jì)的，即該濾波器的誤差為：即要求最小方法來(lái)設(shè)計(jì)濾波器，這樣的濾波器更忠實(shí)于理想濾波器（即濾波系數(shù)更接近于理想濾波器）。證明如下：因此，幅度頻譜差值越小，實(shí)際濾波器就越接近理想濾波器

20、。而等波紋濾波器是通過(guò)最大加權(quán)誤差最小化來(lái)實(shí)現(xiàn)，其誤差為：要求該誤差最小來(lái)實(shí)現(xiàn)濾波器，得出來(lái)的濾波系數(shù)較窗函數(shù)設(shè)計(jì)相差較遠(yuǎn)。以下通過(guò)對(duì)例3中的h1及h2作比較。%sigsum是用來(lái)對(duì)數(shù)組各元素進(jìn)行求和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;  結(jié)果如下：y1 =   1.9099e-004 y2 =     0.0278