8.用双线性变换法设计IIR数字滤波器---数字信号处理实验报告

计算机与信息工程学院验证性实验报告专业：通信工程年级/班级：2011级2013—2014学年第一学期课程名称数字信号处理指导教师段新涛本组成员学号姓名实验地点计科楼111实验时间周五7-8节项目名称用双线性变换法设计IIR数字滤波器实验类型验证性一、实验目的(1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法。(2)掌握数字滤波器的计算机仿真方法。(3)通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。二、实验内容和要求1、用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为：在通带内频率低于0.2π时，最大衰减小于1dB；在阻带内[0.3π，π]频率区间上，最小衰减大与15dB。2、以0.02π为采样间隔，打印出数字滤波器在频率区间[0，π/2]的幅频响应特性曲线。3、用所设计的滤波器对实际心电图信号采样序列（在本实验后面给出）进行仿真滤波处理，并打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。三、实验主要仪器设备和材料计算机，MATLAB四、实验方法复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容，用双线性变换法设计数字滤波器系统函数zH。其中满足本实验要求的数字滤波器系统函数为：212121612155.09044.013583.00106.117051.02686.1110007378.0zzzzzzzzHzHkk31（4.1）式中：3211212121，，，kzCzBzzAzHkkk（4.2）2155.09044.03583.00106.17051.02686.109036.0332211CBCBCBA，，，根据设计指标，调用MATLAB信号处理工具箱buttord和butter，也可以得到zH。由公式（3.1）和（3.2）可见，滤波器zH由三个二阶滤波器zH1、zH2和zH3级联而成，如图3-1所示。编写滤波器仿真程序，计算zH对心电图采样序列x(n)的响应序列y(n)。设yk(n)为第k级二阶滤波器Hk(z)的输出序列，yk-1(n)为输入序列，如图3-1所示。由（4.2）式可得到差分方程为：21212111nyCnyBnAynAynAynykkkkkkkk（4.3）当k=1时，yk-1(n)=x(n)。所以H(z)对x(n)的总响应序列y(n)可以用顺序迭代算法得到。即依次对k=1，2，3，求解差分方程（4.3），最后得到y3(n)=y(n)。仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序。也可以直接调用MATLABfilter函数实现仿真。zH1zH2zH3nxny1ny2nyny3图3-1滤波器zH的组成在通用计算机上运行仿真滤波程序，并调用通用绘图子程序，完成实验内容（2）和（3）。本实验中用到的心电图信号采用序列x(n)人体心电图信号在测量过程往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。在实验中，以x(n)作为输入序列，滤除其中的干扰成分。五、实验程序及结果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(1,2,1);plot(fk,20*log10(abs(Hk)));gridon;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(1,2,2);plot(wk/pi,20*log10(abs(Hz)));gridon;xlabel('\omega/\pi');ylabel('幅度(dB)');axis([0,1,-100,5]);title('(b)');00.51-100-500/幅度(dB)(b)00.51-100-500/幅度(dB)(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];k=1;closeall;figure(1);subplot(2,2,1);n=0:55;%更正stem(n,x,'.');axis([0,56,-100,50]);%更正holdon;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);%Thefunctionistofiltethesingalxx=y;ifk==2A=A1;endifk==3A=A2;endk=k+1;endsubplot(2,2,3)n=0:55;%更正stem(n,y,'.');axis([0,56,-15,5]);holdon;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]);gridon;%更正title('滤波器的幅频响应曲线');02040-100-50050nx(n)心电图信号采集序列x(n)02040-15-10-505ny(n)三级滤波后的心电图信号00.10.20.30.4-40-200滤波器的幅频响应曲线教师签名：年月日