MATLAB—FIR数字滤波器设计

  本文是基于MATLAB的数字滤波器设计，所有数据基于计算机内部处理，因而都是离散信号，所以采用的是数字滤波器。从实现的网络结构或者单位脉冲响应来看，数字滤波器可以分为FIR(Finite Impulse Response , 有限脉冲响应)滤波器和IIR(Infinite Impulse Response , 无线脉冲响应)滤波器。本文主要从两个方面介绍FIR滤波器的MATLAB设计：集成函数和工具箱。

1 FIR滤波器的原理

  从时域看，FIR的一般表达式H(z)如下：

${\sum }_{n=0}^{N-1}h(n)z^{-n}=h(0)+h(1)z^{-1}+...+h(N-1)z^{-(N-1)}$

  从表达式至少可以看出两点：系统只在原点处有极点，可以用抽头法构建模型，每一个乘法器的系数即为抽头系数。
由于N是有限值，因此抽头系数项是有限的，可以得到h(n)：　

${\sum }_{n=0}^{N-1}h(0)\delta (n)+h(1)\delta (n-1)+...+h(N-1)\delta (n-(N-1))$

  可以看出h(n)，即系统的脉冲响应的项数是有限的，因此成为有限脉冲响应。

2 FIR滤波器的特点

  FIR一个突出的优点就是具有严格的线性相位特性，但并不是所有结构的FIR都具备此特性，只有当FIR滤波器的单位脉冲响应满足对称条件时，FIR才具有线性相位特性。

2.1 相位特性

  研究FIR的相位特性时，将结构分为奇对称和偶对称两种。

2.1.1 偶对称

当FIR单位脉冲响应满足偶对称，即有：

$h(n)=h(M-n),0\le n\le M$

代入H(z)的表达式：

${\sum }_{n=0}^{M}h(n)z^{-n}={\sum }_{n=0}^{M}h(M-n)z^{-n}$

令M-n=k , 对等号右边部分做变换：

$H(z)={\sum }_{k=0}^{M}h(k)z^{-(M-k)}=z^{-M}{\sum }_{n=0}^{M}h(k)z^k=z^{-M}H(z^{-1})$

所以H(z)有两种表达方式，那么H(z)可以表示为：

$H(z)=\frac{1}{2}[H(z)+z^{-M}H(z)]=\frac{1}{2}{\sum }_{k=0}^{M}h(n)[z^{-n}+z^{-M}{z}^n]$

提取一项$z^{-\frac{M}{2}}$，可得到滤波器频响：

$H(e^{jw})=e^{-jw\frac{M}{2}}{\sum }_{n=0}^{M}h(n)cos[w(\frac{M}{2}-n)]=A(w)e^{-jw\frac{M}{2}}$

根据相频定义，可知FIR相频为：

$\phi (w)=-\frac{M}{2}w$

故满足线性相位特性。

2.1.2 奇对称

当FIR单位脉冲响应满足偶对称，即有：

$h(n)=-h(n-M),0\le n\le M$

由2.1.1推导思路可得：

$H(e^{jw})=e^{-jw\frac{M}{2}+\frac{\pi }{2}}{\sum }_{n=0}^{M}h(n)sin[w(\frac{M}{2}-n)]=A(w)e^{-jw\frac{M}{2}+\frac{\pi }{2}}$

根据相频定义，可知FIR相频为：

$\phi (w)=-\frac{M}{2}w+\frac{\pi }{2}$

故满足线性相位特性。

2.2 幅度特性

  研究FIR的幅度特性时，将结构分别分为偶数和奇数的偶对称和奇对称四种。推导方式和相位特性类似，这里直接给出结论。

单位脉冲响应特征	相位特性	幅度特性	滤波器种类
偶对称，偶整数	线性相位	对于w=0，π ，2π为偶对称	适合各种滤波器
偶对称，奇整数	线性相位	对于w=π为奇对称，对于w=0、2π为偶对称，w=π处为0	不适合高通，带阻
奇对称，偶整数	线性相位，附加90°相移	对于w=0、π、2π均为奇对称，在w=0、π、2π处都为0	只适合带通
偶对称，奇整数	线性相位，附加90°相移	对于w=0、2π均为奇对称，在w=π处为偶对称，在w=0、2π处为0	适合高通、带通

3 几种滤波器函数

3.1 fir1()

• fir1函数语法形式

b=fir1(n,wn,'ftype',window,'noscale');
% n     滤波器阶数
% wn    类型和意义与ftype有关，wn的取值范围为(0,1)，1代表fs的1/2。
% 当wn为单值，表示截止频率。若ftype为'low'，表示低通；若ftype为'high'，表示高通。
% 当wn为[wn1 wn2]，即两个元素组成的向量，则表示带通或带阻，ftype对应'bandpass'或'stop'。
% 当wn由多个数组成的向量[w1 w2...wn]，其中w1<w2<…<wn，则fir1返回带0<ω<w1，w1<ω<w2，…，wn<ω<1的n阶多带滤波器。
% ftype 滤波器类型
% 'low'       -低通
% 'high'      -高通
% 'bandpass'  -带通
% 'stop'      -带阻
% 'DC-0'      -第一带为带阻
% 'DC-1'	  -第一带为带通
% window  默认为海明窗(Hamming),还有汉宁窗(Hanning),布拉克曼窗(Blackman)和凯塞窗(Kaiser)
% noscale 指定归一化滤波器幅度

3.2 fir2()

  fir2在fir1的基础上，还可以指定频段内的理想幅值，即任意响应滤波器。

• fir2函数语法形式

b=fir2(n,fm,m,npt,lap,window);
% n        滤波器阶数
% f和m     f是一个频段向量，在(0,1)内分布,对应归一化频率;m是对应频段理想幅值
% npt      正整数，指定对幅频进行插值的点数，默认为512
% lap      插值时，将非连续点转化成连续点的点数，默认为25
% window   窗函数，默认为海明窗

3.3 kaiserord()

  凯塞窗，根据一些期望的参数，得到设计滤波器所需要的参数，可以用于设计最优滤波器。

• kaiserord函数语法形式

[n,wn,beta,ftype]=kaiserord(f,a,dev,fs);
% 输入
% f   偶数个元素的向量，表示过渡带的起点到终点，过渡带的增益。
% eg:[1200 1300 1500 1600]表示过渡带1200~1300Hz,1500~1600Hz。fs为采样频率。如果没有fs，f元素取值在(0,1)，即自动归一化;如果有fs，f就为实际频率。
% a   f确定了过渡带，a确定各频段的理想幅值。通带设置为1，阻带设置为0。
% 假设f=[0 0.3 0.4 0.6 0.7 1.0] , a=[0 1 0 0 0.5 0.5]
% dev 向量，各频段纹波。
% 输出
% n     返回滤波器阶数
% wn    向量，返回滤波器截止频率点
% beta  返回凯塞窗的beta值
% ftype 返回滤波器类型 

3.4 firpm()

  既能设计任意响应，又能设计最优滤波器，还能附加90°相移

• firpm函数语法形式

[b,delta]=firpm(n,f,a,w,'ftype');
% n,f,a与kaiserord函数一样
% w,对应频段幅值的权值,权值越高，越接近理想状态
% ftype 滤波器结构。默认为偶对称脉冲响应，'hilbert'为奇对称结构，即具有90°相移特性
% delta 最大纹波

4 仿真实例

  需求：利用凯塞窗设计一个低通FIR滤波器，过渡带为1000~1500Hz，采样频率为8000Hz，通带纹波最大为0.01，阻带纹波最大为0.05。利用海明窗及firpm函数设计相同的低通滤波器，截止频率为1500Hz，滤波器阶数为凯赛窗函数求取的值。绘出三种方法设计的幅度频率响应曲线。

4.1 代码

%% 低通；过渡带为1~1.5KHz，采样频率8KHz，通带最大纹波0.01，阻带最大纹波0.05

% 获取凯塞窗参数,后续fir1要设计凯赛窗的低通滤波器
fs=8000;                                   % 采样频率
fc=[1000 1500];                            % 过渡带；1000~1500过渡带.通带:0~(1000*2/fs),阻带:(1500*2/fs)~1
mag=[1 0];                                 % 理想幅值.0~1000*2/fs是1，1500*2/fs~1是0
dev=[0.01 0.05];                           % 纹波
[n,wn,beta,ftype]=kaiserord(fc,mag,dev,fs);% 获取凯塞窗参数,根据过渡带幅值返回ftype类型

% fir1设计凯塞窗和海明窗滤波器
h_kaiser=fir1(n,wn,ftype,kaiser(n+1,beta));
h_hamm=fir1(n,fc(2)*2/fs);                 %fc(2)=1500Hz,LPF截止频率

% 设计最优滤波器
fpm=[0 fc(1)*2/fs fc(2)*2/fs 1];           % firpm频段向量,归一化[0 1]
magpm=[1 1 0 0];% firpm幅值向量
h_pm=firpm(n,fpm,magpm);                   %设计低通滤波器

% 幅频
m_kaiser=20*log(abs(fft(h_kaiser,1024)))/log(10);
m_hamm=20*log(abs(fft(h_hamm,1024)))/log(10);
m_pm=20*log(abs(fft(h_pm,1024)));

% 设置幅频响应的横坐标为Hz
x_f=[0:1:length(m_kaiser)/2]*fs/length(m_kaiser);

% 只显示正频率部分
m1=m_kaiser(1:length(x_f));
m2=m_hamm(1:length(x_f));
m3=m_pm(1:length(x_f));

% 绘制幅频曲线
plot(x_f,m1,'-',x_f,m2,'-.',x_f,m3,'--');
xlabel('频率(Hz)');
ylabel('幅度(dB)');
legend('凯塞窗','海明窗','最优滤波器');
grid;

4.2 仿真分析

  如图所示，绘制了三种滤波器。如果要用凯塞窗设计fir1，那么先要利用kaiserord函数获取凯塞窗的相关参数，即第一种滤波器。如果设计海明窗，则输入截止频率即可(要归一化后的频率，否则需要输入采样频率)。设计最优滤波器输入频段和相应幅度即可，这里需要和前两种做比较，因此都是输入n阶。