1、實驗四 IIR數字濾波器的設計及網絡結構一、實驗目的 1了解IIR數字濾波器的網絡結構。2.掌握模擬濾波器、IIR數字濾波器的設計原理和步驟。 3學習編寫數字濾波器的設計程序的方法。 二、實驗內容數字濾波器：是數字信號處理技術的重要內容。它的主要功能是對數字信號進行處理，保留數字信號中的有用成分，去除信號中的無用成分。1數字濾波器的分類濾波器的種類很多，分類方法也不同。（1）按處理的信號劃分：模擬濾波器、數字濾波器（2）按頻域特性劃分；低通、高通、帶通、帶阻。（3）按時域特性劃分：FIR、IIR 2IIR數字濾波器的傳遞函數及特點數字濾波器是具有一定傳輸特性的數字信號處理裝置。它的輸入和輸出均

2、為離散的數字信號，借助數字器件或一定的數值計算方法，對輸入信號進行處理，改變輸入信號的波形或頻譜，達到保留信號中有用成分去除無用成分的目的。如果加上A/D、D/A轉換，則可以用于處理模擬信號。設IIR濾波器的輸入序列為x(n)，則IIR濾波器的輸入序列x(n)與輸出序列y(n)之間的關系可以用下面的方程式表示： （5-1）其中，和是濾波器的系數，其中中至少有一個非零。與之相對應的差分方程為： （5-2）由傳遞函數可以發現無限長單位沖激響應濾波器有如下特點：（1） 單位沖激響應h(n)是無限長的。 （2） 系統傳遞函數H(z)在有限z平面上有極點存在。（3） 結構上存在著輸出到輸入的反饋，也就是

3、結構上是遞歸型的。3IIR濾波器的結構IIR濾波器包括直接型、級聯型和并聯型三種結構： 直接型：優點是簡單、直觀。但由于系數bm 、ak與零、極點對應關系不明顯，一個bm 或ak的改變會影響H(z)所有零點或極點的分布，所以一方面，bm 、ak對濾波器性能的控制關系不直接，調整困難；另一方面，零、極點分布對系數變化的靈敏度高，對有限字長效應敏感，易引起不穩定現象和較大誤差。Matlab實現：filter( )函數實現IIR數字濾波器直接形式。格式為： y=filter(b,a,x)b，a為差分方程輸入、輸出系數向量（或系統函數的分子、分母多項式，降冪），x為輸入序列，y為輸出序列。其中，傳遞函

4、數（tf）形式則a=1 a1 a2 aNb=b0 b1 b2 bM 級聯型：基于因式分解，將系統函數H(z)分解為因子乘積的形式。 (5-3） 級聯型結構：Matlab實現：tf2zp( )函數用于求系統函數的零、極點和增益常數，zp2sos ( )函數則根據tf2zp( )函數結果求出各基本節系數。格式為：z,p,K=tf2zp(b,a)；sos=zp2sos(z,p,K)；b，a為差分方程輸入、輸出系數向量（系統函數的分子、分母多項式，降冪）。其中，零極點增益形式（zp）：則零點向量z=z1 z2 zM-1;極點向量p=z1,z2,zN-1k為系統增益。二階分式形式（sos）為：把H(z)

5、劃成二階因式則其二階因式為： 并聯型：基于部分分式展開，將系統函數H(z)分解為部分分式和的形式。 （5-4）并聯型結構：Matlab實現：residue( )函數可以實現并聯型結構，有兩種格式：K,r,p=residue(b,a); b,a=residue(b,a);其中，部分分式形式：則極點向量 p=p(1) p(2) p(n)其對應系數向量r=r(1) r(2) r(n)余數多項式系數向量k=k(1) k(2) k(M-N+1)【實例5-1】已知三階IIR數字濾波器的系統函數 求：直接形式的單位采樣響應h(n)；級聯型結構的各基本節系數；并聯型結構的部分分式系數。解：MATLAB源程序為

6、 b=3,5/3,2/3;a=1,1/6,1/3,-1/6;x=1,zeros(1,50);y=filter(b,a,x);n=0:50;plot(n,y);b=3,5/3,2/3,0;a=1,1/6,1/3,-1/6; z,p,K=tf2zp(b,a);sos=zp2sos(z,p,K);b=3,5/3,2/3;a=1,1/6,1/3,-1/6;K,r,p=residue(b,a);KK1=K(1),K(2);zz1=z(1),z(2);b2,a2=residue(KK1,zz1,0);5IIR數字濾波器的具體設計（1）巴氏模擬原型濾波器的設計巴氏模擬低通濾波器幅度平方函數為MATLAB工具

7、箱函數buttap,buttord和butter是巴氏濾波器設計函數。1）Z,P,K=buttap(N)該格式用于計算N階巴氏歸一化（3dB截止頻率c=1）模擬低通原型濾波器系統函數的零、極點和增益。得到的系統函數為：如果要從計算得到的零、極點得到系統函數的分子和分母向量B和A，可以調用結構轉換函數B,A=zp2tf(Z,P,K)。2）N,wc=buttord(wp,ws,Rp,Rs,s)該格式用于計算巴氏濾波器的階數N和3dB截止頻率wc。其中：N濾波器的階數。Wc3dB截止頻率的歸一化值。wp、ws通帶、阻帶邊界頻率的歸一化值。要求：，1表示數字頻率。Rp、Rs通帶最大、阻帶最小衰減。s可

8、選項，直接設計模擬濾波器，此時wp和ws均為實際模擬角頻率。3）B,A=butter(N,wc,ftype，s)計算N階巴氏數字濾波器系統函數分子和分母多項式的系數向量B和A。 【實例5-2】已知通帶截止頻率，通帶最大衰減，阻帶截止頻率，阻帶最小衰減，設計巴特沃斯型模擬低通濾波器。解：MATLAB源程序為Wp=2*pi*5000;Ws=2*pi*12000;Rp=5;As=30;N,Wc=buttord(Wp,Ws,Rp,Rs,s);B,A=butter(N,Wc,s);freqs(B,A);【實例5-3】設計階數為3，5，10，15的巴氏模擬原型濾波器。并畫出幅頻響應曲線。解：MATLAB源

9、程序為for i=1:4; switch i case 1 N=3; case 2 N=5; case 3 N=10; case 4; N=15; end;z,p,k=buttap(N);b,a=zp2tf(z,p,k);h,w=freqs(b,a,n);Ah=abs(h);subplot(2,2,i),plot(w,Ah);axis(0 2 0 1);xlabel(w/wc);ylabel(|H(jw)|.2);title(filer N=,num2str(N);grid;end;【實例5-4】設通帶、阻帶截止頻率fp=0.5kHz、fs=1.2kHz，通帶、阻帶最大衰減Rp=1dB,Rs=

10、30dB，要求設計巴氏低通濾波器。解：MATLAB源程序為: OmegaP=2*pi*500;OmegaS=2*pi*1200; Rp=1;Rs=30; N,OmegaC=buttord(OmegaP,OmegaS,Rp,Rs,s); %確定階數N z0,p0,k0=buttap(N);%確定傳遞函數 z=z0*OmegaC;%去歸一化 k=k0*OmegaCN; p=p0*OmegaC; bs=k*real(poly(z); as=real(poly(p); freqs(bs,as);（2）IIR數字濾波器的設計模擬濾波器Ha(s) 轉換成數字濾波器H(z)應滿足要求： (1) 因果穩定的模

11、擬濾波器轉換成數字濾波器，仍是因果穩定的；(2)數字濾波器的頻率響應模仿模擬濾波器的頻響，s平面的虛軸映射z平面的單位圓，相應的頻率之間成線性關系。脈沖響應不變法和雙線性變換法都滿足如上要求。脈沖響應不變法用數字濾波器的單位脈沖響應序列h(n)模仿模擬濾波器的沖激響應ha(t),讓h(n)正好等于ha(t)的采樣值，即h(n)=ha(nT)，其中T為采樣間隔。雙線性變換法s平面與z平面之間滿足以下映射關系：s平面的虛軸單值地映射于z平面的單位圓上，s平面的左半平面完全映射到z平面的單位圓內。雙線性變換不存在混疊問題。雙線性變換時一種非線性變換，這種非線性引起的幅頻特性畸變可通過預畸而得到校正。

12、以低通數字濾波器為例，將設計步驟歸納如下：(1)確定數字濾波器的性能指標：通帶臨界頻率fp、阻帶臨界頻率fs；通帶內的最大衰減Rp；阻帶內的最小衰減Rs； (2)確定相應的數字角頻率，p=2fp；s=2fs； (3)計算經過預畸的相應模擬低通原型的頻率，； (4)根據p和s計算模擬低通原型濾波器的階數N，并求得低通原型的傳遞函數Ha(s)； (5)用上面的雙線性變換公式代入Ha(s)，求出所設計的傳遞函數H(z)； (6)分析濾波器特性，檢查其指標是否滿足要求。1)沖激響應不變法MATLAB實現：采用impinvar函數bz,az=impinvar(b,a,fs,tol) bz,az轉換的分子

13、分母向量 fs采樣頻率，默認值為1Hz tol誤差容限，表示轉換后的離散系統函數是否有重復的極點。2）雙線性變換法MATLAB實現采用bilinear函數：bz, az=bilinear(bs,as,fs)把模擬濾波器的傳遞函數轉換為數字的其中：bs,as模擬濾波器的分子、分母向量bz,az數字濾波器的分子、分母向量【實例5-5】基于Butterworth模擬濾波器原型，使用雙線性轉換設計數字濾波器，其中參數指標為：通帶截止頻率： 通帶波動值：阻帶截止頻率： 阻帶波動值：解：首先確定濾波器的階數N，同時根據確定0.5。接著使用bilinear進行雙線性轉換，最后繪制在頻域上的各種圖像，其源代碼如下：數字濾波器指標wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;%將數字濾波器指標反轉變化為模擬濾波器的參數T=1;fs=1/T;omegap=(2/T)*tan(wp/T);omegas=(2/T)*tan(ws/T);%butterworth原型模擬濾波器的設計N,Wc=buttord(omegap,omegas,Rp,As,s);B,A=butter(N,Wc,s);%雙線性變換b,a=