利用MATLAB仿真軟件系統(tǒng)結(jié)合雙線(xiàn)性變換法設(shè)計(jì)一個(gè)數(shù)字巴特沃斯高通IIR濾波器要點(diǎn)

1、 摘要Matlab是一個(gè)矩陣設(shè)計(jì)平臺(tái)，傳統(tǒng)數(shù)字濾波器設(shè)計(jì)需要大量的計(jì)算，但是利用Matlab可以快速實(shí)現(xiàn)濾波器的設(shè)計(jì)與仿真，而且頻譜分析功能強(qiáng)大，在數(shù)字信號(hào)處理中發(fā)揮了巨大的作用。本次實(shí)驗(yàn)中，用雙線(xiàn)性不變法設(shè)計(jì)高通巴特沃斯IIR數(shù)字濾波器，介紹了設(shè)計(jì)步驟，然后在Matlab環(huán)境下進(jìn)行了仿真與調(diào)試，實(shí)現(xiàn)了設(shè)計(jì)目標(biāo)。關(guān)鍵詞：Matlab 數(shù)字濾波器 雙線(xiàn)性變換法 IIR摘要Matlab是一個(gè)矩陣設(shè)計(jì)平臺(tái)，傳統(tǒng)數(shù)字濾波器設(shè)計(jì)需要大量的計(jì)算，但是利用Matlab可以快速實(shí)現(xiàn)濾波器的設(shè)計(jì)與仿真，而且頻譜分析功能強(qiáng)大，在數(shù)字信號(hào)處理中發(fā)揮了巨大的作用。本次實(shí)驗(yàn)中，用雙線(xiàn)性不變法設(shè)計(jì)高通巴特沃斯IIR數(shù)字濾

2、波器，介紹了設(shè)計(jì)步驟，然后在Matlab環(huán)境下進(jìn)行了仿真與調(diào)試，實(shí)現(xiàn)了設(shè)計(jì)目標(biāo)。關(guān)鍵詞：Matlab 數(shù)字濾波器 雙線(xiàn)性變換法 IIR摘要Matlab是一個(gè)矩陣設(shè)計(jì)平臺(tái)，傳統(tǒng)數(shù)字濾波器設(shè)計(jì)需要大量的計(jì)算，但是利用Matlab可以快速實(shí)現(xiàn)濾波器的設(shè)計(jì)與仿真，而且頻譜分析功能強(qiáng)大，在數(shù)字信號(hào)處理中發(fā)揮了巨大的作用。本次實(shí)驗(yàn)中，用雙線(xiàn)性不變法設(shè)計(jì)高通巴特沃斯IIR數(shù)字濾波器，介紹了設(shè)計(jì)步驟，然后在Matlab環(huán)境下進(jìn)行了仿真。關(guān)鍵詞：Matlab 數(shù)字濾波器 雙線(xiàn)性變換法 1設(shè)計(jì)要求和說(shuō)明利用MATLAB仿真軟件系統(tǒng)結(jié)合雙線(xiàn)性變換法設(shè)計(jì)一個(gè)數(shù)字巴特沃斯高通IIR濾波器。MATLAB工具箱為濾波器的設(shè)

3、計(jì)應(yīng)用提供了豐富而簡(jiǎn)便的方法，使原來(lái)的非常繁瑣復(fù)雜的程序設(shè)計(jì)變成簡(jiǎn)單的程序調(diào)用。1.1 設(shè)計(jì)原理濾波器，顧名思義，就是對(duì)系統(tǒng)輸入信號(hào)進(jìn)行濾波。那個(gè)數(shù)字濾波器的數(shù)學(xué)運(yùn)算通常用兩種方法來(lái)表示。一種是頻域法，即利用FFT快速運(yùn)算辦法對(duì)輸入信號(hào)進(jìn)行離散傅里葉變換，分析其頻譜，然后根據(jù)所希望的頻率特性進(jìn)行濾波，再利用傅里葉反變換來(lái)輸出出時(shí)域信號(hào)。N階低通巴特沃斯濾波器的特性為： 其中，c為通帶寬度，即截止頻率。當(dāng)階次N增大時(shí)，濾波器的特性曲線(xiàn)變的更加陡峭，其特性就越接近理想的矩形幅頻特性。巴特沃斯濾波器屬于全極點(diǎn)設(shè)置，他的極點(diǎn)由以下公式來(lái)確定。式中，。所以，在s平面上有2N個(gè)極點(diǎn)。這些極點(diǎn)是等間隔的分布

4、在半徑為c的圓周上，并且這些極點(diǎn)都是成復(fù)共軛對(duì)出現(xiàn)的，極點(diǎn)位置與虛軸是對(duì)稱(chēng)的，且不在虛軸上。雙線(xiàn)性變換法是將平面壓縮變換到某一中介平面的一條橫帶里，再通過(guò)標(biāo)準(zhǔn)變換關(guān)系將此帶變換到整個(gè)z平面上去，這樣就使平面與平面之間建立一一對(duì)應(yīng)的單值關(guān)系，消除了多值變換性。為了將平面的軸壓縮到平面的軸上的到一段上，可以通過(guò)以下的正切變換來(lái)實(shí)現(xiàn)：這樣當(dāng)由經(jīng)變化到時(shí)，由經(jīng)過(guò)0變化到，也映射到了整個(gè) 軸。將這個(gè)關(guān)系延拓到整個(gè)平面和平面，則可以得到再將平面通過(guò)標(biāo)準(zhǔn)變換關(guān)系映射到平面，即令得到同樣對(duì)z求解，得到雙線(xiàn)性變換平面的映射關(guān)系如圖一：圖一 雙線(xiàn)性變換平面映射圖這樣的變換叫做雙線(xiàn)性變換。1.2 設(shè)計(jì)思路數(shù)字濾波器

5、一般設(shè)計(jì)過(guò)程為： 按照實(shí)際需要，確定濾波器性能要求； 用一個(gè)因果、穩(wěn)定的離散線(xiàn)性時(shí)不變系統(tǒng)去逼近這一個(gè)性能指標(biāo)； 用有限精度的計(jì)算去實(shí)現(xiàn)所設(shè)計(jì)的系統(tǒng)； 通過(guò)模擬驗(yàn)證所設(shè)計(jì)的系統(tǒng)是否符合給定的性能要求。 IIR濾波器設(shè)計(jì)的主要方法是先設(shè)計(jì)低通模擬濾波器，然后轉(zhuǎn)換為高通、帶通或帶阻數(shù)字濾波器。對(duì)于其他如高通，帶通，則通過(guò)頻率變換轉(zhuǎn)換為設(shè)計(jì)相應(yīng)的高通，帶通等。在設(shè)計(jì)的全過(guò)程的各個(gè)步驟，matlab都提供相應(yīng)的工具箱函數(shù)，使得IIR數(shù)字濾波器設(shè)計(jì)變得非常簡(jiǎn)單。總的來(lái)說(shuō)，我的設(shè)計(jì)思路主要有以下兩種：思路一：從歸一化模擬低通原型出發(fā)，先在模擬域內(nèi)經(jīng)頻率變換成為所需類(lèi)型的模擬濾波器；然后進(jìn)行雙線(xiàn)性變換，由S

6、域變換到Z域，而得到所需類(lèi)型的數(shù)字濾波器。歸一化模擬低通原型數(shù)字高,帶通或帶阻模擬高,帶通或帶阻 模擬域 沖激響應(yīng)不變法頻率變換 雙線(xiàn)性變換法圖二 先頻率變換再離散思路二：先進(jìn)行雙線(xiàn)性變換，將模擬低通原型濾波器變換成數(shù)字低通濾波器；然后在Z域內(nèi)經(jīng)數(shù)字頻率變換為所需類(lèi)型的數(shù)字濾波器。 數(shù)字高,帶通或帶阻數(shù)字原型低通歸一化模擬低通原型 數(shù)字域 雙線(xiàn)性變換法 頻率變換圖三 先離散再頻率變換以上兩種思路都可以，我最后選擇了第一種思路進(jìn)行設(shè)計(jì)，即先在模擬域內(nèi)經(jīng)頻率變換成為所需類(lèi)型的模擬濾波器；然后進(jìn)行雙線(xiàn)性變換，由S域變換到Z域，而得到所需類(lèi)型的數(shù)字濾波器。1.3 設(shè)計(jì)過(guò)程下面結(jié)合一個(gè)實(shí)例來(lái)說(shuō)明巴特沃斯

7、濾波器的設(shè)計(jì)過(guò)程。 例：設(shè)計(jì)一個(gè)巴特沃斯濾波器函數(shù)。該濾波器的性能指標(biāo)為：通帶的截止頻率為c=10000rad/s， 通帶的最大衰減Ap=3 dB ，阻帶的截止頻率s=40000rad/s ，阻帶的最小衰減 As=35 dB 。 解： （1）求相關(guān)參數(shù)（2） 確定參數(shù)N取N=3 ， 根據(jù) 得到所以 ，極點(diǎn)形式可以表示為 即滿(mǎn)足系統(tǒng)性能指標(biāo)的函數(shù)2 Matlab軟件的程序設(shè)計(jì)程序設(shè)計(jì)中用到主要函數(shù)：z,p,k=buttap(N); %創(chuàng)建Buttord低通濾波器原型Bap,Aap=zp2tf(z,p,k); %由零極點(diǎn)轉(zhuǎn)換為傳遞函數(shù)的形式Bbs,Abs=lp2hp(Bap,Aap,Wn); %模

8、擬低通變高通Bbz,Abz=bilinear(Bbs,Abs,FS); %雙線(xiàn)性變換2.1 數(shù)字域指標(biāo)變換成模擬域指標(biāo) 其程序?yàn)椋篺p = 400 ; Fs = 300; Rp = 1; Rs = 20; wp =fp*2*pi; ws =fs*2*pi;FS=1000;T=1/FS;程序執(zhí)行結(jié)果為：2.2 數(shù)字域頻率進(jìn)行預(yù)畸變令 Wp=wp/(FS); Ws=ws/(FS); 其畸變程序?yàn)椋簑p2=2*tan(Wp/2)/T; ws2=2*tan(Ws/2)/T; 經(jīng)過(guò)預(yù)畸變，可以發(fā)現(xiàn)頻率變?yōu)?3 模擬濾波器的設(shè)計(jì)其程序?yàn)镹,Wn = buttord(wp2,ws2,Rp,Rs,s) %設(shè)計(jì)

9、模擬濾波器z,p,k=buttap(N); %創(chuàng)建Buttord低通濾波器原型Bap,Aap=zp2tf(z,p,k); %由零極點(diǎn)轉(zhuǎn)換為傳遞函數(shù)的形式figure(1)freqs(Bap,Aap);%模擬低通濾波器的頻率響應(yīng)title(模擬濾波器（低通原型）的頻率響應(yīng))Bbs,Abs=lp2hp(Bap,Aap,Wn); %模擬低通變高通figure(2)freqs(Bbs,Abs);title(模擬濾波器(高通原型)的頻率響應(yīng))程序執(zhí)行后可以發(fā)現(xiàn)其頻率響應(yīng)為：模擬濾波器的低通原型波形如圖四所示：圖四 模擬濾波器（低通原型）的頻率響應(yīng) 模擬濾波器高通原型的波形如圖五所示：圖五 模擬濾波器（高

10、通原型）的頻率響應(yīng)2.4 模擬濾波器變成數(shù)字濾波器其程序?yàn)椋築bz,Abz=bilinear(Bbs,Abs,FS); %用雙線(xiàn)性變換法設(shè)計(jì)數(shù)字濾波器 freqz(Bbz,Abz,512,FS); 程序運(yùn)行的結(jié)果如圖六：圖六 數(shù)字濾波器的頻率響應(yīng)由于使用的是雙線(xiàn)性不變法設(shè)計(jì)的，其相位為非線(xiàn)性。此處主要是基于要獲得嚴(yán)格的頻率響應(yīng)，以及較準(zhǔn)確地控制截止頻率的位置，故畫(huà)出了詳細(xì)的幅頻響應(yīng)。得到數(shù)字濾波器的詳細(xì)頻率響應(yīng)圖形圖七 詳細(xì)的幅頻響應(yīng)（通帶衰減為5dB）當(dāng)通帶變?yōu)?dB（之前為1dB）時(shí)，阻帶不變時(shí)，其波形如下圖十。對(duì)比圖七可知，其在通帶處的衰減變?yōu)榱?dB，曲線(xiàn)平滑了一些。從調(diào)試的過(guò)程中發(fā)現(xiàn)

11、：通帶衰減越小，濾波器的性能越好。 阻帶衰減越大，濾波器的性能越好。3.程序調(diào)試中出現(xiàn)的問(wèn)題 freqs(Bbs,Abs)與freqz(Bbz,Abz,512,FS)剛開(kāi)始時(shí)分得不是很清楚，后來(lái)發(fā)現(xiàn)前者是求模擬域頻率響應(yīng)的，而后者是求數(shù)字域頻率響應(yīng)的。后者的用法如下：h,w=freqz(b,a,n,whole);或h,f=freqz(b,a,n,whole,Fs);式中，b,a為數(shù)字濾波器分子和分母多項(xiàng)式的系數(shù)，n為復(fù)數(shù)頻率的響應(yīng)點(diǎn)數(shù)，為整數(shù)，最好為2的冪，缺省時(shí)為512；Fs為采樣頻率，單位Hz。如果給定該值，則f位置輸出為頻率Hz，若沒(méi)有給定，則按角頻率（Angular frequency

12、）給定f的頻率矢量；whole表示返回的頻率f或w值包含z平面整個(gè)單位圓頻率矢量，即02；缺省時(shí)，頻率f或w值包含z平面上半單位圓（0）之間等間距n個(gè)點(diǎn)頻率矢量。h為復(fù)頻率響應(yīng)；w為n點(diǎn)頻率向量（單位rad）；f為n點(diǎn)頻率向量(Hz)。函數(shù)返回值缺省時(shí)，繪制幅頻響應(yīng)和相頻響應(yīng)圖4小結(jié) 總的來(lái)說(shuō)，這次課程設(shè)計(jì)讓我對(duì) MATLAB 有了更深刻的了解，對(duì)數(shù)字濾波器的設(shè)計(jì)流程有了大致的了解，掌握了一些設(shè)計(jì)濾波器的基本方法，提高了理論用于實(shí)踐的能力，掌握了更多專(zhuān)業(yè)相關(guān)的使用知識(shí)與技能。同時(shí)，也暴露了我很多的不足，在以后的學(xué)習(xí)中， 將進(jìn)一步發(fā)現(xiàn)并克服缺點(diǎn)。參考文獻(xiàn)1 陳懷琛 MATLAB在電子信息課程中的

13、應(yīng)用（第二版） 電子工業(yè)出版社 20062 郭仕劍 MATLAB7.X數(shù)字信號(hào)處理 人民郵電出版社 20063 陳貴明 用MATLAB語(yǔ)言處理數(shù)字信號(hào)與數(shù)字圖像 科學(xué)出版社 20004 王家文 MATLAB7.0圖形圖像處理 國(guó)防工業(yè)出版社 2006附錄：實(shí)驗(yàn)總程序：fp = 400 ; fs = 300; Rp = 1; %通帶最大衰減Rp=1dBRs = 20; %阻帶最小衰減Rs=20dBwp =fp*2*pi; ws =fs*2*pi;FS=1000;T=1/FSWp=wp/(FS);Ws=ws/(FS);wp2=2*tan(Wp/2)/T; %預(yù)畸變求濾波器通帶臨界頻率ws2=2*t

14、an(Ws/2)/T; %預(yù)畸變求濾波器阻帶臨界頻率 N,Wn = buttord(wp2,ws2,Rp,Rs,s) z,p,k=buttap(N); %創(chuàng)建Buttord低通濾波器原型Bap,Aap=zp2tf(z,p,k); %由零極點(diǎn)轉(zhuǎn)換為傳遞函數(shù)的形式figure(1)freqs(Bap,Aap); %模擬低通濾波器的頻率響應(yīng)title(模擬濾波器（低通原型）的頻率響應(yīng))Bbs,Abs=lp2hp(Bap,Aap,Wn); %模擬低通變高通figure(2)freqs(Bbs,Abs);title(模擬濾波器的頻率響應(yīng)) Bbz,Abz=bilinear(Bbs,Abs,FS); %雙

15、線(xiàn)性變換figure(3)freqz(Bbz,Abz,512,FS);title(數(shù)字濾波器的頻率響應(yīng))hw,w=freqz(Bbz,Abz,512); figure(4)plot(w/pi,20*log10(abs(hw);gridaxis(0,1,-200,10)title( 數(shù)字濾波器詳細(xì)的頻率響應(yīng))xlabel(w/pi);ylabel(幅度(dB);wp=0.8*pi;ws=0.6*pi;OmegaP=2*1000*tan(wp/2);OmegaS=2*1000*tan(ws/2);lamdas=OmegaP/OmegaS;N=0.5*log10(10.(20/10)-1)/(10.(1/10)-1)/log10(lamdas);Wn= 4.8890e+003az=0 0 0 0 1;bz=1 2.613 3.414,2.613,1;Bbs,Abs=lp2hp(az,bz,Wn) Bbz,Abz=bilinear(Bbs,Abs,1000); hw,w=freqz(Bbz,Abz,512); figure(5)plot(w/pi,20*log10(abs(hw);gridaxis(0,1,-200,10)set(gca,XTickMode,manual,XTick,0,Ws/(pi),Wp/