MATLAB常用的FIR滤波器设计方法之窗函数法

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滤波器设计方法在另一章中介绍。