FIR滤波器窗函数设计.doc

FIR滤波器窗函数设计

课题名称：FIR滤波器窗函数设计

FIR滤波器窗函数设计
引言：
数字滤波器(Digital Filter)是指输入、输出都是离散时间信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。在许多数字信号处理系统中，如图像信号处理等，有限冲激响应（FIR）滤波器是最常用的组件之一，它完成信号预调、频带选择和滤波等功能。FIR滤波器虽然在截止频率的边沿陡峭性能上不及无限冲激响应（IIR）滤波器，但是却具有严格的线性相位特性，稳定性好，能设计成多通带（或多阻带）滤波器组，所以能够在数字信号处理领域得到广泛的应用。
一、 数字滤波器的分类
（1） 根据系统响应函数的时间特性分为两类
1. FIR（Finite Impulse Response）数字滤波器网络
特点：不存在反馈支路，其单位冲激响应为有限长。
2. IIR（Infinite Impulse Response）数字滤波器网络
特点：存在反馈支路，即信号流图中存在环路，其单位冲激响应为无限长。
（2） FIR数字滤波器IIR数字滤波器的区别
1. 从性能上来说，IIR滤波器传递函数包括零点和极点两组可调因素，对极点的惟一限制是在单位圆内。因此可用较低的阶数获得高的选择性，所用的存储单元少，计算量小，效率高。但是这个高效率是以相位的非线性为代价的。选择性越好，则相位非线性越严重。FIR滤波器传递函数的极点固定在原点，是不能动的，它只能靠改变零点位置来改变它的性能。所以要达到高的选择性，必须用较高的阶数；对于同样的滤波器设计指标，FIR滤波器所要求的阶数可能比IIR滤波器高5-10倍，但是 FIR滤波器可以得到严格的线性相位。
2. 从结构上看，IIR滤波器必须采用递归结构，极点位置必须在单位圆内，否则系统将不稳定。相反，FIR滤波器只要采用非递归结构，不论在理论上还是在实际的有限精度运算中都不存在稳定性问题，因此造成的频率特性误差也较小。此外FIR滤波器可以采用快速傅里叶变换算法，在相同阶数的条件下，运算速度可以快得多。
3. 从设计工具看，IIR滤波器可以借助于模拟滤波器的成果，因此一般都有有效的封闭形式的设计公式可供准确计算，计算工作量比较小，对计算工具的要求不高。
从上面的简单比较可以看到IIR与FIR滤波器各有所长，所以在实际应用时应该从多方面考虑来加以选择。从使用要求上来看，在对相位要求不敏感的场合，如语言通信等，选用IIR较为合适，这样可以充分发挥其经济高效的特点；对于图像信号处理，数据传输等以波形携带信息的系统，则对线性相位要求较高。如果有条件，采用FIR滤波器较好。另外，不论IIR和FIR，阶数越高，信号延迟越大；同时在IIR滤波器中，阶数越高，系数的精度要求越高，否则很容易造成有限字长的误差使极点移到单位园外。因此在阶数选择上是综合考虑的。
二、 FIR数字滤波器的设计
(一) FIR滤波器原理
FIR滤波器的系统输入输出差分方程为：
所以FIR滤波器的系统函数为：

由于FIR滤波器的单位脉冲响应h(n)是一个有限长序列，H (z)是的(N−1)次多项式，它在Z平面上有(N−1)个零点，同时在原点有(N−1)阶重极点。因此，H(z)永远稳定。FIR滤波器设计的任务是选择有限长度的h(n)，使传输函数满足一定的幅度特性和线性相位要求。由于FIR滤波器很容易实现严格的线性相位，所以FIR 数字滤波器设计的核心思想是求出有限的脉冲响应来逼近给定的频率响应。
(二) FIR滤波器的窗函数设计原理
窗函数法的设计思想是按照所要求的理想滤波器频率响应，设计一个FIR滤波器，使之频率响应来逼近。先由的傅里叶反变换导出理想滤波器的冲激响应序列，即：

由于是矩形频率特性，所以是一无限长的序列，且是非因果的，而要计的FIR滤波器的冲激响应序列是有限长的，所以要用有限长的序列h(n)来逼近无限长的序列，最有效的方法是截断，或者说用一个有限长度的窗口函数w(n)序列来截取，即: 。
按照复卷积公式，在时域中的乘积关系可表示成在频域中的周期性卷积关系，即可得所设计的FIR滤波器的频率响应:

其中，为截断窗函数的频率特性。由此可见，实际的FIR数字滤波器的频率响应逼近理想滤波器频率响应的好坏，完全取决于窗函数的频率特性。
如果w(n)具有下列形式：

w(n)相当于一个矩形，我们称之为矩形窗。即我们可采用矩形窗函数w(n)将无限脉冲响应截取一段来近似为。经过加矩形窗后所得的滤波器实际频率响应能否很好地逼近理想频率响应呢？下图给出了理想滤波器加矩形窗后的情况。理想低通滤波器的频率响应如图中左上角图，矩形窗的频率响应为左下角图。根据卷积定理，即得实际滤波器的频率响应图形为图中右图。

由图可看出，加矩形窗后使实际频率响应偏离理想频率响应，主要影响有三个方面：
（1）理想幅频特性陡直边缘处形成过渡带，过渡带宽取决于矩形窗函数频率响应的主瓣宽度。
（2）过渡带两侧形成肩峰和波纹，这是矩形窗函数频率响应的旁瓣引起的，旁瓣相对值越大，旁瓣越多，波纹越多。
（3）随窗函数宽度N的增大，矩形窗函数频率响应的主瓣宽度减小，但不改变旁瓣的相对值。
为了改善滤波器的性能，需使窗函数谱满足：主瓣尽可能窄，以使设计出来的滤波器有较陡的过渡带；第一副瓣面积相对主瓣面积尽可能小，即能量尽可能集中在主瓣，外泄少，使设计出来的滤波器的肩峰和余振小逼近于理想滤波器。但是这两个条件是相互矛盾的，实际应用中，折衷处理，兼顾各项指标。
(三) MATLAB信号处理中提供的窗函数
上边只考虑了矩形窗，如果我们使窗的主瓣宽度尽可能地窄，旁瓣尽可能地小，可以获得性能更好的滤波器，通过改变窗的形状来达到这个目的。在数字信号处理的发展过程中形成了不同于矩形窗的很多窗函数，这些窗函数在主瓣和旁瓣特性方面各有特点，可满足不同的要求。为此，用窗函数法设计FIR数字滤波器时，要根据给定的滤波器性能指标选择窗口宽度N和窗函数w(n)。下面具体介绍几类类窗函数及其特性。
1. 矩形窗
矩形窗函数的时域形式可以表示为：

它的频域特性为：

2. 汉宁窗函数
汉宁窗函数的时域形式可以表示为：

它的频域特性为：

其中，为矩形窗函数的幅度频率特性函数。
汉宁窗函数的最大旁瓣值比主瓣值低31dB，但是主瓣宽度比矩形窗函数的主瓣宽度增加了1倍，为8π/N。
3. 海明窗函数
海明窗函数的时域形式可以表示为:

它的频域特性为:

其中，为矩形窗函数的幅度频率特性函数。
海明窗函数的最大旁瓣值比主瓣值低41dB，但它和汉宁窗函数的主瓣宽度是一样大的。
4. 三角窗函数
三角窗是最简单的频谱函数为非负的一种窗函数。三角窗函数的时域形式可以表示为：
当n为奇数时：

当n为偶数时：

它的频域特性为：

三角窗函数的主瓣宽度为8π/N，比矩形窗函数的主瓣宽度增加了一倍，但是它的旁瓣宽度却小得多。
(四) 各窗函数的图形及幅频特性
1. 用MATLAB编程绘制各种窗函数的形状，窗函数的长度为21。程序如下：
clf;Nwin=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); %海明窗
stext='海明窗';
case 4
w=triang(Nwin); %三角窗
stext='三角窗';
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 on;
end
程序运行结果如下图：

2. 用MATLAB编程，采用512个频率点绘制各窗函数的幅频特性。程序如下：
clf;Nf=512; %窗函数复数频率特性的数据点数
Nwin=20; %窗函数数据长度
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
[y,f]=freqz(w,1,Nf); %求解窗函数的幅频特性
mag=abs(y); %求得窗函数幅频特性
posplot=['2,2,' int2str(ii)];
subplot(posplot);
plot(f/pi,20* log10(mag/max(mag))); %绘制窗函数的幅频特性
xlabel('归一化频率');ylabel('振幅/dB');
title(stext);grid on;
end
程序运行结果如下图：

由上图各窗函数的幅频特性可以看到各种窗函数都具有明显的主瓣和旁瓣。主瓣频宽和旁瓣的幅值衰减特性决定了窗函数的应用场合。矩形窗具有最窄的主瓣，但也有最大的旁瓣峰值（第一旁瓣衰减为13dB）；不同窗函数在这两方面的特点是不同的，因此应根据具体的问题进行选择。通常来讲，海明窗和汉宁窗的主瓣，具有较小的旁瓣和较大的衰减速度，是较为常用的窗函数。下表列出了各种窗函数主瓣和旁瓣的特征：
各窗函数的特征表
窗函数
主瓣频宽
第一旁瓣相对主瓣衰减（dB）
矩形窗

-13
汉宁窗

-31
海明窗

-41
三角窗

-25
此外，主旁瓣频率宽度还与窗函数长度N有关。增加窗函数长度N将减小窗函数的主瓣宽度，但不能减小旁瓣幅值衰减的相对值（分贝数），这个值是由窗函数决定的。例如：绘制矩形窗函数的幅频响应，窗长度分别为：(1)N=10;(2)N=20; (3)N=50;(4)N=100时的图形如下：

由上图可以看出，随着N的增大，主瓣和旁瓣都变窄，但第一旁瓣相对主瓣的幅值下降分贝数相同，第二旁瓣相对第一旁瓣幅值下降的分贝数也相同。然而，随着N的增大，旁瓣数也增多，减少主瓣宽度和抑制旁瓣是一对矛盾，不可兼得，只能根据不同用途折衷处理。
三、 运用窗函数设计数字滤波器
用于信号分析中的窗函数可根据滤波器的指标进行不同选择。用于滤波器的窗函数，一般要求窗函数的主瓣宽度窄，以获得较好的过渡带；旁瓣相对值尽可能少，增加通带的平稳度和增大阻带的衰减。
窗函数设计法的基本原理是用一定宽度窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列，主要设计步骤为：
（1）通过傅里叶逆变换获得理想滤波器的单位脉冲响应；
（2）由性能指标（阻带衰减的分贝数）根据窗函数的特征表的值确定满足阻带衰减的窗函数类型，和窗口长度N ；
（2）根据，求实际滤波器的单位脉冲响应；
（3）检验滤波器的性能。
下面根据上面步骤设计FIR滤波器：
一、例1要求用窗函数设计一个线性相位FIR低通滤波器，并满足性能指标：通带边界的归一化频率=0.5,阻带边界的归一化频率=0.66，阻带衰减不小于30dB，通带波纹不大于3dB。假设一个信号，其中f1=5Hz,f2=20Hz。信号的采样频率为50Hz。试将原信号与通过滤波器的信号进行比较。
根据要求：阻带衰减不小于30dB, 根据窗函数的特征表选取汉宁窗满足滤波要求；在窗函数设计法中，要求设计的频率归一化到区间内，Nyquist频率对应于，因此通带和阻带边界频率为和。
程序如下：
wp=0.5*pi;ws=0.66*pi; %滤波器边界频率
wdelta=ws-wp; %过渡带宽
N=ceil(8*pi/wdelta) %根据过渡带宽得滤波器所用窗函数的最小长度
Nw=N;
wc=(wp+ws)/2; %理想低通滤波器截止频率在通带和阻带边界频率的中点
n=0:N-1;
alpha=(N-1)/2; %求滤波器的相位延迟
m=n-alpha+eps; %eps为MATLAB系统的精度
hd=sin(wc*m)./(pi*m); %求理想滤波器脉冲响应
win=hanning(Nw); %采用汉宁窗
h=hd.*win'; %在时间域乘积对应于频率域的卷积
b=h;
figure(1)
[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;
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); %可采用此函数给出滤波器的输出
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('输出信号')
程序运行结果如下图：

图一：低通滤波器幅频和相频响应

图二：输入输出信号
结果分析：该设计的要求通带边界=0.5,阻带边界频率=0.66，对应于50Hz的采样频率通带边界频率为=/2*Fnormal=50/2*0.5=12.5Hz， =50/2*0.66=16.5Hz，其中为采样频率，Fnormal为归一化频率。由图一中滤波器幅频响应可以看到，在小于12.5Hz的频段上，几乎看不到下降，即满足通带波纹不大于3dB的要求。在大于16.5Hz的频段上，阻带衰减大于30dB，满足设计要求。由图一中相频响应可见，在通带范围内，相位频率为一条直线，表明该滤波器为线性相位。图二给出了滤波器的输入信号和输出信号，输入信号包括5Hz和20Hz的信号，由图一可知，20Hz的信号不能通过该滤波器，通过滤波器后只剩下5Hz的信号，输出结果也证明了这一点。但是由于FIR滤波器所需的阶数较高，信号延迟(N-1)/2也较大，输出信号前面有一段直线就是延迟造成的。上述程序显示的N取50才能满足设计要求。这样相位延迟为(N-1)/2*1/Fs=0.49s，可以看到输出信号前面一段直线的距离大约为0.49s。验证了FIR滤波器相位延迟的理论。
二、标准型FIR滤波器
上面的设计中运用理想脉冲响应与窗函数乘积的方法给出了滤波器传递函数的设计方法。其实MATLAB已将上述方法复合成一个函数，提供基于上述原理设计标准型FIR滤波器的工具函数。fir1就是采用经典窗函数法设计线性相位FIR数字滤波器的函数，且具有标准低通、带通、高通、带阻等类型。函数调用格式为：
b=fir1(n,wn[,'ftype',window])
式中，n为FIR滤波器的阶数，对于高通,带阻滤波器，n需取偶数；wn为滤波器截止频率，范围为0~1（归一化频率）。对于带通,带阻滤波器，wn=[w1,w2](w1<w2)； ‘ftype'为滤波器的类型：缺省时为低通或带通滤波器；'high'为高通滤波器；‘stop'为带阻滤波器，'DC-1'为第一频带为通带的多带滤波器；'DC-0'为第一频带为阻带的多带滤波器。window为窗函数列向量，其长度为n+1。缺省时，自动取海明窗。MATLAB提供的窗函数有boxcar、hanning、hamming、bartlett、blackman、kaiser、chebwin。b为FIR滤波器系数向量，长度为n+1。FIR滤波器的传递函数具有下列形式：

用函数fir1设计的FIR滤波器的群延迟为n/2。考虑到n阶滤波器系数个数为N，即n+1，这里的延迟与前面所讲的(N-1)/2的延迟一致。注意这里的滤波器的最小阶数比窗函数的长度少1。
运用MATLAB工具函数设计上面FIR滤波器的程序如下：
wp=0.5*pi;ws=0.66*pi; %滤波器的边界频率
wdelta=ws-wp; %过渡带宽度
N=ceil(8* pi/wdelta); %求解滤波器的最小阶数，根据汉宁窗求主瓣宽
Wn=(0.5+0.66)*pi/2; %截止频率取通带和阻带边界频率的中点
b=fir1(N,Wn/pi,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;
程序运行结果与图一中的滤波器幅频和相频响应是一致的。
例二，设计一个48阶FIR带通滤波器，通带边界的归一化频率为0.35和0.65。假设一个信号，其中含有f1=1Hz，f2=10Hz,f3=20Hz，三种频率成分信号的采样频率为50Hz。试将原信号与通过滤波器的信号进行比较。
编写程序如下：
wp=[0.35 0.65];N=48; %通带边界频率（归一化频率）和滤波器阶数
Fs=50; %采样频率50HZ
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');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*sin(2*pi*f3*t);%输入信号
%y=filter(b,1,x); %可以采用过滤器进行滤波
y=fftfilt(b,x); %采用fftfilt对输入信号滤波
figure(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')
图三：带通滤波器幅频和相频响应

图四：输入输出信号
结果分析：程序运行结果见图三和图四。通带归一化频率0.35和0.65对应于采样频率为50Hz的通带范围为8.75Hz和16.25Hz。由图三可见满足这一设计要求。在这个频带范围内的相位满足线性相位，符合FIR滤波器的一般特点。图四为检测滤波器的输入信号和输出信号。输入信号中含有1Hz,10Hz和20Hz的信号。根据图三上图可知，1Hz和20Hz的频率在阻带范围内，不能通过该滤波器，只有10Hz的信号可以通过该滤波器，输出信号表明了这一点。滤波器的相位延迟根据N/2*0.02s=0.48s得到，输出信号前面的直线部分大体为这个时间延迟,另外滤波后周期为10Hz的信号相位（红线开始部分），跟滤波前的相位（输入信号开始部分）也一致，说明通过该滤波器滤波后没有改变信号的相位。
四、 总结与体会
由于FIR滤波器就有较好的线性相位特性，使其在图像信号处理，数据传输等以波形携带信息的系统领域中有着广泛的应用。此外，FIR滤波器的冲激响应为有限长序列，其系统函数为一多项式，它所含的极点均为原点，所以FIR滤波器永远是稳定的。本次课题主要学习了FIR数字滤波器的设计原理、窗函数在FIR滤波器设计中的应用并用MATLAB实现了FIR低通带通滤波器的设计。通过这次学习，我不但掌握了FIR数字滤波器窗函数的基本知识及其实际应用的技巧，还提高了自己MATLAB软件应用方面的能力，收获很多。