实验四 IIR数字滤波器的设计
1、实验目的
（1） 掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的具体方法和原理，熟悉双线性变换法和脉冲响应不变法设计低通、带通IIR数字滤波器的计算机编程；
（2） 观察双线性变换法和脉冲响应不变法设计的数字滤波器的频域特性，了解双线性变换法和脉冲响应不变法的特点和区别；
（3） 熟悉Butterworth滤波器、Chebyshev滤波器和椭圆滤波器的频率特性。
2、实验原理与方法
IIR数字滤波器的设计方法可以概括为如图所示，本实验主要掌握IIR滤波器的第一种方法，即利用模拟滤波器设计IIR数字滤波器，这是IIR数字滤波器设计最常用的方法。利用模拟滤波器设计，需要将模拟域的Ha(s)转换为数字域H(z)，最常用的转换方法为脉冲响应不变法和双线性变换法。
（1）脉冲响应不变法
用数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应ha(t)，让h(n)正好等于ha(t)的采样值，即
其中T为采样间隔。如果以Ha(s)及H(z)分别表示ha(t)的拉氏变换及h(n)的Z变换，则
在MATLAB中，可用函数impinvar实现从模拟滤波器到数字滤波器的脉冲响应不变映射。
（2）双线性变换法
S平面与z平面之间满足下列映射关系
或 
S平面的虚轴单值地映射于z平面的单位圆上，s平面的左半平面完全映射到z平面的单位圆内。双线性变换不存在频率混叠问题。
在MATLAB中，可用函数bilinear实现从模拟滤波器到数字滤波器的双线性变换映射。
双线性变换是一种非线性变换，即
，这种非线性引起的幅频特性畸变可通过预畸变得到校正。
（3）设计步骤
IIR数字滤波器的设计过程中，模拟滤波器的设计是关键。模拟滤波器的设计一般是采用分布设计的方式，这样设计原理非常清楚，具体步骤如前文所述。MATLAB信号处理工具箱也提供了模拟滤波器设计的完全工具函数：butter、cheby1、cheby2、ellip、besself。用户只需一次调用就可完成模拟滤波器的设计，这样虽简化了模拟滤波器的设计过程，但设计原理却被屏蔽了。
模拟滤波器设计完成之后，利用impinvar或bilinear函数将模拟滤波器映射为数字滤波器，即完成了所需数字滤波器的设计。
下图给出了实际低通、高通、带通和带阻滤波器的幅频特性和各截止频率的含义。另外，为了描述过渡带的形状，还引入了通带衰减和阻带衰减的概念。
图 实际滤波器的幅频特性和各截止频率的含义
通带衰减：
dB
阻带衰减：
dB
在MATLAB信号处理工具箱中，通常用Rp和Rs来表示αp和αs。
3、实验内容
（1）参照教材5.5节所述滤波器设计步骤，利用双线性变换法设计一个Chebyshev I型数字高通滤波器，观察通带损耗和阻带衰减是否满足要求。已知滤波器的指标为fp=0.3kHz，αp=1.2dB，fs=0.2kHz，αs=20dB，T=1ms。
（2）已知fp=0.2kHz，αp=1dB，fs=0.3kHz，αs=25dB，T=1ms，分别用脉冲响应不变法和双线性变换法设计一个Butterworth数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。
（3）设计一个数字带通滤波器，通带范围为0.25π~0.45π，通带内最大衰减3dB，0.15π以下和0.55π以上为阻带，阻带内最小衰减为15dB，试采用Butterworth或ellip（椭圆）模拟低通滤波器设计。
（4）利用双线性变换法设计一个带宽为0.08π的10阶椭圆带阻滤波器以滤除数字频率为0.44π的信号，选择合适的阻带衰减值，画出幅度响应。产生下面序列的201个样本
， n=0，2，…，200
并将它通过这个带阻滤波器进行处理（filter函数），讨论所得到的结果。
4、实验报告
（1）简述实验目的和实验原理。
（2）按实验步骤附上所设计的滤波器传递函数H(z)及相应的幅频特性曲线，定性分析所得到的图形，判断设计是否满足要求。
（3）总结脉冲响应不变法和双线性变换法的特点及设计全过程。
（4）收获与建议。
5、实验源程序
%用双线性变换法设计一个Chebyshev型高通滤波器程序如下
Rp=1.2;Rs=20;T=0.001;fp=300;fs=200;
%求出待设计的数字滤波器的边界频率
wp=2*pi*fp*T;ws=2*pi*fs*T;
%预畸变
wp1=(2/T)*tan(wp/2);ws1=(2/T)*tan(ws/2);
%设计模拟滤波器
[n,wn]=cheb1ord(wp1,ws1,Rp,Rs,'s');
[b,a]=cheby1(n,Rp,wn,'high','s');
%双线性变换
[bz,az]=bilinear(b,a,1/T);
[db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi,db);
axis([0,1,-30,2]);
%用双线性变换法设计一个Butterworth型数字低通滤波器程序如下
Rp=1;Rs=25;T=0.001;fp=300;fs=200;
%求出待设计的数字滤波器的边界频率
wp=2*pi*fp*T;ws=2*pi*fs*T;
%预畸变
wp1=(2/T)*tan(wp/2);ws1=(2/T)*tan(ws/2);
%设计模拟滤波器
[n,wn]=buttord(wp1,ws1,Rp,Rs,'s');
[b,a]=butter(n,wn,'low','s');
%双线性变换
[bz,az]=bilinear(b,a,1/T);
[db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi,db);
axis([0,1,-30,2]);
%用脉冲响应不变法设计一个Butterworth数字低通滤波器的程序如下：
%模拟滤波器的技术要求
wp=400*pi;ws=600*pi;Rp=1;Rs=25;
%求模拟滤波器的系统函数
[n,wn]=buttord(wp,ws,Rp,Rs,'s')
[b,a]=butter(n,wn,'s')
%求模拟滤波器的频率响应，w取（0~1000pi）rad/s
[db,mag,pha,w]=freqs_m(b,a,500*2*pi);
%绘图，为了使模坐标显示频率f （单位Hz），将原变量w（模拟角频率，单位为rad/s）进行了处理
plot(w/(2*pi),db,'LineWidth',2,'Color','b');
axis([0,500,-20,1]);hold on
%脉冲响应不变法
fs=1000;[bz,az]=impinvar(b,a,fs);
%求数字滤波器的频率响应
[db,mag,pha,w]=freqz_m(bz,az);
%绘图，为了与模拟滤波器的频响在同一坐标中绘出，需要将数字频率w转换成模拟频率f，转换公式为f=w*fs/2*pi
plot(0.5*fs*w/pi,db,'LineWidth',2,'Color','r');
axis([0,599,-20,1]);hold off
%采用ellip（椭圆）模拟低通滤波器设计，其程序如下：
%确定所需类型数字滤波器的技术指标
Rp=3;Rs=15;T=0.001;
wp1=0.25*pi;wp2=0.45*pi;ws1=0.15*pi;ws2=0.55*pi;
%将所需类型数字滤波器的技术指标转换成模拟滤波器的技术指标
wp3=(2/T)*tan(wp1/2);wp4=(2/T)*tan(wp2/2);
ws3=(2/T)*tan(ws1/2);ws4=(2/T)*tan(ws2/2);
%将所需类型数字滤波器的技术指标转换成模拟滤波器的技术指标,设计模拟滤波器
wp=[wp3,wp4];ws=[ws3,ws4];
[n,wn]=ellipord(wp,ws,Rp,Rs,'s');[z,p,k]=ellipap(n,Rp,Rs);[b,a]=zp2tf(z,p,k);
%频率更换
w0=sqrt(wp3*wp4);Bw=wp4-wp3;
[b1,a1]=lp2bp(b,a,w0,Bw);
%双线性变换法
[bz,az]=bilinear(b1,a1,1/T);
[db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi,db);
axis([0,1,-50,2]);
%采用Butterworth模拟低通滤波器设计，其程序如下：
%确定所需类型数字滤波器的技术指标
Rp=3;Rs=15;T=0.001;
wp1=0.25*pi;wp2=0.45*pi;ws1=0.15*pi;ws2=0.55*pi;
%将所需类型数字滤波器的技术指标转换成模拟滤波器的技术指标
wp3=(2/T)*tan(wp1/2);wp4=(2/T)*tan(wp2/2);
ws3=(2/T)*tan(ws1/2);ws4=(2/T)*tan(ws2/2);
%将所需类型数字滤波器的技术指标转换成模拟滤波器的技术指标,设计模拟滤波器
wp=[wp3,wp4];ws=[ws3,ws4];
[n,wn]=buttord(wp,ws,Rp,Rs,'s');[z,p,k]=buttap(n);[b,a]=zp2tf(z,p,k);
%频率更换
w0=sqrt(wp3*wp4);Bw=wp4-wp3;
[b1,a1]=lp2bp(b,a,w0,Bw);
%双线性变换法
[bz,az]=bilinear(b1,a1,1/T);
[db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi,db);
axis([0,1,-50,2]);
第二篇：实验四 用双线性变换法设计IIR数字滤波器
实验四 用双线性变换法设计IIR数字滤波器
一、 实验目的
1、熟悉用双线性变换法设计IIR设计数字滤波器的原理与方法。
2、掌握数字滤波器的计算机仿真方法。
3、通过观察对实际心电图信号的滤波作用,观察数字滤波的感性知识。
二、 实验内容及步骤
(1) 用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为:在通带内频率低于0.2π时,最大衰减小于1dB；在阻带内[0.3π , π]频率区间上,最小衰减大于15dB。
(2) 打印出数字滤波器在频率区间[ 0, 0.5π]上的幅频衰减曲线，和[?4π,4π]上的幅频响应曲线。
(3) 用所设计的滤波器对实际心电图信号采样序列(在本实验后面给出)进行仿真滤波处理,并分别打印出滤波前后的心电图波形图,观察总结滤波作用与效果。
3．实验步骤
（1）复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计 IIR 数字滤波器 的内容，用双线性变换法设计满足设计指标的数字滤波器系统函数 H(z)。
0.0007378(1+z?1)6
H(Z)=(1?1.2686z?1+0.705z?2)(1?1.0106z?1+0.3583z?2)(1?0.904z?1+0.2155z?2)
=∏HK(z) (1.1)k=13
A(1+2z?1+z?2), κ=1,2,3 (1.2) Hk(z)=?1?21?Bkz?Ckz
A=0.09036
B1=1.2686, C1=?0.7051 B2=1.0106, C3=?0.3583
Β3=0.9044, C1=?0.2155
由（1.1）式和（1.2）式可见，滤波器H(z)由三个二阶滤波器H1(z)、H2(z)和H3(z)级联组成，如图1,1 所示。
图 1.1
（2）打印出数字滤波器在频率区间[ 0, 0.5π]上的幅频衰减曲线，和[?4π,4π]上的幅频响应曲线。主要使用的MATLAB函数：freqz() ( 参考电脑桌面上数字信号处理文件夹内的文件名为matlab的pdf文档第4、5页的内容) 。
下面的M文件举例说明如何绘制第一个二阶滤波器的幅频衰减曲线和幅频响应曲线。
部分函数功能简介：
max(x): 向量x的元素的最大值
eps：系统的浮点（Floating-point）精确度
axis: 坐标轴的控制函数，调用格式如下
axis([xmin,xmax,ymin,ymax])
grid on 在图形中绘制坐标网格
(3）编写滤波器仿真程序，计算H(z)对心电图信号采样序列x(n)的响应序列y(n)。设 yk(n)为第 k 级二阶滤波器Hk(z)的输出序列，yk?1(n) 为输入序列，
如图 4.1所示。由（1.2）式可得到差分方程：
yk(n)=Ayk?1(n)+2Ayk?1(n?1)+Ayk?1(n?2)+Bkyk(n?1)+Ckyk(n?2) (1.3)
当κ=1时，yk-1(n)= x(n) 。所以H(z)对x(n)的总响应序列y(n)可以用顺序选代
算法得到。即依次对k=1,2,3，求解差分方程(1.3)，最后得到 y3(n)=y(n)。
仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序。也可以直接调用 MATLAB filter函数实现仿真。
调用MATLAB filter()函数对实际心电图信号滤波( 参考电脑桌面上数字信号处理文件夹内的文件名为matlab的pdf文档第3、4页的内容) 。并调用DFT()函数绘制滤波前、后信号的波形图，和幅频特性曲线。
下面的M文件举例说明如何绘制调用MATLAB filter()函数，完成第一个二阶滤波器对实际心电图信号滤波。
附录：心电图信号采样序列 x(n)
人体心电图信号在测量过程中往往受到工业高频干扰,所以必须经过低通滤波处理后,才能作为判断心脏功能的有用信息.下面给出一实际心电图信号采样序列样式本x(n),其中存在高频干扰.在实验中,以x(n)作为输入序列,滤除其中的干扰成分。
xn=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,
0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,
4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];
请大家注意：
本学期数字信号处理实验只写前三次的实验报告。请大家在本学期15周周末在开放实验中心（202.202.43.114）网站上查询、核对自己的实验成绩。