1、基于mat lab 的低通濾波器摘要:調(diào)用MA TLAB 信號(hào)處理工具箱中濾波通過觀察濾波器輸入輸出信號(hào)的時(shí)域波形及其頻譜，建立數(shù)字濾波的概念。應(yīng)用最廣泛的是雙線性變換法?；驹O(shè)計(jì)過程是：先將給定的數(shù)字濾波器的指標(biāo)轉(zhuǎn)換成過渡模擬濾波器的指標(biāo)； 設(shè)計(jì)過渡模擬濾波器；將過渡模擬濾波器系統(tǒng)函數(shù)轉(zhuǎn)換成數(shù)字濾波器的系統(tǒng)函數(shù)。MA TLAB 信號(hào)處理工具箱中的各種IIR 數(shù)字濾波器設(shè)計(jì)函數(shù)都是采用雙線性變換法。關(guān)鍵詞:濾波器，matlab ，c 語言，聲音Abstract : call MATLAB signal processing toolbox filtering through the obser

2、vation filter of input and output signalstime domain waveform and spectrum, establish the concept of digital filter. One of the most widely applied is double linear transformation method. The basic design process is: (1 to a given digital filters index converted into analog filter transition index;

3、(2 the design transition simulation filter; (3 transition simulation filter system function will be converted into digital filter system function. MATLAB signal processing toolbox digital filter function design of IIR is bilinear transformation method.Keywords: filter, matlab, the c language, the vo

4、ice一. 任務(wù)：用matlab 軟件設(shè)計(jì)IIR 模擬、數(shù)字以及各種窗函數(shù)的FIR 低通濾波器二. 設(shè)計(jì)目的：（1）了解matlab 軟件的用途以及用法；（2）了解用沖激響應(yīng)不變法設(shè)計(jì)模擬低通濾波器；（3）了解用脈沖響應(yīng)不變法設(shè)計(jì)的巴特沃思數(shù)字低通濾波器； （4）了解基于漢明窗函數(shù)的FIR 低通濾波器的設(shè)計(jì)。三. 設(shè)計(jì)內(nèi)容：低通濾波器的技術(shù)要求用圖形表示如下： 1 通帶 過渡帶 阻帶如上圖表示了一個(gè)頻域設(shè)計(jì)、一維低通濾波器的技術(shù)要求圖。ws 和wp 分別稱為通帶截止頻率和阻帶截止頻率。通帶頻率范圍為0wwp，阻帶頻率范圍wsw。從wp 到ws 稱為過渡帶，用w表示，在過渡帶里，幅頻特性單調(diào)下降

5、。在通帶和阻帶內(nèi)的衰減一般用數(shù) dB 表示，通帶內(nèi)允許最大衰減是s，阻帶內(nèi)允許最小衰減是p。定義分別為：Matlab信號(hào)處理工具箱中提供了設(shè)計(jì)巴特沃思模擬濾波器的函數(shù)buttord 、 buttap 和butter ，格式如下：, (, , , , C P S P S N W buttord W W R R s =用于計(jì)算巴特沃思模擬低通濾波器的階N 和3dB 截止頻率Wc (即本書中的符號(hào)c 。其中，Wp 和Ws 分別是濾波器的通帶截止頻率p 和阻止截止頻率s ，單位為rad/s；Rp 和Rs 分別是通帶最大衰減系數(shù)p 和阻帶最小衰減系數(shù)s ，單位為dB 。, , ( z p G butta

6、p N =用于計(jì)算N 階巴特沃思?xì)w一化（c =1）模擬低通濾波器系統(tǒng)函數(shù)的零、極點(diǎn)和增益因子，返回長度為N 的向量z 和p 分別給出N 個(gè)零點(diǎn)和極點(diǎn)，G 是濾波器增益。得到的濾波器系統(tǒng)函數(shù)形式如下：1212( ( ( ( ( ( ( a N a a N Q s s z s z s z H s GP s s p s p s p -=-其中，k z 和k p 分別是向量z 和p 的第k 個(gè)元素。如果要從零、極點(diǎn)得到系統(tǒng)函數(shù)的分子和分母多項(xiàng)式系數(shù)向量B 和A ，可以調(diào)用結(jié)構(gòu)轉(zhuǎn)換函數(shù)：, 2(, , B A zp tf z p G = 結(jié)構(gòu)轉(zhuǎn)換后系統(tǒng)函數(shù)的形式為：111111( ( ( M M Ma

7、N N N（1）巴特沃思模擬低通濾波器設(shè)計(jì)指標(biāo)如下： wp=0.4*pi, ws=0.45*pi, Rp=5,Rs=10; 實(shí)驗(yàn)所用到的程序如下p j0p j (e 20lg (e H H =sj0s j (e20lg (e H H =wp=0.4*pi;ws=0.45*pi;Rp=5;Rs=10; N,wc=buttord(wp,ws,Rp,Rs,s ; B,A=butter(N,wc,s ; w=0:0.001:6;H,W=freqs(B,A,w; H=20*log10(abs(H; plot(w,H,grid on ;xlabel(頻率/(rad/s;ylabel( 幅度/dB經(jīng)matl

8、ab 仿真以后得到的濾波器如下圖所示： 圖1.1 巴特沃思模擬低通濾波器運(yùn)行后在matlab 控制窗口上顯示如下結(jié)果： （2）為了方便觀察故設(shè)計(jì)了不同階次N 的濾波器，即Butterworth 模擬原型低通濾波器，其程序代碼如下：n=0:0.01:2; for ii=1:4 switch iicase 1,N=2; case 2,N=5; case 3,N=10; case 4,N=20;end z,p,k=buttap(N; b,a=zp2tf(z,p,k; H,w=freqs(b,a,n; magH2=(abs(H.2; hold on ;0123456-60-50-40-30-20-10

9、010頻率/(rad/s幅度/d B巴特沃什模擬低通濾波器plot(w,magH2; endxlabel(w/wc; ylabel(|H(jw|2;title(Butterworth 模擬原型低通濾波器 ; text(1.5,0.18,n=2 text(1.3,0.08,n=5 text(1.16,0.08,n=10 text(0.93,0.98,n=20 grid on ; 仿真后結(jié)果圖形如圖：圖1.2 Butterworth模擬原型低通濾波器（3）利用沖激響應(yīng)不變法，把系統(tǒng)函數(shù)為651(2+=s s s s H a 的模擬濾波器變換成等價(jià)的數(shù)字濾波器，采樣間隔T=1S。先經(jīng)過計(jì)算得到：21

10、32651 (2+=+=s s s s s s H a ，則有如下式子 21112130067. 01851. 012209. 00000. 1112 (-+-=-=zz z z e T z e T z H T T Matlab 實(shí)驗(yàn)參考程序如下:B=1,1; A=1,5,6; T=1; Fs=1/T;Bz,Az=impinvar(B,A,Fs; %用沖激響應(yīng)不變法將模擬濾波器變換成數(shù)字濾波器0.20.40.60.800.20.40.60.811.21.4w/wc|H (j w |2Butterworth 模擬原型濾波器運(yùn)行結(jié)果：Bz =1.0000 -0.2209Az = 1.0000 -0

11、.1851 0.0067得到的數(shù)據(jù)和先前計(jì)算得到的數(shù)據(jù)相同，故利用沖激響應(yīng)不變化法把模擬濾波器設(shè)計(jì)成數(shù)字濾波器正確。（1）步驟1 中產(chǎn)生的低通模擬濾波器的系統(tǒng)函數(shù)，經(jīng)計(jì)算為：經(jīng)過仿真后的結(jié)果，與上式對(duì)比，所得結(jié)果進(jìn)分析相同。此處為N=7的低通模擬濾波器（2）在Butterworth 模擬原型低通濾波器中，可以看到，濾波器的幅頻特性隨著濾波器階次N 的增加而變得越來越好。（3）步驟3中，利用沖激響應(yīng)不變法，把系統(tǒng)函數(shù)的模擬濾波器變換成等價(jià)的數(shù)字濾波器，在matlab 仿真軟件中驗(yàn)證了沖激響應(yīng)不變法的正確性。用MA TLAB 進(jìn)行典型的數(shù)字濾波器的設(shè)計(jì)，一般步驟如下:（1）將設(shè)計(jì)指標(biāo)歸一化處理。如

12、果采用雙線性變換法，還需進(jìn)行預(yù)畸變。（2）根據(jù)歸一化頻率，確定最小階數(shù)N 和頻率參數(shù)Wn 。可供選用的階數(shù)擇函數(shù)有:buttord，cheblord ，cheb2ord ，ellipord 等。（3）運(yùn)用最小階數(shù)N 設(shè)計(jì)模擬低通濾波器原型。模擬低通濾波器的創(chuàng)建函數(shù)有:buttap,cheblap, cheb2ap,ellipap 和besselap ，這些函數(shù)輸出的是零極點(diǎn)式形式，還要用zp2tf 函數(shù)轉(zhuǎn)換成分子分母多項(xiàng)式形式。如果想根據(jù)最小階數(shù)直接設(shè)計(jì)模擬低通濾波器原型，可用butter,chebyl,cheby2,ellip,bessel 等函數(shù)，只是注意要將函數(shù)中的Wn 設(shè)為1。（4）根

13、據(jù)第2步的頻率參數(shù)Wn ，模擬低通濾波原型轉(zhuǎn)換模擬低通、高通、帶通、帶阻濾波器，可用函數(shù)分別是：lp21p,lp2hp,lp2bp,lp2bs 。（5）運(yùn)用脈沖響應(yīng)不變法或雙線性變法把模濾波器轉(zhuǎn)數(shù)字濾波器，調(diào)用的函數(shù)是impinvar 和bilinear 。脈沖響應(yīng)不變法適用于采樣頻率大于4倍截止頻率的銳截止低通帶通濾波器，而雙線性變換法適合于相位特性要求不高的各型濾波器。fp=2100; fs=8000; Fs=20000; Rp=0.5; Rs=30;T=1/Fs; %設(shè)計(jì)指標(biāo) W1p=fp/Fs*2;W1s=fs/Fs*2;%求歸一化頻率N,Wn=buttord(W1p,W1s,Rp,R

14、s,s;%確定butterworth 的最小介數(shù)N 和頻率參數(shù)Wn z,p,k=buttap(N; %設(shè)計(jì)模擬低通原型的零極點(diǎn)增益參數(shù) bp,ap=zp2tf(z,p,k; %將零極點(diǎn)增益轉(zhuǎn)換成分子分母參數(shù)）（）（1893. 1s 4144. 152. 024144. 1483. 12s 1.4144 2.143S s23625. 3 (+=s s s s H abs,as=lp2lp(bp,ap,Wn*pi*Fs; %將低通原型轉(zhuǎn)換為模擬低通bz,az=impinvar(bs,as,Fs; %用脈沖響應(yīng)不變法進(jìn)行模數(shù)變換sys=tf(bz,az,T; %給出傳輸函數(shù)H(ZH,W=freqz(

15、bz,az,512,Fs; %生成頻率響應(yīng)參數(shù)subplot(2,1,1;plot(W,20*log10(abs(H; %繪制幅頻響應(yīng)grid on; %加坐標(biāo)網(wǎng)格xlabel(頻率/Hz;ylabel(振幅/dB;subplot(2,1,2;plot(W,abs(H; grid on;xlabel(頻率/Hz;ylabel(振幅/H;運(yùn)行后的波形如下： 圖2.1 巴特沃思和橢圓數(shù)字低通濾波器N =4bz = 0.0000 0.0999 0.1914 0.0252az= 1.0000 -1.4336 1.0984 -0.4115 0.06273.3基于漢明窗函數(shù)的FIR 低通濾波器的設(shè)計(jì)設(shè)低通

16、濾波器通帶截止頻率為0.2p =，阻帶截止頻率為0.4s =，通帶紋波p A =0.5dB，最小阻帶衰減r A =50dB。對(duì)于線性相位FIR 濾波器通常采用窗函數(shù)法設(shè)計(jì)。窗函數(shù)設(shè)計(jì)FIR 濾波器基本思想，就是從時(shí)域出發(fā), 把理想的無限長的 (n h d 用一定形狀的窗函數(shù)截取成有限長的( h n ，以此( h n 來逼近( d h n ，從而使所得到的頻率響應(yīng)( jw H e 與所要求的理想頻率響應(yīng)( jw d H e 相接近。理想的數(shù)字濾波器頻率特性( jw d H e 是無法實(shí)現(xiàn)的，F(xiàn)IR 的設(shè)計(jì)就是要尋找一個(gè)可以得到的頻率特性( jw H e 來逼近( jw d H e ，這相當(dāng)于用一

17、個(gè)可實(shí)現(xiàn)的單位脈沖響應(yīng)( h n 去逼近一個(gè)理想單位脈沖響( d h n 。要想用一個(gè)有限長的因果序列去逼近它，最簡單的方法是截取n 從0N-1的一段來表示它，即：( ( d h n h n = （0n N-1）；當(dāng)n 為其他的數(shù)時(shí)，( h n =0。同時(shí)，為了保證線性相位，還要滿足偶對(duì)稱：( (1 h n h N n =- ，這就好像通過一個(gè)窗口觀看到的一段( d h n ，因此( h n 就表示成( d h n 和一個(gè)“ 窗口函數(shù)” 的乘積。這樣，對(duì)( h n 的求解就變?yōu)椋? ( ( d h n h n W n =*這里的( W n 就稱為窗口函數(shù)。窗口函數(shù)相當(dāng)于對(duì)( d h n 在時(shí)域

18、上截取一段，勢必造成頻域的矩形窗口的失真，結(jié)果就是截取出的信號(hào)也相應(yīng)失真。為了補(bǔ)償這種失真，只有改變原來窗口的形狀，修正經(jīng)過時(shí)域截取后的窗口失真。改善窗口函數(shù)的標(biāo)準(zhǔn)：（1）盡量減少窗口頻譜的旁瓣，使能量集中在主瓣中，以減少波動(dòng)的幅度和次數(shù)，提高阻帶衰減。這樣，主瓣就不能太窄；（2）主瓣的寬度應(yīng)盡量窄，以得到盡量窄的過渡帶。工程中常用的窗函數(shù)有以下六種：矩形窗，三角形窗，漢寧窗，漢明窗，布萊克曼窗，凱澤窗(可變窗函數(shù) 。表1為六種窗函數(shù)的性能指標(biāo)。窗函數(shù) 窗普性能指標(biāo) 加窗后濾波器性能指標(biāo)- - 旁瓣峰值 dB 主瓣寬度 2/n 過濾寬度2/n 阻帶最小衰減 矩形窗 -13 2 0.9 -21三

19、角窗 -25 4 2.1 -25漢寧窗 -31 4 3.1 -44漢明窗 - 41 4 3.3 -53布萊克曼窗 -57 6 5.5 -74凱澤窗 -57 5 5 -80表3.1 常用傳函數(shù)的性能指標(biāo)wp=0.2*pi;wr=0.4*pi;trwidth=wr-wp; %過渡帶寬度N=ceil(6.64*pi/trwidth+1; %濾波器的長度wc=(wr+wp/2;hd=ideal_lp(wc,N;w_ham=(hamming(N;h=hd.*w_ham;db,mag,pha,grd,w=freqz_m(h,1;delta_w=2*pi/1000;Ap=-(min(db(1:1:wp/de

20、lta_w+1;Ar=-round(max(db(wr/delta_w+1:1:501;subplot(221;stem(n,hd;title(理想單位脈沖響應(yīng)hd(nsubplot(222;plot(n,w_ham;title(海明窗w(nsubplot(223;stem(n,h;title(實(shí)際單位脈沖響應(yīng)h(nsubplot(224;plot(w/pi,db;title(幅度響應(yīng)(dB輸出結(jié)果如下： 圖3.10.0220404383682397 0.00935693617842729 -0.0232325534213008下面我把濾波器的階數(shù)人為的設(shè)置成15，我們得到如下的輸出結(jié)果： 圖3.2由上圖可以看到：通帶紋波為1.6322dB ，最小阻帶衰