IIR数字滤波器设计和应用.doc

IIR数字滤波器设计和应用一、实验目的1、熟悉IIR数字滤波器的设计方法。2、掌握模拟滤波器的matlab实现。3、熟悉用脉冲响应不变法设计IIR数字滤波器的原理与方法。4、熟悉用双线性变换法设计IIR数字滤波器的原理与方法。5、掌握数字滤波器的计算机仿真方法。6、通过观察对实际心电图信号的滤波作用， 获得数字滤波的感性知识。二、实验内容和要求1、用Matlab语言分别设计巴特奥斯低通滤波器和切比雪夫低通滤波器，其技术指标为：通带截止频率，通带最大衰减；阻带起始频率，阻带最小衰减；要求：求出他们的零点、极点、阶数、增益等，并画出图形作比较。（1）巴特奥斯低通滤波器：%Design a Butterworth Analog bandpass filter%Desired performents of the filterwp=5000*2*pi; ws=10000*2*pi; Rp=3; Rs=30; %compute oder and cutoff frequency N,Wn=buttord(wp,ws,Rp,Rs,s); N Fc=Wn/(2*pi) %Compute the filter b,a=butter(N,Wn,s);%Outputw=linspace(1,3000,1000)*2*pi;H=freqs(b,a,w);magH=abs(H);phaH=unwrap(angle(H);plot(w/(2*pi),20*log10(magH);xlabel(Frequency(Hz);ylabel(Magnidute(dB);title(Butterworth Analog bandpass filter)grid on运算结果：N =5Fc = 5.0124e+003（2）切比雪夫低通滤波器：%Chebyshev lowpass filter %Desired performentswp=5000*2*pi;ws=10000*2*pi;Rp=3;Rs=30;%Compute Order Nebs=sqrt(10(Rp/10)-1);A=10(Rs/20);Wc=wpWr=ws/wp;g=sqrt(A*A-1)/ebs;N1=log10(g+sqrt(g*g-1)/log10(Wr+sqrt(Wr*Wr-1);N=ceil(N1)运算结果：Wc = 3.1416e+004N = 42、切比雪夫低通滤波器，其技术指标为：通带截止频率，通带最大衰减；阻带起始频率，阻带最小衰减；要求：求出他们的零点、极点、阶数、增益等，并画出图形作比较。3、 用脉冲响应不变法把模拟滤波器 转换为数字滤波器H(z)，其中采样周期为T=0.1s。%Impulse invariance method %From analog to digital filter conversion%Model of analog filter b=1 1;a=1 5 6;Ts=0.1;R,Ps, K=residue(b,a);Pz=exp(Ps*Ts);disp(Use direct principle)bz, az=residue(R,Pz,K) disp(Use function IMPINVAR)bz1, az1=impinvar(b,a,1/Ts)运算结果Use direct principle：bz = 1.0000 -0.8966az = 1.0000 -1.5595 0.6065Use function “IMPINVAR”：bz1 = 0.1000 -0.0897az1 = 1.0000 -1.5595 0.60654、利用双线性变换法设计一个巴特奥斯低通IIR数字滤波器。(1) 设计的指标参数为：在通带内频率低于0.2时，最大衰减小于1dB，在阻带内0.3，频率区间上，最小衰减大于15dB. (2) 以 0.02为采样间隔， 打印出数字滤波器在频率区间0, /2上的幅频响应特性曲线。 (3) 用所设计的滤波器对实际心电图信号采样序列(在本实验后面给出)进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图， 观察总结滤波作用与效果。例如假设最后求出的数字滤波器系统函数H(z)为：式中A=0.090 36 B1=1.2686,C1=-0.7051 B2=1.0106,C2=-0.3583 B3=0.9044,C3=-0.2155 由(1式)和(2式)可见， 滤波器H(z)由三个二阶滤波器H1(z), H2(z)和H3(z)级联组成， 如图 1 所示。图1 滤波器H(z)的组成(2) 编写滤波器仿真程序, 计算H(z)对心电图信号采样序列x(n)的响应序列y(n)。设yk(n)为第k级二阶滤波器Hk(z)的输出序列， y k-1(n)为输入序列， 如图1 所示。 由(2式)可得到差分方程： yk(n)=Ayk-1(n)+2Ayk-1(n-1)+Ayk-1(n-2)+Bkyk(n-1)+Ckyk(n-2) (3式) (3) 在通用计算机上运行仿真滤波程序， 并调用通用绘图子程序， 完成实验内容(2)和(3)。人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。在实验中，以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;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 xx=y;if k=2;A=A1;endif k=3A=A2;endk=k+1;endsubplot(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.18072, 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); title(滤波器的幅频响应曲线) ;三、思考题1. 答：二者一致。2.答：（1）优点：一定是递归结构、对频率分量的选择性好（零极点可同时起作用）、相同性能下阶次较低、运算误差大，有可能出现极限环振荡、设计有封闭形式的公式，一次完成