IIR數字濾波器設計(I)數字濾波器簡介數字濾波器是對數字信號處理實現濾波的線性時不變系統。數字濾波實質上是一種運算過程，實現對信號的運算處理。輸入數字信號(數字序列)通過特定的運算轉變為輸出的數字序列，因此，數字濾波器本質上是一個完成特定運算的數字計算過程，也可以理解為是一臺計算機。描述離散系統輸出與輸入關系的卷積和差分方程只是給數字信號濾波器提供運算規則，使其按照這個規則完成對輸入數據的處理。時域離散系統的頻域特性：Y(ejw)=X(ejw)H(ejw)，其中Y(ejw),X(ejw)分別是數字濾波器的輸出序列和輸入序列的頻域特性(或稱為頻譜特性)，H(ejw)是數字濾波器的單位取樣響應的頻譜，又稱為數字濾波器的頻域響應。可以看出，輸入序列的頻譜X(ejw)經過濾波后變為X(ejw)H(e>)，因此，只要按照輸入信號頻譜的特點和處理信號的目的，適當選擇H(ejw)，使得濾波后的X(ejw)H(ejw)滿足設計的要求。IIR濾波器簡介數字濾波器在數字信號處理的各種應用中發揮著十分重要的作用。它是通過對采樣數據信號進行數學運算處理來達到濾波的目的。其中無限沖擊響應數字濾波器也稱IIR是采用對離散采樣數據作差分方程運算來進行濾波。IIRDF的優點在于可以利用模擬濾波器設計的結果，然后采用雙線性變換法或沖激響應不變法將模擬濾波器轉換成數字濾波器，而模擬濾波器的設計方便簡單、有大量的圖表可查。IIRDF的設計還可利用二階子系統的串聯來有效實現。其主要缺點在于Nyquist間隔內無法實現線性相位,若需線性相位，比如在圖象處理和數據傳輸中都要求信道具有線性相位特性，則需要用全通網絡進行相位校正。IIR濾波器的傳遞函數模型IIR數字濾波器的基本網絡結構可用差分方程、單位脈沖響應及系統函數來描述，其差分方程一般表示為y(n)=送bx(n-i)_工ay(n-i) 式(2.1)iii—0 i二1式中x(n)和y(n)分別表示輸入和輸出信號序列，a和b是濾波器系數。設ii式(2.1)中輸入信號x(n)和輸出信號y(n)在n-0以前處于零起始狀態，則有

H(Z)= =i=0 =i=11＋yaz-iyaz-1＋yaz-iyaz-i式(2.2)對于線性時不變(LTI)系統來說，=上式中a和b式(2.2)IIR濾波器的狀態方程模型狀態方程是描述系統的一種常用的方式。在系統有不可見的狀態變量時可使用濾波器狀態方程模型。連續LTI系統的狀態方程可以寫成x(t)=Ax(t)+Bu(t)y(t)=Cx(t)+Du(t)式(2.3)式(2.4)其中A,B，C,D分別為常數矩陣，在MATLAB中，一般情況下，系統的狀態方程可以簡記為(A,B,C,D),如果D三0,則系統的狀態方程模型可以簡記為(A,B,式(2.3)式(2.4)對于離散系統來說,狀態方程模型可以寫成x[(i+1)T]=Ax(iT)+Bu(iT)y[(i+1)T]=Cx(i+1)T+Du(i+1)T式(2.5)式(2.6)這樣的離散時間系統的狀態方程模型在MATLAB中寫成(A,B,C,D)或(A,B,C)式(2.5)式(2.6)IIR濾波器的零極點模型零極點模型實際是傳遞函數的另一種表達,它的原理是分別對原系統傳遞函數的分子和分母進行分解因式,以獲得系統的零極點表示形式,對單輸入單輸出系統來說,零極點模型可以簡單的表示為：my(s+z)i式(2.7)工(s+p) (s+pi)(s+p2)"式(2.7)工(s+p) (s+pi)(s+p2)"+叮ii它們既可以是實數也可以是復數，而K稱為系統的增益。在MATLAB中簡記為(Z,P,K)。但為了保證濾波器的穩定性，系統的所有極點都位于左半S平面。如果穩定系統所有的零點都位于左半S平面，則稱該系統為最小相位系統，否則稱為非最小相位系統。如果系統中的某個零點值恰好等于其中一個極點的值,則它們之間可以對消以直接獲得一個完全等效的低價系統。IIR數字濾波器的實現方法沖激響應不變法沖激響應不變法［5］的變換原理是使數字濾波器的單位沖激響應序列h(n)模仿模擬濾波器的單位沖激響應序列h(t)。將模擬濾波器的單位沖激響應加以等a簡隔抽樣，使h(n)正好等于h(t)的抽樣值，即滿足ah(n)=h(nT) 式(3.1)a其中T是抽樣周期。如果令H(s)是h(t)的拉普拉斯變換，H(z)為h(n)的Z變換，根據抽樣序aa列的Z變換和模擬信號的拉普拉斯變換的關系，可得H(z)I =￡H(s-jk) 式(3.2)z=esTT a Tk=—g可以看出，利用沖激響應不變法將模擬濾波器變換數字濾波器，實際上是先對模擬濾波器的系統函數H(s)做周期延拓，再經過z=esT的映射變換，從而得a到數字濾波器的系統函數h(z)。假設在S平面上，S在jq軸上取值，Z在Z平面內的單位圓周e?上取值，則可得到數字濾波器的頻率響應H(e。)和模擬濾波器的頻率響應H(jQ)間的對應關系：1占 CO—2兀k /oo\H(ejO)=￡H(j ) 式(3.3)TaTk=—g如圖3.1所示，在S平面上每一條寬為2兀/T的橫條都將重疊地映射到整個Z平面上，而每一條的左半邊映射到Z平面單位圓內，右半邊映射到Z平面單位圓外，S平面虛軸映射Z平面單位圓上，虛軸上每一段長為2兀/T的線段都映射到Z平面單位圓上一周。由于S平面每一橫條都要重疊地映射到Z平面上，這正好反映H(z)和H(s)的周期延拓函數的變換關系z=esT，故有沖激響應不變法a并不是簡單的從S平面映射到Z平面。' J^[z]—“江圖3.1沖激響應不變法的映射關系因為不是簡單的一一映射關系，且對于任何一個實際的模擬濾波器，它的頻率響應是不可能真正帶限的。因而將不可避免的出現頻率的重疊，即混疊失真。數字濾波器的頻率響應不能重現模擬濾波器的頻率響應，只有當模擬濾波器的頻率響應在超過重疊頻率后的衰減很大時，混疊失真才會很小，此時才能夠滿足設計要求。綜上所述，沖激響應不變法具有以下特點：模擬頻率和數字頻率的轉換是線性的，并保持了模擬濾波器的時域瞬態特性。當模擬濾波器的頻率響應不是嚴格限帶時，則用沖激不變法設計出的數字濾波器在頻率出現混疊失真。由于(2)而使這種設計方法受到限制，即當g(j)不嚴格限帶或g(t)變化不穩定。巴特沃斯低通濾波器的設計巴特沃斯濾波器的特性巴特沃斯濾波器的特點是同頻帶內的頻率響應曲線最大限度平坦，沒有起伏，而在阻頻帶則逐漸下降為零。在振幅的對數對角頻率的波特圖上

從某一邊界角頻率開始，振幅隨著角頻率的增加而逐步減少，趨向負無窮大。一階巴特沃斯濾波器的衰減率為每倍頻6分貝，每十倍頻20分貝。二階巴特沃斯濾波器的衰減率為每倍頻12分貝、三階巴特沃斯濾波器的衰減率為每倍頻18分貝、如此類推。巴特沃斯濾波器的振幅對角頻率單調下降，并且也是唯一的無論階數，振幅對角頻率曲線都保持同樣的形狀的濾波器。只不過濾波器階數越高，在阻頻帶振幅衰減速度越快。其他濾波器高階的振幅對角頻率圖和低級數的振幅對角頻率有不同的形狀。模擬低通巴特沃斯濾波器以巴特沃斯函數為濾波器的系統函數，它的幅度平方函數表示為：IH(IH(j0)|2a式中的N0=0時：c1— j01+￡2( )2Nj0c為正整數，表示濾波器的階數。0為通帶的截止頻率，因為當式(3.10)|H(j0)|2=丄(歸一化后的巴特沃斯濾波器，￡=1)a2|H(j0)|=丄，進而有20log(|Ha(j0)|)=3dB，可a 2 10|H(j0)|a c分貝。巴特沃斯濾波器在通帶有最為平坦的幅度特性，即N階低通濾波器在。二0處幅度平方函數|h(jQ)|2的前(2N-1)階為零。阻帶內隨著頻率的升高而單調下降。濾波器特征完全由階數N決定，當N增加時，通帶更加平坦，也更接近理想的低通濾波器的特性。低通巴特沃斯濾波器的設計過程包括以下兩個方面［7］：sA，0=0p根據給定的通帶和阻帶指標確定階數N。通常在設計時會給定技術指標,包括通帶的截止頻率0、通帶內最大衰減A(sA，0=0p最小衰減A(dB)。假設給定的，0=0時通帶的最大衰減為A，0=0時阻帶s c p s的最小衰減為A。其幅度平方函數圖如圖3.5所示。這時：S通帶容限為：|H(j0)|>-(|0|<0)；1+￡2c阻帶容限為：IH(jQ)l<1(l0l>Q)。a 1+九2 c其中S,九均為與衰減有關的參數。根據式(3.10)可得到:圖3.5巴特沃斯濾波器幅度平方函數1A=-10log( )=10log(1+圖3.5巴特沃斯濾波器幅度平方函數1A=-10log( )=10log(1+￡2)p 1+￡21A=—10log( )=10log(1+尢2)s 1+尢2由式(3.11a)、(3.11b)兩式可得：￡=(10。」Ap—1)0.5 ，X=(100.1As—1)0.5對于同一濾波器，其幅度平方函數相等，于是可得式(3.11a)式(3.11b)式(3.12)/jQ、

￡2( -)2NjQcX2進而可得N為正整數的取值要求：N>Xlog()10￡式(3.13)logjQ當給定的參數指標是歸一化的，即Q=1,￡=1時，上式表示為:log(九)N> to-log(0)10s⑵從幅度平方函數確定系統函數H(s)。令s=j0代入(2.8)式得到:a1H(s)H(-s)=—aa s1+￡2( )2Nj0c巴特沃斯濾波器的全部零點都在s=8處，沒有有限零點，只求解上式的分母多項式的根。當設計歸一化的巴特沃斯濾波器時，它的極點均勻分布在單位圓上。得到H(s)的極點為：a./2k一1 . 2k一1 k一{1，2…，(3)N+為奇數2 式(314)s,=-sin( 兀)+jcos( 兀)k={123 N 式(3?14丿pk 2N 2N 1，，2…3 為偶數 N令歸一化的系統函數用H(s令歸一化的系統函數用H(s)，則在左半平面上的極點組成H(s)，H(s)aN和H(s)的關系為H(s)=H(s/0)。a aN c的極點表達式形式可表示為：aNaNH(s)aH(s)

a1H(s)1H(s)=a(s一s)(s一s)(s一s)???(s-s)123pk式(3.15)⑶利用MATLAB設計巴特沃斯濾波器的實例。一、設計內容設計一個數字巴特沃斯低通濾波器，設計指標如下：① =① =0.2兀pR=1dBpA=15dB

s米樣時間間隔T=1S。二、設計要求1)用單位沖激響應不變變換法進行設計。（2）給出詳細的濾波器設計說明書。（3）給出經過運行是正確的程序清單并加上詳細的注釋。（4）畫出所設計濾波器的幅度特性和相位特性。三、設計說明（1） 把數字頻率轉換為模擬頻率：Q=WT，0=wT。ppss（2） 計算巴特沃斯模擬濾波器的截止頻率0和階數N。c（3） 設計巴特沃斯模擬低通濾波器，給出參數b和a（此處使用了MATLAB中的buttap（N）函數。（4） 把模擬濾波器用單位沖激響應不變變換法轉換成數字濾波器（此處使用了MATLAB中的residuez函數）。（5） 變直接形式為并聯形式，并給出結構圖。（6）畫出幅度特性和相位特性。四、設計步驟框圖數字濾波器的設計步驟如圖1所示參數變換巴特沃斯模擬低通濾波器巴特沃斯模擬濾波器離散化處理直接型轉變成為并聯型五、Matlab程序代碼%主程序wp=0.2*pi;ws=0.3*pi;Rp=1；%通帶波動；As=15;%阻帶衰減T=1;%采樣周期Omigrp=wp*T;Omigrs=ws*T;t1=[1,2*0.2588,0.2588人2+0.9659人2];b=1;a=conv(conv(t1,t1),t1);[N,wc]=buttord(wp,ws,Rp,As,'s')%計算巴特沃斯數字濾波器的階數N和截止頻率wc[B,A]=butter(N,wc,'s')[z,p,k]=buttap(N)%模擬低通原型零、極點系數和增益因子[bz,az]=impinvar(b,a)%用沖激響應不變法將模擬濾波器轉化為數字濾波器，采樣頻率默認1Hzwz=[0:pi/512:pi];hz1=freqz(bz,az,wz);%巴特沃斯模擬低通濾波器頻率響應[C,B,A]=dir2par(b,a)%直接型轉換成并聯型%繪圖subplot(1,2,1);plot(wz/pi,abs(hz1)/hz1(1));title('幅度響應')xlabel('');ylabel('|H|');axis([0,1,0,1]);grid;subplot(1,2,2);plot(wz/pi,hz1/pi);title('相位響應')xlabel('');ylabel('單位：pi');axis([0,1,-1,2]);grid;%直接型轉換成并聯型子程序function[C,B,A]=dir2par(b,a);M=length(b);N=length(a);[r1,p1,C]=residuez(b,a);p=cplxpair(p1,10000000*eps);I=cplxcomp(p1,p);r=r1(I);K=floor(N/2);B=zeros(K,2);A=zeros(K,3);ifK*2==N;fori=1:2:(N-2)Brow=r(i:1:(i+1),:);Arow=p(i:1:(i+1),:);[Brow,Arow]=residuez(Brow,Arow,[]);B(fix((i+1)/2),:)=real(Brow);A(fix((i+1)/2),:)=real(Arow);end[Brow,Arow]=residuez(r(N-1),p(N-1),[]);B(K,:)=[real(Brow)0];A(K,:)=[real(Arow)0];elsefori=1:2:(N-1)Brow=r(i:1:(i+1),:);Arow=p(i:1:(i+1),:);[Brow,Arow]=residuez(Brow,Arow,[]);B(fix