FIR濾波器設(shè)計

1、第7章FIR濾波器設(shè)計第六章我們介紹了無限沖激響應(yīng)（IIR）濾波器的設(shè)計方法。其中最常用的由模擬濾波器轉(zhuǎn)換為數(shù)字濾波器的方法為雙線性變換法,因為這種方法無混疊效應(yīng)，效果較好。但通過6-11、圖 6-13、圖前面的例子我們看到，IIR數(shù)字濾波器相位特性不好（非線性，如圖6-15 ）,也不易控制。然而在現(xiàn)代信號處理中，例如圖像處理、數(shù)據(jù)傳輸、雷達接收以及一些要求較高的系統(tǒng)中對相位特性要求較為嚴格,這種濾波器就無能為力了。 改善相位特性的方法是采用有限沖激響應(yīng)濾波器。本章首先對FIR濾波器原理及其使用函數(shù)作基本介紹,然后重點介紹窗函數(shù)法設(shè)計FIR濾波器，并對最優(yōu)濾波器設(shè)計函數(shù)進行介紹。7.1 FIR

2、濾波器原理概述及濾波函數(shù)7.1.1 FIR濾波器原理及設(shè)計方法分類根據(jù)第6章對數(shù)字濾波器的介紹，我們知道FIR濾波器的傳遞函數(shù)為:Y(z)X(z)y( n) b(0)x( n)N 1b(m)x(nm 0b(1)x(nm)因此，F(xiàn)IR濾波器又稱為卷積濾波器。1)b(N 1)x( n N 1)b(n) x(n)根據(jù)第4章中所描述的系統(tǒng)頻率響應(yīng),FIR濾波器的頻率響應(yīng)表達式為:N 1H ejb(n)en 0jn(7-2)信號通過FIR濾波器不失真條件與（6-6 ）式所描述的相同，即濾波器在通帶內(nèi)具有恒1(7-1)h(n)z n0可得FIR濾波器的系統(tǒng)差分方程為:定的幅頻特性和線性相位特性。理論上可以

3、證明（這里從略）:當FIR濾波器的系數(shù)滿足下列中心對稱條件:b(n) b(N 1 n) 或 b(n) b( N 1 n)(7-3)時，濾波器設(shè)計在逼近平直幅頻特性的同時，還能獲得嚴格的線性相位特性。線性相位FIR濾波器的相位滯后和群延遲在整個頻帶上是相等且不變的。對于一個 濾波器，群延遲為常數(shù)，即濾波后的信號簡單地延遲常數(shù)個時間步長。N階的線性相位FIR這一特性使通帶頻率內(nèi)信號通過濾波器后仍保持原有波形形狀而無相位失真。本章主要介紹的 FIR數(shù)字濾波器設(shè)計方法及MATLAB 信號處理工具箱提供的FIR數(shù)字濾波器設(shè)計函數(shù)，見表7-1。由于篇幅所限，本章我們主要介紹窗函數(shù)法和最優(yōu)化設(shè)計方法。表7-

4、1 FIR濾波器設(shè)計的主要方法函數(shù)設(shè)計方法說明工具函數(shù)窗函數(shù)法理想濾波器加窗處理fir1（單頻帶），fir2（多頻帶）,kaiserord最優(yōu)化設(shè)計平方誤差最小化逼近理firls,想幅頻響應(yīng)或Park-McClellan算法產(chǎn)生等波紋濾波器remez,remezord約束最小一乘逼近在滿足最大誤差限制條 件下使整個頻帶平方誤差最 小化fircls,fircls1升余弦函數(shù)具有光滑、正弦過渡帶Fircos的低通濾波器設(shè)計7.1.2 FIR數(shù)字濾波器濾波函數(shù)相對于IIR濾波器的濾波函數(shù)，F(xiàn)IR數(shù)字濾波器濾波函數(shù)除了dim pulse 和dste p僅適用于IIR濾波器外，其他各種函數(shù)可直接應(yīng)用于F

5、IR濾波器，只是輸入的分母多項式向量a=1。另外，MATLAB還提供了一個函數(shù) fftfilt,該函數(shù)利用效率高的基于FFT算法實現(xiàn)對數(shù)據(jù)的濾波，該函數(shù)只適用于FIR濾波器，調(diào)用形式為:y=fftfilt(b,x, n)式中，b為FIR濾波器的系數(shù)向量；x為輸入數(shù)據(jù)；n為FFT長度，缺省時，函數(shù)選用最佳 的FFT長度，y為濾波器的輸出。該函數(shù)執(zhí)行下面的操作:n=len gth(x);y=ifft(fft(x).*fft(b, n) ./ft(a, n);應(yīng)注意，y=fftfilt(b,x)等價于 y=filter(b,a,x)。7.2 FIR濾波器的窗函數(shù)設(shè)計7.2.1 窗函數(shù)的基本原理b，即

6、系統(tǒng)單位脈FIR濾波器設(shè)計的主要任務(wù)是根據(jù)給定的性能指標確定濾波器的系數(shù)沖序列h(n)，它是一個有限長序列。FIR濾波器的理想頻率響應(yīng)，可寫成復(fù)數(shù)形式的Fourier級數(shù)形式:Hd ejhdn(7-4)式中，hd(n)是對應(yīng)的單位脈沖響應(yīng)序列。這說明濾波器的頻率響應(yīng)和單位脈沖響應(yīng)互為Fourier變換對。因此其單位脈沖響應(yīng)可由下式求得,hd n求得序列hd n后，通過z變換，可得到Hd(z)nhd(n)z n(7-6)注意，這里hd n為無限長序列，因此 Hd Z是物理上不可實現(xiàn)的。如何變成物理上可實現(xiàn)呢？ 一個自然的想法是只取其中的某些項，即只截取hd n中的一部分，比如n=0,N-1 ,

7、N為正整數(shù)。這種處理相當于將hd n ,n=-s與函數(shù)w(n)相乘，w(n)具有下列形式:w(n)0, n 0，n N1,0 n Nw(n)相當于一個矩形，我們稱之為矩形窗。即我們可采用矩形窗函數(shù)w(n)將無限脈沖響應(yīng)hd n截取一段h(n)來近似為hd n，這種截取在數(shù)學(xué)上表示為:h(n)=hd n w(n)(7-7)這里應(yīng)該強調(diào)的是，加窗函數(shù)不是可有可無的，而是將設(shè)計變?yōu)槲锢砜蓪崿F(xiàn)所必須的。截取之后的濾波器傳遞函數(shù)變?yōu)?N 1H(z)h( n)znn 0(7-8)式中，N為窗口寬度，H(z)是物理可實現(xiàn)系統(tǒng)。為了獲得線性相位，F(xiàn)IR濾波器h(n)必須滿足中心對稱條件(即 7-3式)，序列h

8、(n)的延遲為N 1/2。這種方法的基本原理是用一定寬度的矩形窗函數(shù)截取無限脈沖響應(yīng)序列獲得有限長的脈沖響應(yīng)序列，從而得到FIR濾波器的脈沖響應(yīng)，故稱為 FIR濾波器的窗函數(shù)設(shè)計法。7-1經(jīng)過加矩形窗后所得的濾波器實際頻率響應(yīng)能否很好地逼近理想頻率響應(yīng)呢？圖示意給出了理想濾波器加矩形窗后的情況。理想低通濾波器的頻率響應(yīng)如圖中左上角圖,形窗的頻率響應(yīng)為左下角圖。時間域內(nèi)的乘積(7-7 )式要求實際頻率響應(yīng)為這兩個頻率響應(yīng)函數(shù)在頻域內(nèi)的卷積(卷積定理)，即得到圖形為圖 7-1 (右圖)。圖7-1 FIR濾波器理想與實際頻率響應(yīng)由圖可看出，加矩形窗后使實際頻率響應(yīng)偏離理想頻率響應(yīng)，主要影響有三個方面

9、:(1 )理想幅頻特性陡直邊緣處形成過渡帶，過渡帶寬取決于矩形窗函數(shù)頻率響應(yīng)的主 瓣寬度。(2)過渡帶兩側(cè)形成肩峰和波紋，這是矩形窗函數(shù)頻率響應(yīng)的旁瓣引起的，旁瓣相對值越大，旁瓣越多，波紋越多。(3)隨窗函數(shù)寬度 N的增大，矩形窗函數(shù)頻率響應(yīng)的主瓣寬度減小，但不改變旁瓣的 相對值。為了改善FIR濾波器性能，要求窗函數(shù)的主瓣寬度盡可能窄，以獲得較窄的過渡帶；旁瓣相對值盡可能小，數(shù)量盡可能少，以獲得通帶波紋小，阻帶衰減大，在通帶和阻帶內(nèi)均平穩(wěn)的特點，這樣可使濾波器實際頻率響應(yīng)更好地逼近理想頻率響應(yīng)。這里我們明確兩個概念： 截斷和頻譜泄漏。 信號是無限長的，而在進行信號處理時只能x(t)乘以有限長的

10、窗函數(shù)采取有限長信號，所以需要將信號“截斷”。在信號處理中，“截斷”被看成是用一個有限 長的“窗口”看無限長的信號，或者從分析的角度是無限長的信號(7-9 )w(t)。由傅立葉變換性質(zhì)可知，時間域內(nèi)的乘積對應(yīng)于頻率域的卷積，即x(t) w(t)X(f) W(f)這里，x(t)是頻寬有限信號，而w(t)是頻寬無限信號,表示互為Fourier變換對。截斷后的信號也必須是頻寬無限信號,這樣就是有限頻帶的信號分散到無限頻帶中去，這樣處理的前驅(qū)者們設(shè)計了不同于矩形窗的很多窗函數(shù),這些窗函數(shù)在主瓣和旁瓣特性方面各有就產(chǎn)生了所謂頻譜泄漏。 從能量的角度來看，頻譜泄漏也是能量的泄漏，因為加窗后使原來頻譜泄漏是

11、不可避免的，但要盡信號集中的窄頻帶內(nèi)的能量分散到無限的頻帶寬度范圍內(nèi)。量減小。上邊只考慮了矩形窗， 如果我們使窗的主瓣寬度盡可能地窄，旁瓣盡可能地小，可以獲其實數(shù)字信號得性能更好的濾波器， 能否改變窗的形狀而達到這個目的呢？回答是肯定的。特點，可滿足不同的要求。為此，用窗函數(shù)法設(shè)計FIR數(shù)字濾波器時，要根據(jù)給定的濾波器性能指標選擇窗口寬度 N和窗函數(shù)w(n)。下面我們介紹窗函數(shù)。7.2.2 MATLAB信號處理中提供的窗函數(shù)(1)矩形窗：前面分析中所用的矩形窗可用下面函數(shù)來實現(xiàn)w=boxcar (N)，N為窗的長度(以下函數(shù)與此同)，w為返回的窗函數(shù)序列。(2 )漢寧窗：w=ha nnin g

12、(N)漢寧窗的表達式為:kw(k) o.51 cos2 廠(7-10)(3 )哈明窗：w=hammi ng(N)哈明窗的表達式為:kw(k 1) 0.54 0.46cos 2 廿0,1,., N 1(7-11)(4)Bartlett 窗：w=bartlett(N)Bartlett窗的表達式為:N為奇數(shù)時,w(k)(7-12)2(k 1)N 1 '2 2(k 1)N 11 k NJN為偶數(shù)時,w(k)2(kN2(N k)1)(7-13)Blackman 窗: w= blackman(N)Blackma n窗的表達式為:點。w(k)Blackma n0.42 0.5COS 2窗比其他相同尺

13、寸窗k 10.08N 11,., N(7-14)（哈明窗,漢寧窗）具有主瓣較寬和旁瓣泄漏較小的特(6 )三角窗：w=tria ng(N)三角窗的表達式為: 當N為奇數(shù)時,當N為偶數(shù)時,w(k)吾12(N k 1)(7-15)(7-16),1 k Nw(k)N 1222(N k 1), kNN2三角窗和Bartlett窗十分類似。三角窗的兩端值不為零，而Bartlett窗則為零，這一點可從例7-1中看出。(7) Kaiser 窗:w=kaiser( n,beta)其中，beta 是 Kaiser窗參數(shù)，影響窗旁瓣幅值的衰減率。Kaiser窗表達式:If2I0 b 1 缶(7-17)w(k)式中，

14、|0.是修正過的零階Bessel函數(shù)。Kaiser窗用于濾波器設(shè)計時，若旁瓣幅值為d ，則0.11028.7 ,500 40.584221 .0.0788621 ,215021(7-18 )0,(8)Chebyshev 窗: w=chebwin(n,r)式中，r是窗口的旁瓣幅值在主瓣以下的分貝數(shù)。Chebyshev窗具有主瓣寬度最小，而旁瓣等咼，高度可調(diào)整的特點。F面我們在MATLAB觀看各種窗函數(shù)的形狀和頻率域圖象來驗證上述所講特點。21 O【例7-1】 用MATLAB編程繪制各種窗函數(shù)的形狀。窗函數(shù)的長度為%Sa mp7_1cifNwi n=21; n=O:Nwi n-1;%數(shù)據(jù)總數(shù)和序列

15、序號figure(1)for ii=1:4switch iicase 1w=boxcar(Nwi n);%矩形窗stext='矩形窗;case 2w=ha nnin g(Nw in);%漢寧窗stext='漢寧窗;case 3w=hammi ng(Nwi n);%哈明窗stext='哈明窗;case 4w=bartlett(Nwi n);%Bartlett 窗stext='Bartelett 窗'endpos pl ot='2,2,' in t2str(ii);%指定繪制窗函數(shù)的圖形位置sub plot(pospl ot);stem (n

16、 ,w);%繪出窗函數(shù)hold onplot (n ,w,'r');%繪制包絡(luò)線xlabel(' n');ylabel('w( n)');title(stext);hold off;grid onxlabel(' n');ylabel('w( n)');title(stext);endfigure (2)for ii=1:4switch iicase 1w=blackma n(Nwin);%Blackma n 窗stext='Blackma n 窗'case 2w=tria ng(Nwi n);%三

17、角窗stext='三角窗'case 3w=kaiser(Nwi n,4);%Kaiser 窗stext='Kaiser 窗(Beta=4)'case 4w=chebwi n( Nwi n,40);%Chebyshev 窗stext='Chebyshev 窗(r=40)'endpos pl ot='2,2,' in t2str(ii);%繪出窗函數(shù)sub plot(pospl ot);stem( n,w);hold onhold oftgrid on;end程序運行結(jié)果見圖7-2 ?？梢钥吹礁鞣N窗函數(shù)的形狀。矩形窗漢寧窗哈明窗Bar

18、telett 窗switch ii三角窗Blackman 窗Kaiser 窗（Beta=4）wChebyshev 窗（r=40）圖7-2各種窗函數(shù)的時間域形狀【例7-2】用MATLAB編程，采用512個頻率點繪制各種窗函數(shù)的幅頻特性。%Sa mp7_2clf;Nf=512;%窗函數(shù)復(fù)數(shù)頻率特性的數(shù)據(jù)點數(shù)Nwi n=20;%窗函數(shù)數(shù)據(jù)長度figured）for ii=1:4case 1w=boxcar(Nwi n);%矩形窗stext='矩形窗;case 2w=ha nnin g(Nw in);%漢寧窗stext='漢寧窗;case 3w=hammi ng (Nwi n);%哈明

19、窗stext='哈明窗;case 4w=bartlett(Nwi n);% Bartlett 窗stext='Bartelett 窗'endy,f=freqz(w,1,Nf); %求解窗函數(shù)的幅頻特性，窗函數(shù)相當于一個數(shù)字濾波器mag=abs(y);% 求得窗函數(shù)幅頻特性pos pl ot='2,2,' in t2str(ii);sub plot(pospl ot);plot(f/pi,20* Iog10(mag/max(mag);%繪制窗函數(shù)的幅頻特性xlabel('歸一化頻率');ylabel('振幅 /dB');ti

20、tle(stext);grid on;endfigure (2)for ii=1:4switch iicase 1w=blackma n(Nwin);%Blackma n 窗%三角窗%Kaiser 窗w=chebwi n( Nwi n,40);%Chebyshev 窗stext='Blackma n 窗'case 2w=tria ng(Nwi n);stext='三角窗'case 3w=kaiser(Nwi n,4);stext='Kaiser 窗(Beta=4)'case 4stext='Chebyshev 窗(r=40)'en

21、dy,f=freqz(w,1,Nf);%以Nf點數(shù)求解窗函數(shù)的幅頻響應(yīng)mag=abs(y);%求得窗函數(shù)幅頻響應(yīng)繪制幅頻響應(yīng)pos pl ot='2,2,' in t2str(ii);sub plot( pos plot); plot(f/pi,20* log10(mag/max(mag); % xlabel('歸一化頻率');ylabel('振幅 /dB');title(stext);grid on;end程序運行結(jié)果見圖7-3 ?？梢钥吹礁鞣N窗函數(shù)的幅頻形狀。對照該圖可知這些窗函數(shù)具有上面所分析的窗函數(shù)的特征。矩形窗漢寧窗歸一化頻率哈明窗-1

22、200'l1嚴0'1 '11. C'1 r,i ''；L1-20-40-60-80-1000.5歸一化頻率.廠'dn 1、廣、1 'II 1A1iII1-20-40-60-80-1000.5歸一化頻率Bartelett 窗歸一化頻率-12001Blackman 窗歸一化頻率-100011| ' 1,11 111，! ' 11Kaiser 窗（Beta=4）-20-40-60-800.5歸一化頻率1三角窗歸一化頻率Chebyshev 窗（r=40）歸一化頻率圖7-3各種窗函數(shù)的幅頻形狀由圖7-3 可見，各種窗函數(shù)都具

23、有明顯的主瓣(Main lobe )和旁瓣(Sidelobe )。主矩形窗具有最窄的主瓣，但也有最瓣頻寬和旁瓣的幅值衰減特性決定了窗函數(shù)的應(yīng)用場合。大的旁瓣峰值(第一旁瓣衰減為13 dB)； Blackman窗具有最大的旁瓣衰減，但也具有最寬的主瓣寬度。不同窗函數(shù)在這兩方面的特點是不同的，因此應(yīng)根據(jù)具體的問題進行選擇。是較為常用的窗函數(shù)。通常來講，哈明窗和漢寧窗的主瓣具有較小的旁瓣和較大的衰減速度,表7-2總結(jié)了各種窗函數(shù)主瓣和旁瓣的特征（理論分析可參考其他的數(shù)字信號處理教材），大家可對照窗函數(shù)的幅頻形狀（圖 7-3 ）認真理解體會。表7-2各種窗函數(shù)的特點窗函數(shù)主瓣寬第一旁瓣相對主瓣衰減（分

24、貝）矩形窗4-13漢寧窗8N-31哈明窗8N-41Bartlett 窗8N-25Blackma n窗12N-57三角窗8N-25Kaiser 窗可調(diào)整可調(diào)整Chebyshev窗可調(diào)整可調(diào)整主旁瓣頻率寬度還與窗函數(shù)長度N有關(guān)。增加窗函數(shù)長度 N將減小窗函數(shù)的主瓣寬度,但不能減小旁瓣幅值衰減的相對值（分貝數(shù)），這個值是由窗函數(shù)決定的。這個特點可由下面的例子清楚地看出。【例7-3】繪制矩形窗函數(shù)的幅頻響應(yīng)，窗長度分別為：（1）N=10;（2）N=20;(3)N=50;(4)N=100.%Sa mp7_3clf;Nf=512;for ii=1:4switch iicase 1Nwi n=10;case

25、 2Nwi n=20;case 3Nwi n=50;case 4Nwin=100;endw=boxcar(Nwi n);%矩形窗y,f=freqz(w,1,Nf);%用不同的窗長度求得復(fù)數(shù)頻率特性mag=abs(y);%求得幅頻特性pos pl ot='2,2,' in t2str(ii);%指定繪圖位置sub plot(pospl ot);plot (f/pi,20*log10(mag/max(mag);%繪出幅頻形狀xlabel('歸一化頻率');ylabel('振幅 /dB');stext='N=' i nt2str(Nwi

26、 n);%給出標題，指出所用的數(shù)據(jù)個數(shù)title(stext);grid on;endN=10歸一化頻率-700產(chǎn) -'11 I' 1'11 .1!1 "N=20-10-20-30-40-50-600.5歸一化頻率1N=50B幅 振歸一化頻率N=100歸一化頻率圖7-4數(shù)據(jù)長度不同的矩形窗的幅頻形狀程序運行結(jié)果見圖7-4。顯然，隨著 N的增大，主瓣和旁瓣都變窄，但第一旁瓣相對這樣，隨著主瓣的幅值下降分貝數(shù)相同，第二旁瓣相對第一旁瓣幅值下降的分貝數(shù)也相同。N的增大，旁瓣也得到抑制，有力地減少了頻譜泄漏，但不能完全消除。減少主瓣寬度和抑制旁瓣是一對矛盾，不可兼得，

27、只能根據(jù)不同用途折衷處理。723 運用窗函數(shù)設(shè)計數(shù)字濾波器用于信號分析中的窗函數(shù)可根據(jù)用戶的不同要求選擇。用于濾波器的窗函數(shù)，一般要求F面采用實例說明如何根據(jù)上面步驟設(shè)計FIR濾波器。窗函數(shù)的主瓣寬度窄， 以獲得較好的過渡帶； 旁瓣相對值盡可能少， 增加通帶的平穩(wěn)度和增 大阻帶的衰減?；诖昂瘮?shù)的FIR數(shù)字濾波器設(shè)計的算法十分簡單，其主要步驟為:(1)對濾波器理想頻域幅值響應(yīng)進行傅立葉逆變換獲得理想濾波器的單位脈沖響應(yīng)hd(n)。一般假定理想低通濾波器的截止頻率為c，其幅頻特性滿足(7-19)H ej則根據(jù)傅立葉逆變換，單位脈沖響應(yīng)為:其中,(2)函數(shù)類型hd( n)丄 cejndsin cn

28、2c為信號延遲。由性能指標(阻帶衰減的分貝數(shù))根據(jù)表w(n)。(7-20)7-2第3列的值確定滿足阻帶衰減的窗IIR濾波器一樣，F(xiàn)IR濾波器也要確定滿足性能指標的濾波器最小階數(shù)。由前面的討論(圖 7-1)濾波器的階數(shù)越高，濾波器的幅頻特性越好，但數(shù)據(jù)處理的費用也越高，因此像7-2可知，濾波器的主瓣寬度相當于過渡帶寬，因此，使過渡帶寬近似于窗函數(shù)主瓣寬(表中的第二列)可求得滿足性能指標的窗口長度N，此時，信號延遲為(N-1)/2。(3)求實際濾波器的單位脈沖響應(yīng)h(n):根據(jù)h(n)=h d(n)*w(n)。(4 )檢驗濾波器的性能?？稍O(shè)定一些信號采用7.1.2節(jié)指出的函數(shù)或 6.3.2節(jié)所給的

29、函數(shù)進行濾波。wdelta=ws-w p;%過渡帶寬【例7-4】 用窗函數(shù)設(shè)計一個線性相位FIR低通濾波器，并滿足性能指標：通帶邊界的歸一化頻率wp=0.5,阻帶邊界的歸一化頻率ws=0.66 ,阻帶衰減不小于 30dB ,通帶波紋不大于3dB。假設(shè)一個信號，其中f1=5Hz,f2=20Hz。信號的采樣頻率為 50Hz。試將原信號與通過濾波器的信號進行比較。由題意，阻帶衰減不小于30dB,根據(jù)表7-2,選取漢寧窗，因為漢寧窗的第一旁瓣相對主0 區(qū)間瓣衰減為31dB，滿足濾波要求。在窗函數(shù)設(shè)計法中，要求設(shè)計的頻率歸一化到內(nèi)，Nyquist頻率對應(yīng)于，因此通帶和阻帶邊界頻率為0.5 和0.66。程

30、序如下%Sa mp7_4wp=0.5* pi;ws=0.66* pi;%濾波器邊界頻率N=ceil(8* pi/wdelta)%根據(jù)過渡帶寬等于表7-2中漢寧窗函數(shù)主瓣寬求得濾波器所用窗函數(shù)的最小長度Nw=N;wc=(w p+ws)/2;%截止頻率在通帶和阻帶邊界頻率的中n=0:N-1;al pha=(N-1)/2;%求濾波器的相位延遲m=n-alp ha+e ps;%eps為MATLAB系統(tǒng)的精度hd=si n( wc*m)./( pi*m);%由(7-20 )式求理想濾波器脈沖響應(yīng)win=ha nnin g(Nw);%采用漢寧窗h=hd.*w in'%在時間域乘積對應(yīng)于頻率域的卷積

31、b=h;figure(1)H,f=freqz(b,1,512,50); % 采用50 Hz的采樣頻率繪出該濾波器的幅頻和相頻響應(yīng)sub plot(2,1,1), plot(f,20*log10(abs(H)xlabel('頻率/Hz');ylabel('振幅 /dB');grid on;sub plot(2,1,2), plot(f,180/pi*u nwrap(an gle(H)xlabel('頻率/Hz');ylabel(' 相位 /o');grid on;%im pz(b,1);%可采用此函數(shù)給出濾波器的脈沖響應(yīng)%zplan

32、 e(b,1);%可采用此語句給出濾波器的零極點圖%grpdelay(b,1);%可采用此函數(shù)給出濾波器的群延遲f1=3;f2=20;%檢測輸入信號含有兩種頻率成分dt=0.02; t=0:dt:3;%采樣間隔和檢測信號的時間序列x=s in (2* pi*f1* t)+cos(2* pi*f2* t);%檢測信號%y=filter(b,1,x);%可采用此函數(shù)給出濾波器的輸出y=fftfilt(b,x);%給出濾波器的輸出figure(2)sub plot(2,1,1), p lot(t,x),title('輸入信號')%繪輸入信號%繪輸出信號繪出延遲到的時刻sub plot

33、(2,1,2) ,p lot(t,y)hold on; p lot(1 1*(N-1)/2*dt,ylim, 'r') %xlabel('時間 /s'),title('輸出信號')程序運行結(jié)果見圖 7-5和圖7-6。該例設(shè)計通帶邊界wp=0.5,阻帶邊界頻率 ws=0.66 ， 對應(yīng)于 50Hz 的采樣頻率通帶邊界頻率為fp=Fs/2* Fn ormal=50/2*0.5=12.5Hz,fs=50/2*0.66=16.5Hz,其中Fs為采樣頻率，F(xiàn)normal為歸一化頻率。由圖 7-5上圖可以看到，在小于12.5Hz的頻段上，幾乎看不到下降，即滿

34、足通帶波紋不大于3dB的要求。在大于16.5Hz的頻段上，阻帶衰減大于30dB，滿足設(shè)計要求。由圖 7-5下圖可見，在通帶范圍內(nèi)，相位頻率為一條直線，表明該濾波器為線性相位。圖7-6給出了濾波器的輸入信號和輸出信號，輸入信號包括3Hz和20Hz的信號，由圖7-5可知，20Hz的信號不能通過該濾波器，通過濾波器后只剩下3Hz的信號，輸出結(jié)果也證明了這一點。但要注意,由于FIR濾波器所需的階數(shù)較高，信號延遲（N-1）/2也較大，輸出信號前面有一段直線就是延遲造成的。上述程序顯示的N取50才能滿足設(shè)計要求。這樣相位延遲為（N-1）/2*1/Fs=0.49s，可以看到輸出信號前面一段直線的距離大約為0

35、.49s。驗證了 FIR濾波器相位延遲的理論。在輸出信號的前部，有一些小信號，這是截斷信號邊界所致,后面的部分就沒有了這種信號。若采用零相位的filtfilt函數(shù)（說明見第六章第三節(jié)）輸出，則可最大限度地減小邊界的影響。50-150-'I 宀 r .'|II1"i t.0-50-1005201015頻率/Hz25圖7-5例7-4所設(shè)計濾波器的幅頻響應(yīng)（上圖）和相頻響應(yīng)（下圖）0-1輸入信號I' I I間:0.511.522.53輸岀信號圖7-6例7-4所設(shè)計濾波器的輸入和輸出信號-2匚07.2.4 標準型FIR濾波器7.2.3節(jié)給出了運用理想脈沖響應(yīng)與窗函數(shù)乘

36、積的方法給出了濾波器傳遞函數(shù)的設(shè)計方法。其實MATLAB已將上述方法復(fù)合成一個函數(shù)，提供基于上述原理設(shè)計標準型FIR濾波器的工具函數(shù)。fir1就是采用經(jīng)典窗函數(shù)法設(shè)計線性相位FIR數(shù)字濾波器的函數(shù)， 且具有標準低通，帶通，高通，帶阻等類型。函數(shù)調(diào)用格式為:b=fir1( n,wn ,ft yp e',w in dow)式中，n為FIR濾波器的階數(shù)，對于高通,帶阻濾波器，n需取偶數(shù)；wn為濾波器截止頻率,范圍為01 (歸一化頻率)。對于帶通，帶阻濾波器，wn=w1,w2(w1<w2);對于多帶濾波 器，女0 wn=w1 ,w2 ,w3,w4,頻率分段為：0<w<w1 ,

37、w1<w<w2 ,w2<w<w3 ,w3<w<w4 。'high'為高通濾波器； stoP' ftype'為濾波器的類型：缺省時為低通或帶通濾波器;為帶阻濾波器，'DC-1'為第一頻帶為通帶的多帶濾波器；'DC-0'為第一頻帶為阻帶的多帶濾波器。wi ndow為窗函數(shù)列向量，其長度為n+1。缺省時，自動取哈明窗。MATLAB提供的窗函數(shù)有 boxcar、hanning、hamming、bartlett 、blackman 、kaiser、chebwin ，調(diào)用方式見上節(jié)。b為FIR濾波器系數(shù)向量

38、,長度為n+1。FIR濾波器的傳遞函數(shù)具有下列形式:1 2b(z) b(1)b(2)zb(3)zb(n 1)z n(7-21)用函數(shù)fir1設(shè)計的FIR濾波器的群延遲為 n/2?？紤]到n階濾波器系數(shù)個數(shù)為N，即n+1，這里的延遲與前面所講的(N-1)/2的延遲一致。注意這里的濾波器的最小階數(shù)比窗函數(shù)的長度少1。7-4?！纠?-5】 用窗函數(shù)設(shè)計一個線性相位FIR低通濾波器，技術(shù)指標同上節(jié)例%Sa mp7_5wp=0.5* pi;ws=0.66* pi;%濾波器的邊界頻率wdelta=ws-w p;%過渡帶寬度N=ceil(8* pi/wdelta);%求解濾波器的最小階數(shù)，根據(jù)漢寧窗主瓣寬Wn

39、=(0.5+0.66)* pi/2;%截止頻率取通帶和阻帶邊界頻率的中點設(shè)計FIR濾波器,注意fir1要求輸入歸一化頻率b=fir1(N,Wn/pi,ha nnin g(N+1);%H,f=freqz(b,1,512,50); % 采用50Hz的采樣頻率求出頻率響應(yīng)sub plot(2,1,1), plot(f,20*log10(abs(H)xlabel('頻率/Hz');ylabel('振幅 /dB');grid on;sub plot(2,1,2), plot(f,180/pi*u nwrap(an gle(H)xlabel('頻率/Hz')

40、;ylabel(' 相位 /o');grid on;程序運行與上例中的圖 7-5 一致?！纠?-6】 設(shè)計一個48階FIR帶通濾波器，通帶邊界的歸一化頻率為0.35和0.65。假設(shè)50Hz。一個信號，其中含有f1=1Hz , f2=10Hz,f3=20Hz，三種頻率成分信號的采樣頻率為試將原信號與通過濾波器的信號進行比較。%Sa mp7_6wp=0.35 0.65;N=48;%通帶邊界頻率（歸一化頻率）和濾波器階數(shù)sub plot(2,1,2) ,p lot(t,y)%繪出輸出信號波形Fs=50;b=fir1(N,w p);%設(shè)計FIR帶通濾波器figured)H,f=freq

41、z(b,1,512,Fs);%以50Hz為采樣頻率求出濾波器頻率響應(yīng)sub plot(2,1,1), plot(f,20*log10(abs(H)xIabelC 頻率/Hz');ylabel('振幅 /dB');grid on;sub plot(2,1,2), plot(f,180/pi*u nwrap(an gle(H)xlabel('頻率/Hz');ylabel(' 相位 /o');grid on;f1=1;f2=10;f3=20;%輸入信號的三種頻率成分t=0:1/50:3;%時間序列x=s in (2* pi *f1*t)+0.5

42、*cos(2* pi*f2*t)+0.5*si n(2*p i*f3*t);%輸入信號%y=filter(b,1,x);%可以采用過濾器進行濾波y=fftfilt(b,x);%采用fftfilt對輸入信號濾波figure (2)sub plot(2,1,1), plot(t,x),title('輸入信號'）%繪出輸入信號波形hold on;p lot(N/2*0.02*o nes(1,2),ylim, 'r')%繪制延遲到的時刻titleC 輸出信號),xlabel('時間 /s')程序運行結(jié)果見圖 7-7和圖7-8。通帶歸一化頻率0.35和0.

43、65對應(yīng)于采樣頻率為50Hz的通帶范圍為 8.75Hz和16.25Hz（采用6-19式計算）。由圖7-7上圖可見滿足這一設(shè)計要求。在這個頻帶范圍內(nèi)的相位滿足線性相位，符合FIR濾波器的一般特點。圖7-8為檢測濾波器的輸入信號和輸出信號。輸入信號中含有伯z,10Hz 和20Hz 的信號。根據(jù)圖7-7上圖可知，1Hz和20Hz的頻率在阻帶范圍內(nèi),不能通過該濾波器，只有10Hz的信號可以N/2*0.02s=0.48s得至 U,通過該濾波器，輸出信號表明了這一點。濾波器的相位延遲根據(jù)輸出信號前面的直線部分大體為這個時間延遲，另外濾波后周期為 10Hz的信號相位（紅線開始部分），跟濾波前的相位（信號開始

44、部分）也一致，說明通過該濾波器濾波后沒有改變 信號的相位。50-100'|1.'J1! !111> 1:IIIIII0-505201015頻率/Hz25圖7-7例7-6設(shè)計濾波器的幅頻特性（上圖）和相頻特性（下圖）o位 相01.5-20.6輸入信號10.50-0.5-1-1.50.51.52.53輸岀信號0-0.400.40.20-0.2Illi0.5II1； !' i. !11II1111?I'll,圖7-8【例7-7】FIR低通濾波器階數(shù)為計此濾波器并對信號濾波器輸出的第81由于采樣頻率為i| I!1 !' 'll1.5時間/s2.5例

45、7-6濾波器的輸入信號和輸出信號40，截止頻率為200Hz，采樣頻率為 Fs=1000Hz。試設(shè)x(t)=si n(2* pi *f1*t)+si n(2* pi *f2*t)濾波，f1=50Hz,f2=250Hz ，選取個采樣點到第1000Hz ，241個采樣點之間的信號并與對應(yīng)的輸入信號進行比較。所以該濾波器的歸一化頻率的1對應(yīng)于 Nyquist 頻率500Hz，因此歸一化截止頻率為200/500（參看（6-19 ）式）。%Sa mp7_7clf;N=1000;Fs=1000;%數(shù)據(jù)總數(shù)和采樣頻率fc=200;n=0:N-1;t =n/Fs;%時間序列f1=50;f2=250;x=s in

46、 (2* pi *f1*t)+si n(2* pi *f2*t);%輸入信號b=fir1(40,fc*2/Fs);%設(shè)計40階的低通濾波器，歸一化截止頻率據(jù)6-19式y(tǒng)fft=fftfilt(b,x,256);%對數(shù)據(jù)進行濾波n1=81:241;t1=t( n1);%選擇采樣點間隔x1=x( n1);%與采樣點對應(yīng)的輸入信號sub plot(2,1,1); plot(t1,x1); grid on;%繪制輸入信號title（'輸入信號'）;n2=n1-40/2;t2=t( n2);%輸出信號，扣除了相位延遲N/2%繪制輸出信號y2=yfft( n2);sub plot(2,1,

47、2); plot(t2,y2);title（'輸出信號'）;grid on; xlabel(' 時間 /s')程序輸出結(jié)果見圖7-9。可見經(jīng)過濾波器的濾波，完全濾去了250Hz的高頻成分，只剩下50Hz的低頻成分。2輸入信號clf;N=300;Fs=1000;%數(shù)據(jù)總數(shù)和采樣頻率-0.08I'i11 rI'：.! rIU'l1IIHii1II1 ii 1!fI N'Ji1；i1)H 1 fNui i! ：i'i1 H! 1 !；. ll山 h!',) 1Il刖 ii! hl'汁口ll 1 1 !劇 '

48、;,i 1 H艸 H ： 11：Hl 抑 MuMtJH'；!'IE' JI1 , , 1 1FJi !'JH'i；,1,1'.!(iiiH H lly；III ' I' ili'卩ii1'n叩IIH !ii!：iI'1Ii1 ：; i ：i叩1II1'| 11.J1 1.1ii10-10.10.120.140.160.180.20.220.24輸岀信號-n丨;1. ,1 ! 1-A1 11 11 11i1 i h1 1'!1 :-fl1 1 ! 1 i :*. ti i-:八1 .f 1 11!1J ! i;11111'丨1:111!i ii1!1!1i11 11 111111 11111111 1;111111111111