赞
踩
目录
数字滤波技术是数字信号处理的一个重要组成部分,滤波器的设计是信号处理的核心问题之一。根据FIR滤波器的原理,提出了FIR滤波器的窗函数设计法,并对常用的几种窗函数进行了比较。给出了在MATLAB环境下,用窗函数法设计FIR滤波器的过程和设计实例。仿真结果表明,设计的FIR滤波器的各项性能指标均达到了指定要求,设计过程简便易行。该方法为快速、高效地设计FIR滤波器提供了一个可靠而有效的途径。
随着信息时代的到来,数字信号处理已经成为当今一门极其重要的学科和技术,并且在通信、语音、图像、自动控制等众多领域得到了广泛的应用。在数字信号处理中,数字滤波器占有极其重要的地位,它具有精度高、可靠性好、灵活性大等特点。现代数字滤波器可以用软件或硬件两种方式来实现。软件方式实现的优点是可以通过滤波器参数的改变去调整滤波器的性能。
FIR滤波器实质上是一个分节的延迟线,把每一节的输出加权累加,便得到滤波器的输出。对于FIR滤波器,幅度上只需满足以下两个条件之一,就能构成线性相位FIR滤波器。根据冲激响应的时域特性,数字滤波器可分为无限长冲激响应滤波器(IIR)和有限长冲激响应滤波器(FIR)。FIR的突出优点是:系统总是稳定的、易于实现线性相位、允许设计多通带(或多阻带)滤波器,但与IIR相比,在满足同样阻带衰减的情况下需要的阶数较高。滤波器的阶数越高,占用的运算时间越多,因此在满足指标要求的情况下应尽量减少滤波器的阶数。
FIR滤波器的基本结构可以理解为一个分节的延时线,把每一节的输出加权累加,可得到滤波器的输出。FIR滤波器的冲激响应h(n)是有限长的,数学上M阶FIR滤波器可以表示为:
FIR滤波器的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数的问题,设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。
窗函数设计法是一种通过截短和计权的方法使无限长非因果序列成为有限长脉冲响应序列的设计方法。通常在设计滤波器之前,应该先根据具体的工程应用确定滤波器的技术指标。在大多数实际应用中,数字滤波器常常被用来实现选频操作,所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。
用窗函数法设计FIR滤波器的步骤如下:
(1)根据过渡带宽及阻带衰减要求,选择窗函数的类型并估计窗口长度N(或阶数M=N-1)。窗函数类型可根据最小阻带衰减AS独立选择,因为窗口长度N对最小阻带衰减AS没有影响。在确定窗函数类型以后,可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N。设待求滤波器的过渡带宽为△ω,它与窗口长度N近似成反比。窗函数类型确定后,其计算公式也确定了,不过这些公式是近似的,得出的窗口长度还要在计算中逐步修正。原则是在保证阻带衰减满足要求的情况下,尽量选择较小的N。在N和窗函数类型确定后,即可调用MATLAB中的窗函数求出窗函数wd(n)。
(2)根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)。如果给出待求滤波器的频率响应为Hd(ejω),则理想的单位脉冲响应可以用下面的傅里叶反
变换式求出:
在一般情况下,hd(n)是不能用封闭公式表示的,需要采用数值方法表示。从ω=0到ω=2π采样N点,采用离散傅里叶反变换(IDFT)即可求出。
(3)计算滤波器的单位脉冲响应h(n)。它是理想单位脉冲响应和窗函数的乘积,即h(n)=hd(n)•wd(n),在MATLAB中用点乘命令表示为h=hd•wd。
(4)验算技术指标是否满足要求。为了计算数字滤波器在频域中的特性,可调用freqz子程序,如果不满足要求,可根据具体情况,调整窗函数类型或长度,直到满足要求为止。
使用窗函数法设计时要满足以下两个条件:
(1)窗谱主瓣尽可能地窄,以获得较陡的过渡带;
(2)尽量减少窗谱的最大旁瓣的相对幅度,也就是使能量尽量集中于主瓣,减小峰肩和纹波,进而增加阻带的衰减。
根据工程经验,给定的滤波器指标参数一般为通带截止频率ωp、阻带截止频率ωs、实际通带波动Rp和最小阻带衰减As。窗函数设计的经验公式为:
在实际工程中常用的窗函数有五种,即矩形窗、三角窗、汉宁窗、海明窗和凯泽窗。这些窗函数在MATLAB中分别用boxcar、triang、hanning、hamming、kaiser实现,它们之间的性能比较如表1所示。
窗类型 旁瓣峰值 主瓣峰值 最小阻带衰减
矩形窗 13dB 4π/M 21dB
三角窗 25dB 8π/M 25dB
汉宁窗 31dB 8π/M 44dB
海明窗 41dB 8π/M 53dB
凯泽窗 57dB 12π/M 74dB
用窗函数设计高通滤波器,性能指标如下:通带截止频率ωs=0.2π,阻带截止频率ωp=0.3π,实际通带波动Rp=0.25dB,最小阻带衰减As=70dB。
- As=70;
- ωs=0.2*π;
- ωp=0.3*π
- tr_width=ωp-ωs; %计算过渡带宽
- M=ceil((As-7.95)*2*π/(14.36*tr_width)+1)+1; 按凯泽窗计算滤波器长度
- disp([’滤波器的长度为’,num2str(M)]);
- beta=0.1102*(As-8.7); %计算凯泽窗的β值
- n=[0:1:M-1];
- disp([’线性相位斜率为’,num2str(beta)]);
- w_kai=(kaiser(M,beta))’; %求凯泽窗函数
- ωc=(ωs+ωp)/2;
- hd=ideal_lp(π,M)-ideal_lp(ωc,M); %求理想脉冲响应
- h=hd*w_kai; %设计的脉冲响应为理想脉冲响应与窗函数乘积
- [db,mag,pha,grd,ω]=freqz_m(h,[1]);
- delta_ω=2*π/1000;
- Rp=-(min(db(ωp/delta_ω+1:1:501)));
- disp([’实际通带波动为’,num2str(Rp)]); %以下为作图程序
- As=-round(max(db(1:1:ωs/delta_ω+1)));
- disp([’最小阻带衰减为’,num2str(As)]);
- subplot(1,1,1);
- subplot(2,2,1);
- stem(n,hd);
- title(’理想脉冲响应’);
- axis([0 M-1 -0.4 0.8]);
- ylabel(’hd(n)’);
- subplot(2,2,2);
- stem(n,w_kai);
- title(’凯泽窗’);
- axis([0 M-1 0 1.1]);
- ylabel(’wd(n)’);
- subplot(2,2,3);
- stem(n,h);
- title(’实际脉冲响应’);
- axis([0 M-1 -0.4 0.8]);
- xlabel(’n’);ylabel(’h(n)’);
- subplot(2,2,4);
- plot(ω/π,db);
- title(’幅度响应/dB’);
- axis([0 1 -100 10]);
- grid;
- xlabel(’以π为单位的频率’);
- ylabel(’分贝数/dB’);
用窗函数设计低通滤波器,性能指标如下:通带截止频率ωp=0.1π,阻带截止频率ωs=0.25π,实际通带波动Rp=0.10dB,最小阻带衰减As=40dB。
分析:从表1可以看出,汉宁窗、海明窗和凯泽窗能提供大于40dB的最小阻带衰减。但汉宁窗的旁瓣峰值较小,而主瓣宽度和海明窗一样。可以使滤波器的阶数较少,所以选用汉宁窗进行设计,程序主要部分如下:
- ωp=0.10*π;
- ωs=0.25*π;
- tr_width=ωs-ωp; %计算过渡带宽
- M=ceil(6.6*/tr_width)+1; %按汉宁窗计算滤波器长度
- disp([’滤波器的长度为’,num2str(M)]);
- n=0:M-1;
- ωc=(ωs+ωp)/2; %截止频率取为两边缘频率的平均值
- hd=ideal_lp(ωc,M); %求理想脉冲响应
- w_han=(hanning(M))’; %求汉宁窗函数
- h=hd*w_han; %设计的脉冲响应为理想脉冲响应与窗函数乘积
- [db,mag,pha,grd,ω]=freqz_m(h,[1]);%以下为作图语句
- delta_ω=2*π/1000;
- Rp=-(min(db(1:1: ωp/delta_ω+1)));
- disp([’实际通带波动为’,num2str(Rp)]); %以下为作图程序
- As=-round(max(db(ωs/delta_ω+1:1:501)));
- disp([’最小阻带衰减为’,num2str(As)]);
- subplot(221)
- stem(n,hd);
- title(’理想冲击响应’),
- axis([0 M-1 -0.1 0.3]);
- ylabel(’hd(n)’);
- subplot(222)
- stem(n,w_han);
- title(’汉宁窗’),
- axis([0 M-1 0 1.1]);
- ylabel(’wd(n)’);
- subplot(223)
- stem(n,h);
- title(’实际冲击响应’),
-
- axis([0 M-1 -0.1 0.3]);
- xlabel(’n’);
- ylabel(′h(n)′);
- subplot(224);
- plot(ω/π,db);
- title(′幅度响应(db)′);
- axis([0 1 -100 10]),
- grid;
- xlabel(′以π为单位的频率′);
- ylabel(′分贝数′);
1.IIR数字滤波器的系统函数可以写成封闭函数的形式。
2.IIR数字滤波器采用递归型结构,即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成,可以组合成直接型、正准型、级联型、并联型四种结构形式,都具有反馈回路。由于运算中的舍入处理,使误差不断累积,有时会产生微弱的寄生振荡。
3.IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果,如巴特沃斯、契比雪夫和椭圆滤波器等,有现成的设计数据或图表可查,其设计工作量比较小,对计算工具的要求不高。在设计一个IIR数字滤波器时,我们根据指标先写出模拟滤波器的公式,然后通过一定的变换,将模拟滤波器的公式转换成数字滤波器的公式。
4.IIR数字滤波器的相位特性不好控制,对相位要求较高时,需加相位校准网络。
IIR 数字滤波器可用一个n阶差分方程
不难看出,数字滤波器与模拟滤波器的设计思路相仿,其设计实质也是寻找一组系数{b,a},去逼近所要求的频率响应,使其在性能上满足预定的技术要求;不同的是模拟滤波器的设计是在S平面上用数学逼近法去寻找近似的所需特性H(S),而数字滤波器则是在Z平面寻找合适的H(z)。
- 通过合适类型模拟滤波器进行数字带通设计程序
- fp= [480,520];fs=[450,550] %模拟通带、阻带频率
- wp=[480,520]*pi*2 %模拟通带角频率
- ws=[450,550]*pi*2 %模拟阻带角频率
- rp=3;rs=20 %通带波动、阻带衰减
- 巴特沃斯型模拟带通滤波器设计
- [n,wn]=buttord (wp,ws,rp,rs,'s')
- [b,a]=butter(n,wn,'s') %模拟带通滤波函数系数
- 巴特沃斯型模拟带通滤波器频率响应
- [ha,w]= freqs(b,a)
- ma=abs(ha);pha=unwrap(angle(ha))
- subplot(421);plot(w/(2*pi),ma) %模拟幅频曲线
- subplot(423);plot(w/(2 pi),pha) %模拟相频曲线
- 冲击响应不变法进行离散化设计
- fo=5000 %采样频率
- [bn,an]=impinvar(b,a,5000) %数字带通滤波函数系数
- 巴特沃斯型数字带通滤波器频率响应
- [hz,w]=freqz(bn,an)
- mz=abs(hz);phz=unwrap(angle(hz))
- subplot(422);plot(w,mz) %数字滤波器幅频曲线
- subplot(424);plot(w,phz) %数字滤波器相频曲线
- hi=impz(bn,an) %数字滤波器冲击响应
- subplot(425),plot(hi) %冲击响应曲线
- n=0:300;t=n/fo
- xl=2*square(2*pi*500*t) %500Hz方波信号
- subplot(426);plot(x1) %500Hz方波波形
- yi=conv(hi,x1) %时域卷积输出
- subplot(427);plot(yi) %卷积输出波形
- y1=filter(bn,an,x1) %数字滤波函数输出
- subplot(428);plot(y1) %数字滤波器输出波形
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。