实验三-数字IIR滤波器的设计.doc

实验三 数字IIR滤波器的设计一. 实验目的1. 加深对选频滤波概念的理解；2. 学会用MATLAB设计IIR数字低通、高通、带通和带阻滤波器。二. 实验原理2.1 数字滤波器的结构2.1.1 直接型传递函数形式用分子和分母多项式的系数两个向量来表示。Num=b(1) b(2) b(nb+1)Den=a(1) a(2) a(na+1)例：num=1 -3 11 -27 18;den=16 12 2 -4 -1;2.1.2 零极点增益形式用零点向量z、极点向量p以及增益K表示num=1 -3 11 -27 18;den=16 12 2 -4 -1;z,p,k=tf2zp(num,den)z = 0.0000 + 3.0000i0.0000 - 3.0000i2.0000 1.0000 p =-0.5000 + 0.5000i-0.5000 - 0.5000i0.5000 -0.2500 k =0.06252.1.3 级联型二阶因子级联形式用的数组SOS表示sos=zp2sos(z,p,k)sos = 0.0625 -0.1875 0.1250 1.0000 -0.2500 -0.1250 1.0000 -0.0000 9.0000 1.0000 1.0000 0.50002.1.4 并联型部分分式展开式形式用函数residuez（）把传递函数形式转换为部分分式形式。例：num=1 -3 11 -27 18;den=16 12 2 -4 -1;r,p,c=residuez(num,den)r = -5.0250 - 1.0750i -5.0250 + 1.0750i 0.9250 27.1875 p = -0.5000 + 0.5000i -0.5000 - 0.5000i 0.5000 -0.2500 c =-18然后将共轭极点两两合并成如下的一阶节和二阶节的并联。MATLAB Signal Processing Toolbox 中的系统模型之间的相互转换函数有：z,p,k=tf2zp(b,a)b,a=zp2tf(z,p,k)sos=zp2sos(z,p,k)z,p,k=sos2zp(sos)b,a=sos2tf(sos)r,p,k=residuez(b,a)b,a=residuz(r,p,k)2.2 IIR数字滤波器设计2.2.1 MATLAB中模拟滤波器设计函数介绍 设计流程图：给定模拟滤波器的技术指标采用真实角频率（弧度/秒）计算模拟滤波器的阶次N和截止频率Wc。(利用buttord,cheb1ord,cheb2ord,ellipord等)设计模拟低通滤波器原型G(p)。(利用buttap,cheb1ap,cheb2ap,ellipap等)模拟滤波器完全设计函数H(s)（butter,cheby1,cheby2,ellip）由模拟低通原型经频率变换获得所需要的低通、高通、带通和带阻滤波器H(s)。（利用lp2lp,lp2hp,lp2bp,lp2bs）1. 阶次计算函数；(1) Butterworth模拟滤波器N, Wn = buttord(Wp, Ws, Rp, Rs, s)；(2) Chebyshev I型模拟滤波器N, Wn = cheb1ord(Wp, Ws, Rp, Rs, s)(3) Chebyshev II型模拟滤波器N, Wn = cheb2ord(Wp, Ws, Rp, Rs, s)(4) 椭圆滤波器 N, Wn = ellipord(Wp, Ws, Rp, Rs, s)其中：Wp, Ws分别是通带和阻带的截止频率，单位为弧度/秒。对低通和高通，Wp, Ws都是标量，对带通和带阻，Wp, Ws是12的向量。Rp, Rs 分别是通带和阻带的衰减(dB)。N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率，它和Wp稍有不同。2模拟低通原型滤波器设计函数 （注意：得到的是归一化的传递函数）（1） Butterworth模拟滤波器z, p, k=buttap(N)；滤波器传递函数具有如下形式： (2) Chebyshev I型模拟滤波器z, p, k=cheb1ap(N,Rp)滤波器传递函数具有如下形式：(3) Chebyshev II型模拟滤波器z, p, k= cheb2ap (N,Rs)滤波器传递函数具有如下形式： (4) 椭圆滤波器z, p, k=ellipap(N,Rp,Rs)滤波器传递函数具有如下形式： 其中，N是欲设计的低通原型滤波器的阶次，z, p, k是设计出的极点、零点及增益。3频率变换函数将模拟低通原型转换为实际的低通、高通、带通及带阻滤波器。 (1) B, A=lp2lp(b, a, wo),(2) B, A=lp2hp(b, a, wo)(3) B, A=lp2bp(b, a, wo, Bw),(4) B, A=lp2bs(b, a, wo, Bw)b, a 是AF LP 的分子、分母的系数向量，B, A是转换后的的分子、分母的系数向量；在（1）、（2）中，wo是低通或高通滤波器的截止频率；在（3）、（4）中，wo是带通或带阻滤波器的中心频率，Bw是其带宽。， 4模拟滤波器完全设计函数（=2+3）(1) Butterworth模拟滤波器b,a=butter(N,Wn,s) % 低通或带通滤波器b,a=butter(N,Wn,high,s) %高通滤波器b,a=butter(N,Wn,stop,s) %带阻滤波器 (2) Chebyshev I型模拟滤波器b,a= cheby1 (N,Rp,Wn,s) % 低通或带通滤波器b,a= cheby1 (N,Rp,Wn,high,s) %高通滤波器b,a= cheby1(N,Rp,Wn,stop,s) %带阻滤波器(3) Chebyshev II型模拟滤波器b,a= cheby2 (N,Rs,Wn,s) % 低通或带通滤波b,a= cheby2 (N,Rs,Wn,high,s) %高通滤波器b,a= cheby2(N,Rs,Wn,stop,s) %带阻滤波器 (4) 椭圆滤波器b,a= ellip(N,Rp,Rs,Wn,s) % 低通或带通滤波器b,a= ellip (N,Rp,Rs,Wn,high,s) %高通滤波器b,a= ellip(N,Rp,Rs,Wn,stop,s) %带阻滤波器例1：设计一个Butterworth模拟低通滤波器，通带截至频率fp=5KHz，阻带截至频率fs=10KHz；通带最大衰减Ap=3，阻带最小衰减As=30。画出该滤波器的频率响应。%Design a Butterworth Analog lowpass filter%SpecificationsWp=5000*2*pi;Ws=10000*2*pi;Ap=3;As=30;%Compute order and Cutoff frequencyN,Wc=buttord(Wp,Ws,Ap,As,s);%design analog lowpass filter prototypez,p,k=buttap(N) %归一化的传递函数G(p)b,a=zp2tf(z,p,k);%frequency transform,得到反归一化传递函数G(s)bt,at=lp2lp(b,a,Wp);w=0:20*2*pi:20000*2*pi;H=freqs(bt,at,w);subplot(211),plot(w/(2*pi),abs(H);xlabel(Frequency(Hz);ylabel(Magnitude);title(Butterworth Analog lowpass filter);grid%Using filter whole design fuctionbb,aa=butter(N,Wc,s);w=0:20*2*pi:20000*2*pi;HH=freqs(bb,aa,w);subplot(212),plot(w/(2*pi),abs(HH);xlabel(Frequency(Hz);ylabel(Magnitude);title(Butterworth Analog lowpass filter);grid例2：设计一个Chebyshev模拟带通滤波器，通带截至频率分别为1000Hz和2000Hz，阻带截至频率分别为500Hz和2500Hz；通带最大衰减Ap=1，阻带最小衰减As=100。画出该滤波器的频率响应。%Design a Chebyshev Analog bandpass filter%SpecificationsWp=1000 2000*2*pi;Ws=500 2500*2*pi;Ap=1;As=100%Compute order and Cutoff frequencyN,Wc=cheb1ord(Wp,Ws,Ap,As,s);%design analog lowpass filter prototypez,p,k=cheb1ap(N,Ap);b,a=zp2tf(z,p,k);%find the center frequency and bandwidthWo=sqrt(Wp(1)*Wp(2);Bw=Wp(2)-Wp(1);%frequency transformbt,at=lp2bp(b,a,Wo,Bw);w=0:20*2*pi:4000*2*pi;H=freqs(bt,at,w);subplot(211),plot(w/(2*pi),abs(H);xlabel(Frequency(Hz);ylabel(Magnitude);title( Chebyshev Analog bandpass filter);grid%Using filter whole design fuctionbb,aa=cheby1(N,Ap,Wc,s);w=0:20*2*pi:4000*2*pi;H=freqs(bb,aa,w);subplot(212),plot(w/(2*pi),abs(H);xlabel(Frequency(Hz);ylabel(Magnitude);title( Chebyshev Analog bandpass filter);grid计算模拟滤波器的阶次N和截止频率Wc。(利用buttord,cheb1ord,cheb2ord,ellipord等)设计模拟低通滤波器原型G(p)。(利用buttap,cheb1ap,cheb2ap,ellipap等)由模拟低通原型经频率变换获得所需要的低通、高通、带通和带阻滤波器H(s)。（利用lp2lp,lp2hp,lp2bp,lp2bs）将模拟滤波器离散化，获得IIR数字滤波器H(z)（利用impinvar,bilinear）转成模拟滤波器的技术指标给定数字滤波器的技术指标给定数字滤波器的技术指标2.2.2 MATLAB中IIR数字滤波器设计函数给采用归一化频率01对应采用真实角频率（弧度/秒）计算数字滤波器的阶次N和截止频率Wc。(利用buttord,cheb1ord,cheb2ord,ellipord等)模拟滤波器完全设计函数H(s)（butter,cheby1,cheby2,ellip）数字滤波器完全设计函数H(z)（butter,cheby1,cheby2,ellip）1. 阶次计算函数；(1) Butterworth数字滤波器N, Wn = buttord(Wp, Ws, Rp, Rs)；(2) Chebyshev I型数字滤波器N, Wn = cheb1ord(Wp, Ws, Rp, Rs)(3) Chebyshev II型数字滤波器N, Wn = cheb2ord(Wp, Ws, Rp, Rs)(4) 椭圆数字滤波器 N, Wn = ellipord(Wp, Ws, Rp, Rs)其中：Wp, Ws分别是通带和阻带的数字标准化频率，取值范围在01之间，标准化频率1对应的数字频率为。对低通和高通，Wp, Ws都是标量，对带通和带阻，Wp, Ws是12的向量。Rp, Rs 分别是通带和阻带的衰减(dB)。N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率，它和Wp稍有不同。2数字滤波器完全设计函数（采用的是bilinear变换）(1) Butterworth数字滤波器b,a=butter(N,Wn) % 低通或带通滤波器b,a=butter(N,Wn,high) %高通滤波器b,a=butter(N,Wn,stop) %带阻滤波器 (2) Chebyshev I型数字滤波器b,a= cheby1 (N,Rp,Wn) % 低通或带通滤波器b,a= cheby1 (N,Rp,Wn,high) %高通滤波器b,a= cheby1(N,Rp,Wn,stop) %带阻滤波器(4) Chebyshev II型数字滤波器b,a= cheby2 (N,Rs,Wn) % 低通或带通滤波器b,a= cheby2 (N,Rs,Wn,high) %高通滤波器b,a= cheby2(N,Rs,Wn,stop) %带阻滤波器 (4) 椭圆滤波器b,a= ellip(N,Rp,Rs,Wn,s) % 低通或带通滤波器b,a= ellip (N,Rp,Rs,Wn,high) %高通滤波器b,a= ellip(N,Rp,Rs,Wn,stop) %带阻滤波器3impinvar 用冲激响应不变法实现模拟滤波器到数字滤波器的变换。Bz, Az=impinvar(B, A, Fs)4bilinear 双线性变换用双线性变换实现模拟滤波器到数字滤波器的变换。Bz, Az=bilinear(B, A, Fs)式中B, A分别是模拟滤波器传递函数G（s）的分子、分母多项式的系数向量， Bz, Az分别是数字滤波器系统函数H(z)的分子、分母多项式的系数向量，Fs是抽样频率。例3：设计一个Butterworth数字低通滤波器，通带截至频率fp=100Hz，阻带截至频率fs=300Hz，采样频率Fs=1000Hz；通带最大衰减Ap=3，阻带最小衰减As=20。画出该滤波器的频率响应。用该滤波器对50Hz和400Hz的两个正弦信号之和x=sin(2*pi*50*n/Fs)+sin(2*pi*400*n/Fs)进行滤波，画出滤波前后的信号波形。%Design a Butterworth Digital lowpass filter%Using bilinear Trasform% DF specificationsfp=100;fs=300;Fs=1000;Ap=3;As=20;T=1/Fs;%Convert digital frequency to analog frequencyWp=(2/T)*tan(fp*pi/Fs);Ws=(2/T)*tan(fs*pi/Fs);%Compute order and cutoff frequencyN,Wc=buttord(Wp,Ws,Ap,As,s)%computer the analog filter prototypez,p,k=buttap(N);B,A=zp2tf(z,p,k)%frequency transformb,a=lp2lp(B,A,Wc);%bilinear Trasformbz,az=bilinear(b,a,1/T)disp(Butterworth lowpass filter);%w=0:0.01*pi:pi;Nn=128;freqz(bz,az,Nn,Fs);%-%Nomalized Angular FrequencyWp=fp*2/Fs;Ws=fs*2/Fs;%Compute order and cutoff frequencyN,Wc=buttord(Wp,Ws,Ap,As); %Desigh digital highpass filterdz,cz=butter(N,Wc)%outputfigure,w=0:0.01*pi:pi;freqz(dz,cz,w); %-n=0:100;x=sin(2*pi*50*n/Fs)+sin(2*pi*400*n