FIR濾波器很多工科出身的人都不會陌生,在我們的學習和工作中,也常常需要設計FIR濾波器。因為FIR濾波器有兩個特點:濾波器是穩定的以及具有線性相位。FIR濾波器在信號處理相關領域當然也包括本人所在的雷達信號處理領域有著廣泛的應用。本文主要介紹MATLAB最常用的FIR濾波器設計方法之窗函數法。其他的方法將在另一章中介紹。
窗函數法是一種基礎且普遍應用的FIR濾波器設計方法。首先需要根據性能指標(如主瓣寬度、旁瓣衰減等)確定適合的窗函數。 主瓣寬度、旁瓣衰減是一對情敵,想要主瓣寬度窄且旁瓣衰減大,那是電視劇裡都不會出現的情況。實際中,需要根據自己的任務指標權衡。此外,還需要確定階數。然後就可以用fir1函數設計濾波器了。
b=fir1(n,wn,'ftype',window)
其中:
b:我們設計的fir濾波器係數,長度為n+1;b跟過渡帶的寬度有關,設計時根據性能要求確定。
n:濾波器的階數。注意,b的長度為n+1。
wn:濾波器的截止頻率,可以是一個標量或者多元素的向量。取值範圍0<wn<1,wn = 1對應於奈奎斯特採樣頻率(採樣頻率/2)。wn是單個值時,為低通/高通濾波器,ftype可以為low/high;當wn為有兩個元素的向量[w1 w2 ],w1 < w2 , 為帶通/帶阻濾波器,ftype可以為bandpass/stop; 如果Wn有兩個以上元素[w1 w2 …… wn],w1< w2 <…… <wn,ftype可以為'DC-0' | 'DC-1'。ftype為'DC-0'代表第一個帶(0~w1 )為阻帶(系統默認),ftype為'DC-1'代表第一個帶(0~w1 )為通帶。wn對應於濾波器歸一化增益-6dB處。
window:表示使用的窗函數,最常用的是漢明窗(Hamming)、漢寧窗(Hanning)、三角窗(bartlett、triang)、矩形窗(boxcar)、布萊克曼窗(Blackman)、chebwin(切比雪夫窗)、凱賽窗(Kaiser);默認是漢明窗(Hamming)。各種窗的差別主要在集中於主瓣的能量和分散在所有旁瓣的能量之比。
例如我們需要設計一個50階,截止頻率ω = 0.3π,使用漢明窗的低通濾波器。
b = fir1(50,0.3,'low',hamming(51));
freqz(b,1,512)
注意:窗函數長度和濾波器係數b的長度應一致。
改一下需求,需要設計一個50階,通帶為 0.3π <ω <0.6π,旁瓣衰減50dB,採用切比雪夫窗的帶通濾波器。
b = fir1(50,[0.3 0.6],'bandpass',chebwin(51,50));
freqz(b,1,512)
對這個濾波器做個小測試。
fs = 200;%採樣頻率
f1 = 10;
f2 = 50;
f3 = 80;
t = linspace(0,1,fs);
x = 2*sin(2*pi*f1*t)+5*sin(2*pi*f2*t)+3*sin(2*pi*f3*t);
plot(abs(fft(x,512)));
title('原始信號頻譜')
b = fir1(50,[0.3 0.6],'bandpass',chebwin(51,50));
y = filter(b,1,x);
figure;
plot(abs(fft(y,512)));
title('濾波後信號頻譜')
可見頻率為10Hz和80Hz的分量被濾除掉了。
用漢明窗(Hamming)、漢寧窗(Hanning)設計呢?如下圖所示。
再來說一下ftype為'DC-0' | 'DC-1'的時候。舉兩個例子。
設計一個46階,阻帶為ω <0.4π,0.6π <ω <0.9π的濾波器。
bnd = [0.4 0.6 0.9];
b_0 = fir1(46,bnd,'DC-0');
fvtool(b_0,1);
如果我們需要把上面的阻帶換為通帶呢?binggo!把ftype:DC-0換成ftype:DC-1 就可以了。
把兩個濾波器畫在一起,如下圖:
clc;clear;
bnd = [0.4 0.6 0.9];
b_0 = fir1(46,bnd,'DC-0');
b_1 = fir1(46,bnd,'DC-1');
h = fvtool(b_0,1,b_1,1);
legend(h,'DC-0','DC-1')
kaiser窗
Kaiser窗是一種最優化窗,具有很好的旁瓣抑制性能。
[n,Wn,beta,ftype] = kaiserord(f,a,dev,fs)
f是一個矢量,表示帶的邊緣頻率。f的長度是a的長度*2-2。關於f和a定義了一個分段常數響應函數。下面會通過例子說明。範圍為0~fs/2。
dev指定le通帶紋波和阻帶衰減。表示每個頻帶輸出濾波器的頻率響應與其期望值之間最大允許的偏差。
fs:採樣頻率Hz為單位。
比如設計一個通帶為0~1KHz,阻帶為1.5~4KHz,5%的通帶紋波和阻帶衰減為50dB的低通濾波器。
fsamp = 8000;
fcuts = [1000 1500];
mags = [1 0];
devs = [0.05 0.0035];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fsamp);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
fvtool(hh,1);
再如設計一個帶阻濾波器。
fsamp = 8000;
fcuts = [1000 1500 2500 3000 ];
mags = [1 0 1];
devs = [0.05 0.0035 0.05];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fsamp);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
fvtool(hh,1);
下面我們分析一下這個濾波器。fsamp ,採樣頻率。
fcuts ,帶邊緣頻率。fcuts 應小於fsamp/2。
mags = [1 0 1];
我們從上圖中可以看到在0<f<1000Hz時,濾波器為通帶。1000<f<1500時是過渡帶,沒有mags 對應。在1500<f<2500時,濾波器為阻帶,在2500<f<3000時,是過渡帶,沒有mags 對應。在3000<f<4000時,濾波器為通帶。mags 的作用就是定義某一段頻率區間是通帶還是阻帶。
其他的FIR濾波器設計方法在另一章中介紹。