姓名：班級：學號：時間：2011年6月設計題目基于MATLAB的巴特沃斯低通濾波器的設計設計要求.通過實驗加深對巴特沃斯低通濾波器基本原理的理解。.學習編寫巴特沃斯低通濾波器的MATLAB仿真程序.濾波器的性能指標如下：通帶截止頻率fp=5kHz,通帶最大衰減a=2dB,阻帶截止頻率fs=12kHz，阻帶最小衰減a=30dBp s

1.設計原理巴特沃斯低通濾波器簡介：巴特沃斯濾波器是電子濾波器的一種，特點是通頻帶內的頻率響應曲線最大限度平坦，沒有起伏，而在阻頻帶則逐漸下降為零。這種濾波器最先由英國工程師斯替芬?巴特沃斯(StephenButterworth)在1930年發表在英國《無線電工程》期刊的一篇論文中提出的，可以構成低通、高通、帶通和帶阻四種組態，是目前最為流行的一類數字濾波器，經過離散化可以作為數字巴特沃思濾波器，較模擬濾波器具有精度高、穩定、靈活、不要求阻抗匹配等眾多優點，因而在自動控制、語音、圖像、通信、雷達等眾多領域得到了廣泛的應用，是一種具有最大平坦幅度響應的低通濾波器。巴特沃斯低通濾波器的設計原理：巴特沃斯低通濾波器的幅度平方函數Ha(j。)2用下式表示：、c 1H(j。)2=-5—a 1+(——)2N。c設計過程其中N為濾波器的階數。當。=0時，巴(jQ)=1；。=。c時，Ha(jQ)|=1八，*。c是3dB截止頻率。。=%時，。逐漸增大，幅度下降非常迅速。。、N同幅度特性關系如圖1.1所示。N決定了幅度下降速度，N越大，通帶就越平坦，過渡帶也隨之變窄，阻帶幅度同過渡帶下降的速度越迅速，總體頻響特性同理想低通濾波器的實際誤差越小。設計過程圖1.1。、N同幅度特性關系用s代替j。，把幅度平方函數|Ha(j。)|2變成s的函數：、 、 1H(s)H(-s)= aa1+(—)2Nj。cs=o+j。，此公式說明了幅度平方函數有2N個極點，極點s可以用下面的公式來表達： kj8(2k+1) .8,2k+1 ./12k+1、s=(j。e2n)=。eJ2eJK2n=。eJK2+2nk c c ck=0,1,2,…，2N-1。2N個極點等間隔分布在半徑為Q的圓上，間隔是兀/Nrad。c如圖1.2所示：圖1.2三階巴特沃斯濾波器極點分布為形成穩定的濾波器，2N個極點中只取s平面左半平面的N個極點構成H（s），而右半平面的N個極點構成H（-s）。H（s）的表示式為aH（s）=—QN—

a H（s-s）kk=0例如N=3,通過下式可以計算出6個極點.2s=QeJ3兀0c當N=3時，.2s=.2s=QeJ3兀0c當N=3時，.2s=QeJ3兀0c.4 .5s=Qej兀，s=QeJ3K，s=Qe%兀，s=Qej2兀1c 2c 3c 4c6個極點中位于左半平面的三個分別為：7s=。ej3k5c取s平面左半平面的極點s，0H（s）=—a，cs1.4s=QeJ3兀2cs組成H（s）：2 aQ3

c

.2.2（s+Q）（s-Qj3兀）（s-Q-j3兀）

c c c一 1將H（s）= 對3dB截止頻率aH（J"）QQ

k=0 cc .、 1H（s）表示為：H（s）= a a H（J"）QQ

k=0ccQc歸一化后的丑=j九，p稱為歸一化拉氏復變量。X=Q/Q，Q cc九稱為歸一化頻率。經過歸一化后巴特沃斯濾波器的傳輸函數為:、 1Ha(p)= H(p-p)kK=0式中，p=s/。為歸一化極點，kkc(1.1)s為位于左半平面的極點用下式表示:k12k+(1.1)s為位于左半平面的極點用下式表示:kp=ej2+2n，k=0,1，..N-1k將極點表示式代入(1.1)式，得到的H(p)的分母是p的N階多項式，用下aH(p)=

a0下面來確定N：由技術指標ap在定義b+bp+bp2d fbpH(p)=

a0下面來確定N：由技術指標ap在定義b+bp+bp2d fbpN-1+pNN-1Qa和Qs確定。=-10lgH(jQ)as=-10lg|H(jQs)、， 1Ng2=L1+(—)2NQc(1.2)(1.3)(1.4)中，將Q=Qp和Q=Q分別代入(1.4)式中，得到H(jQ)2和H(jQs)再將H(jQ)2和H(jQs)2代入(1.2)和(1.3)式中，得到:apI1a?a=-10lg[pa=-10lg[1Q1+(—Q1——],p)2N整理得:Q1+(一)Qc一q、一1+(_p_)2N=10ac1+(,)2N=10as/10c由(1.5)和(1.6)式得到：Q ,'10a ,1(4)N=J*"。1Q\110as/10-1s令p/10(1.5)(1.6)- ：1S/10-1(1.7)九=Q/Q，k=Ip(1.7)spspsp110/10—1則N由下式表示:

zlgk(1.8)N二——.sp(1.8)lg.sp取大于等于N的最小整數。關于3dB截止頻率Q式求出，關于3dB截止頻率Q式求出，Q1+(p)2N=10ap/10Qc1+( )2N=10a/10Qc由(1.9)式得到：，如果技術指標中沒有給出，可以按照(1.7)式或(1.8)(1.9)(1.10)Q=Q(100.1ap_A_-1)-2N由(1.10)式得到:1Q=Q(100.1as-1)-2Ncs2設計方案方案一：用沖激響應不變法設計巴特沃斯低通濾波器沖激響應不變法是使數字濾波器的單位沖激序列h(n)模仿模擬濾波器的單位沖激響應ha(t)。將模擬濾波器的單位沖激響應加以等間隔抽樣，使h(n)正好等于ha(t)的抽樣值，即滿足：h(n)=ha(nT)其中T是抽樣周期。如果令Ha(s)是ha⑴的拉普拉斯變換，H(z)為h(n)的z變換，利用抽樣序列的z變換與模擬信號的拉普拉斯變換的關系，得：TOC\o"1-5"\h\z\o"CurrentDocument"V 1V( 2兀、X(z) =-工X(s-jkQ)=-工Xs-j-kz=esTT a sT aI Tk=-^ k=-^可看出，脈沖響應不變法將模擬濾波器的S平面變換成數字濾波器的Z平面，這個從s到z的變換z=eST是從S平面變換到Z平面的標準變換關系式。沖激響應不變法使得數字濾波器的單位脈沖響應完全模仿模擬濾波器的單位沖激響應，也就是時域逼近良好，而且模擬頻率Q和數字頻率①之間呈線性關系①二QT。但是此方法有一個明顯的缺點就是有頻率響應的混疊效應，所以，脈沖響應不變法只適用于限帶的模擬濾波器(例如，衰減特性很好的低通或帶通濾波器)，而且高頻衰減越快，混疊效應越小。方案二：雙線性變換法設計IIR數字濾波器雙線性變換法是使數字濾波器的頻率響應與模擬濾波器的頻率響應相似的一種變換方法。為了克服多值映射的缺點，采用把整個s平面頻率壓縮方法，將整個頻率軸上的頻率范圍壓縮到n/T?n/T之間，再用z=est轉換到Z平面上。也就是說，第一步先將整個S平面壓縮映射到S1平面的-n/T?n/T一條橫帶里；第二步再通過標準變換關系z=es1t將此橫帶變換到整個Z平面上去。這樣就使S平面與Z平面建立了一一對應的單值關系，消除了多值變換性，也就消除了頻譜混疊現象，映射關系如圖1.3所示。

圖1.3雙線性變換的映射關系為了將S平面的整個虛軸j。壓縮到S1平面j。1軸上的-n/T到n/T段上，可以通過以下的正切變換實現八2T是采樣間隔。O=tanT是采樣間隔。當。1由-n/T經過0變化到n/T時，。由-s經過0變化到+處也即映射了整個j。軸。將式（1-5）寫成e41T/2將此關系解析延拓到整個S平面和S1平面，令j。=s，j。1=s1,則得esT/2一e-sT/2 2 JsT）21一e-ys二一— 1一二—tanh一二一 Tes^T/2+e-s；T/2TI2）T1+e-s1T再將S1平面通過以下標準變換關系映射到Z平面z=es1T從而得到S平面和Z平面的單值映射關系為：21-z-1s- TOC\o"1-5"\h\zT1+Z-1 （2-1）tT 2——s—sT （2-2）式（2-1）與式（2-2）是S平面與Z平面之間的單值映射關系，這種變換都是兩個線性函數之比，因此稱為雙線性變換。雙線性變換法與沖激響應不變法相比，其主要的優點是避免了頻率響應的混疊現象，雖然在線性方面有些欠缺，但是可以通過頻率的預畸來加以校正且計算比沖激響應不變法方便，實現起來比較容易，所以，本設計選擇用雙線性變換法設計巴特沃斯低通濾波器。3設計步驟由于沖擊響應不變法的不足，下面以雙線性變換法設計一個巴特沃思低通濾波器。通過仔細研究MATLAB軟件自帶的信號處理工具箱下％“八五八8口IR%\toolbox\signal\signal\buttord.m和％MATLABDIR%\toolbox\signal\signal\butter.m，可以清楚看到:MATLAB默認采用雙線性變換法設計，其實完全可以采用更簡潔的方法。兇濾波器的性能指標如下：通帶截止頻率fp=5kHz，通帶最大衰減。廣2dB，

阻帶截止頻率fs=12kHz,阻帶最小衰減a130dB3.1MATLAB中所需函數ATLAB的信號處理工具箱提供了濾波器的函數buttap、buttord、butter。H($)=由[z,p,k]=buttap(n)函數可設計出nH($)=(S-p(1))(S-p(2))A(s-p(n))kSn+bSnkSn+bSn-1+A+bs+bn-i 1 0=1rad/s，得到巴特沃斯濾波器歸一化結果，如表1所示。H($)=其中bnb7b6b5b4b3b2b1b011.000021.41421.000032.00002.00001.000042.61313.41422.61311.000053.23615.23615.23613.23611.000063.86377.46419.14167.46413.86371.000074.494010.097814.591814.591810.09784.49401.000085.125813.137121,846225.846221,846213.13715.12581.00000=3n，令3c c表1n=1~8階的巴特沃斯濾波器系數buttord函數可在給定濾波器性能的情況下，選巴特沃斯濾波器的階數n和截止頻率3，從而可用butter函數設計巴特沃斯濾波器的傳遞函數。⑼c[n,3]=buttord(3p,3,Rp,Rs對)可得到足性能的模擬巴特沃斯濾波器c c的最小階數n及截止頻率3，其中3P為通帶的拐角頻率，3s為阻帶的拐角頻率，3P和3s的單位均為rad/s;Rs為通帶區的最大波動系數,Rp為Rs阻帶區的最小衰減系數，Rs和Rp的單位都為dB。[b,a]=butter(n,3,‘s’)可設計截止頻率為3的n階低通模擬巴特沃斯濾波c c器，其傳遞函數為:[10]〃/、 B(s) b(1)sn+b(2)sn-1+A+b(n+1)H(s)= = A(s) sn+a(2)sn-1+A+a(n+1)3.2具體設計步驟經過總結，巴特沃斯低通濾波器的設計步驟大致為:(1)通過Qp,ap(2)根據公式p(1)通過Qp,ap(2)根據公式pks lg九 0sp.，2k+1、=￡-2+2N)，k=0,1，…N-1，求出歸一化極點p，將p4.14.1圖形結果顯示:..、 i代入H(p)= 中，得出歸一化傳輸函數H(p)。aH(p-p) akk=0(3)將H(p)去歸一化。將p=s/Qc代入H(p)之中，從而得到實際的濾波a器傳輸函數H(s)。a3.3MATLAB程序MATLAB程序如下：(1).輸入信號的時域波形和頻譜的設計：closeall；N=256;t=linspace(0,1,N);dt=t(2)-t(1);xt=cos(2*pi*4*t)+cos(2*pi*10*t)+cos(2*pi*20*t);f=(0:(N/2-1))/(dt*N);Xt=fft(xt,N);subplot(2,1,1),plot(t(1:128),xt(1:128))；subplot(2,1,2),plot(f(1:64),abs(Xt(1:64))),xlabel('f(kHz)')(2)輸出的巴特沃斯低通濾波器的波特圖的設計Qc=5.2775;b0=1;b1=3.2361;b2=5.2361;b3=b2;b4=b1;Q=f%Q=linspace(0,25,N);Ha=QcA5./(0*Q).A5+b4*Qc*(j*Q).A4+b3*QcA2*0*Q).A3+b2*QcA3*(j*Q).A2+b1*QcA4*(j*Q)+b0*QcA5);L=length(Ha)Has=20*log10(abs(Ha));figure,plot(Q(1:64),Has(1:64),Q,-30,'r*',12,Has,'*',5,Has,'*'),axis([030-702]),xlabel('f(kHz)'),ylabel('201g(abs(H_{a}(j{\Omega})))(dB)');(3)巴特沃斯低通濾波后的時域波形和頻譜的設計Yt=Xt(1:L).*Ha;yt=ifft(Yt);figure,subp1ot(2,1,1),p1ot(t(1:128),abs(yt(1:128))),subp1ot(2,1,2),p1ot(f(1:64),abs(Yt(1:64)));x1abe1('f(kHz)')4.運行結果及分析4.運行結果及分析-40 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5150100-■??-40 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5150100-■??50_00 10 20 30 40 50 60 70f(kHz)圖4.1輸入信號的時域波形和頻譜-700 5 10 15 20 25 30f(kHz)圖4.2輸出的巴特沃斯低通濾波器的波特圖0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5150100500 0 10 20 30 40 50 60 70f(kHz)圖4.3進行巴特沃斯低通濾波后的時域波形和頻譜4.2結果分析(1)由圖4.1與圖4.3的比較可知：經過巴特沃斯低通濾波器后，當信號頻率大于10Hz時可以被濾除，起到良好的低通濾波效果，而且可以將輸入信號的波形經過濾波后的波形比原波形的振幅有所減小，去除了很多由于噪聲或其他因素所產生的