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

2、濾波器原理概述及濾波函數(shù)7.1.1 FIR濾波器原理及設(shè)計(jì)方法分類根據(jù)第 6 章對數(shù)字濾波器的介紹，我們知道FIR濾波器的傳遞函數(shù)為： (7-1)可得FIR濾波器的系統(tǒng)差分方程為： 因此，F(xiàn)IR濾波器又稱為卷積濾波器。根據(jù)第 4 章中所描述的系統(tǒng)頻率響應(yīng)，F(xiàn)IR濾波器的頻率響應(yīng)表達(dá)式為： （7-2）信號通過FIR濾波器不失真條件與（6-6）式所描述的相同，即濾波器在通帶內(nèi)具有恒定的幅頻特性和線性相位特性。理論上可以證明（這里從略）：當(dāng)FIR濾波器的系數(shù)滿足下列中心對稱條件： （7-3）時，濾波器設(shè)計(jì)在逼近平直幅頻特性的同時，還能獲得嚴(yán)格的線性相位特性。線性相位FIR濾波器的相位滯后和群延遲在整

3、個頻帶上是相等且不變的。對于一個 N 階的線性相位FIR濾波器，群延遲為常數(shù)，即濾波后的信號簡單地延遲常數(shù)個時間步長。這一特性使通帶頻率內(nèi)信號通過濾波器后仍保持原有波形形狀而無相位失真。本章主要介紹的FIR數(shù)字濾波器設(shè)計(jì)方法及 MATLAB 信號處理工具箱提供的了FIR數(shù)字濾波器設(shè)計(jì)函數(shù)，見表7-1。由于篇幅所限，本章我們主要介紹窗函數(shù)法和最優(yōu)化設(shè)計(jì)方法。表7-1 FIR濾波器設(shè)計(jì)的主要方法函數(shù)設(shè)計(jì)方法說明工具函數(shù)窗函數(shù)法理想濾波器加窗處理fir1(單頻帶) ， fir2(多頻帶) ， kaiserord最優(yōu)化設(shè)計(jì)平方誤差最小化逼近理想幅頻響應(yīng)或Park-McClellan 算法產(chǎn)生等波紋濾波

4、器firls ， remez,remezord約束最小二乘逼近在滿足最大誤差限制條件下使整個頻帶平方誤差最小化fircls,fircls1升余弦函數(shù)具有光滑、正弦過渡帶的低通濾波器設(shè)計(jì)Fircos7.1.2 FIR數(shù)字濾波器濾波函數(shù)相對于IIR 濾波器的濾波函數(shù)，F(xiàn)IR數(shù)字濾波器濾波函數(shù)除了dimpulse和dstep僅適用于IIR濾波器外，其他各種函數(shù)可直接應(yīng)用于FIR濾波器，只是輸入的分母多項(xiàng)式向量a=1。另外，MATLAB還提供了一個函數(shù)fftfilt,該函數(shù)利用效率高的基于FFT算法實(shí)現(xiàn)對數(shù)據(jù)的濾波，該函數(shù)只適用于FIR濾波器，調(diào)用形式為：y=fftfilt(b,x,n)式中，b為FI

5、R濾波器的系數(shù)向量；x為輸入數(shù)據(jù)；n為FFT長度，缺省時，函數(shù)選用最佳的FFT長度，y為濾波器的輸出。該函數(shù)執(zhí)行下面的操作：n=length(x);y=ifft(fft(x).*fft(b,n)./fft(a,n);應(yīng)注意，y=fftfilt(b,x)等價于y=filter(b,a,x)。7.2 FIR濾波器的窗函數(shù)設(shè)計(jì)7.2.1 窗函數(shù)的基本原理FIR濾波器設(shè)計(jì)的主要任務(wù)是根據(jù)給定的性能指標(biāo)確定濾波器的系數(shù)b，即系統(tǒng)單位脈沖序列h(n)，它是一個有限長序列。FIR濾波器的理想頻率響應(yīng)，可寫成復(fù)數(shù)形式的Fourier級數(shù)形式： (7-4)式中，hd(n)是對應(yīng)的單位脈沖響應(yīng)序列。這說明濾波器的

6、頻率響應(yīng)和單位脈沖響應(yīng)互為Fourier變換對。因此其單位脈沖響應(yīng)可由下式求得， (7-5)求得序列后，通過z變換，可得到 (7-6)注意，這里為無限長序列，因此是物理上不可實(shí)現(xiàn)的。如何變成物理上可實(shí)現(xiàn)呢？一個自然的想法是只取其中的某些項(xiàng)，即只截取中的一部分，比如n=0,N-1，N為正整數(shù)。這種處理相當(dāng)于將,n=-與函數(shù)w(n)相乘，w(n)具有下列形式：w(n)相當(dāng)于一個矩形，我們稱之為矩形窗。即我們可采用矩形窗函數(shù)w(n)將無限脈沖響應(yīng)截取一段h(n)來近似為，這種截取在數(shù)學(xué)上表示為： h(n)= w(n) (7-7)這里應(yīng)該強(qiáng)調(diào)的是，加窗函數(shù)不是可有可無的，而是將設(shè)計(jì)變?yōu)槲锢砜蓪?shí)現(xiàn)所必須

7、的。截取之后的濾波器傳遞函數(shù)變?yōu)椋?(7-8)式中，N為窗口寬度，H(z)是物理可實(shí)現(xiàn)系統(tǒng)。為了獲得線性相位，F(xiàn)IR濾波器h(n)必須滿足中心對稱條件（即7-3式），序列h(n)的延遲為。這種方法的基本原理是用一定寬度的矩形窗函數(shù)截取無限脈沖響應(yīng)序列獲得有限長的脈沖響應(yīng)序列，從而得到FIR濾波器的脈沖響應(yīng)，故稱為FIR濾波器的窗函數(shù)設(shè)計(jì)法。經(jīng)過加矩形窗后所得的濾波器實(shí)際頻率響應(yīng)能否很好地逼近理想頻率響應(yīng)呢？圖 7-1 示意給出了理想濾波器加矩形窗后的情況。理想低通濾波器的頻率響應(yīng)如圖中左上角圖，矩形窗的頻率響應(yīng)為左下角圖。時間域內(nèi)的乘積（7-7）式要求實(shí)際頻率響應(yīng)為這兩個頻率響應(yīng)函數(shù)在頻域內(nèi)的

8、卷積（卷積定理），即得到圖形為圖7-1（右圖）。圖 7-1 FIR濾波器理想與實(shí)際頻率響應(yīng)由圖可看出，加矩形窗后使實(shí)際頻率響應(yīng)偏離理想頻率響應(yīng)，主要影響有三個方面：（1）理想幅頻特性陡直邊緣處形成過渡帶，過渡帶寬取決于矩形窗函數(shù)頻率響應(yīng)的主瓣寬度。（2）過渡帶兩側(cè)形成肩峰和波紋，這是矩形窗函數(shù)頻率響應(yīng)的旁瓣引起的，旁瓣相對值越大，旁瓣越多，波紋越多。（3）隨窗函數(shù)寬度N的增大，矩形窗函數(shù)頻率響應(yīng)的主瓣寬度減小，但不改變旁瓣的相對值。為了改善FIR濾波器性能，要求窗函數(shù)的主瓣寬度盡可能窄，以獲得較窄的過渡帶；旁瓣相對值盡可能小，數(shù)量盡可能少，以獲得通帶波紋小,阻帶衰減大，在通帶和阻帶內(nèi)均平穩(wěn)的特

9、點(diǎn)，這樣可使濾波器實(shí)際頻率響應(yīng)更好地逼近理想頻率響應(yīng)。這里我們明確兩個概念：截?cái)嗪皖l譜泄漏。信號是無限長的，而在進(jìn)行信號處理時只能采取有限長信號，所以需要將信號“截?cái)唷薄Ｔ谛盘柼幚碇校?“截?cái)唷北豢闯墒怯靡粋€有限長的“窗口”看無限長的信號，或者從分析的角度是無限長的信號x(t)乘以有限長的窗函數(shù)w(t)。由傅立葉變換性質(zhì)可知,時間域內(nèi)的乘積對應(yīng)于頻率域的卷積,即 （7-9）這里，x(t)是頻寬有限信號，而w(t)是頻寬無限信號，表示互為Fourier變換對。截?cái)嗪蟮男盘栆脖仨毷穷l寬無限信號，這樣就是有限頻帶的信號分散到無限頻帶中去，這樣就產(chǎn)生了所謂頻譜泄漏。從能量的角度來看，頻譜泄漏也是能量的

10、泄漏，因?yàn)榧哟昂笫乖瓉硇盘柤械恼l帶內(nèi)的能量分散到無限的頻帶寬度范圍內(nèi)。頻譜泄漏是不可避免的，但要盡量減小。上邊只考慮了矩形窗，如果我們使窗的主瓣寬度盡可能地窄，旁瓣盡可能地小，可以獲得性能更好的濾波器，能否改變窗的形狀而達(dá)到這個目的呢？回答是肯定的。其實(shí)數(shù)字信號處理的前驅(qū)者們設(shè)計(jì)了不同于矩形窗的很多窗函數(shù)，這些窗函數(shù)在主瓣和旁瓣特性方面各有特點(diǎn)，可滿足不同的要求。為此，用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器時，要根據(jù)給定的濾波器性能指標(biāo)選擇窗口寬度N和窗函數(shù)w(n)。下面我們介紹窗函數(shù)。7.2.2 MATLAB信號處理中提供的窗函數(shù)（1）矩形窗：前面分析中所用的矩形窗可用下面函數(shù)來實(shí)現(xiàn)w=boxc

11、ar (N)，N 為窗的長度（以下函數(shù)與此同），w為返回的窗函數(shù)序列。（2）漢寧窗：w=hanning(N)漢寧窗的表達(dá)式為： (7-10)（3）哈明窗：w=hamming(N)哈明Hamming窗的表達(dá)式為： (7-11)（4）Bartlett窗：w=bartlett(N)Bartlett 窗的表達(dá)式為：當(dāng) N 為奇數(shù)時， (7-12)當(dāng) N 為偶數(shù)時， (7-13)（5） Blackman 窗：w= blackman(N)Blackman 窗的表達(dá)式為： (7-14)Blackman 窗比其他相同尺寸窗 (哈明Hamming窗,漢寧Hanning窗) 具有主瓣較寬和旁瓣泄漏較小的特點(diǎn)。（6

12、）三角窗：w=triang(N)三角窗的表達(dá)式為：當(dāng) N 為奇數(shù)時， (7-15)當(dāng) N 為偶數(shù)時， (7-16)三角窗和Bartlett窗十分類似。三角窗的兩端值不為零，而Bartlett窗則為零，這一點(diǎn)可從例7-1中看出。（7）Kaiser窗：w=kaiser(n,beta)其中，beta是Kaiser窗參數(shù),影響窗旁瓣幅值的衰減率。Kaiser窗表達(dá)式： (7-17)式中， I0.是修正過的零階 Bessel 函數(shù)。Kaiser窗用于濾波器設(shè)計(jì)時，若旁瓣幅值為，則 （ 7-18 ）（8） Chebyshev窗：w=chebwin(n,r)式中, r 是窗口的旁瓣幅值在主瓣以下的分貝數(shù)。C

13、hebyshev窗具有主瓣寬度最小，而旁瓣等高,高度可調(diào)整的特點(diǎn)。下面我們在MATLAB觀看各種窗函數(shù)的形狀和頻率域圖象來驗(yàn)證上述所講特點(diǎn)?！纠?-1】 用MATLAB編程繪制各種窗函數(shù)的形狀。窗函數(shù)的長度為21。%Samp7_1clfNwin=21;n=0:Nwin-1; %數(shù)據(jù)總數(shù)和序列序號figure(1)for ii=1:4 switch ii case 1 w=boxcar(Nwin); %矩形窗 stext='矩形窗' case 2 w=hanning(Nwin); %漢寧窗 stext='漢寧窗' case 3 w=hamming(Nwin); %

14、哈明窗 stext='哈明窗' case 4 w=bartlett(Nwin); %Bbartlett窗 stext='Bartelett窗' end posplot='2,2,' int2str(ii); %指定繪制窗函數(shù)的圖形位置 subplot(posplot); stem(n,w); %繪出窗函數(shù) hold on plot (n ,w,'r'); %繪制包絡(luò)線 xlabel('n'); ylabel('w(n)'); title(stext); hold off; grid onendfig

15、ure(2)for ii=1:4 switch ii case 1 w=blackman(Nwin); %Blackman 窗 stext='Blackman窗' case 2 w=triang(Nwin); %三角窗 stext='三角窗' case 3 w=kaiser(Nwin,4); %Kaiser窗 stext='Kaiser窗(Beta=4)' case 4 w=chebwin(Nwin,40); %Chebyshev 窗 stext='Chebyshev窗(r=40)' end posplot='2,2,&#

16、39; int2str(ii); subplot(posplot); stem(n,w); %繪出窗函數(shù) hold on plot (n,w,'r'); %繪出包絡(luò)線 xlabel('n');ylabel('w(n)');title(stext); hold off;grid on;end程序運(yùn)行結(jié)果見圖 7-2 ?？梢钥吹礁鞣N窗函數(shù)的形狀。圖 7-2 各種窗函數(shù)的時間域形狀【例 7-2】 用 MATLAB 編程，采用512個頻率點(diǎn)繪制各種窗函數(shù)的幅頻特性。%Samp7_2clf;Nf=512; %窗函數(shù)復(fù)數(shù)頻率特性的數(shù)據(jù)點(diǎn)數(shù)Nwin=20; %

17、窗函數(shù)數(shù)據(jù)長度figure(1)for ii=1:4 switch ii case 1 w=boxcar(Nwin); %矩形窗 stext='矩形窗' case 2 w=hanning(Nwin); %漢寧窗 stext='漢寧窗' case 3 w=hamming (Nwin); %哈明窗 stext='哈明窗' case 4 w=bartlett(Nwin); % Bartlett窗 stext='Bartelett窗' end y,f=freqz(w,1,Nf); %求解窗函數(shù)的幅頻特性，窗函數(shù)相當(dāng)于一個數(shù)字濾波器 mag

18、=abs(y);%求得窗函數(shù)幅頻特性 posplot='2,2,' int2str(ii); subplot(posplot); plot(f/pi,20* log10(mag/max(mag); %繪制窗函數(shù)的幅頻特性 xlabel('歸一化頻率');ylabel('振幅/dB'); title(stext);grid on;endfigure(2)for ii=1:4 switch ii case 1 w=blackman(Nwin); %Blackman 窗 stext='Blackman窗' case 2 w=triang

19、(Nwin); %三角窗 stext='三角窗' case 3 w=kaiser(Nwin,4); %Kaiser窗 stext='Kaiser窗(Beta=4)' case 4 w=chebwin(Nwin,40); %Chebyshev 窗 stext='Chebyshev窗(r=40)' end y,f=freqz(w,1,Nf); %以 Nf點(diǎn)數(shù)求解窗函數(shù)的幅頻響應(yīng) mag=abs(y); %求得窗函數(shù)幅頻響應(yīng) posplot='2,2,' int2str(ii); subplot(posplot);plot(f/pi,2

20、0* log10(mag/max(mag); %繪制幅頻響應(yīng) xlabel('歸一化頻率');ylabel('振幅/dB'); title(stext);grid on;end程序運(yùn)行結(jié)果見圖 7-3 ?？梢钥吹礁鞣N窗函數(shù)的幅頻形狀。對照該圖可知這些窗函數(shù)具有上面所分析的窗函數(shù)的特征。圖 7-3 各種窗函數(shù)的幅頻形狀由圖 7-3 可見，各種窗函數(shù)都具有明顯的主瓣（ Mainlobe ）和旁瓣（ Sidelobe ）。主瓣頻寬和旁瓣的幅值衰減特性決定了窗函數(shù)的應(yīng)用場合。矩形窗具有最窄的主瓣，但也有最大的旁瓣峰值（第一旁瓣衰減為 13 dB）；Blackman窗具有

21、最大的旁瓣衰減，但也具有最寬的主瓣寬度。不同窗函數(shù)在這兩方面的特點(diǎn)是不同的，因此應(yīng)根據(jù)具體的問題進(jìn)行選擇。通常來講，哈明Hamming窗和漢寧Hanning窗的主瓣，具有較小的旁瓣和較大的衰減速度，是較為常用的窗函數(shù)。表7-2總結(jié)了各種窗函數(shù)主瓣和旁瓣的特征（理論分析可參考其他的數(shù)字信號處理教材），大家可對照窗函數(shù)的幅頻形狀（圖7-3）認(rèn)真理解體會。表7-2 各種窗函數(shù)的特點(diǎn)窗函數(shù)主瓣寬第一旁瓣相對主瓣衰減（分貝）矩形窗-13漢寧Hanning窗-31哈明Hamming窗-41Bartlett窗-25Blackman 窗-57三角窗-25Kaiser窗可調(diào)整可調(diào)整Chebyshev 窗可調(diào)整可

22、調(diào)整主旁瓣頻率寬度還與窗函數(shù)長度N有關(guān)。增加窗函數(shù)長度N將減小窗函數(shù)的主瓣寬度，但不能減小旁瓣幅值衰減的相對值（分貝數(shù)），這個值是由窗函數(shù)決定的。這個特點(diǎn)可由下面的例子清楚地看出。【例7-3】繪制矩形窗函數(shù)的幅頻響應(yīng)，窗長度分別為：(1)N=10;(2)N=20; (3)N=50;(4)N=100.%Samp7_3clf;Nf=512;for ii=1:4 switch ii case 1 Nwin=10; case 2 Nwin=20; case 3 Nwin=50; case 4 Nwin=100; end w=boxcar(Nwin); %矩形窗 y,f=freqz(w,1,Nf); %

23、用不同的窗長度求得復(fù)數(shù)頻率特性 mag=abs(y); %求得幅頻特性 posplot='2,2,' int2str(ii); %指定繪圖位置 subplot(posplot); plot (f/pi,20*log10(mag/max(mag); %繪出幅頻形狀 xlabel('歸一化頻率');ylabel('振幅/dB'); stext='N=' int2str(Nwin); %給出標(biāo)題，指出所用的數(shù)據(jù)個數(shù) title(stext);grid on;end圖 7-4 數(shù)據(jù)長度不同的矩形窗的幅頻形狀程序運(yùn)行結(jié)果見圖7-4。顯然，隨

24、著N的增大，主瓣和旁瓣都變窄，但第一旁瓣相對主瓣的幅值下降分貝數(shù)相同，第二旁瓣相對第一旁瓣幅值下降的分貝數(shù)也相同。這樣，隨著N的增大，旁瓣也得到抑制，有力地減少了頻譜泄漏，但不能完全消除。減少主瓣寬度和抑制旁瓣是一對矛盾，不可兼得，只能根據(jù)不同用途折衷處理。7.2.3 運(yùn)用窗函數(shù)設(shè)計(jì)數(shù)字濾波器用于信號分析中的窗函數(shù)可根據(jù)用戶的不同要求選擇。用于濾波器的窗函數(shù)，一般要求窗函數(shù)的主瓣寬度窄，以獲得較好的過渡帶；旁瓣相對值盡可能少，增加通帶的平穩(wěn)度和增大阻帶的衰減。基于窗函數(shù)的FIR數(shù)字濾波器設(shè)計(jì)的算法十分簡單，其主要步驟為：(1)對濾波器理想頻域幅值響應(yīng)進(jìn)行傅立葉逆變換獲得理想濾波器的單位脈沖響應(yīng)

25、hd(n)。一般假定理想低通濾波器的截止頻率為，其幅頻特性滿足 (7-19)則根據(jù)傅立葉逆變換，單位脈沖響應(yīng)為： (7-20)其中，為信號延遲。（2）由性能指標(biāo)（阻帶衰減的分貝數(shù)）根據(jù)表7-2第3列的值確定滿足阻帶衰減的窗函數(shù)類型w(n)。濾波器的階數(shù)越高，濾波器的幅頻特性越好，但數(shù)據(jù)處理的費(fèi)用也越高，因此像IIR濾波器一樣，F(xiàn)IR濾波器也要確定滿足性能指標(biāo)的濾波器最小階數(shù)。由前面的討論（圖7-1）可知，濾波器的主瓣寬度相當(dāng)于過渡帶寬，因此，使過渡帶寬近似于窗函數(shù)主瓣寬（表7-2中的第二列）可求得滿足性能指標(biāo)的窗口長度N，此時，信號延遲為(N-1)/2。（3）求實(shí)際濾波器的單位脈沖響應(yīng)h(n)

26、：根據(jù)h(n)=hd(n)*w(n)。（4）檢驗(yàn)濾波器的性能?？稍O(shè)定一些信號采用進(jìn)行濾波。下面采用實(shí)例說明如何根據(jù)上面步驟設(shè)計(jì)FIR濾波器?！纠?7-4】 用窗函數(shù)設(shè)計(jì)一個線性相位FIR低通濾波器，并滿足性能指標(biāo)：通帶邊界的歸一化頻率wp=0.5,阻帶邊界的歸一化頻率ws=0.66，阻帶衰減不小于30dB，通帶波紋不大于3dB。假設(shè)一個信號，其中f1=5Hz,f2=20Hz。信號的采樣頻率為50Hz。試將原信號與通過濾波器的信號進(jìn)行比較。由題意，阻帶衰減不小于30dB,根據(jù)表7-2,選取漢寧hanning窗，因?yàn)闈h寧Hanning窗的第一旁瓣相對主瓣衰減為31dB，滿足濾波要求。在窗函數(shù)設(shè)計(jì)法

27、中，要求設(shè)計(jì)的頻率歸一化到0區(qū)間內(nèi)，Nyquist頻率對應(yīng)于，因此通帶和阻帶邊界頻率為0.5和0.66。程序如下%Samp7_4wp=0.5*pi;ws=0.66*pi; %濾波器邊界頻率wdelta=ws-wp; %過渡帶寬N=ceil(8*pi/wdelta) %根據(jù)過渡帶寬等于表 7-2中漢寧Hanning窗函數(shù)主瓣寬求得濾波器所用窗函數(shù)的最小長度Nw=N;wc=(wp+ws)/2; %截止頻率在通帶和阻帶邊界頻率的中點(diǎn)n=0:N-1;alpha=(N-1)/2; %求濾波器的相位延遲m=n-alpha+eps; %eps為MATLAB系統(tǒng)的精度hd=sin(wc*m)./(pi*m);

28、 %由（7-20）式求理想濾波器脈沖響應(yīng)win=hanning(Nw); %采用漢寧窗h=hd.*win' %在時間域乘積對應(yīng)于頻率域的卷積b=h; figure(1)H,f=freqz(b,1,512,50); %采用 50 Hz 的采樣頻率繪出該濾波器的幅頻和相頻響應(yīng)subplot(2,1,1),plot(f,20*log10(abs(H)xlabel('頻率/Hz');ylabel('振幅/dB');grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel('頻率/Hz')

29、;ylabel('相位/o');grid on;%impz(b,1); %可采用此函數(shù)給出濾波器的脈沖響應(yīng)%zplane(b,1); %可采用此語句給出濾波器的零極點(diǎn)圖%grpdelay(b,1); %可采用此函數(shù)給出濾波器的群延遲f1=3;f2=20; %檢測輸入信號含有兩種頻率成分dt=0.02; t=0:dt:3; %采樣間隔和檢測信號的時間序列x=sin(2*pi*f1* t)+cos(2* pi*f2* t); %檢測信號%y=filter(b,1,x); %可采用此函數(shù)給出濾波器的輸出y=fftfilt(b,x); %給出濾波器的輸出figure(2)subplot

30、(2,1,1), plot(t,x),title('輸入信號') %繪輸入信號subplot(2,1,2),plot(t,y) % 繪輸出信號hold on; plot(1 1*(N-1)/2*dt,ylim, 'r') %繪出延遲到的時刻xlabel('時間/s'),title('輸出信號')程序運(yùn)行結(jié)果見圖7-5和圖7-6。該例設(shè)計(jì)通帶邊界wp=0.5,阻帶邊界頻率ws=0.66，對應(yīng)于50Hz的采樣頻率通帶邊界頻率為fp=Fs/2*Fnormal=50/2*0.5=12.5Hz, fs=50/2*0.66=16.5Hz, 其

31、中Fs為采樣頻率，F(xiàn)normal為歸一化頻率。由圖7-5上圖可以看到，在小于12.5Hz的頻段上，幾乎看不到下降，即滿足通帶波紋不大于3dB的要求。在大于16.5Hz的頻段上，阻帶衰減大于30dB，滿足設(shè)計(jì)要求。由圖7-5下圖可見，在通帶范圍內(nèi)，相位頻率為一條直線，表明該濾波器為線性相位。圖7-6給出了濾波器的輸入信號和輸出信號，輸入信號包括3Hz和20Hz的信號，由圖7-5可知，20Hz的信號不能通過該濾波器，通過濾波器后只剩下3Hz的信號，輸出結(jié)果也證明了這一點(diǎn)。但要注意，由于FIR濾波器所需的階數(shù)較高，信號延遲(N-1)/2也較大，輸出信號前面有一段直線就是延遲造成的。上述程序顯示的N取

32、50才能滿足設(shè)計(jì)要求。這樣相位延遲為(N-1)/2*1/Fs=0.49s，可以看到輸出信號前面一段直線的距離大約為0.49s。驗(yàn)證了FIR濾波器相位延遲的理論。在輸出信號的前部，有一些小信號，這是截?cái)嘈盘栠吔缢?，后面的部分就沒有了這種信號。若采用零相位的filtfilt函數(shù)（說明見第六章第三節(jié)）輸出，則可最大限度地減小邊界的影響。圖 7-5 例7-4所設(shè)計(jì)濾波器的幅頻響應(yīng)（上圖）和相頻響應(yīng)（下圖）圖7-6 例7-4所設(shè)計(jì)濾波器的輸入和輸出信號7.2.4 標(biāo)準(zhǔn)型FIR濾波器節(jié)給出了運(yùn)用理想脈沖響應(yīng)與窗函數(shù)乘積的方法給出了濾波器傳遞函數(shù)的設(shè)計(jì)方法。其實(shí)MATLAB已將上述方法復(fù)合成一個函數(shù)，提供

33、基于上述原理設(shè)計(jì)標(biāo)準(zhǔn)型FIR濾波器的工具函數(shù)。fir1就是采用經(jīng)典窗函數(shù)法設(shè)計(jì)線性相位FIR數(shù)字濾波器的函數(shù)，且具有標(biāo)準(zhǔn)低通，帶通，高通,帶阻等類型。函數(shù)調(diào)用格式為：b=fir1(n,wn,'ftype',window)式中，n為FIR濾波器的階數(shù)，對于高通,帶阻濾波器，n需取偶數(shù)；wn為濾波器截止頻率，范圍為01（歸一化頻率）。對于帶通,帶阻濾波器，wn=w1,w2(w1<w2)；對于多帶濾波器，如wn=w1，w2，w3,w4,頻率分段為：0<w<w1，w1<w<w2，w2<w<w3，w3<w<w4。 ftype'

34、為濾波器的類型：缺省時為低通或帶通濾波器；'high'為高通濾波器；stop'為帶阻濾波器，'DC-1'為第一頻帶為通帶的多帶濾波器；'DC-0'為第一頻帶為阻帶的多帶濾波器。window為窗函數(shù)列向量，其長度為n+1。缺省時，自動取Hamming哈明窗。MATLAB提供的窗函數(shù)有boxcar、hanning、hamming、bartlett、blackman、kaiser、chebwin，調(diào)用方式見上節(jié)。b為FIR濾波器系數(shù)向量，長度為n+1。FIR濾波器的傳遞函數(shù)具有下列形式： (7-21)用函數(shù)fir1設(shè)計(jì)的FIR濾波器的群延遲為n

35、/2?？紤]到n階濾波器系數(shù)個數(shù)為N，即n+1，這里的延遲與前面所講的(N-1)/2的延遲一致。注意這里的濾波器的最小階數(shù)比窗函數(shù)的長度少1?！纠?-5】 用窗函數(shù)設(shè)計(jì)一個線性相位FIR低通濾波器，技術(shù)指標(biāo)同上節(jié)例7-4。%Samp7_5wp=0.5*pi;ws=0.66*pi; %濾波器的邊界頻率wdelta=ws-wp; %過渡帶寬度N=ceil(8* pi/wdelta);%求解濾波器的最小階數(shù)，根據(jù)漢寧Hanning 窗主瓣寬Wn=(0.5+0.66)*pi/2;%截止頻率取通帶和阻帶邊界頻率的中點(diǎn)b=fir1(N,Wn/pi,hanning(N+1);%設(shè)計(jì)FIR濾波器,注意fir1要

36、求輸入歸一化頻率H,f=freqz(b,1,512,50); %采用50Hz的采樣頻率求出頻率響應(yīng)subplot(2,1,1),plot(f,20*log10(abs(H)xlabel('頻率/Hz');ylabel('振幅/dB');grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel('頻率/Hz');ylabel('相位/o');grid on;程序運(yùn)行與上例中的圖7-5一致。【例 7-6】 設(shè)計(jì)一個48階FIR帶通濾波器，通帶邊界的歸一化頻率假設(shè)一個信號，其中

37、含有f1=1Hz，f2=10Hz,f3=20Hz，三種頻率成分信號的采樣頻率為50Hz。試將原信號與通過濾波器的信號進(jìn)行比較。%Samp7_6wp=0.35 0.65;N=48; %通帶邊界頻率（歸一化頻率）和濾波器階數(shù)Fs=50;b=fir1(N,wp); %設(shè)計(jì)FIR帶通濾波器figure(1)H,f=freqz(b,1,512,Fs); %以50Hz為采樣頻率求出濾波器頻率響應(yīng)subplot(2,1,1),plot(f,20*log10(abs(H)xlabel('頻率/Hz');ylabel('振幅/dB');grid on;subplot(2,1,2)

38、,plot(f,180/pi*unwrap(angle(H)xlabel('頻率/Hz');ylabel('相位/o');grid on;f1=1;f2=10;f3=20; %輸入信號的三種頻率成分t=0:1/50:3; %時間序列x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t);%輸入信號%y=filter(b,1,x); %可以采用過濾器進(jìn)行濾波y=fftfilt(b,x); %采用 fftfilt 對輸入信號濾波figure(2)subplot(2,1,1), plot(t,x),title(&#

39、39;輸入信號')%繪出輸入信號波形subplot(2,1,2),plot(t,y) %繪出輸出信號波形hold on;plot(N/2*0.02*ones(1,2),ylim, 'r') %繪制延遲到的時刻title('輸出信號'),xlabel('時間/s')程序運(yùn)行結(jié)果見圖7-7和圖7-8。通帶歸一化頻率對應(yīng)于采樣頻率為50Hz的通帶范圍為8.75Hz和16.25Hz(采用6-19式計(jì)算)。由圖7-7上圖可見滿足這一設(shè)計(jì)要求。在這個頻帶范圍內(nèi)的相位滿足線性相位，符合FIR濾波器的一般特點(diǎn)。圖7-8為檢測濾波器的輸入信號和輸出信號。輸

40、入信號中含有1Hz,10Hz和20Hz的信號。根據(jù)圖7-7上圖可知，1Hz和20Hz的頻率在阻帶范圍內(nèi)，不能通過該濾波器，只有10Hz的信號可以通過該濾波器，輸出信號表明了這一點(diǎn)。濾波器的相位延遲根據(jù)N/2*0.02s=0.48s得到，輸出信號前面的直線部分大體為這個時間延遲,另外濾波后周期為10Hz的信號相位（紅線開始部分），跟濾波前的相位（信號開始部分）也一致，說明通過該濾波器濾波后沒有改變信號的相位。圖 7-7 例 7-6 設(shè)計(jì)濾波器的幅頻特性（上圖）和相頻特性（下圖）圖 7-8 例 7-6 濾波器的輸入信號和輸出信號【例7-7】FIR低通濾波器階數(shù)為40，截止頻率為200Hz，采樣頻率

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

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

43、出信號');grid on; xlabel('時間/s')程序輸出結(jié)果見圖7-9?？梢娊?jīng)過濾波器的濾波，完全濾去了250Hz的高頻成分，只剩下50Hz的低頻成分。圖 7-9 例7-7設(shè)計(jì)濾波器的輸入信號和輸出信號【例7-8】設(shè)計(jì)采樣頻率為1000Hz,阻帶頻率從100Hz200Hz的100階的帶阻FIR濾波器。對信號x(t)=sin(2*pi*f1*t)+sin(2*pi*f2*t)濾波，f1=50Hz，f2=150Hz，并與對應(yīng)的輸入信號進(jìn)行比較。%Samp7_8clf;N=300;Fs=1000; %數(shù)據(jù)總數(shù)和采樣頻率Order=100; %濾波器階數(shù)n=0:N-1

44、;t=n/Fs; %時間序列wn=100 200/(Fs/2); %邊界頻率轉(zhuǎn)換為歸一化頻率，據(jù)6-19式b=fir1(Order,wn,'stop'); % 設(shè)計(jì) 100 階的帶阻濾波器figure(1)H,f=freqz(b,1,512,Fs);subplot(2,1,1),plot(f,20*log10(abs(H)xlabel('頻率/Hz');ylabel('振幅/dB');grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel('頻率/Hz');ylabel

45、('相位/o');grid on;f1=50;f2=150; %輸入信號頻率x=sin(2*pi*f1*t)+sin(2*pi*f2*t); %輸入信號y=fftfilt(b,x,256);% 對數(shù)據(jù)進(jìn)行濾波figure(2)subplot(2,1,1);plot(t,x); grid on; %繪制輸入信號title('輸入信號');subplot(2,1,2);plot(t,y); %繪制輸出信號hold on;plot(Order/2/Fs*ones(1,2),ylim, 'r') %繪制延遲到的時刻title('輸出信號'

46、);grid on; xlabel('時間/s')程序輸出結(jié)果見圖7-10和圖7-11。由圖7-10上圖可見，阻帶范圍為100200Hz，完全符合設(shè)計(jì)要求。在通帶的相位是線性的。由圖7-11可見，濾波器濾除了150Hz（在阻帶范圍內(nèi)）的信號，保留了50Hz的信號。相位延遲100/2/Fs=0.05s，與圖形一致。圖 7-10 例7-8設(shè)計(jì)帶阻濾波器的幅頻（上圖）和相頻特性（下圖）圖 7-11例7-8設(shè)計(jì)濾波器的輸入和輸出信號該小節(jié)只給出了FIR低通，帶通和帶阻濾波器的例子，請大家在課下自己設(shè)計(jì)高通，第一頻帶為通頻帶和第一頻帶為阻頻帶的多帶濾波器的例子，以加深對該函數(shù)的理解。7.

47、2.5 多頻帶FIR濾波器除了設(shè)計(jì)標(biāo)準(zhǔn)型FIR濾波器外，MATLAB信號處理工具箱還提供另一種基于窗函數(shù)濾波器設(shè)計(jì)的工具函數(shù)fir2，用于設(shè)計(jì)具有任意形狀頻率響應(yīng)的FIR濾波器，其調(diào)用格式為：b=fir2（n，f,m,npt,window）式中，n為濾波器的階數(shù)；f和m分別為濾波器期望幅頻響應(yīng)的頻率向量和幅值向量，取值范圍為01(歸一化頻率)。m，f具有相同的長度,window為窗函數(shù)，得到列向量，長度必須為n+1。缺省時自動取hamming哈明窗；npt為對頻率響應(yīng)進(jìn)行內(nèi)插的點(diǎn)數(shù)，缺省時為512。b為FIR濾波器系數(shù)向量，長度為n+1，濾波器具有與（7-21）式相同的形式?！纠?7-9】 用

48、窗函數(shù)設(shè)計(jì)一個多頻帶的FIR濾波器，濾波器階數(shù)分別為10和100，濾波器的特性同上一章例6-12，即f=0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0,m=0 0 1 1 0 0 1 1 1 0 0,比較理想和實(shí)際濾波器的幅頻響應(yīng)。假設(shè)一個信號，其中f1=12Hz,f2=36Hz。信號的采樣頻率為100Hz。試將原信號與通過濾波器的信號進(jìn)行比較。%Samp7_9clff=0:0.1:1; %歸一化頻率點(diǎn)數(shù)m=0 0 1 1 0 0 1 1 1 0 0; %幅頻特性值Order=10; % 濾波器的階數(shù)b=fir2(Order,f,m,hamming(Orde

49、r+1); %設(shè)計(jì)濾波器h,w=freqz(b,1,128); %計(jì)算濾波器的頻率響應(yīng)subplot(2,1,1)plot(f,m,w/pi,abs(h),'r:') %繪制理想幅頻響應(yīng)和設(shè)計(jì)的濾波器幅頻響應(yīng)legend('理想特性', '實(shí)際設(shè)計(jì)') %給出圖例title('Order=10');xlabel('歸一化頻率');ylabel('振幅');Order=100;b=fir2(Order,f,m,hamming(Order+1); %設(shè)計(jì)階數(shù)為100的濾波器h,w=freqz(b,1,1

50、28); %計(jì)算濾波器的頻率響應(yīng)subplot(2,1,2),plot(f,m,w/pi,abs(h),'r:'); %繪制理想幅頻響應(yīng)和設(shè)計(jì)的幅頻響應(yīng)ylim(0 1)legend('理想特性', '實(shí)際設(shè)計(jì)') %給出圖例title('Order=100');xlabel('歸一化頻率');ylabel('振幅');f1=12;f2=36; % 輸入信號的兩種頻率成分t=0:1/100:2; % 時間序列x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %輸入信號y=ff

51、tfilt(b,x); %對輸入信號進(jìn)行濾波figure(2)subplot(2,1,1), plot(t,x),title('輸入信號') %繪制輸入信號subplot(2,1,2),plot(t,y) %繪制輸出信號hold on;plot(Order/2/100*ones(1,2),ylim, 'r') %繪制延遲到的時刻title('輸出信號'),xlabel('時間/s')程序輸出結(jié)果見圖7-12和圖7-13。由該例可知，只有取100階時，實(shí)際濾波器的幅頻響應(yīng)才逼近理想濾波器的幅頻響應(yīng)。與第6章例6-12的輸出比較可知，相同性能的FIR濾波器階數(shù)比IIR濾波器要高得多。另外，該例中，信號含有兩種頻率成分，12Hz和36Hz，由于信號的采樣頻率為100Hz，因此Nyquist頻率為50Hz，則歸一化頻率的0.2對應(yīng)于10Hz，0.3對應(yīng)于15Hz，因此，0.20.3的通帶對應(yīng)于1015Hz。同理，歸一化頻率中的0.60.8的通帶對應(yīng)于3040Hz。可見12Hz和36Hz的波均在通帶的范圍內(nèi)，因此均可通過。該例的輸出結(jié)果與例6-12的輸出結(jié)果比較可知，F(xiàn)IR濾波器的相位延遲比IIR濾波器的相位延遲大得多。圖 7-12 例