1、實驗三 數字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 =

2、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 級聯(lián)型二階因子級聯(lián)形式用的數組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 并聯(lián)型部分分式展開式形式用函數residuez（）把傳遞函數形式轉換為部分分式形式。例：num=1 -

3、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然后將共軛極點兩兩合并成如下的一階節(jié)和二階節(jié)的并聯(lián)。MATLAB Signal Processing Toolbox 中的系統(tǒng)模型之間的相互轉換函數有：z,p,k=tf2zp(b,a)b,a=zp2tf(z,p,k)sos=zp2sos(z,p,k)z,p,k=

4、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）由模擬低通原型經頻率變換獲得所

5、需要的低通、高通、帶通和帶阻濾波器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分別是通帶和阻帶的截止頻率，單位為弧度

6、/秒。對低通和高通，Wp, Ws都是標量，對帶通和帶阻，Wp, Ws是1×2的向量。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)濾波器傳遞函數

7、具有如下形式： (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是低通或

8、高通濾波器的截止頻率；在（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) %帶阻濾波器(

9、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；通帶最大

10、衰減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 tran

11、sform,得到反歸一化傳遞函數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

12、*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

13、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(

14、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

15、,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,che

16、b2ap,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,che

17、by2,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之

18、間，標準化頻率1對應的數字頻率為。對低通和高通，Wp, Ws都是標量，對帶通和帶阻，Wp, Ws是1×2的向量。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) % 低通

19、或帶通濾波器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) %帶

20、阻濾波器3impinvar 用沖激響應不變法實現(xiàn)模擬濾波器到數字濾波器的變換。Bz, Az=impinvar(B, A, Fs)4bilinear 雙線性變換用雙線性變換實現(xiàn)模擬濾波器到數字濾波器的變換。Bz, Az=bilinear(B, A, Fs)式中B, A分別是模擬濾波器傳遞函數G（s）的分子、分母多項式的系數向量， Bz, Az分別是數字濾波器系統(tǒng)函數H(z)的分子、分母多項式的系數向量，F(xiàn)s是抽樣頻率。例3：設計一個Butterworth數字低通濾波器，通帶截至頻率fp=100Hz，阻帶截至頻率fs=300Hz，采樣頻率Fs=1000Hz；通帶最大衰減Ap=3，阻帶最小衰減As=

21、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)*

22、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/Fs);imp=1;zeros