IIR数字带通滤波器设计
《数字信号处理课程设计报告》
题 目: IIR 数字带通滤波器设计
学 院:
专 业:
班 级:
姓 名:
指导教师:
2012年 6月24日
目 录
1数字滤波器设计原理…………………………………………………3
1.1数字滤波器简介 ………………………………………………3
1.2 IIR滤波器的设计原理 ………………………………………3 2 IIR数字滤波器设计方法……………………………………………4
2.1用脉冲相应不变法设计IIR 数字滤波器 ……………………4
2.2用双线性变换法设计IIR 数字滤波器 ………………………7 3 IIR数字带通滤波器设计过程………………………………………9
3.1设计步骤 ………………………………………………………9
3.2程序流程框图…………………………………………………10
3.3 MATLAB 程序……………………………………………………11 4运行结果及分析 ……………………………………………………12 5总结 …………………………………………………………………13 6参考书目 ……………………………………………………………14
基于MATLAB 的IIR 数字带通滤波器设计
一、数字滤波器设计原理
1.1 数字滤波器简介
数字滤波器是一种用来过滤时间离散信号的数字系统,通过对抽样数据进行数学处理来达到频域滤波的目的。可以设计系统的频率响应,让它满足一定的要求,从而对通过该系统的信号的某些特定的频率成分进行过滤,这就是滤波器的基本原理。如果系统是一个连续系统,则滤波器称为模拟滤波器。如果系统是一个离散系统,则滤波器称为数字滤波器。
信号 通过线性系统后,其输出 就是输入信号 和系统冲激响应 的卷积。除了 外, 的波形将不同于输入波形 。从频域分析来看,信号通过线性系统后,输出信号的频谱将是输入信号的频谱与系统传递函数的乘积。除非 为常数,否则输出信号的频谱将不同于输入信号的频谱,某些频率成分 较大的模,因此, 中这些频率成分将得到加强,而另外一些频率成分 的模很小甚至为零, 中这部分频率分量将被削弱或消失。因此,系统的作用相当于对输入信号的频谱进行加权。
1.2 IIR 滤波器的设计原理
IIR 数字滤波器的设计一般是利用目前已经很成熟的模拟滤波器的设计方法来进行设计,通常采用模拟滤波器原型有butterworth 函数、chebyshev 函数、bessel 函数、椭圆滤波器函数等。
IIR 数字滤波器的设计步骤:
(1) 按照一定规则把给定的滤波器技术指标转换为模拟低通滤波器的技术指标;
(2) 根据模拟滤波器技术指标设计为响应的模拟低通滤波器;
(3) 很据脉冲响应不变法和双线性不变法把模拟滤波器转换为数字滤波器;
(4) 如果要设计的滤波器是高通、带通或带阻滤波器,则首先把它们的技术指标转化为模拟低通滤波器的技术指标,设计为数字低通滤波器,最后通过频率转换的方法来得到所要的滤波器。
二、IIR 数字滤波器设计方法
IIR 数字滤波器是一种离散时间系统,其系统函数为
假设M ≤N ,当M >N 时, 系统函数可以看作一个IIR 的子系统和一个(M-N)的FIR 子系统的级联。IIR 数字滤波器的设计实际上是求解滤波器的系数和 ,它是数学上的一种逼近问题,即在规定意义上(通常采用最小均方误差准则)去逼近系统的特性。如果在S 平面上去逼近,就得到模拟滤波器;如果在z 平面上去逼近,就得到数字滤波器。
1. 用脉冲相应不变法设计IIR 数字滤波器
利用模拟滤波器来设计数字滤波器,也就是使数字滤波器能模仿模拟滤波器的特性,这种模仿可以从不同的角度出发。脉冲响应不变法是从滤波器的脉冲响应出发,使数字滤波器的单位脉冲响应序列h (n ) 模仿模拟滤波器的冲激响应h a (t ) ,即将h a (t ) 进行等间隔采样,使h (n ) 正好等于h a (t ) 的采样值,满足
h (n )=h a (nT )
式中, T 是采样周期。
如果令H a (s ) 是h a(t ) 的拉普拉斯变换,H (z ) 为h (n ) 的Z 变换,利用采样序列的
Z 变换与模拟信号的拉普拉斯变换的关系得
X (
z ) z =e sT 1∞1∞2π⎫⎛=∑X a (s -jk Ωs ) =∑X a s -j k ⎪T k =-∞T k =-∞⎝T ⎭(1-1)
则可看出,脉冲响应不变法将模拟滤波器的S 平面变换成数字滤波器的Z 平面,这个从s 到z 的变换z =esT 是从S 平面变换到Z 平面的标准变换关系式。
图1-1脉冲响应不变法的映射关系
由(1-1)式,数字滤波器的频率响应和模拟滤波器的频率响应间的关系为
1∞⎛ω-2πk ⎫j ω(1-2) H (e ) =∑H a j ⎪T k =-∞⎝T ⎭
这就是说,数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。正如采样定理所讨论的,只有当模拟滤波器的频率响应是限带的,且带限于折叠频率以内时,即
Ωs (1-3) H a (j Ω) =0T 2
才能使数字滤波器的频率响应在折叠频率以内重现模拟滤波器的频率响应,|Ω|≥= π而不产生混叠失真,即
1⎛ω⎫|ω|
产生周期延拓分量的频谱交叠,即产生频率响应的混叠失真,如图7-4所示。这时数字滤波器的频响就不同于原模拟滤波器的频响,而带有一定的失真。当模拟滤波器的频率响应在折叠频率以上处衰减越大、越快时,变换后频率响应混叠失真就越小。这时,采用脉冲响应不变法设计的数字滤波器才能得到良好的效果。
对某一模拟滤波器的单位冲激响应h a (t ) 进行采样,采样频率为f s ,若使f s
增加,即令采样时间间隔(T =1/f s )减小,则系统频率响应各周期延拓分量之间
相距更远,因而可减小频率响应的混叠效应。
脉冲响应不变法优缺点:
从以上讨论可以看出,脉冲响应不变法使得数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应,也就是时域逼近良好,而且模拟频率Ω和数字频率ω之间呈线性关系ω=ΩT 。因而,一个线性相位的模拟滤波器(例如贝塞尔滤波器)通过脉冲响应不变法得到的仍然是一个线性相位的数字滤波器。
脉冲响应不变法的最大缺点是有频率响应的混叠效应。所以,脉冲响应不变法只适用于限带的模拟滤波器(例如,衰减特性很好的低通或带通滤波器) ,而且高频衰减越快,混叠效应越小。至于高通和带阻滤波器,由于它们在高频部分不衰减,因此将完全混淆在低频响应中。如果要对高通和带阻滤波器采用脉冲响应不变法,就必须先对高通和带阻滤波器加一保护滤波器,滤掉高于折叠频率以上的频率,然后再使用脉冲响应不变法转换为数字滤波器。当然这样会进一步增加设计复杂性和滤波器的阶数。
2. 用双线性变换法设计IIR 数字滤波器
脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S 平面到Z平面是多值的映射关系所造成的。为了克服这一缺点,可以采用非线性频率压缩方法,将整个频率轴上的频率范围压缩到-π/T ~π/T 之间,再用z =esT 转换到Z 平面上。也就是说,第一步先将整个S 平面压缩映射到S 1平面的-π/T ~
π/T 一条横带里;第二步再通过标准变换关系z =es 1T 将此横带变换到整个Z 平面上去。这样就使S 平面与Z 平面建立了一一对应的单值关系,消除了多值变换性,也就消除了频谱混叠现象,映射关系如图1-3所示。
S 平面S 1平面Z 平面
图1-3双线性变换的映射关系
为了将S 平面的整个虚轴j Ω压缩到S1平面j Ω1轴上的-π/T 到π/T 段上,可以通过以下的正切变换实现
2⎛ΩT ⎫Ω=tan 1⎪T ⎝2⎭式中, T 仍是采样间隔。 (1-5)
当Ω1由-π/T 经过0变化到π/T 时,Ω由-∞经过0变化到+∞,也即映射了整个j Ω轴。将式(1-5)写成
2e j Ω1T /2-e j Ω1T /2j Ω=⋅j Ω1T /2-j Ω1T /2T e +e
将此关系解析延拓到整个S 平面和S1平面,令j Ω=s ,j Ω1=s 1,则得
-s T 2e s 1T /2-e -s 1T /22⎛s 1T ⎫21-e 1s =⋅s 1T /2-s 1T /2=tanh ⎪=⋅-s T T e +e T ⎝2⎭T 1+e 1
再将S1平面通过以下标准变换关系映射到Z 平面
z =es 1T
从而得到S 平面和Z 平面的单值映射关系为:
21-z -1(1-6) s =-1T 1+z T 21+s +s (1-7) =z =T 21-s -s 2T
式(1-6)与式(1-7)是S 平面与Z 平面之间的单值映射关系,这种变换都是两个线性函数之比,因此称为双线性变换
式(1-5)与式(1-6)的双线性变换符合映射变换应满足的两点要求。 首先, 把z =ej ω,可得
21-e -j ω2⎛ω⎫s ==j tan ⎪=j Ω-j ωT 1+e T ⎝2⎭即S 平面的虚轴映射到Z 平面的单位圆。
其次,将s =σ+jΩ代入式(1-8),得
2 +σ+j Ωz = -σ-j ΩT 因此 2⎛2⎫2 +σ⎪+Ω ⎝T ⎭|z |=2 2⎫2 -σ⎪+Ω⎝T ⎭
由此看出,当σ0时,|z |>1。也就是说,S 平面的左半平面映射到Z 平面的单位圆内,S 平面的右半平面映射到Z 平面的单位圆外,S 平面的虚轴映射到Z 平面的单位圆上。因此,稳定的模拟滤波器经双线性变换(1-8)
后所得的数字滤波器也一定是稳定的。
双线性变换法优缺点
双线性变换法与脉冲响应不变法相比,其主要的优点是避免了频率响应的混叠现象。这是因为S 平面与Z 平面是单值的一一对应关系。S 平面整个j Ω轴单值地对应于Z 平面单位圆一周,即频率轴是单值变换关系。这个关系如式(1-8)所示,重写如下:
2⎛ω⎫tan ⎪T
⎝2⎭上式表明,S 平面上Ω与Z 平面的ω成非线性的正切关系,如图7-7所示。 Ω=
由图7-7看出,在零频率附近,模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系;但当Ω进一步增加时,ω增长得越来越慢,最后当Ω→∞时,ω终止在折叠频率ω=π处,因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象,从而消除了频率混叠现象。
图1-4双线性变换法的频率变换关系
但是双线性变换的这个特点是靠频率的严重非线性关系而得到的,如式(1-8)及图1-4所示。由于这种频率之间的非线性变换关系,就产生了新的问题。首先,一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器,不再保持原有的线性相位了;其次,这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的,即某一频率段的幅频响应近似等于某一常数(这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性),不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变,如图1-5所示。
图1-5双线性变换法幅度和相位特性的非线性映射
对于分段常数的滤波器,双线性变换后,仍得到幅频特性为分段常数的滤波器,但是各个分段边缘的临界频率点产生了畸变,这种频率的畸变,可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变,然后经变换后正好映射到所需要的数字频率上。
三、IIR 数字带通滤波器设计过程:
根据以上IIR 数字滤波器设计方法,下面运用双线性变换法基于MATLAB 设计一个IIR 带通滤波器,其中带通的中心频率为ωp0=0.5π,; 通带截止频率ωp1=0.4π, ωp2=0.6π; 通带最大衰减αp=3dB;阻带最小衰减αs=15dB;阻带截止频率ωs2=0.7π
1. 设计步骤:
(1)根据任务, 确定性能指标:在设计带通滤波器之前, 首先根据工程实际的需要
确定滤波器的技术指标:
带通滤波器的阻带边界频率关于中心频率ωp0几何对称,因此ws1=wp0- (ws2-wp0)=0.3π
通带截止频率wc1=0.4π,wc2=0.6π; 阻带截止频率wr1=0.3π,wr2=0.7π;阻带最小衰减αs=3dB和通带最大衰减αp=15dB;
(2)用Ω=2/T*tan(w/2)对带通数字滤波器H(z)的数字边界频率预畸变,得到带
通模拟滤波器H(s)的边界频率主要是通带截止频率ωp1, ωp2; 阻带截止频率ωs1,ωs2的转换。
为了计算简便,对双线性变换法一般T=2s
通带截止频率wc1=(2/T)*tan(wp1/2)=tan(0.4π/2)=0.7265
wc2=(2/T)*tan(wp2/2)=tan(0.6π/2)=1.3764
阻带截止频率wr1=(2/T)*tan(ws1/2)=tan(0.3π/2)=0.5095
wr2=(2/T)*tan(ws2/2)=tan(0.7π/2)=1.9626
阻带最小衰减αs=3dB和通带最大衰减αp=15dB;
(3)运用低通到带通频率变换公式λ=(((Ω^2)-(Ω0^2))/(B*Ω)) 将模拟带通滤
波器指标转换为模拟低通滤波器指标。
B=wc2-wc1=0.6499
normwr1=(((wr1^2)-(w0^2))/(B*wr1))=2.236
normwr2=(((wr2^2)-(w0^2))/(B*wr2))=2.236
normwc1=(((wc1^2)-(w0^2))/(B*wc1))=1
normwc2=(((wc2^2)-(w0^2))/(B*wc2))=1
得出,normwc=1,normwr=2.236
模拟低通滤波器指标:normwc=1,normwr=2.236,αp=3dB,αs=15dB
(4)设计模拟低通原型滤波器。用模拟低通滤波器设计方法得到模拟低通滤波器的传输函数Ha(s);借助巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer)滤波器、贝塞尔(Bessel)滤波器等。
(5)调用lp2bp 函数将模拟低通滤波器转化为模拟带通滤波器。
(6)利用双线性变换法将模拟带通滤波器Ha(s)转换成数字带通滤波器H(z).
2. 程序流程框图:
↓
↓
↓
↓
↓
↓
↓
3.MATLAB 程序:
MATLAB 程序如下:
clear
wp0=0.5*pi;wp1=0.4*pi;wp2=0.6*pi; Ap=3;ws2=0.7*pi;As=15;T=2; %数字带通滤波器技术指标
ws1=wp0-(ws2-wp0); %计算带通滤波器的阻带下截止频率
wc1=(2/T)*tan(wp1/2);wc2=(2/T)*tan(wp2/2);
wr1=(2/T)*tan(ws1/2);wr2=(2/T)*tan(ws2/2);
w0=(2/T)*tan(wp0/2); %频率预畸变
B=wc2-wc1; %带通滤波器的通带宽度
normwr1=(((wr1^2)-(w0^2))/(B*wr1));
normwr2=(((wr2^2)-(w0^2))/(B*wr2));
normwc1=(((wc1^2)-(w0^2))/(B*wc1));
normwc2=(((wc2^2)-(w0^2))/(B*wc2)); %带通到低通的频率变换 if abs(normwr1)>abs(normwr2)
normwr=abs(normwr2)
else normwr=abs(normwr1)
end
normwc=1; %将指标转换成归一化模拟低通滤波器的指标
N=buttord(normwc,normwr,Ap,As,'s'); %设计归一化的模拟低通滤波器阶数N 和3db 截止频率
[bLP,aLP]=butter(N,normwc,'s'); %计算相应的模拟滤波器系统函数G(p)
[bBP,aBP]=lp2bp(bLP,aLP,w0,B); %模拟域频率变换,将G(P)变换成模拟带通滤波器H(s)
[b,a]=bilinear(bBP,aBP,0.5); %用双线性变换法将H(s)转换成数字带
通滤波器H(z)
w=linspace (0,2*pi,500);
h=freqz(b,a,w);
subplot(2,1,2);
plot(w,abs(h));
grid on
xlabel('w(rad)')
ylabel('|H(jw)|')
title('频谱函数')
subplot(2,2,1);
plot(w,20*log10(abs(h)));
axis([0,2*pi,-120,20]);
grid on
xlabel('w(rad)')
ylabel('20*lg|H(jw)|(db)')
title('20*lg|H(jw)|--w')
四.运行结果及分析:
图
程序运行结果:normwr=2.2361
由设计流程计算得normwr=2.236与运行结果相同。
低通原型的每一个边界频率都映射为带通滤波器两个相应的边界频率。根据通带截至频率和阻带截至频率与频谱函数曲线比较,满足设计要求。
五.总结:
通过这个实验,对设计带通数字滤波器的整个过程有了很好的掌握。其中双线性变换法,巴特沃斯设计模拟滤波器的运用,也比较熟悉了。
通过对数字带通滤波器的设计,熟悉了MATLAB 的运行环境,初步掌握了MATLAB 语言在数字信号处理中一些基本库函数的调用和编写基本程序等应用;熟悉了滤波器设计的一般原理,对滤波器有了一个感性的认识;学会了数字高通滤波器设计的一般步骤;加深了对滤波器设计中产生误差的原因以及双线性变换法优缺点的理解和认识。总之,使理论联系了实际,巩固并深化了对课本基本知识的认识和理解,使理论得以升华。
参考书目:
1. 程佩青《数字信号处理教程》北京清华大学出版社2007年2月.
2. 赵知劲、刘顺兰《数字信号处理实验》. 浙江大学出版社.
3. S.K.Mitra.Digital Signal Processing:A Computer-Based Approach. NewYork,NewYork:McGraw-Hill,thirded,2006.
4. 肖伟、刘忠等《 MATLAB程序设计与应用》清华大学出版社、北京交通大学出版社.
5. 胡良剑、孙晓君 《 MATLAB数学实验》. 高等教育出版社.
6. 刘兴钊《数字信号处理 电子工业出版社 2010年
7. 黄大伟 《数字滤波器》 中国铁道出版社 1991年
8. 薛山 《MATLAB 基础教程》清华大学出版社 2012年
《计算机控制技术》课程设计指导教师评语