用脉冲响应不变法和双线性变换法设计IIR数字滤波器分析

吉林建筑大学
电气与电子信息工程学院
数字信号处理课程设计报告
设计题目：IIR数字滤波器的设计
专业班级：信工111
学生姓名：**
学号：********
指导教师：高晓红王超
设计时间：2014.01.06－2014.01.10
目录
一、设计目的 (1)
二、设计内容 (1)
三、设计原理 (1)
3.1 数字低通滤波器的设计原理 (1)
3.2 变换方法的原理 (2)
四、设计步骤 (8)
五、数字低通滤波器 MATLAB 编程及幅频特性曲线 (10)
5.1 MATLAB语言编程 (10)
5.2 幅频特性曲线 (12)
六、总结 (13)
七、参考文献 (13)
一、设计目的
课程设计是理论学习的延伸，是掌握所学知识的一种重要手段，对于贯彻理论联系实际、提高学习质量、塑造自身能力等于有特殊作用。

本次课程设计一方面通过MATLAB仿真设计内容，使我们加深对理论知识的理解，同时增强其逻辑思维能力，另一方面对课堂所学理论知识作一个总结和补充。

二、设计内容
已知通带截止频率f
p =0.2kHz,通带最大衰减α
p
=1dB,阻带截止频率
f s =0.3kHz，阻带最小衰减α
s
=25dB，T=1ms，按照以上技术要求，用脉冲响应不
变法和双线性变换法设计巴特沃斯数字低通滤波器，并观察所设计数字滤波器的幅频特性曲线。

三、设计原理
3.1 数字低通滤波器的设计原理
滤波器是自动控制、信号处理和通信领域的重要组成部分,广泛地应用于各种系统中。

数字滤波器是指输入、输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。

与模拟滤波器相比具有很多突出的优点，例如它可以满足滤波器对幅度和相位特性的严格要求，可以避免模拟滤波器所无法克服的电压漂移和噪声问题。

设计数字滤波器，首先要按一定规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标。

根据转换后的技术指标设计模拟低通滤波器G(s)，再按一定的规则将G(s)转换成H(z)。

若是高通、带通或带阻数字滤波器则将它们的技术指标先转化为低通模拟滤波器的技术指标，然后设计出低通G(s)，再将G(s)转换为所需的H(z)将系统函数Ha(s)从s平面转换到z平面的方法有多种，但工程上常用的是脉冲响应不变法和双线性变换法。

这两种方法有着各自的优缺点，在设计巴特沃斯数字低通滤波器时需要根据相应要求或想要达成的效果从中选择。

3.2 变换方法的原理
脉冲响应不变法原理：
Ha(s) 拉氏逆变换 ha(t) 等间隔采样 ha(nT)=h(n) Z 变换 H(z)
脉冲响应不变法是一种将模拟滤波器转化为数字滤波器的基本方法。

它利用模拟滤波器理论设计数字滤波器，也就是使数字滤波器能模仿模拟滤波器的特性，这种模仿可从不同的角度出发。

脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列h(n) 模仿模拟滤波器的冲击响应ha(t), 使h(n)正好等于ha(t)的采样值，即
)()(nT h n h a =
T 为采样周期。

利用模拟滤波器理论设计数字滤波器，也就是使数字滤波器能模仿模拟滤波器的特性，这种模仿可从不同的角度出发。

脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列h(n) 模仿模拟滤波器的冲击响应ha(t), 使h(n)正好等于ha(t)的采样值，即
)()(nT h n h a =
T 为采样周期。

如以Ha(s)及H （z ）分别表示ha(t)的拉氏变换及h(n)的z 变换，即
)]([)(t h L s H a a =
)]([)(n h Z z H =
则根据采样序列z 变换与模拟信号拉氏变换的关系，得:
上式表明，采用脉冲响应不变法将模拟滤波器变换为数字滤波器时，它所完成的S 平面到Z 平面的变换，正是以前讨论的拉氏变换到Z 变换的标准变换关系，即首先对Ha(s)作周期延拓，然后再经过z=e st 的映射关系映射到Z 平面上。

z=e st 的映射关系表明，S 平面上每一条宽为2π/T 的横带部分，都将重叠地映射到Z 平面的整个全部平面上。

每一横带的左半部分映射到Z 平面单位圆以内，每一横带的右半部分映射到Z 平面单位圆以外，jΩ轴映射在单位圆上，但jΩ轴上的每一段2π/T 都对应于绕单位圆一周，如图1所示。

图1 脉冲响应不变法的映射关系
应当指出，z=e st 的映射关系反映的是Ha(s)的周期延拓与H （z ）的关系，而不是Ha(s)本身与H （z ）的关系，因此，使用脉冲响应不变法时，从Ha(s)到H(z)并没有一个由S 平面到Z 平面的简单代数映射关系，即没有一个s=f(z)的代数关系式。

另外，数字滤波器的频响也不是简单的重现模拟滤波器的频响，而是模拟滤波器频响的周期延拓，周期为Ωs =2π/T=2πf s ，即
∑∑∞-∞=⎪⎪⎪⎭⎫ ⎝⎛∞-∞=⎪⎪⎪⎭⎫ ⎝⎛+=+Ω=m a m a j T m j H T T m j j H T e H πωπω2121)(
如果模拟滤波器的频响带限于折叠频率Ωs /2以内，即
0)(=Ωj H a |Ω|≥π/T
这时数字滤波器的频响才能不失真地重现模拟滤波器的频响
)(1)(T
j H T e H a j ωω= |ω|＜π 但任何一个实际的模拟滤波器，其频响都不可能是真正带限的，因此不可避免地存在频谱的交叠，即混淆，这时，数字滤波器的频响将不同于原模拟滤波器的频响而带有一定的失真。

模拟滤波器频响在折叠频率以上衰减越大，失真则越小，这时，采用脉冲响应不变法设计的数字滤波器才能得到良好的效果。

虽然脉冲响应不变法能保证S 平面与Z 平面的极点位置有一一对应的代数关系，但这并不是说整个S 平面与Z 平面就存在这种一一对应的关系，特别是数字滤波器的零点位置与S 平面上的零点就没有一一对应关系，而是随着Ha(s)的极点si 与系数Ai 的不同而不同。

H(ejω) 是Ha(jΩ)的周期延拓（周期为f s ），因Ha(jΩ)并不是带限，即在
超过f s 频率部分并不为0，所以就产生了混迭。

当为低通或带通滤波器时，f s 越大，则Ha(jΩ)的下一周期相隔越远，混迭也就越小。

当为带阻或高通滤波器时，Ha(jΩ)在超过f s /2频率部分全为通带，这样就不满足抽样定理，发生了完全的
混迭，所以脉冲响应不变法不能设计带阻或高通滤波器。

双线性变换法原理：
脉冲响应不变法的主要缺点是产生频率响应的混叠失真。

这是因为从S 平面到Ｚ平面是多值的映射关系所造成的。

为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T ～π/T 之间，再用z e st
=转换到Z 平面上。

也就是说，第一步先将整个S 平面压缩映射到S1平面的-π/T ～π/T 一条横带里；第二步再通过标准变换关系1z e s t =将此横带变换到整个Z 平面上
去。

这样就使S 平面与Z 平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图2所示。

图2双线性变换的映射关系
为了将S 平面的整个虚轴j Ω压缩到S1平面j Ω1轴上的-π/T 到π/T 段上，可以通过以下的正切变换实现 ⎪⎭
⎫ ⎝⎛Ω=
Ω2tan 21T T （1） 式中,T 仍是采样间隔。

当Ω1由-π/T 经过0变化到π/T 时，Ω由-∞经过0变化到+∞，也即映射了整个j Ω轴。

将式（1）写成 2/2/2/2/11112T j T j T j T j e e e e T j ΩΩΩΩ+-=Ω
将此关系解析延拓到整个S 平面和S1平面，令j Ω=s ，j Ω1=s1，则得 T T s s T j T j T j T j e e T T s T e e e e T s 1111111122tan 2212/2/2/2/+-=⎪⎭⎫ ⎝⎛=+-=ΩΩΩΩ
再将S1平面通过以下标准变换关系映射到Z 平面z=e s1T ，从而得到S 平面和Z 平面的单值映射关系为： 1
1
112--+-=z z T s （2）
s T s T T s T z -+=-+
=222121 （3）
式（2）与式（3）是S 平面与Z 平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换。

式（1）与式（2）的双线性变换符合映射变换应满足的两点要求。

首先,把z=e j Ω，可得 Ω=⎪⎭
⎫ ⎝⎛=+-=--j T j e e T s j j 2tan 2112ωωω （4） 即S 平面的虚轴映射到Z 平面的单位圆。

其次，将s=σ+j Ω代入式（4），得
因此
由此看出，当σ<0时，|z|<1；当σ>0时，|z|>1。

也就是说，S 平面的左半平面映射到Z 平面的单位圆内，S 平面的右半平面映射到Z 平面的单位圆外，S 平面的虚轴映射到Z 平面的单位圆上。

因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。

双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。

这是因为S 平面与Z 平面是单值的一一对应关系。

S 平面整个j Ω轴单值地对应于Z 平面单位圆一周，即频率轴是单值变换关系。

这个关系如式（4）所示，重写如下：
上式表明，S平面上Ω与Z平面的ω成非线性的正切关系，如图3所示。

由图3看出，在零频率附近，模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系；但当Ω进一步增加时，ω增长得越来越慢，最后当Ω→∞时，ω终止在折叠频率ω=π处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。

但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如前边式（4）及图3所示。

图3双线性变换法的频率变换关系
由于这种频率之间的非线性变换关系，就产生了新的问题。

首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；
还有一点，就是这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图4所示。

图4双线性变换法幅度和相位特性的非线性映射
对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。

也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。

由于双线性变换法获得的数字滤波器频率响应特性中不会出现混叠现象，因此可以适用于高通、带通和带阻滤波器的设计。

四、设计步骤
脉冲响应不变法设计IIR 数字低通滤波器步骤：
（1）将数字滤波器设计指标转换为相应的模拟滤波器指标。

ππ51042⨯==
ΩT
f p rad/s 1=p α dB ππ51062⨯==ΩT f s rad/s 25=s α dB （2）设计相应的模拟滤波器，得到模拟系统函数H a (s)。

根据单调下降要求，选择巴特沃斯滤波器。

8919.341
101101.01.0=--=p s sp K αα 5.122==p
s sp f f ππλ 76.8lg lg ==
sp sp k N λ 取N=9 （3）按照下面公式，将模拟滤波器系统函数H a (s)转换成数字滤波器系统函数H(z)。

∑=-=9
1)(k k a k s s A s H ，∑=--=9111)(k T S k z e TA z H k 按照上边的式子求解计算相当复杂，我们可以利用MATLAB 软件来简化计算过程，调用MATLAB 信号处理工具箱函数进行设计。

双线性法设计IIR 数字低通滤波器步骤：
（1）列出数字低通技术指标。

ππω4002==p p f rad 1=p α dB
ππω6002==s s f rad 25=s α dB
（2）这里T=0.001s ，预畸变校正计算相应模拟低通的技术指标为。

2
tan 2p p T ω=Ω 1=p α dB 2tan 2s s T ω=
Ω 25=s α dB （3）设计巴特沃斯低通模拟滤波器。

据算阶数N 。

p
s sp ΩΩ=λ 11011010/10/--=
p s sp K αα sp
sp k N λlg lg = 求出Ωc ,，保证阻带技术指标满足要求，通带指标有富余。

根据N 的值查表
得到归一化方程G a (p)，将p=s/Ωc 带入G a (p)，去归一化得到实际的H a (s)。

（4）用双线性变换法将H a (s)转换成数字滤波器H(z),即 11
112|)()(--+-==z z s a s H z H
按照上式计算的过程也是非常复杂，我们可以通过MATLAB 软件使计算得到
简化。

根据所给出的技术指标并加以转换，转换之后，调用MATLAB 软件，然后应
用MATLAB 软件设计MATLAB 程序来简化计算过程，而且应用MATLAB 软件可以得
到相应的滤波器图像，从而更直观的观察所设计的滤波器的相关性质，程序中采
用双线性变换法，调用MATLAB 工具箱函数buttord 和butter 直接设计数字滤波
器
五、数字低通滤波器MATLAB编程及幅频特性曲线
5.1 MATLAB语言编程
（1）模拟滤波器设计程序：
Fs=1000; %采样频率为周期倒数
Wp=400*pi;
Ws=600*pi; %设置归一化通带和阻带截止频率 Ap=1;
As=25; %设置通带最大和最小衰减
N,Wc]=buttord(Wp,Ws,Ap,As,'s'); %调用buttord函数确定巴特沃斯
滤波器阶数
[B,A]=butter(N,Wc,'s'); %调用butter函数设计巴特沃斯滤
波器
W=(0:pi:1000*pi); %指定一段频率值
hs=freqs(B,A,W); %计算模拟滤波器的幅频响应
plot(W/pi/2,abs(hs)/abs(hs(1))); %绘出巴特沃斯模拟滤波器的幅频
特性曲线
grid on;
title('巴特沃斯模拟滤波器');
xlabel('频率 /Hz');
ylabel('归一化幅值 ');
(2) 用脉冲响应不变法设计数字滤波器的程序：
Fs=1000; %采样频率为周期倒数
Wp=400*pi;
Ws=600*pi; %设置归一化通带和阻带截止频率 Ap=1;
As=25; %设置通带最大和最小衰减
[N,Wc]=buttord(Wp,Ws,Ap,As,'s'); %调用buttord函数确定巴特沃斯滤
波器阶数
[B,A]=butter(N,Wc,'s'); %调用butter函数设计巴特沃斯滤
波器
W=(0:pi:1000*pi); %指定一段频率值
[D,C]=impinvar(B,A,Fs); %调用脉冲响应不变法
W=(0:0.001*pi:pi);
Hz=freqz(D,C,W); %返回频率响应
plot(W/pi,abs(Hz)/abs(Hz(1))); %绘出巴特沃斯数字低通滤波器的
幅频特性曲线
grid on;
title('巴特沃斯数字滤波器（脉冲响应不变法）');
xlabel('频率 /Hz');
ylabel('归一化幅值');
(3)用双线性变换法设计数字滤波器的程序：
Fs=1000; %采样频率为周期倒数
Wp=400*pi;
Ws=600*pi; %设置归一化通带和阻带截
止频率
Ap=1;
As=25; %设置通带最大和最小衰减
[N,Wc]=buttord(Wp,Ws,Ap,As,'s'); %调用buttord函数确定巴特沃斯滤
波器阶数
[B,A]=butter(N,Wc,'s'); %调用butter函数设计巴特沃斯滤
波器
W=(0:pi:1000*pi); %指定一段频率值
[D,C]=bilinear(B,A,Fs); %调用双线性变换法
W=(0:0.001*pi:pi);
Hz=freqz(D,C,W); %返回频率响应
plot(W/pi,abs(Hz)/abs(Hz(1))); %绘出巴特沃斯数字低通滤波器的
幅频特性曲线
grid on;
title('巴特沃斯数字滤波器（双线性变换法）');
xlabel('频率 /Hz');
ylabel('归一化幅值');
5.2 幅频特性曲线
六、总结
低通滤波器的应用十分广泛，通过MATLAB，可以很容易地设计出巴特沃斯数字低通滤波器，而且通过MATLAB设计滤波器更方便、快捷,节省了很多实践,提高了编程效率,并且参数的修改也十分方便. 还可以进一步进行优化设计。

通过这几天的设计过程，让我感触最深的就是MATLAB功能的强大性与掌握编程各种函数和语句的重要性。

当然首先要了解所要编程运行的对象的原理。

在课程设计的过程中，我深深的感受到我们所掌握的知识还不够，需要学习的东西太多了，真的是学无止境。

学习的过程很辛苦，但是我们也学到了很多知识，最后感谢老师给我们这次学习的机会，让我发现自己的不足，让我知道以后应该更投入的学习知识，提高自己的学习能力。

七、参考文献
[1] 《数字信号处理——（第2版）》，高西全丁玉美编著，西安电子科技大学
出版社，2001年
[2] 《MATLAB信号处理详解》，陈亚勇主编.，人民邮电出版社，2001年
[3]《数字信号处理教程——MATLAB释疑与实现》陈怀琛编著，电子工业出
版社，2004年
[4] 《信号处理原理及应用》，谢平王娜林洪彪编著，机械工业出版社，2008
年
[5] 《数字信号处理——（第3版）》，高西全丁玉美编著，清华大学出版社，
2010年。