实验八IIR滤波器的算法.doc

实验八IIR滤波器的算法一、实验目的和任务1理解并掌握IIR数字滤波器的机理； 2. 掌握设计IIR数字滤波器的原理和方法；3. 观察IIR滤波器的频域特性。二、实验原理1、设计原理（1）脉冲响应不变法所谓脉冲相应不变法就是使数字滤波器的单位脉冲响应序列h(n)等于模拟滤波器的单位冲击响应和(t)的采样值。（2）双向性变换法：由于s平面和z平面的单值双线性映射关系为s其中T为采样周期。因此，若已知模拟滤波器的传递函数，将上式代入即可得到数字滤波器的系统函数H(z)。在双线性变换中，模拟角频率和数字角频率的变换关系为 可见，和w之间的变换关系为非线性的。3三、IIR滤波器的设计参考实例一：IIR低通滤波（Butterworth）% IIR Lowpass Use Butterworthclear;fs=20;fpass=4;fstop=5;Ap=0.5;As=10;wp=2*pi*fpass/fs;ws=2*pi*fstop/fs;omegap=tan(wp/2);omegas=tan(ws/2);ep=sqrt(10(Ap/10)-1);es=sqrt(10(As/10)-1);N=ceil(log(es/ep)/log(omegas/omegap);omega0=omegap/ep(1/N);K=floor(N/2);for i=1:K theta(i)=pi*(N-1+2*i)/(2*N);endfor i=1:K G(i)=omega02/(1-2*omega0*cos(theta(i)+omega02);endfor i=1:K a1(i)=2*(omega02-1)/(1-2*omega0*cos(theta(i)+omega02);endfor i=1:K a2(i)=(1+2*omega0*cos(theta(i)+omega02)/(1-2*omega0*cos(theta(i)+omega02);endif K(N/2) G0=omega0/(omega0+1);a0=(omega0-1)/(omega0+1);endw=0:pi/300:pi;Hw2=1./(1+(tan(w/2)/omega0).(2*N);plot(w/pi,Hw2);grid;参考实例二：IIR高通滤波（Butterworth）% IIR Hightpass Use Butterworthclear;fs=20;fpass=5;fstop=4;Ap=0.5;As=10;wp=2*pi*fpass/fs;ws=2*pi*fstop/fs;omegap=cot(wp/2);omegas=cot(ws/2);ep=sqrt(10(Ap/10)-1);es=sqrt(10(As/10)-1);N=ceil(log(es/ep)/log(omegas/omegap);omega0=omegap/ep(1/N);K=floor(N/2);for i=1:K theta(i)=pi*(N-1+2*i)/(2*N);endfor i=1:K G(i)=omega02/(1-2*omega0*cos(theta(i)+omega02);endfor i=1:K a1(i)=-2*(omega02-1)/(1-2*omega0*cos(theta(i)+omega02);endfor i=1:K a2(i)=(1+2*omega0*cos(theta(i)+omega02)/(1-2*omega0*cos(theta(i)+omega02);endif K(N/2) G0=omega0/(omega0+1);a0=-(omega0-1)/(omega0+1);endw=(0+eps):pi/300:pi;Hw2=1./(1+(cot(w/2)/omega0).(2*N);plot(w/pi,Hw2);grid;参考实例三：IIR带通滤波（Butterworth）% IIR Bandpass Use Butterworthclear;fs=20;fpa=2;fpb=4;fsa=1.5;fsb=4.5;Ap=0.0877;As=16.9897;wpa=2*pi*fpa/fs;wpb=2*pi*fpb/fs;wsa=2*pi*fsa/fs;wsb=2*pi*fsb/fs;c=sin(wpa+wpb)/(sin(wpa)+sin(wpb);omegap=abs(c-cos(wpb)/sin(wpb);omegasa=(c-cos(wsa)/sin(wsa);omegasb=(c-cos(wsb)/sin(wsb);omegas=min(abs(omegasa),abs(omegasb);ep=sqrt(10(Ap/10)-1);es=sqrt(10(As/10)-1);N=ceil(log(es/ep)/log(omegas/omegap);omega0=omegap/ep(1/N);K=floor(N/2);for i=1:K theta(i)=pi*(N-1+2*i)/(2*N);endfor i=1:K G(i)=omega02/(1-2*omega0*cos(theta(i)+omega02);endfor i=1:K a1(i)=4*c*(omega0*cos(theta(i)-1)/(1-2*omega0*cos(theta(i)+omega02);endfor i=1:K a2(i)=2*(2*c2+1-omega02)/(1-2*omega0*cos(theta(i)+omega02);endfor i=1:K a3(i)=-(4*c*(omega0*cos(theta(i)+1)/(1-2*omega0*cos(theta(i)+omega02);endfor i=1:K a4(i)=(1+2*omega0*cos(theta(i)+omega02)/(1-2*omega0*cos(theta(i)+omega02);endif K(N/2) G0=omega0/(1+omega0);a0(1)=-2*c/(1+omega0);a0(2)=(1-omega0)/(1+omega0);endw=(0+eps):pi/300:pi;Hw2=1./(1+(c-cos(w)./(omega0*sin(w).(2*N);plot(w/pi,Hw2);grid;参考实例四：IIR带阻滤波（Butterworth）% IIR Bandstop Use Butterworthclear;fs=20;fpa=1.5;fpb=4.5;fsa=2;fsb=4;Ap=0.5;As=10;wpa=2*pi*fpa/fs;wpb=2*pi*fpb/fs;wsa=2*pi*fsa/fs;wsb=2*pi*fsb/fs;c=sin(wpa+wpb)/(sin(wpa)+sin(wpb);omegap=abs(sin(wpb)/(c-cos(wpb);omegasa=sin(wsa)/(cos(wsa)-c);omegasb=sin(wsb)/(cos(wsb)-c);omegas=min(abs(omegasa),abs(omegasb);ep=sqrt(10(Ap/10)-1);es=sqrt(10(As/10)-1);N=ceil(log(es/ep)/log(omegas/omegap);omega0=omegap/ep(1/N);K=floor(N/2);theta=zeros(1,K);for i=1:K theta(i)=pi*(N-1+2*i)/(2*N);endG=zeros(1,K);a1=zeros(1,K);a2=zeros(1,K);for i=1:K G(i)=omega02/(1-2*omega0*cos(theta(i)+omega02);endfor i=1:K a1(i)=2*(omega02-1)/(1-2*omega0*cos(theta(i)+omega02);endfor i=1:K a2(i)=(1+2*omega0*cos(theta(i)+omega02)/(1-2*omega0*cos(theta(i)+omega02);endif K(N/2) G0=omega0/(omega0+1);a0=(omega0-1)/(omega0+1);endw=(0+eps):pi/300:pi;Hw2=1./(1+(sin(w)./(omega0*(c-cos(w).(2*N);plot(w/pi,Hw2);grid;参考实例五：IIR低通滤波（chebyshev 1）% IIR Lowpass Use Chebyshev Type 1clear;fs=20;fpass=4;fstop=5;Ap=0.5;As=10;wp=2*pi*fpass/fs;ws=2*pi*fstop/fs;omegap=tan(wp/2);omegas=tan(ws/2);ep=sqrt(10(Ap/10)-1);es=sqrt(10(As/10)-1);e=es/ep;w=omegas/omegap;N=ceil(log(e+sqrt(e2-1)/log(w+sqrt(w2-1);a=log(1/ep+sqrt(1/ep2+1)/N;omega0=omegap*sinh(a);K=floor(N/2);theta=zeros(1,K);omega=zeros(1,K);for i=1:K theta(i)=pi*(N-1+2*i)/(2*N);endfor i=1:K omega(i)=omegap*sin(theta(i);endG=zeros(1,K);a1=zeros(1,K);a2=zeros(1,K);for i=1:K G(i)=(omega02+omega(i)2)/(1-2*omega0*cos(theta(i)+omega02+omega(i)2);endfor i=1:K a1(i)=2*(omega02+omega(i)2-1)/(1-2*omega0*cos(theta(i)+omega02+omega(i)2);endfor i=1:K a2(i)=(1+2*omega0*cos(theta(i)+omega02+omega(i)2)/(1-2*omega0*cos(theta(i)+omega02+omega(i)2);endif K(N/2) G0=omega0/(omega0+1);a0=(omega0-1)/(omega0+1); else H0=sqrt(1/(1+ep2);endf=0:1/300:10;Hf2=1./(1+ep2*(cheby(N,tan(pi*f/fs)/omegap).2);plot(f,abs(Hf2);grid;参考实例六：IIR低通滤波（chebyshev 1）% IIR Lowpass Use Chebyshev Type 2clear;fs=20;fpass=4;fstop=5;Ap=0.5;As=10;wp=2*pi*fpass/fs;ws=2*pi*fstop/fs;omegap=tan(wp/2);omegas=tan(ws/2);ep=sqrt(10(Ap/10)-1);es=sqrt(10(As/10)-1);e=es/ep;w=omegas/omegap;N=ceil(log(e+sqrt(e2-1)/log(w+sqrt(w2-1);a=log(1/ep+sqrt(1/ep2+1)/N;omega0=omegap*sinh(a);K=floor(N/2);theta=zeros(1,K);omega=zeros(1,K);for i=1:K theta(i)=pi*(N-1+2*i)/(2*N);endfor i=1:K omega(i)=omegap*sin(theta(i);endG=zeros(1,K);a1=zeros(1,K);a2=zeros(1,K);for i=1:K G(i)=(omega02+omega(i)2)/(1-2*omega0*cos(theta(i)+omega02+omega(i)2);endfor i=1:K a1(i)=2*(omega02+omega(i)2-1)/(1-2*omega0*cos(theta(i)+omega02+omega(i)2);en