实验三 用双线性变换法设计IIR数字滤波器

实验项目名称：用双线性变换法设计IIR数字滤波器 实验项目性质：验证性实验 所属课程名称：数字信号处理 实验计划学时：2 一. 实验目的

(1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法。 (2)掌握数字滤波器的计算机仿真方法。

(3)通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。 二. 实验内容和要求

(1)用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数

为：在通带内频率低于0.2π时，最大衰减小于1dB；在阻带内[0.3π，π]频率区间上，最小衰减大与15dB。 (2)以0.02π为采样间隔，打印出数字滤波器在频率区间[0，π/2]的幅频

响应特性曲线。 (3)用所设计的滤波器对实际心电图信号采样序列（在本实验后面给出）进

行仿真滤波处理，并打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。 三. 实验主要仪器设备和材料 计算机，MATLAB6.5或以上版本 四. 实验方法、步骤及结果测试

(1)复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波

器的内容，用双线性变换法设计数字滤波器系统函数Hz。其中满足本实验要求的数字滤波器系统函数为：

Hz11.2686z10.7051z20.00073781z1 121211.0106z0.3583z10.9044z0.2155z6Hkz

k13 （3.1）

A12z1z2式中： Hkz，k1，2，3 121BkzCkzA0.09036B11.2686，C10.7051 （3.2）

B21.0106，C20.3583B30.9044，C30.2155根据设计指标，调用MATLAB信号处理工具箱buttord和butter，也可以得到Hz。

由公式（3.1）和（3.2）可见，滤波器Hz由三个二阶滤波器H1z、H2z和H3z级联而成，如图3-1所示。

(2)编写滤波器仿真程序，计算Hz对心电图采样序列x(n)的响应序列y(n)。

xn z(n) 为输入序列，如图3-1Hykz(z) 的输出序列， 级二阶滤波器H2H设yk(n)H为第3k-11zk所示。由（3.2）式可得到差分方程为： n1yknAyk1n图2Ay2Bkykn1Ckykn2（3.3）3-1 zk的组成HAyk1滤波器1n 当k=1时，yk-1(n)=x(n)。所以H(z)对x(n)的总响应序列y(n)可以用顺序迭代

算法得到。即依次对k=1，2，3，求解差分方程（3.3），最后得到y3(n)=y(n)。仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序。也可以直接调用MATLAB filter函数实现仿真。

(3)在通用计算机上运行仿真滤波程序，并调用通用绘图子程序，完成实验

内容（2）和（3）。 (4)本实验中用到的心电图信号采用序列x(n)

y1n y2n y3nyn

人体心电图信号在测量过程往往受到工业高频干扰，所以必须经过低通滤波

处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。在实验中，以x(n)作为输入序列，滤除其中的干扰成分。

{x(n)}={-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}

T=1;Fs=1/T;

wpz=0.2;wsz=0.3;

wp=2*tan(wpz*pi/2);ws=2*tan(wsz*pi/2);rp=1;rs=15; [N,wc]=buttord(wp,ws,rp,rs,'s'); [B,A]=butter(N,wc,'s'); fk=0:1/512:1;wk=2*pi*fk; Hk=freqs(B,A,wk); subplot(2,2,1);

plot(fk,20*log10(abs(Hk)));grid on; xlabel('\\omega/\\pi');ylabel('幅度(dB)'); axis([0,1,-100,5]);title('(b)'); [N,wc]=buttord(wpz,wsz,rp,rs); [Bz,Az]=butter(N,wc); wk=0:pi/512:pi;

Hz=freqz(Bz,Az,wk); subplot(2,2,4);

plot(wk/pi,20*log10(abs(Hz)));grid on; xlabel('\\omega/\\pi');ylabel('幅度(dB)'); axis([0,1,-100,5]);title('(b)');

二；

x=[-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]; subplot(2,2,1); n=0:55; stem(n,x,'.'); xlabel('n'); ylabel('x(n)');

title('x(n)的脉冲响应'); A=0.09036; b1=[A,2*A,A];

a1=[1,-1.2686,0.7051]; h1=filter(b1,a1,x);

[H1,w]=freqz(b1,a1,100); b2=[A,2*A,A];

a2=[1,-1.0106,0.3583]; h2=filter(b2,a2,h1);

[H2,w]=freqz(b2,a2,100); b3=[A,2*A,A];

a3=[1,-0.9044,0.2155]; h3=filter(b3,a3,h2);

[H3,w]=freqz(b3,a3,100); subplot(2,2,2); stem(n,h3,'.'); xlabel('n'); ylabel('y(n)');

title('通过滤波器H1(z),H2(z),H3(z)后的y3(n)函数'); subplot(2,2,3); H4=H1.*(H2); H=H4.*(H3); mag=abs(H);

db=20*log10((mag+eps)/max(mag)); plot(w/pi,db); xlabel('ω/π');

ylabel('20log[Ha3(ejw)]');

title('通过滤波器H1(z),H2(z),H3(z)后的对数频率响应20log[Ha3(ejw)]函数'); grid; figure(2);

N=1024;n=0:N/2-1;

Xk=fft(x,N);AXk=abs(Xk(1:N/2)); f=(0:N/2-1)*Fs/N; f=f/Fs;

subplot(211);plot(f,AXk);title('x(n)的频谱'); xlabel('f');ylabel('| X(k) |');axis([0,0.5,0,400]);

Yk=fft(y,N);AYk=abs(Yk(1:N/2));

subplot(212);plot(f,AYk);title('y(n)的频谱'); xlabel('f');ylabel('| Y(k) |');axis([0,0.5,0,400])

%x(n)的心电脉冲函数

x=[-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] subplot(2,2,1); n=0:55; stem(n,x,'.'); xlabel('n'); ylabel('x(n)');

title('x(n)的心电脉冲函数');

%通过滤波器H1(z)后的y1(n)函数 A=0.09036; b1=[A,2*A,A];

a1=[1,-1.2686,0.7051]; h1=filter(b1,a1,x);

[H1,w]=freqz(b1,a1,100);

%通过滤波器H1(z),H2(Z)后的y2(n)函数 b2=[A,2*A,A];

a2=[1,-1.0106,0.3583]; h2=filter(b2,a2,h1);

[H2,w]=freqz(b2,a2,100);

%通过滤波器H1(z),H2(Z),H3(Z)后的y3(n)函数 b3=[A,2*A,A];

a3=[1,-0.9044,0.2155]; h3=filter(b3,a3,h2);

[H3,w]=freqz(b3,a3,100); subplot(2,2,2); stem(n,h3,'.'); xlabel('n'); ylabel('y(n)');

title('通过滤波器H1(Z),H2(Z),H3(Z)后的y3(n)函数'); subplot(2,2,3); H4=H1.*(H2); H=H4.*(H3); mag=abs(H);

db=20*log10((mag+eps)/max(mag)); plot(w/pi,db); xlabel('w/pi');

ylabel('20log[Ha3(ejw)]');

title('通过滤波器H1(z)、H2(z)、H3(z)后的对数频率响应20log[Ha3(ejw)]函数'); grid;

MATLAB程序清单：

%实验三，用双线性变换法设计IIR数字滤波器

x=[-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]; k=1; close all; figure(1); subplot(2,2,1); n=0:55; %更正 stem(n,x,'.');

axis([0,56,-100,50]); %更正 hold on; n=0:60; m=zeros(61); plot(n,m); xlabel('n'); ylabel('x(n)');

title('心电图信号采集序列x(n)'); B=[0.09036,2*0.09036,0.09036]; A=[1.2686,-0.7051]; A1=[1.0106,-0.3583]; A2=[0.9044,-0.2155]; while(k<=3)

y=filter(B,A,x); %The function is to filte the singal x x=y; if k==2 A=A1; end if k==3 A=A2; end k=k+1;

end

subplot(2,2,3) n=0:55; %更正 stem(n,y,'.'); axis([0,56,-15,5]); hold on; n=0:60; m=zeros(61); plot(n,m); xlabel('n'); ylabel('y(n)');

title('三级滤波后的心电图信号'); %求数字滤波器的幅频特性 A=[0.09036,0.1872,0.09036]; B1=[1,-1.2686,0.7051]; B2=[1,-1.0106,0.3583]; B3=[1,-0.9044,0.2155]; [H1,w]=freqz(A,B1,100); [H2,w]=freqz(A,B2,100); [H3,w]=freqz(A,B3,100); H4=H1.*(H2); H=H4.*(H3); mag=abs(H);

db=20*log10((mag+eps)/max(mag)); subplot(2,2,2) plot(w/pi,db); axis([0,0.5,-50,10]); grid on; %更正

title('滤波器的幅频响应曲线');

程序运行结果：

直接运行程序，结果输出滤波器幅频特性曲线图，有噪声的心电图采集信号

波形图和经过三级二阶滤波器滤波后的心电图信号波形图，见图3-2，对比分析就可以看出低通滤波器滤除信号中高频噪声的滤波效果。

五. 实验报告要求 (1) 简述实验原理及目的。

(2) 由所打印的| H(ejw) |特性曲线及设计过程简述双线性变换法的特点。 (3) 对比滤波前后的心电图信号波形，说明数字滤波器的滤波过程与滤波作用。 (4) 简要回答思考题。 六. 思考题