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

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

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

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

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

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

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

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

9、點，這樣可使濾波器實際頻率響應更好地逼近理想頻率響應。這里我們明確兩個概念：截斷和頻譜泄漏。信號是無限長的，而在進行信號處理時只能采取有限長信號，所以需要將信號“截斷”。在信號處理中， “截斷”被看成是用一個有限長的“窗口”看無限長的信號，或者從分析的角度是無限長的信號x(t)乘以有限長的窗函數w(t)。由傅立葉變換性質可知,時間域內的乘積對應于頻率域的卷積,即 （7-9）這里，x(t)是頻寬有限信號，而w(t)是頻寬無限信號，表示互為Fourier變換對。截斷后的信號也必須是頻寬無限信號，這樣就是有限頻帶的信號分散到無限頻帶中去，這樣就產生了所謂頻譜泄漏。從能量的角度來看，頻譜泄漏也是能量的

10、泄漏，因為加窗后使原來信號集中的窄頻帶內的能量分散到無限的頻帶寬度范圍內。頻譜泄漏是不可避免的，但要盡量減小。上邊只考慮了矩形窗，如果我們使窗的主瓣寬度盡可能地窄，旁瓣盡可能地小，可以獲得性能更好的濾波器，能否改變窗的形狀而達到這個目的呢？回答是肯定的。其實數字信號處理的前驅者們設計了不同于矩形窗的很多窗函數，這些窗函數在主瓣和旁瓣特性方面各有特點，可滿足不同的要求。為此，用窗函數法設計FIR數字濾波器時，要根據給定的濾波器性能指標選擇窗口寬度N和窗函數w(n)。下面我們介紹窗函數。7.2.2 MATLAB信號處理中提供的窗函數（1）矩形窗：前面分析中所用的矩形窗可用下面函數來實現w=boxc

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

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

13、hebyshev窗具有主瓣寬度最小，而旁瓣等高,高度可調整的特點。下面我們在MATLAB觀看各種窗函數的形狀和頻率域圖象來驗證上述所講特點。【例7-1】 用MATLAB編程繪制各種窗函數的形狀。窗函數的長度為21。%Samp7_1clfNwin=21;n=0:Nwin-1; %數據總數和序列序號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); %指定繪制窗函數的圖形位置 subplot(posplot); stem(n,w); %繪出窗函數 hold on plot (n ,w,'r'); %繪制包絡線 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); %繪出窗函數 hold on plot (n,w,'r'); %繪出包絡線 xlabel('n');ylabel('w(n)');title(stext); hold off;grid on;end程序運行結果見圖 7-2 。可以看到各種窗函數的形狀。圖 7-2 各種窗函數的時間域形狀【例 7-2】 用 MATLAB 編程，采用512個頻率點繪制各種窗函數的幅頻特性。%Samp7_2clf;Nf=512; %窗函數復數頻率特性的數據點數Nwin=20; %

17、窗函數數據長度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); %求解窗函數的幅頻特性，窗函數相當于一個數字濾波器 mag

18、=abs(y);%求得窗函數幅頻特性 posplot='2,2,' int2str(ii); subplot(posplot); plot(f/pi,20* log10(mag/max(mag); %繪制窗函數的幅頻特性 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點數求解窗函數的幅頻響應 mag=abs(y); %求得窗函數幅頻響應 posplot='2,2,' int2str(ii); subplot(posplot);plot(f/pi,2

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

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

22、調整主旁瓣頻率寬度還與窗函數長度N有關。增加窗函數長度N將減小窗函數的主瓣寬度，但不能減小旁瓣幅值衰減的相對值（分貝數），這個值是由窗函數決定的。這個特點可由下面的例子清楚地看出。【例7-3】繪制矩形窗函數的幅頻響應，窗長度分別為：(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、用不同的窗長度求得復數頻率特性 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); %給出標題，指出所用的數據個數 title(stext);grid on;end圖 7-4 數據長度不同的矩形窗的幅頻形狀程序運行結果見圖7-4。顯然，隨

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

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

26、：根據h(n)=hd(n)*w(n)。（4）檢驗濾波器的性能。可設定一些信號采用 7.1.2 節指出的函數或6.3.2節所給的函數進行濾波。下面采用實例說明如何根據上面步驟設計FIR濾波器。【例 7-4】 用窗函數設計一個線性相位FIR低通濾波器，并滿足性能指標：通帶邊界的歸一化頻率wp=0.5,阻帶邊界的歸一化頻率ws=0.66，阻帶衰減不小于30dB，通帶波紋不大于3dB。假設一個信號，其中f1=5Hz,f2=20Hz。信號的采樣頻率為50Hz。試將原信號與通過濾波器的信號進行比較。由題意，阻帶衰減不小于30dB,根據表7-2,選取漢寧hanning窗，因為漢寧Hanning窗的第一旁瓣相

27、對主瓣衰減為31dB，滿足濾波要求。在窗函數設計法中，要求設計的頻率歸一化到0區間內，Nyquist頻率對應于，因此通帶和阻帶邊界頻率為0.5和0.66。程序如下%Samp7_4wp=0.5*pi;ws=0.66*pi; %濾波器邊界頻率wdelta=ws-wp; %過渡帶寬N=ceil(8*pi/wdelta) %根據過渡帶寬等于表 7-2中漢寧Hanning窗函數主瓣寬求得濾波器所用窗函數的最小長度Nw=N;wc=(wp+ws)/2; %截止頻率在通帶和阻帶邊界頻率的中點n=0:N-1;alpha=(N-1)/2; %求濾波器的相位延遲m=n-alpha+eps; %eps為MATLAB系

28、統的精度hd=sin(wc*m)./(pi*m); %由（7-20）式求理想濾波器脈沖響應win=hanning(Nw); %采用漢寧窗h=hd.*win' %在時間域乘積對應于頻率域的卷積b=h; figure(1)H,f=freqz(b,1,512,50); %采用 50 Hz 的采樣頻率繪出該濾波器的幅頻和相頻響應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(

29、H)xlabel('頻率/Hz');ylabel('相位/o');grid on;%impz(b,1); %可采用此函數給出濾波器的脈沖響應%zplane(b,1); %可采用此語句給出濾波器的零極點圖%grpdelay(b,1); %可采用此函數給出濾波器的群延遲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); %可采用此函數給出濾波器的輸出y=fftfilt(b,x);

30、%給出濾波器的輸出figure(2)subplot(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('輸出信號')程序運行結果見圖7-5和圖7-6。該例設計通帶邊界wp=0.5,阻帶邊界頻率ws=0.66，對應于50Hz的采樣頻率通帶邊界頻率為fp=Fs/2*Fnormal=50/2*0.5=12.5H

31、z, fs=50/2*0.66=16.5Hz, 其中Fs為采樣頻率，Fnormal為歸一化頻率。由圖7-5上圖可以看到，在小于12.5Hz的頻段上，幾乎看不到下降，即滿足通帶波紋不大于3dB的要求。在大于16.5Hz的頻段上，阻帶衰減大于30dB，滿足設計要求。由圖7-5下圖可見，在通帶范圍內，相位頻率為一條直線，表明該濾波器為線性相位。圖7-6給出了濾波器的輸入信號和輸出信號，輸入信號包括3Hz和20Hz的信號，由圖7-5可知，20Hz的信號不能通過該濾波器，通過濾波器后只剩下3Hz的信號，輸出結果也證明了這一點。但要注意，由于FIR濾波器所需的階數較高，信號延遲(N-1)/2也較大，輸出信

32、號前面有一段直線就是延遲造成的。上述程序顯示的N取50才能滿足設計要求。這樣相位延遲為(N-1)/2*1/Fs=0.49s，可以看到輸出信號前面一段直線的距離大約為0.49s。驗證了FIR濾波器相位延遲的理論。在輸出信號的前部，有一些小信號，這是截斷信號邊界所致，后面的部分就沒有了這種信號。若采用零相位的filtfilt函數（說明見第六章第三節）輸出，則可最大限度地減小邊界的影響。圖 7-5 例7-4所設計濾波器的幅頻響應（上圖）和相頻響應（下圖）圖7-6 例7-4所設計濾波器的輸入和輸出信號7.2.4 標準型FIR濾波器7.2.3節給出了運用理想脈沖響應與窗函數乘積的方法給出了濾波器傳遞函數

33、的設計方法。其實MATLAB已將上述方法復合成一個函數，提供基于上述原理設計標準型FIR濾波器的工具函數。fir1就是采用經典窗函數法設計線性相位FIR數字濾波器的函數，且具有標準低通，帶通，高通,帶阻等類型。函數調用格式為：b=fir1(n,wn,'ftype',window)式中，n為FIR濾波器的階數，對于高通,帶阻濾波器，n需取偶數；wn為濾波器截止頻率，范圍為01（歸一化頻率）。對于帶通,帶阻濾波器，wn=w1,w2(w1<w2)；對于多帶濾波器，如wn=w1，w2，w3,w4,頻率分段為：0<w<w1，w1<w<w2，w2<w&l

34、t;w3，w3<w<w4。 ftype'為濾波器的類型：缺省時為低通或帶通濾波器；'high'為高通濾波器；stop'為帶阻濾波器，'DC-1'為第一頻帶為通帶的多帶濾波器；'DC-0'為第一頻帶為阻帶的多帶濾波器。window為窗函數列向量，其長度為n+1。缺省時，自動取Hamming哈明窗。MATLAB提供的窗函數有boxcar、hanning、hamming、bartlett、blackman、kaiser、chebwin，調用方式見上節。b為FIR濾波器系數向量，長度為n+1。FIR濾波器的傳遞函數具有下列形式

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

36、hanning(N+1);%設計FIR濾波器,注意fir1要求輸入歸一化頻率H,f=freqz(b,1,512,50); %采用50Hz的采樣頻率求出頻率響應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;程序運行與上例中的圖7-5一致。【例 7-6】 設計一個4

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

38、l('振幅/dB');grid on;subplot(2,1,2),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); %可以采用過濾器進行濾波y=fftfilt(b,x); %采用 fftfilt 對輸入信號濾波fi

39、gure(2)subplot(2,1,1), plot(t,x),title('輸入信號')%繪出輸入信號波形subplot(2,1,2),plot(t,y) %繪出輸出信號波形hold on;plot(N/2*0.02*ones(1,2),ylim, 'r') %繪制延遲到的時刻title('輸出信號'),xlabel('時間/s')程序運行結果見圖7-7和圖7-8。通帶歸一化頻率0.35和0.65對應于采樣頻率為50Hz的通帶范圍為8.75Hz和16.25Hz(采用6-19式計算)。由圖7-7上圖可見滿足這一設計要求。在這個頻

40、帶范圍內的相位滿足線性相位，符合FIR濾波器的一般特點。圖7-8為檢測濾波器的輸入信號和輸出信號。輸入信號中含有1Hz,10Hz和20Hz的信號。根據圖7-7上圖可知，1Hz和20Hz的頻率在阻帶范圍內，不能通過該濾波器，只有10Hz的信號可以通過該濾波器，輸出信號表明了這一點。濾波器的相位延遲根據N/2*0.02s=0.48s得到，輸出信號前面的直線部分大體為這個時間延遲,另外濾波后周期為10Hz的信號相位（紅線開始部分），跟濾波前的相位（信號開始部分）也一致，說明通過該濾波器濾波后沒有改變信號的相位。圖 7-7 例 7-6 設計濾波器的幅頻特性（上圖）和相頻特性（下圖）圖 7-8 例 7-

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

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

43、2);subplot(2,1,2);plot(t2,y2); %繪制輸出信號title('輸出信號');grid on; xlabel('時間/s')程序輸出結果見圖7-9。可見經過濾波器的濾波，完全濾去了250Hz的高頻成分，只剩下50Hz的低頻成分。圖 7-9 例7-7設計濾波器的輸入信號和輸出信號【例7-8】設計采樣頻率為1000Hz,阻帶頻率從100Hz200Hz的100階的帶阻FIR濾波器。對信號x(t)=sin(2*pi*f1*t)+sin(2*pi*f2*t)濾波，f1=50Hz，f2=150Hz，并與對應的輸入信號進行比較。%Samp7_8clf

44、;N=300;Fs=1000; %數據總數和采樣頻率Order=100; %濾波器階數n=0:N-1;t=n/Fs; %時間序列wn=100 200/(Fs/2); %邊界頻率轉換為歸一化頻率，據6-19式b=fir1(Order,wn,'stop'); % 設計 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,18

45、0/pi*unwrap(angle(H)xlabel('頻率/Hz');ylabel('相位/o');grid on;f1=50;f2=150; %輸入信號頻率x=sin(2*pi*f1*t)+sin(2*pi*f2*t); %輸入信號y=fftfilt(b,x,256);% 對數據進行濾波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

46、),ylim, 'r') %繪制延遲到的時刻title('輸出信號');grid on; xlabel('時間/s')程序輸出結果見圖7-10和圖7-11。由圖7-10上圖可見，阻帶范圍為100200Hz，完全符合設計要求。在通帶的相位是線性的。由圖7-11可見，濾波器濾除了150Hz（在阻帶范圍內）的信號，保留了50Hz的信號。相位延遲100/2/Fs=0.05s，與圖形一致。圖 7-10 例7-8設計帶阻濾波器的幅頻（上圖）和相頻特性（下圖）圖 7-11例7-8設計濾波器的輸入和輸出信號該小節只給出了FIR低通，帶通和帶阻濾波器的例子，請大家

47、在課下自己設計高通，第一頻帶為通頻帶和第一頻帶為阻頻帶的多帶濾波器的例子，以加深對該函數的理解。7.2.5 多頻帶FIR濾波器除了設計標準型FIR濾波器外，MATLAB信號處理工具箱還提供另一種基于窗函數濾波器設計的工具函數fir2，用于設計具有任意形狀頻率響應的FIR濾波器，其調用格式為：b=fir2（n，f,m,npt,window）式中，n為濾波器的階數；f和m分別為濾波器期望幅頻響應的頻率向量和幅值向量，取值范圍為01(歸一化頻率)。m，f具有相同的長度,window為窗函數，得到列向量，長度必須為n+1。缺省時自動取hamming哈明窗；npt為對頻率響應進行內插的點數，缺省時為51

48、2。b為FIR濾波器系數向量，長度為n+1，濾波器具有與（7-21）式相同的形式。【例 7-9】 用窗函數設計一個多頻帶的FIR濾波器，濾波器階數分別為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,比較理想和實際濾波器的幅頻響應。假設一個信號，其中f1=12Hz,f2=36Hz。信號的采樣頻率為100Hz。試將原信號與通過濾波器的信號進行比較。%Samp7_9clff=0:0.1:1; %歸一化頻率點數m=0 0 1 1 0 0 1 1 1 0 0; %幅頻

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

50、f,m,hamming(Order+1); %設計階數為100的濾波器h,w=freqz(b,1,128); %計算濾波器的頻率響應subplot(2,1,2),plot(f,m,w/pi,abs(h),'r:'); %繪制理想幅頻響應和設計的幅頻響應ylim(0 1)legend('理想特性', '實際設計') %給出圖例title('Order=100');xlabel('歸一化頻率');ylabel('振幅');f1=12;f2=36; % 輸入信號的兩種頻率成分t=0:1/100:2; %

51、時間序列x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %輸入信號y=fftfilt(b,x); %對輸入信號進行濾波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')程序輸出結果見圖7-12和圖7-13。由該例可知，只有取100階時，實際濾波器的幅頻響應才逼近理想濾波器的幅頻響應。與第6章例6-12的輸出比較可知，相同性能的FIR濾波器階數比IIR濾波器要高得多。另外，該例中，信號含有兩種頻率成分，12Hz和36Hz，由于信號的采樣頻率為100Hz，因此Nyquist頻率為50Hz，則歸一化頻率的0.2對應于10Hz，0.3對應于15Hz，因此，0.20.3的通帶對應于1015Hz。同理，歸一化頻率中的0.60.8的通帶對應于3040Hz。可見12Hz和36Hz的波均在通帶的范圍內，因此均可通過。該例的輸出結果與