沖激響應濾波器課程設計

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

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

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

￡2( -)2NjQcX2進而可得N為正整數(shù)的取值要求：N>Xlog()10￡式(3.13)logjQ當給定的參數(shù)指標是歸一化的，即Q=1,￡=1時，上式表示為:log(九)N> to-log(0)10s⑵從幅度平方函數(shù)確定系統(tǒng)函數(shù)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+為奇數(shù)2 式(314)s,=-sin( 兀)+jcos( 兀)k={123 N 式(3?14丿pk 2N 2N 1，，2…3 為偶數(shù) N令歸一化的系統(tǒng)函數(shù)用H(s令歸一化的系統(tǒng)函數(shù)用H(s)，則在左半平面上的極點組成H(s)，H(s)aN和H(s)的關(guān)系為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設計巴特沃斯濾波器的實例。一、設計內(nèi)容設計一個數(shù)字巴特沃斯低通濾波器，設計指標如下：① =① =0.2兀pR=1dBpA=15dB

s米樣時間間隔T=1S。二、設計要求1)用單位沖激響應不變變換法進行設計。（2）給出詳細的濾波器設計說明書。（3）給出經(jīng)過運行是正確的程序清單并加上詳細的注釋。（4）畫出所設計濾波器的幅度特性和相位特性。三、設計說明（1） 把數(shù)字頻率轉(zhuǎn)換為模擬頻率：Q=WT，0=wT。ppss（2） 計算巴特沃斯模擬濾波器的截止頻率0和階數(shù)N。c（3） 設計巴特沃斯模擬低通濾波器，給出參數(shù)b和a（此處使用了MATLAB中的buttap（N）函數(shù)。（4） 把模擬濾波器用單位沖激響應不變變換法轉(zhuǎn)換成數(shù)字濾波器（此處使用了MATLAB中的residuez函數(shù)）。（5） 變直接形式為并聯(lián)形式，并給出結(jié)構(gòu)圖。（6）畫出幅度特性和相位特性。四、設計步驟框圖數(shù)字濾波器的設計步驟如圖1所示參數(shù)變換巴特沃斯模擬低通濾波器巴特沃斯模擬濾波器離散化處理直接型轉(zhuǎn)變成為并聯(lián)型五、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')%計算巴特沃斯數(shù)字濾波器的階數(shù)N和截止頻率wc[B,A]=butter(N,wc,'s')[z,p,k]=buttap(N)%模擬低通原型零、極點系數(shù)和增益因子[bz,az]=impinvar(b,a)%用沖激響應不變法將模擬濾波器轉(zhuǎn)化為數(shù)字濾波器，采樣頻率默認1Hzwz=[0:pi/512:pi];hz1=freqz(bz,az,wz);%巴特沃斯模擬低通濾波器頻率響應[C,B,A]=dir2par(b,a)%直接型轉(zhuǎn)換成并聯(lián)型%繪圖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;%直接型轉(zhuǎn)換成并聯(lián)型子程序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