FIR滤波器设计Word下载.docx
《FIR滤波器设计Word下载.docx》由会员分享，可在线阅读，更多相关《FIR滤波器设计Word下载.docx（64页珍藏版）》请在冰豆网上搜索。
具有光滑、正弦过渡带的低通滤波器设计
Fircos
7.1.2FIR数字滤波器滤波函数
相对于IIR滤波器的滤波函数，FIR数字滤波器滤波函数除了dimpulse和dstep仅适用于IIR滤波器外，其他各种函数可直接应用于FIR滤波器，只是输入的分母多项式向量a=1。
另外，MATLAB还提供了一个函数fftfilt,该函数利用效率高的基于FFT算法实现对数据的滤波，该函数只适用于FIR滤波器，调用形式为：
y=fftfilt（b,x[,n]）
式中，b为FIR滤波器的系数向量；
x为输入数据；
n为FFT长度，缺省时，函数选用最佳的FFT长度，y为滤波器的输出。
该函数执行下面的操作：
n=length（x）;
y=ifft（fft（x）.*fft（b,n）./fft（a,n））;
应注意，y=fftfilt（b,x）等价于y=filter（b,a,x）。
7.2FIR滤波器的窗函数设计
7.2.1窗函数的基本原理
FIR滤波器设计的主要任务是根据给定的性能指标确定滤波器的系数b，即系统单位脉冲序列h（n），它是一个有限长序列。
FIR滤波器的理想频率响应，可写成复数形式的Fourier级数形式：
（7-4）
式中，hd（n）是对应的单位脉冲响应序列。
这说明滤波器的频率响应和单位脉冲响应互为Fourier变换对。
因此其单位脉冲响应可由下式求得，
（7-5）
求得序列
后，通过z变换，可得到
（7-6）
注意，这里
为无限长序列，因此
是物理上不可实现的。
如何变成物理上可实现呢？
一个自然的想法是只取其中的某些项，即只截取
中的一部分，比如n=0,…,N-1，N为正整数。
这种处理相当于将
n=-∞~∞与函数w（n）相乘，w（n）具有下列形式：
w（n）相当于一个矩形，我们称之为矩形窗。
即我们可采用矩形窗函数w（n）将无限脉冲响应
截取一段h（n）来近似为
，这种截取在数学上表示为：
h（n）=
w（n）（7-7）
这里应该强调的是，加窗函数不是可有可无的，而是将设计变为物理可实现所必须的。
截取之后的滤波器传递函数变为：
（7-8）
式中，N为窗口宽度，H（z）是物理可实现系统。
为了获得线性相位，FIR滤波器h（n）必须满足中心对称条件（即7-3式），序列h（n）的延迟为
。
这种方法的基本原理是用一定宽度的矩形窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列，从而得到FIR滤波器的脉冲响应，故称为FIR滤波器的窗函数设计法。
经过加矩形窗后所得的滤波器实际频率响应能否很好地逼近理想频率响应呢？
图7-1示意给出了理想滤波器加矩形窗后的情况。
理想低通滤波器的频率响应如图中左上角图，矩形窗的频率响应为左下角图。
时间域内的乘积（7-7）式要求实际频率响应为这两个频率响应函数在频域内的卷积（卷积定理），即得到图形为图7-1（右图）。
图7-1FIR滤波器理想与实际频率响应
由图可看出，加矩形窗后使实际频率响应偏离理想频率响应，主要影响有三个方面：
（1）理想幅频特性陡直边缘处形成过渡带，过渡带宽取决于矩形窗函数频率响应的主瓣宽度。
（2）过渡带两侧形成肩峰和波纹，这是矩形窗函数频率响应的旁瓣引起的，旁瓣相对值越大，旁瓣越多，波纹越多。
（3）随窗函数宽度N的增大，矩形窗函数频率响应的主瓣宽度减小，但不改变旁瓣的相对值。
为了改善FIR滤波器性能，要求窗函数的主瓣宽度尽可能窄，以获得较窄的过渡带；
旁瓣相对值尽可能小，数量尽可能少，以获得通带波纹小,阻带衰减大，在通带和阻带内均平稳的特点，这样可使滤波器实际频率响应更好地逼近理想频率响应。
这里我们明确两个概念：
截断和频谱泄漏。
信号是无限长的，而在进行信号处理时只能采取有限长信号，所以需要将信号“截断”。
在信号处理中，“截断”被看成是用一个有限长的“窗口”看无限长的信号，或者从分析的角度是无限长的信号x（t）乘以有限长的窗函数w（t）。
由傅立叶变换性质可知,时间域内的乘积对应于频率域的卷积,即
（7-9）
这里，x（t）是频宽有限信号，而w（t）是频宽无限信号，
表示互为Fourier变换对。
截断后的信号也必须是频宽无限信号，这样就是有限频带的信号分散到无限频带中去，这样就产生了所谓频谱泄漏。
从能量的角度来看，频谱泄漏也是能量的泄漏，因为加窗后使原来信号集中的窄频带内的能量分散到无限的频带宽度范围内。
频谱泄漏是不可避免的，但要尽量减小。
上边只考虑了矩形窗，如果我们使窗的主瓣宽度尽可能地窄，旁瓣尽可能地小，可以获得性能更好的滤波器，能否改变窗的形状而达到这个目的呢？
回答是肯定的。
其实数字信号处理的前驱者们设计了不同于矩形窗的很多窗函数，这些窗函数在主瓣和旁瓣特性方面各有特点，可满足不同的要求。
为此，用窗函数法设计FIR数字滤波器时，要根据给定的滤波器性能指标选择窗口宽度N和窗函数w（n）。
下面我们介绍窗函数。
7.2.2MATLAB信号处理中提供的窗函数
（1）矩形窗：
前面分析中所用的矩形窗可用下面函数来实现w=boxcar（N），N为窗的长度（以下函数与此同），w为返回的窗函数序列。
（2）汉宁窗：
w=hanning（N）
汉宁窗的表达式为：
（7-10）
（3）哈明窗：
w=hamming（N）
哈明窗的表达式为：
（7-11）
（4）Bartlett窗：
w=bartlett（N）
Bartlett窗的表达式为：
当N为奇数时，
（7-12）
当N为偶数时，
（7-13）
（5）Blackman窗：
w=blackman（N）
Blackman窗的表达式为：
（7-14）
Blackman窗比其他相同尺寸窗（哈明窗,汉宁窗）具有主瓣较宽和旁瓣泄漏较小的特点。
（6）三角窗：
w=triang（N）
三角窗的表达式为：
（7-15）
（7-16）
三角窗和Bartlett窗十分类似。
三角窗的两端值不为零，而Bartlett窗则为零，这一点可从例7-1中看出。
（7）Kaiser窗：
w=kaiser（n,beta）
其中，beta是Kaiser窗参数,影响窗旁瓣幅值的衰减率。
Kaiser窗表达式：
（7-17）
式中，I0[.]是修正过的零阶Bessel函数。
Kaiser窗用于滤波器设计时，若旁瓣幅值为
，则
（7-18）
（8）Chebyshev窗：
w=chebwin（n,r）
式中,r是窗口的旁瓣幅值在主瓣以下的分贝数。
Chebyshev窗具有主瓣宽度最小，而旁瓣等高,高度可调整的特点。
下面我们在MATLAB观看各种窗函数的形状和频率域图象来验证上述所讲特点。
【例7-1】用MATLAB编程绘制各种窗函数的形状。
窗函数的长度为21。
%Samp7_1
clf
Nwin=21;
n=0:
Nwin-1;
%数据总数和序列序号
figure
（1）
forii=1:
4
switchii
case1
w=boxcar（Nwin）;
%矩形窗
stext='
矩形窗'
;
case2
w=hanning（Nwin）;
%汉宁窗
汉宁窗'
case3
w=hamming（Nwin）;
%哈明窗
哈明窗'
case4
w=bartlett（Nwin）;
%Bartlett窗
Bartelett窗'
end
posplot=['
2,2,'
int2str（ii）];
%指定绘制窗函数的图形位置
subplot（posplot）;
stem（n,w）;
%绘出窗函数
holdon
plot（n,w,'
r'
）;
%绘制包络线
xlabel（'
n'
ylabel（'
w（n）'
title（stext）;
holdoff;
gridon
end
figure
（2）
w=blackman（Nwin）;
%Blackman窗
Blackman窗'
w=triang（Nwin）;
%三角窗
三角窗'
w=kaiser（Nwin,4）;
%Kaiser窗
Kaiser窗（Beta=4）'
w=chebwin（Nwin,40）;
%Chebyshev窗
Chebyshev窗（r=40）'
plot（n,w,'
%绘出包络线
ylabel（'
title（stext）;
gridon;
程序运行结果见图7-2。
可以看到各种窗函数的形状。
图7-2各种窗函数的时间域形状
【例7-2】用MATLAB编程，采用512个频率点绘制各种窗函数的幅频特性。
%Samp7_2
clf;
Nf=512;
%窗函数复数频率特性的数据点数
Nwin=20;
%窗函数数据长度
w=hamming（Nwin）;
%Bartlett窗
[y,f]=freqz（w,1,Nf）;
%求解窗函数的幅频特性，窗函数相当于一个数字滤波器
mag=abs（y）;
%求得窗函数幅频特性
plot（f/pi,20*log10（mag/max（mag）））;
%绘制窗函数的幅频特性
归一化频率'
振幅/dB'
%以Nf点数求解窗函数的幅频响应
%求得窗函数幅频响应
plot（f/pi,20*log10（mag/max（mag）））;
%绘制幅频响应
程序运行结果见图7-3。
可以看到各种窗函数的幅频形状。
对照该图可知这些窗函数具有上面所分析的窗函数的特征。
图7-3各种窗函数的幅频形状
由图7-3可见，各种窗函数都具有明显的主瓣（Mainlobe）和旁瓣（Sidelobe）。
主瓣频宽和旁瓣的幅值衰减特性决定了窗函数的应用场合。
矩形窗具有最窄的主瓣，但也有最大的旁瓣峰值（第一旁瓣衰减为13dB）；
Blackman窗具有最大的旁瓣衰减，但也具有最宽的主瓣宽度。
不同窗函数在这两方面的特点是不同的，因此应根据具体的问题进行选择。
通常来讲，哈明窗和汉宁窗的主瓣具有较小的旁瓣和较大的衰减速度，是较为常用的窗函数。
表7-2总结了各种窗函数主瓣和旁瓣的特征（理论分析可参考其他的数字信号处理教材），大家可对照窗函数的幅频形状（图7-3）认真理解体会。
表7-2各种窗函数的特点
窗函数
主瓣宽
第一旁瓣相对主瓣衰减（分贝）
矩形窗
-13
汉宁窗
-31
哈明窗
-41
Bartlett窗
-25
Blackman窗
-57
三角窗
Kaiser窗
可调整
Chebyshev窗
主旁瓣频率宽度还与窗函数长度N有关。
增加窗函数长度N将减小窗函数的主瓣宽度，但不能减小旁瓣幅值衰减的相对值（分贝数），这个值是由窗函数决定的。
这个特点可由下面的例子清楚地看出。
【例7-3】绘制矩形窗函数的幅频响应，窗长度分别为：
（1）N=10;
（2）N=20;
（3）N=50;
（4）N=100.
%Samp7_3
Nwin=10;
Nwin=20;
Nwin=50;
case4
Nwin=100;
%用不同的窗长度求得复数频率特性
%求得幅频特性
%指定绘图位置
plot（f/pi,20*log10（mag/max（mag）））;
%绘出幅频形状
stext=['
N='
int2str（Nwin）];
%给出标题，指出所用的数据个数
图7-4数据长度不同的矩形窗的幅频形状
程序运行结果见图7-4。
显然，随着N的增大，主瓣和旁瓣都变窄，但第一旁瓣相对主瓣的幅值下降分贝数相同，第二旁瓣相对第一旁瓣幅值下降的分贝数也相同。
这样，随着N的增大，旁瓣也得到抑制，有力地减少了频谱泄漏，但不能完全消除。
减少主瓣宽度和抑制旁瓣是一对矛盾，不可兼得，只能根据不同用途折衷处理。
7.2.3运用窗函数设计数字滤波器
用于信号分析中的窗函数可根据用户的不同要求选择。
用于滤波器的窗函数，一般要求窗函数的主瓣宽度窄，以获得较好的过渡带；
旁瓣相对值尽可能少，增加通带的平稳度和增大阻带的衰减。
基于窗函数的FIR数字滤波器设计的算法十分简单，其主要步骤为：
（1）对滤波器理想频域幅值响应进行傅立叶逆变换获得理想滤波器的单位脉冲响应hd（n）。
一般假定理想低通滤波器的截止频率为
，其幅频特性满足
（7-19）
则根据傅立叶逆变换，单位脉冲响应为：
（7-20）
其中，
为信号延迟。
（2）由性能指标（阻带衰减的分贝数）根据表7-2第3列的值确定满足阻带衰减的窗函数类型w（n）。
滤波器的阶数越高，滤波器的幅频特性越好，但数据处理的费用也越高，因此像IIR滤波器一样，FIR滤波器也要确定满足性能指标的滤波器最小阶数。
由前面的讨论（图7-1）可知，滤波器的主瓣宽度相当于过渡带宽，因此，使过渡带宽近似于窗函数主瓣宽（表7-2中的第二列）可求得满足性能指标的窗口长度N，此时，信号延迟
为（N-1）/2。
（3）求实际滤波器的单位脉冲响应h（n）：
根据h（n）=hd（n）*w（n）。
（4）检验滤波器的性能。
可设定一些信号采用7.1.2节指出的函数或6.3.2节所给的函数进行滤波。
下面采用实例说明如何根据上面步骤设计FIR滤波器。
【例7-4】用窗函数设计一个线性相位FIR低通滤波器，并满足性能指标：
通带边界的归一化频率wp=0.5,阻带边界的归一化频率ws=0.66，阻带衰减不小于30dB，通带波纹不大于3dB。
假设一个信号，其中f1=5Hz,f2=20Hz。
信号的采样频率为50Hz。
试将原信号与通过滤波器的信号进行比较。
由题意，阻带衰减不小于30dB,根据表7-2,选取汉宁窗，因为汉宁窗的第一旁瓣相对主瓣衰减为31dB，满足滤波要求。
在窗函数设计法中，要求设计的频率归一化到0~
区间内，Nyquist频率对应于
，因此通带和阻带边界频率为0.5
和0.66
程序如下
%Samp7_4
wp=0.5*pi;
ws=0.66*pi;
%滤波器边界频率
wdelta=ws-wp;
%过渡带宽
N=ceil（8*pi/wdelta）%根据过渡带宽等于表7-2中汉宁窗函数主瓣宽求得滤波器所用窗函数的最小长度
Nw=N;
wc=（wp+ws）/2;
%截止频率在通带和阻带边界频率的中点
N-1;
alpha=（N-1）/2;
%求滤波器的相位延迟
m=n-alpha+eps;
%eps为MATLAB系统的精度
hd=sin（wc*m）./（pi*m）;
%由（7-20）式求理想滤波器脉冲响应
win=hanning（Nw）;
%采用汉宁窗
h=hd.*win'
%在时间域乘积对应于频率域的卷积
b=h;
[H,f]=freqz（b,1,512,50）;
%采用50Hz的采样频率绘出该滤波器的幅频和相频响应
subplot（2,1,1）,plot（f,20*log10（abs（H）））
xlabel（'
频率/Hz'
subplot（2,1,2）,plot（f,180/pi*unwrap（angle（H）））
相位/^o'
%impz（b,1）;
%可采用此函数给出滤波器的脉冲响应
%zplane（b,1）;
%可采用此语句给出滤波器的零极点图
%grpdelay（b,1）;
%可采用此函数给出滤波器的群延迟
f1=3;
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）;
%给出滤波器的输出
subplot（2,1,1）,plot（t,x）,title（'
输入信号'
）%绘输入信号
subplot（2,1,2）,plot（t,y）%绘输出信号
holdon;
plot（[11]*（N-1）/2*dt,ylim,'
）%绘出延迟到的时刻
时间/s'
）,title（'
输出信号'
）
程序运行结果见图7-5和图7-6。
该例设计通带边界wp=0.5,阻带边界频率ws=0.66，对应于50Hz的采样频率通带边界频率为fp=Fs/2*Fnormal=50/2*0.5=12.5Hz,fs=50/2*0.66=16.5Hz,其中Fs为采样频率，Fnormal为归一化频率。
由图7-5上图可以看到，在小于12.5Hz的频段上，几乎看不到下降，即满足通带波纹不大于3dB的要求。
在大于16.5Hz的频段上，阻带衰减大于30dB，满足设计要求。
由图7-5下图可见，在通带范围内，相位频率为一条直线，表明该滤波器为线性相位。
图7-6给出了滤波器的输入信号和输出信号，输入信号包括3Hz和20Hz的信号，由图7-5可知，20Hz的信号不能通过该滤波器，通过滤波器后只剩下3Hz的信号，输出结果也证明了这一点。
但要注意，由于FIR滤波器所需的阶数较高，信号延迟（N-1）/2也较大，输出信号前面有一段直线就是延迟造成的。
上述程序显示的N取50才能满足设计要求。
这样相位延迟为（N-1）/2*1/Fs=0.49s，可以看到输出信号前面一段直线的距离大约为0.49s。
验证了FIR滤波器相位延迟的理论。
在输出信号的前部，有一些小信号，这是截断信号边界所致，后面的部分就没有了这种信号。
若采用零相位的filtfilt函数（说明见第六章第三节）输出，则可最大限度地减小边界的影响。
图7-5例7-4所设计滤波器的幅频响应（上图）和相频响应（下图）
图7-6例7-4所设计滤波器的输入和输出信号
升级会员 