基于MATLAB的巴特沃斯低通濾波器的設(shè)計

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

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

a H（s-s）kk=0例如N=3,通過下式可以計算出6個極點(diǎn).2s=QeJ3兀0c當(dāng)N=3時，.2s=.2s=QeJ3兀0c當(dāng)N=3時，.2s=QeJ3兀0c.4 .5s=Qej兀，s=QeJ3K，s=Qe%兀，s=Qej2兀1c 2c 3c 4c6個極點(diǎn)中位于左半平面的三個分別為：7s=。ej3k5c取s平面左半平面的極點(diǎn)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稱為歸一化拉氏復(fù)變量。X=Q/Q，Q cc九稱為歸一化頻率。經(jīng)過歸一化后巴特沃斯濾波器的傳輸函數(shù)為:、 1Ha(p)= H(p-p)kK=0式中，p=s/。為歸一化極點(diǎn)，kkc(1.1)s為位于左半平面的極點(diǎn)用下式表示:k12k+(1.1)s為位于左半平面的極點(diǎn)用下式表示:kp=ej2+2n，k=0,1，..N-1k將極點(diǎn)表示式代入(1.1)式，得到的H(p)的分母是p的N階多項式，用下aH(p)=

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

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

圖1.3雙線性變換的映射關(guān)系為了將S平面的整個虛軸j。壓縮到S1平面j。1軸上的-n/T到n/T段上，可以通過以下的正切變換實(shí)現(xiàn)八2T是采樣間隔。O=tanT是采樣間隔。當(dāng)。1由-n/T經(jīng)過0變化到n/T時，。由-s經(jīng)過0變化到+處也即映射了整個j。軸。將式（1-5）寫成e41T/2將此關(guān)系解析延拓到整個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平面通過以下標(biāo)準(zhǔn)變換關(guān)系映射到Z平面z=es1T從而得到S平面和Z平面的單值映射關(guān)系為：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平面之間的單值映射關(guān)系，這種變換都是兩個線性函數(shù)之比，因此稱為雙線性變換。雙線性變換法與沖激響應(yīng)不變法相比，其主要的優(yōu)點(diǎn)是避免了頻率響應(yīng)的混疊現(xiàn)象，雖然在線性方面有些欠缺，但是可以通過頻率的預(yù)畸來加以校正且計算比沖激響應(yīng)不變法方便，實(shí)現(xiàn)起來比較容易，所以，本設(shè)計選擇用雙線性變換法設(shè)計巴特沃斯低通濾波器。3設(shè)計步驟由于沖擊響應(yīng)不變法的不足，下面以雙線性變換法設(shè)計一個巴特沃思低通濾波器。通過仔細(xì)研究MATLAB軟件自帶的信號處理工具箱下％“八五八8口IR%\toolbox\signal\signal\buttord.m和％MATLABDIR%\toolbox\signal\signal\butter.m，可以清楚看到:MATLAB默認(rèn)采用雙線性變換法設(shè)計，其實(shí)完全可以采用更簡潔的方法。兇濾波器的性能指標(biāo)如下：通帶截止頻率fp=5kHz，通帶最大衰減。廣2dB，

阻帶截止頻率fs=12kHz,阻帶最小衰減a130dB3.1MATLAB中所需函數(shù)ATLAB的信號處理工具箱提供了濾波器的函數(shù)buttap、buttord、butter。H($)=由[z,p,k]=buttap(n)函數(shù)可設(shè)計出nH($)=(S-p(1))(S-p(2))A(s-p(n))kSn+bSnkSn+bSn-1+A+bs+bn-i 1 0=1rad/s，得到巴特沃斯濾波器歸一化結(jié)果，如表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階的巴特沃斯濾波器系數(shù)buttord函數(shù)可在給定濾波器性能的情況下，選巴特沃斯濾波器的階數(shù)n和截止頻率3，從而可用butter函數(shù)設(shè)計巴特沃斯濾波器的傳遞函數(shù)。⑼c[n,3]=buttord(3p,3,Rp,Rs對)可得到足性能的模擬巴特沃斯濾波器c c的最小階數(shù)n及截止頻率3，其中3P為通帶的拐角頻率，3s為阻帶的拐角頻率，3P和3s的單位均為rad/s;Rs為通帶區(qū)的最大波動系數(shù),Rp為Rs阻帶區(qū)的最小衰減系數(shù)，Rs和Rp的單位都為dB。[b,a]=butter(n,3,‘s’)可設(shè)計截止頻率為3的n階低通模擬巴特沃斯濾波c c器，其傳遞函數(shù)為:[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具體設(shè)計步驟經(jīng)過總結(jié)，巴特沃斯低通濾波器的設(shè)計步驟大致為:(1)通過Qp,ap(2)根據(jù)公式p(1)通過Qp,ap(2)根據(jù)公式pks lg九 0sp.，2k+1、=￡-2+2N)，k=0,1，…N-1，求出歸一化極點(diǎn)p，將p4.14.1圖形結(jié)果顯示:..、 i代入H(p)= 中，得出歸一化傳輸函數(shù)H(p)。aH(p-p) akk=0(3)將H(p)去歸一化。將p=s/Qc代入H(p)之中，從而得到實(shí)際的濾波a器傳輸函數(shù)H(s)。a3.3MATLAB程序MATLAB程序如下：(1).輸入信號的時域波形和頻譜的設(shè)計：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)輸出的巴特沃斯低通濾波器的波特圖的設(shè)計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)巴特沃斯低通濾波后的時域波形和頻譜的設(shè)計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.運(yùn)行結(jié)果及分析4.運(yùn)行結(jié)果及分析-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進(jìn)行巴特沃斯低通濾波后的時域波形和頻譜4.2結(jié)果分析(1)由圖4.1與圖4.3的比較可知：經(jīng)過巴特沃斯低通濾波器后，當(dāng)信號頻率大于10Hz時可以被濾除，起到良好的低通濾波效果，而且可以將輸入信號的波形經(jīng)過濾波后的波形比原波形的振幅有所減小，去除了很多由于噪聲或其他因素所產(chǎn)生的