FIR濾波器窗函數(shù)設(shè)計

1、課題名稱：FIR濾波器窗函數(shù)設(shè)計FIR濾波器窗函數(shù)設(shè)計引言：數(shù)字濾波器(Digital Filter)是指輸入、輸出都是離散時間信號，通過一定運算關(guān)系改變輸入信號所含頻率成分的相對比例或者濾除某些頻率成分的器件。在許多數(shù)字信號處理系統(tǒng)中，如圖像信號處理等，有限沖激響應(yīng)（FIR）濾波器是最常用的組件之一，它完成信號預(yù)調(diào)、頻帶選擇和濾波等功能。FIR濾波器雖然在截止頻率的邊沿陡峭性能上不及無限沖激響應(yīng)（IIR）濾波器，但是卻具有嚴格的線性相位特性，穩(wěn)定性好，能設(shè)計成多通帶（或多阻帶）濾波器組，所以能夠在數(shù)字信號處理領(lǐng)域得到廣泛的應(yīng)用。一、 數(shù)字濾波器的分類（1） 根據(jù)系統(tǒng)響應(yīng)函數(shù)的時間特性分為兩類

2、1. FIR（Finite Impulse Response）數(shù)字濾波器網(wǎng)絡(luò)特點：不存在反饋支路，其單位沖激響應(yīng)為有限長。2. IIR（Infinite Impulse Response）數(shù)字濾波器網(wǎng)絡(luò)特點：存在反饋支路，即信號流圖中存在環(huán)路，其單位沖激響應(yīng)為無限長。（2） FIR數(shù)字濾波器IIR數(shù)字濾波器的區(qū)別1. 從性能上來說，IIR濾波器傳遞函數(shù)包括零點和極點兩組可調(diào)因素，對極點的惟一限制是在單位圓內(nèi)。因此可用較低的階數(shù)獲得高的選擇性，所用的存儲單元少，計算量小，效率高。但是這個高效率是以相位的非線性為代價的。選擇性越好，則相位非線性越嚴重。FIR濾波器傳遞函數(shù)的極點固定在原點，是不能動

3、的，它只能靠改變零點位置來改變它的性能。所以要達到高的選擇性，必須用較高的階數(shù)；對于同樣的濾波器設(shè)計指標，F(xiàn)IR濾波器所要求的階數(shù)可能比IIR濾波器高5-10倍，但是 FIR濾波器可以得到嚴格的線性相位。2. 從結(jié)構(gòu)上看，IIR濾波器必須采用遞歸結(jié)構(gòu)，極點位置必須在單位圓內(nèi)，否則系統(tǒng)將不穩(wěn)定。相反，F(xiàn)IR濾波器只要采用非遞歸結(jié)構(gòu)，不論在理論上還是在實際的有限精度運算中都不存在穩(wěn)定性問題，因此造成的頻率特性誤差也較小。此外FIR濾波器可以采用快速傅里葉變換算法，在相同階數(shù)的條件下，運算速度可以快得多。3. 從設(shè)計工具看，IIR濾波器可以借助于模擬濾波器的成果，因此一般都有有效的封閉形式的設(shè)計公式

4、可供準確計算，計算工作量比較小，對計算工具的要求不高。從上面的簡單比較可以看到IIR與FIR濾波器各有所長，所以在實際應(yīng)用時應(yīng)該從多方面考慮來加以選擇。從使用要求上來看，在對相位要求不敏感的場合，如語言通信等，選用IIR較為合適，這樣可以充分發(fā)揮其經(jīng)濟高效的特點；對于圖像信號處理，數(shù)據(jù)傳輸?shù)纫圆ㄐ螖y帶信息的系統(tǒng)，則對線性相位要求較高。如果有條件，采用FIR濾波器較好。另外，不論IIR和FIR，階數(shù)越高，信號延遲越大；同時在IIR濾波器中，階數(shù)越高，系數(shù)的精度要求越高，否則很容易造成有限字長的誤差使極點移到單位園外。因此在階數(shù)選擇上是綜合考慮的。二、 FIR數(shù)字濾波器的設(shè)計(一) FIR濾波器原

5、理FIR濾波器的系統(tǒng)輸入輸出差分方程為：所以FIR濾波器的系統(tǒng)函數(shù)為：由于FIR濾波器的單位脈沖響應(yīng)h(n)是一個有限長序列，H (z)是Z-1的(N1)次多項式，它在Z平面上有(N1)個零點，同時在原點有(N1)階重極點。因此，H(z)永遠穩(wěn)定。FIR濾波器設(shè)計的任務(wù)是選擇有限長度的h(n)，使傳輸函數(shù)H(ej)滿足一定的幅度特性和線性相位要求。由于FIR濾波器很容易實現(xiàn)嚴格的線性相位，所以FIR 數(shù)字濾波器設(shè)計的核心思想是求出有限的脈沖響應(yīng)來逼近給定的頻率響應(yīng)。(二) FIR濾波器的窗函數(shù)設(shè)計原理窗函數(shù)法的設(shè)計思想是按照所要求的理想濾波器頻率響應(yīng)Hd(ej)，設(shè)計一個FIR濾波器，使之頻率

6、響應(yīng)Hej=來逼近Hd(ej)。先由Hd(ej)的傅里葉反變換導出理想濾波器的沖激響應(yīng)序列hd(n)，即：hdn=12-Hd(ej)ejnd由于Hd(ej)是矩形頻率特性，所以hdn是一無限長的序列，且是非因果的，而要計的FIR濾波器的沖激響應(yīng)序列是有限長的，所以要用有限長的序列h(n)來逼近無限長的序列hdn，最有效的方法是截斷hdn，或者說用一個有限長度的窗口函數(shù)w(n)序列來截取hdn，即: hn=w(n)hdn。按照復卷積公式，在時域中的乘積關(guān)系可表示成在頻域中的周期性卷積關(guān)系，即可得所設(shè)計的FIR濾波器的頻率響應(yīng):H(ej)=12-HdejW(ej（-）)d其中，W(ej)為截斷窗函

7、數(shù)的頻率特性。由此可見，實際的FIR數(shù)字濾波器的頻率響應(yīng)H(ej)逼近理想濾波器頻率響應(yīng)Hd(ej)的好壞，完全取決于窗函數(shù)的頻率特性W(ej)。如果w(n)具有下列形式：w(n)相當于一個矩形，我們稱之為矩形窗。即我們可采用矩形窗函數(shù)w(n)將無限脈沖響應(yīng)hdn截取一段hn來近似為hdn。經(jīng)過加矩形窗后所得的濾波器實際頻率響應(yīng)能否很好地逼近理想頻率響應(yīng)呢？下圖給出了理想濾波器加矩形窗后的情況。理想低通濾波器的頻率響應(yīng)Hd(ej)如圖中左上角圖，矩形窗的頻率響應(yīng)W(ej)為左下角圖。根據(jù)卷積定理，即得實際濾波器的頻率響應(yīng)H(ej)圖形為圖中右圖。由圖可看出，加矩形窗后使實際頻率響應(yīng)偏離理想頻率

8、響應(yīng)，主要影響有三個方面：（1）理想幅頻特性陡直邊緣處形成過渡帶，過渡帶寬取決于矩形窗函數(shù)頻率響應(yīng)的主瓣寬度。（2）過渡帶兩側(cè)形成肩峰和波紋，這是矩形窗函數(shù)頻率響應(yīng)的旁瓣引起的，旁瓣相對值越大，旁瓣越多，波紋越多。（3）隨窗函數(shù)寬度N的增大，矩形窗函數(shù)頻率響應(yīng)的主瓣寬度減小，但不改變旁瓣的相對值。為了改善濾波器的性能，需使窗函數(shù)譜滿足：主瓣盡可能窄，以使設(shè)計出來的濾波器有較陡的過渡帶；第一副瓣面積相對主瓣面積盡可能小，即能量盡可能集中在主瓣，外泄少，使設(shè)計出來的濾波器的肩峰和余振小逼近于理想濾波器。但是這兩個條件是相互矛盾的，實際應(yīng)用中，折衷處理，兼顧各項指標。(三) MATLAB信號處理中提

9、供的窗函數(shù)上邊只考慮了矩形窗，如果我們使窗的主瓣寬度盡可能地窄，旁瓣盡可能地小，可以獲得性能更好的濾波器，通過改變窗的形狀來達到這個目的。在數(shù)字信號處理的發(fā)展過程中形成了不同于矩形窗的很多窗函數(shù)，這些窗函數(shù)在主瓣和旁瓣特性方面各有特點，可滿足不同的要求。為此，用窗函數(shù)法設(shè)計FIR數(shù)字濾波器時，要根據(jù)給定的濾波器性能指標選擇窗口寬度N和窗函數(shù)w(n)。下面具體介紹幾類類窗函數(shù)及其特性。1. 矩形窗矩形窗函數(shù)的時域形式可以表示為：它的頻域特性為：2. 漢寧窗函數(shù)漢寧窗函數(shù)的時域形式可以表示為： 它的頻域特性為： 其中，為矩形窗函數(shù)的幅度頻率特性函數(shù)。漢寧窗函數(shù)的最大旁瓣值比主瓣值低31dB，但是主

10、瓣寬度比矩形窗函數(shù)的主瓣寬度增加了1倍，為8/N。3. 海明窗函數(shù)海明窗函數(shù)的時域形式可以表示為: 它的頻域特性為:其中，為矩形窗函數(shù)的幅度頻率特性函數(shù)。海明窗函數(shù)的最大旁瓣值比主瓣值低41dB，但它和漢寧窗函數(shù)的主瓣寬度是一樣大的。4. 三角窗函數(shù)三角窗是最簡單的頻譜函數(shù)W(ej)為非負的一種窗函數(shù)。三角窗函數(shù)的時域形式可以表示為：當n為奇數(shù)時：當n為偶數(shù)時：它的頻域特性為：三角窗函數(shù)的主瓣寬度為8/N，比矩形窗函數(shù)的主瓣寬度增加了一倍，但是它的旁瓣寬度卻小得多。(四) 各窗函數(shù)的圖形及幅頻特性1. 用MATLAB編程繪制各種窗函數(shù)的形狀，窗函數(shù)的長度為21。程序如下：clf;Nwin=21

11、;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); %海明窗 stext='海明窗'case 4 w=triang(Nwin); %三角窗 stext='三角窗' end posplot='2,2,' int2str(ii); subplot(posplot); stem(n,

12、w); %繪出窗函數(shù) hold on plot (n,w,'r'); %繪出包絡(luò)線 xlabel('n');ylabel('w(n)');title(stext); hold off;grid on;end程序運行結(jié)果如下圖：2. 用MATLAB編程，采用512個頻率點繪制各窗函數(shù)的幅頻特性。程序如下：clf;Nf=512; %窗函數(shù)復數(shù)頻率特性的數(shù)據(jù)點數(shù)Nwin=20; %窗函數(shù)數(shù)據(jù)長度figure(1)for ii=1:4 switch ii case 1 w=boxcar(Nwin); %矩形窗 stext='矩形窗' ca

13、se 2 w=hanning(Nwin); %漢寧窗 stext='漢寧窗'case 3 w=hamming(Nwin); %海明窗 stext='海明窗'case 4 w=triang(Nwin); %三角窗 stext='三角窗' end y,f=freqz(w,1,Nf); %求解窗函數(shù)的幅頻特性 mag=abs(y); %求得窗函數(shù)幅頻特性 posplot='2,2,' int2str(ii); subplot(posplot); plot(f/pi,20* log10(mag/max(mag); %繪制窗函數(shù)的幅頻特性

14、xlabel('歸一化頻率');ylabel('振幅/dB'); title(stext);grid on;end程序運行結(jié)果如下圖：由上圖各窗函數(shù)的幅頻特性可以看到各種窗函數(shù)都具有明顯的主瓣和旁瓣。主瓣頻寬和旁瓣的幅值衰減特性決定了窗函數(shù)的應(yīng)用場合。矩形窗具有最窄的主瓣，但也有最大的旁瓣峰值（第一旁瓣衰減為13dB）；不同窗函數(shù)在這兩方面的特點是不同的，因此應(yīng)根據(jù)具體的問題進行選擇。通常來講，海明窗和漢寧窗的主瓣，具有較小的旁瓣和較大的衰減速度，是較為常用的窗函數(shù)。下表列出了各種窗函數(shù)主瓣和旁瓣的特征：各窗函數(shù)的特征表窗函數(shù)主瓣頻寬第一旁瓣相對主瓣衰減（dB）

15、矩形窗4N-13漢寧窗8N-31海明窗8N-41三角窗8N-25此外，主旁瓣頻率寬度還與窗函數(shù)長度N有關(guān)。增加窗函數(shù)長度N將減小窗函數(shù)的主瓣寬度，但不能減小旁瓣幅值衰減的相對值（分貝數(shù)），這個值是由窗函數(shù)決定的。例如：繪制矩形窗函數(shù)的幅頻響應(yīng)，窗長度分別為：(1)N=10;(2)N=20; (3)N=50;(4)N=100時的圖形如下：由上圖可以看出，隨著N的增大，主瓣和旁瓣都變窄，但第一旁瓣相對主瓣的幅值下降分貝數(shù)相同，第二旁瓣相對第一旁瓣幅值下降的分貝數(shù)也相同。然而，隨著N的增大，旁瓣數(shù)也增多，減少主瓣寬度和抑制旁瓣是一對矛盾，不可兼得，只能根據(jù)不同用途折衷處理。三、 運用窗函數(shù)設(shè)計數(shù)字濾

16、波器用于信號分析中的窗函數(shù)可根據(jù)濾波器的指標進行不同選擇。用于濾波器的窗函數(shù)，一般要求窗函數(shù)的主瓣寬度窄，以獲得較好的過渡帶；旁瓣相對值盡可能少，增加通帶的平穩(wěn)度和增大阻帶的衰減。窗函數(shù)設(shè)計法的基本原理是用一定寬度窗函數(shù)截取無限脈沖響應(yīng)序列獲得有限長的脈沖響應(yīng)序列，主要設(shè)計步驟為：（1）通過傅里葉逆變換獲得理想濾波器的單位脈沖響應(yīng)hdn；（2）由性能指標（阻帶衰減的分貝數(shù)）根據(jù)窗函數(shù)的特征表的值確定滿足阻帶衰減的窗函數(shù)類型w(n)，和窗口長度N；（2）根據(jù)hn=w(n)hdn，求實際濾波器的單位脈沖響應(yīng)hn；（3）檢驗濾波器的性能。下面根據(jù)上面步驟設(shè)計FIR濾波器：一、例1要求用窗函數(shù)設(shè)計一個

17、線性相位FIR低通濾波器，并滿足性能指標：通帶邊界的歸一化頻率wp=0.5,阻帶邊界的歸一化頻率ws=0.66，阻帶衰減不小于30dB，通帶波紋不大于3dB。假設(shè)一個信號，其中f1=5Hz,f2=20Hz。信號的采樣頻率為50Hz。試將原信號與通過濾波器的信號進行比較。根據(jù)要求：阻帶衰減不小于30dB, 根據(jù)窗函數(shù)的特征表選取漢寧窗滿足濾波要求；在窗函數(shù)設(shè)計法中，要求設(shè)計的頻率歸一化到0區(qū)間內(nèi)，Nyquist頻率對應(yīng)于，因此通帶和阻帶邊界頻率為0.5和0.66。程序如下：wp=0.5*pi;ws=0.66*pi; %濾波器邊界頻率wdelta=ws-wp; %過渡帶寬N=ceil(8*pi/w

18、delta) %根據(jù)過渡帶寬得濾波器所用窗函數(shù)的最小長度Nw=N;wc=(wp+ws)/2; %理想低通濾波器截止頻率在通帶和阻帶邊界頻率的中點n=0:N-1;alpha=(N-1)/2; %求濾波器的相位延遲m=n-alpha+eps; %eps為MATLAB系統(tǒng)的精度hd=sin(wc*m)./(pi*m); %求理想濾波器脈沖響應(yīng)win=hanning(Nw); %采用漢寧窗h=hd.*win' %在時間域乘積對應(yīng)于頻率域的卷積b=h; figure(1)H,f=freqz(b,1,512,50);%采用50Hz的采樣頻率繪出濾波器的幅頻和相頻響應(yīng)subplot(2,1,1),p

19、lot(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;f1=5;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); %可采用

20、此函數(shù)給出濾波器的輸出y=fftfilt(b,x); %對信號進行濾波器的輸出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('輸出信號')程序運行結(jié)果如下圖：圖一：低通濾波器幅頻和相頻響應(yīng)圖二：輸入輸出信號結(jié)果分析：該設(shè)計的要求通帶邊界wp=0.5,阻帶邊界頻率ws=

21、0.66，對應(yīng)于50Hz的采樣頻率通帶邊界頻率為fp=Fs/2*Fnormal=50/2*0.5=12.5Hz， fs=50/2*0.66=16.5Hz，其中Fs為采樣頻率，F(xiàn)normal為歸一化頻率。由圖一中濾波器幅頻響應(yīng)可以看到，在小于12.5Hz的頻段上，幾乎看不到下降，即滿足通帶波紋不大于3dB的要求。在大于16.5Hz的頻段上，阻帶衰減大于30dB，滿足設(shè)計要求。由圖一中相頻響應(yīng)可見，在通帶范圍內(nèi)，相位頻率為一條直線，表明該濾波器為線性相位。圖二給出了濾波器的輸入信號和輸出信號，輸入信號包括5Hz和20Hz的信號，由圖一可知，20Hz的信號不能通過該濾波器，通過濾波器后只剩下5Hz的

22、信號，輸出結(jié)果也證明了這一點。但是由于FIR濾波器所需的階數(shù)較高，信號延遲(N-1)/2也較大，輸出信號前面有一段直線就是延遲造成的。上述程序顯示的N取50才能滿足設(shè)計要求。這樣相位延遲為(N-1)/2*1/Fs=0.49s，可以看到輸出信號前面一段直線的距離大約為0.49s。驗證了FIR濾波器相位延遲的理論。二、標準型FIR濾波器上面的設(shè)計中運用理想脈沖響應(yīng)與窗函數(shù)乘積的方法給出了濾波器傳遞函數(shù)的設(shè)計方法。其實MATLAB已將上述方法復合成一個函數(shù)，提供基于上述原理設(shè)計標準型FIR濾波器的工具函數(shù)。fir1就是采用經(jīng)典窗函數(shù)法設(shè)計線性相位FIR數(shù)字濾波器的函數(shù)，且具有標準低通、帶通、高通、帶

23、阻等類型。函數(shù)調(diào)用格式為：b=fir1(n,wn,'ftype',window)式中，n為FIR濾波器的階數(shù)，對于高通,帶阻濾波器，n需取偶數(shù)；wn為濾波器截止頻率，范圍為01（歸一化頻率）。對于帶通,帶阻濾波器，wn=w1,w2(w1<w2)； ftype'為濾波器的類型：缺省時為低通或帶通濾波器；'high'為高通濾波器；stop'為帶阻濾波器，'DC-1'為第一頻帶為通帶的多帶濾波器；'DC-0'為第一頻帶為阻帶的多帶濾波器。window為窗函數(shù)列向量，其長度為n+1。缺省時，自動取海明窗。MATLAB

24、提供的窗函數(shù)有boxcar、hanning、hamming、bartlett、blackman、kaiser、chebwin。b為FIR濾波器系數(shù)向量，長度為n+1。FIR濾波器的傳遞函數(shù)具有下列形式：用函數(shù)fir1設(shè)計的FIR濾波器的群延遲為n/2?？紤]到n階濾波器系數(shù)個數(shù)為N，即n+1，這里的延遲與前面所講的(N-1)/2的延遲一致。注意這里的濾波器的最小階數(shù)比窗函數(shù)的長度少1。運用MATLAB工具函數(shù)設(shè)計上面FIR濾波器的程序如下：wp=0.5*pi;ws=0.66*pi;%濾波器的邊界頻率wdelta=ws-wp; %過渡帶寬度N=ceil(8* pi/wdelta); %求解濾波器的

25、最小階數(shù)，根據(jù)漢寧窗求主瓣寬Wn=(0.5+0.66)*pi/2; %截止頻率取通帶和阻帶邊界頻率的中點b=fir1(N,Wn/pi,hanning(N+1); %設(shè)計FIR濾波器,注意fir1要求輸入歸一化頻率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('頻率/H

26、z');ylabel('相位/o');grid on;程序運行結(jié)果與圖一中的濾波器幅頻和相頻響應(yīng)是一致的。例二，設(shè)計一個48階FIR帶通濾波器，通帶邊界的歸一化頻率為0.35和0.65。假設(shè)一個信號，其中含有f1=1Hz，f2=10Hz,f3=20Hz，三種頻率成分信號的采樣頻率為50Hz。試將原信號與通過濾波器的信號進行比較。編寫程序如下：wp=0.35 0.65;N=48; %通帶邊界頻率（歸一化頻率）和濾波器階數(shù)Fs=50; %采樣頻率50HZb=fir1(N,wp); %設(shè)計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),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*s