實(shí)驗(yàn)4基于MATLAB的FIR數(shù)字濾波器設(shè)計(jì)

1、實(shí)驗(yàn)4基于MATLAB的FIR數(shù)字濾波器設(shè)計(jì)實(shí)驗(yàn)?zāi)康模杭由顚?duì)數(shù)字濾波器的常用指標(biāo)和設(shè)計(jì)過(guò)程的理解.實(shí)驗(yàn)原理：低通濾波器的常用指標(biāo):1pH()1p,forH()S,forS通帶邊緣頻率P,阻帶邊緣頻率S通帶起伏P,通帶峰值起伏p2021g10(1p)dBsppSPassbandStopbandTransitionbandFig1TypicalmagnitudespecificationforadigitalLPF阻帶起伏s,最小阻帶衰減S2021g10(s)dB數(shù)字濾波器有IIR和FIR兩種類(lèi)型,它們的特點(diǎn)和設(shè)計(jì)方法不同在MATLAB中,可以用b=fir1(N,Wn,'ftypeaper

2、)等函數(shù)輔助設(shè)計(jì)FIR數(shù)字濾波器.N代表濾波器階數(shù)；Wn代表濾波器的截止頻率(歸一化頻率),當(dāng)設(shè)計(jì)帶通和帶阻濾波器時(shí),Wn為雙元素相量；ftype代表濾波器類(lèi)型,如high'高通,'stop'帶阻等；taper為窗函數(shù)類(lèi)型,默認(rèn)為海明窗,窗系數(shù)需要實(shí)現(xiàn)用窗函數(shù)blackman,hamming,hanningchebwin,kaiser產(chǎn)生.例1用凱塞窗設(shè)計(jì)一FIR低通濾波器,通帶邊界頻率p0.3,阻帶邊界頻率s50dBo0.5,阻帶衰減不小于解首先由過(guò)渡帶寬和阻帶衰減來(lái)決定凱塞窗的N和0.2上圖給出了以上設(shè)計(jì)的頻率特性,(a)為N=30直接截取的頻率特性(b)為凱塞窗設(shè)

3、計(jì)的頻率特性.凱塞窗設(shè)計(jì)對(duì)應(yīng)的MATLAB程序?yàn)椋簑n=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);axis(0,1,-80,10);grid;xlabel('歸一化頻率/');ylabel('幅度/dB');還可以使用n,Wn,beta,fty

4、pe=kaiserord(f,a,dev)函數(shù)來(lái)估計(jì)濾波器階數(shù)等,得到凱塞窗濾波器：fcuts=0.30.5;%歸一化頻率omega/pimags=10;devs=0.0510A(-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)用MATLAB信號(hào)處理工具箱函數(shù)remezord來(lái)計(jì)算等波紋濾波器階數(shù)N和加權(quán)函數(shù)W(o),調(diào)用函數(shù)remez可進(jìn)行等波紋濾波器的設(shè)計(jì),直接求出濾波

5、器系數(shù).函數(shù)remezord中的數(shù)組fedge為通帶和阻帶邊界頻率,數(shù)組mval是兩個(gè)邊界處的幅值,而數(shù)組dev是通帶和阻帶的波動(dòng),fs是采樣頻率單位為Hz.例2利用雷米茲交替算法設(shè)計(jì)等波紋濾波器,設(shè)計(jì)一個(gè)線(xiàn)性相位低通FIR數(shù)字濾波器,其指標(biāo)為：通帶邊界頻率fc=800Hz,阻帶邊界fr=1000Hz,通帶波動(dòng)阻帶最小衰減At=40dB,采樣頻率fs=4000Hz.在MATLAB中可以用remezord和remez兩個(gè)函數(shù)設(shè)計(jì),其結(jié)果如圖2,MATLAB程序如下：fedge=8001000;mval=10;dev=0.05590.01;fs=4000;N,fpts,mag,wt=remezor

6、d(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'一、實(shí)驗(yàn)內(nèi)容：利用MATLAB編程設(shè)計(jì)一個(gè)數(shù)字帶通濾波器,指標(biāo)要求如下：通帶邊緣頻率：pi0.45,P20.65,通帶峰值起伏：p1dBo阻帶邊緣頻率：S10.3,S20.75,最小阻帶衰減：s40dBO分別用窗函數(shù)法和等波紋濾波器法設(shè)計(jì)兩種FIR數(shù)字濾波器.實(shí)驗(yàn)要求：給出FIR數(shù)字濾波器的沖激響應(yīng),繪出它們的幅度和

7、相位頻響曲線(xiàn),討論它們各自的實(shí)現(xiàn)形式和特點(diǎn).1-1用窗函數(shù)法實(shí)現(xiàn)：調(diào)用函數(shù)n,wn,bta,ftype=kaiserordf,a,dev,fs參數(shù)：f=0.30.450.650.8為對(duì)應(yīng)數(shù)字頻率si0.3,pi0.45,P20.65,S20.75a=010為由f指定的各個(gè)頻帶上的幅值向量,一般只有0和1表示；和f長(zhǎng)度關(guān)系為2*a的長(zhǎng)度一2=f的長(zhǎng)度devs=0.010.10870.01用于指定各個(gè)頻帶輸出濾波器的頻率響應(yīng)與其期望幅值之間的最大輸出誤差或偏差,長(zhǎng)度與a相等,計(jì)算公,通帶波動(dòng)衰減式：阻帶衰減誤差=誤差=fs缺省值為2HZ>>n,wn,bta,ftype=kaiseror

8、d(0.30.450.650,8,010,0,010.10870.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)>>grid> >xlabel('歸一化頻率w');ylabel('幅度/db

9、9;);> >subplot(2,1,2)> >plot(w1/pi,angle(hh1)> >grid> >xlabel('歸一化頻率w');ylabel('相位/rad');h1=Columns1through80.00410.0055-0.0091-0.0018-0.0056-0.00000.0391-0.0152Columns9through16-0.03810.0077-0.02930.09400.0907-0.2630-0.05170.3500Columns17through24-0.0517-0.26

10、300.09070.0940-0.02930.0077-0.0381-0.0152Columns25through310.0391-0.0000-0.0056-0.0018-0.00910.00550.0041圖4-1如果直接用freqz(h1,1,256);畫(huà)圖得:0石口00.10.20.30.4050.60.7080.91NDrmali2edFrequencyktira(i/sample)s<dp=ce(cw-1000-500-1000-1300-0Q.10.2D,30.40.50.60.708口.91MormalizedFrequency(xjiracj/samplej5C0(Sa

11、EOJP)3gfCL1-2用等波紋法設(shè)計(jì):調(diào)用函數(shù)n,fpts,mag,wt=remezord(f,a,dev)f=0.30.450.650.8a=010dev=0.010.10870.01其含義同函數(shù)n,wn,bta,ftype=kaiserord(f,a,dev,fs)中的參數(shù)相同.>>n,fpts,mag,wt=remezord(0.30.450.650.8,010,0.010.10870.01);%用remezord函數(shù)估算出remez函數(shù)要用到的階n、歸一化頻帶邊緣矢量fpts、頻帶內(nèi)幅值響應(yīng)矢量mag及加權(quán)矢量w,使remez函數(shù)設(shè)計(jì)出的濾波器滿(mǎn)足f、a及dev指定的性

12、能要求.> >h2=remez(n,fpts,mag,wt)%設(shè)計(jì)出等波紋濾波器> >hh2,w2=freqz(h2,1,256);> >figure(2)> >subplot(2,1,1)>>plot(w2/pi,20*log10(abs(hh2)> >grid> >xlabel('歸一化頻率w');ylabel('幅度/db');> >subplot(2,1,2)> >plot(w2/pi,angle(hh2)> >grid> >

13、;xlabel('歸一化頻率w');ylabel('相位/rad');> >h2h2=Columns1through9-0.00130.0092-0.0255-0.06420.11770.0922-0.2466-0.04660.3116Columns10through17-0.0466-0.24660.09220.1177-0.0642-0.02550.0092-0.00137016o.o,化o204-4PE坦m00.10.20,304050.60.7080.91NDrmahiedFrequency(ktira(J/satriple)(Llp)<

14、;dp=ch(cw圖4-2用freqz(h2,1,256)直接得圖：50o巧口-1CO5CO0,10.2D,3040.50.60.708口.9NormalizedFrequencyrad/sample)、對(duì)課本作業(yè)9.23畫(huà)圖2-1)用漢寧窗實(shí)現(xiàn)：1)手動(dòng)計(jì)算>>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=freq

15、z(h1,1,256);> >figure(1)> >subplot(2,1,1)>>plot(w1,20*log10(abs(hh1)> >grid> >xlabel('數(shù)字頻率w/rad');ylabel('幅度/db');> >subplot(2,1,2)> >plot(w1,angle(hh1)> >grid> >xlabel('數(shù)字頻率w/rad');ylabel('相位/rad');圖4-32-2)用自帶的fir

16、l函數(shù)：> >n=59;> >wn=3/85/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('數(shù)字頻率w/rad');ylabel('幅度/db');> >grid> >subplot(2,1,2)> >

17、plot(w2,angle(hh2)> >xlabel('數(shù)字頻率w/rad');ylabel('相角/rad');> >grid圖4-42-3)用等波紋法設(shè)計(jì)：調(diào)用函數(shù)n,fpts,mag,wt=remezord(f,a,dev)f=0.26250.3750.6250.7375a=010>>dev=0.010.10870.01>>n,fpts,mag,wt=remezord(0.26250.3750.6250.7375,010,0.010.10870.01);h2=remez(n,fpts,mag,wt);>

18、;>hh2,w2=freqz(h2,1,256);>>figure(2)>>subplot(2,1,1)>>plot(w2/pi,20*log10(abs(hh2)> >grid> >xlabel('歸一化頻率w');ylabel('幅度/db');> >subplot(2,1,2)> >plot(w2/pi,angle(hh2)> >grid> >xlabel('歸一化頻率w');ylabel('相位/rad');&

19、gt; >h2畫(huà)圖為：歸一化頻率加另外帶通濾波器還可以用低通和高通級(jí)聯(lián)的方法實(shí)現(xiàn),由于步驟過(guò)多,這里不作討論.三、實(shí)驗(yàn)分析：1總結(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)稍微增加或降低階次.如果截止頻率在Nyquist頻率附近,或者設(shè)定的dev值較大,那么得不到正確結(jié)果.2)濾波器實(shí)現(xiàn)形式及特點(diǎn)：由于一般的濾波器在利用窗函數(shù)是其通帶波紋和阻帶波

20、紋不同(一般為第一個(gè)阻帶波紋最大)因此,在滿(mǎn)足第一個(gè)阻帶衰減旁瓣時(shí),比其頻率高的旁瓣,它們的衰減都大大超出fs要求.而根據(jù)阻帶最減與項(xiàng)數(shù)的近似關(guān)系NP(At)可得當(dāng)阻帶衰減越大,所需項(xiàng)數(shù)TW越多.2)等波紋設(shè)計(jì)和窗函數(shù)設(shè)計(jì)不同之處在于實(shí)現(xiàn)形式及特點(diǎn)：窗函數(shù)設(shè)計(jì)是通過(guò)最小平方積分方法來(lái)設(shè)計(jì)的即該濾波器的誤差為:2E窗1H理想()H實(shí)際()Id即要求E窗最小方法來(lái)設(shè)計(jì)濾波器,這樣的濾波器更忠實(shí)于理想濾波器(即濾波系數(shù)更接近于理想濾波器.證實(shí)如下：2.E窗1H理想()H實(shí)際()|d*(H理想()H實(shí)際()(H理想()H實(shí)際()d2*(|H理想()|H理想()H實(shí)P()2H壬!想()H實(shí)際()|H實(shí)際

21、()|2)d=22*(|h理想n|h理想nh實(shí)際n*h理想nh實(shí)際n|h實(shí)際n|)2=2(|h理想nh實(shí)際n|)因此,幅度頻譜差值越小,實(shí)際濾波器就越接近理想濾波器.而等波紋濾波器是通過(guò)最大加權(quán)誤差最小化來(lái)實(shí)現(xiàn),其誤差為：E等波紋W()(H理想()H實(shí)際()要求該誤差最小來(lái)實(shí)現(xiàn)濾波器,得出來(lái)的濾波系數(shù)較窗函數(shù)設(shè)計(jì)相差較遠(yuǎn).以下通過(guò)對(duì)題目中的h1及h2作比擬.%sigsum是用來(lái)對(duì)數(shù)組各元素進(jìn)行求和functiony=sigsum(n1,n2,n,x);y=0;fori=n1+1-min(n):n2+1-min(n)y=y+x(i);end>>n=0.001:30.001;>&

22、gt;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).A2)/31;> >y2=sigsum(0,16,0:16,(abs(delta2).A2)/17;> >y1y1=1.9099e-004> >y2y2=0.0278由此得到用窗函數(shù)實(shí)現(xiàn)的濾波系數(shù)比用等波紋濾波器系數(shù)的每一項(xiàng)更接近于理想濾波器(y1為用窗函數(shù)實(shí)現(xiàn)的與理想濾波器的差值,y2為用等波紋濾波器實(shí)現(xiàn)的與理想濾波器的差值)；比照?qǐng)D4-1及圖4-2可見(jiàn)在幅度頻譜上等波紋濾波器阻帶邊緣比用窗函數(shù)實(shí)現(xiàn)的更平滑(理想濾波器為垂直下降的).從設(shè)計(jì)的角度考慮,由于窗函數(shù)設(shè)計(jì)法都是通過(guò)已有的