基于MATLAB的IIR數(shù)字濾波器的設(shè)計(jì)

1、基于Matlab 的IIR 數(shù)字濾波器的設(shè)計(jì)一 緒論數(shù)字濾波是數(shù)字信號(hào)處理的重要內(nèi)容，是由乘法器、加法器和單位延時(shí)器組成的一種運(yùn)算過(guò)程，其功能是對(duì)輸入離散信號(hào)進(jìn)行運(yùn)算處理，以達(dá)到改變信號(hào)頻譜的目的。數(shù)字濾波器根據(jù)頻域特性可分為低、高通、帶通和帶阻四個(gè)基本類型；根據(jù)時(shí)域特性可分為無(wú)限長(zhǎng)單位沖激響應(yīng)IIR(Infinite Impulse Response濾波器和有限長(zhǎng)單位沖激響應(yīng)FIR(Finite Impulse Response濾波器。數(shù)字濾波在通信、圖像編碼、語(yǔ)音編碼、雷達(dá)等許多領(lǐng)域中有著十分廣泛的應(yīng)用。目前，數(shù)字信號(hào)濾波器的設(shè)計(jì)圖像處理、數(shù)據(jù)壓縮等方面的應(yīng)用取得了令人矚目的進(jìn)展和成就。鑒

2、于此，數(shù)字濾波器的設(shè)計(jì)就顯得尤為重要。MATLAB是美國(guó)MathWorks 公司推出的一套用于工程計(jì)算的可視化高性能語(yǔ)言與軟件環(huán)境。MATLAB 為數(shù)字濾波的研究和應(yīng)用提供了一個(gè)直觀、高效、便捷的利器。它以矩陣運(yùn)算為基礎(chǔ)，把計(jì)算、可視化、程序設(shè)計(jì)融合到了一個(gè)交互式的工作環(huán)境中。MATLAB 推出的工具箱使各個(gè)領(lǐng)域的研究人員可以直觀方便地進(jìn)行科學(xué)研究、工程應(yīng)用，其中的信號(hào)處理(signal processing 、圖像處理(image processing, 小波(wavelet等工具箱為數(shù)字濾波研究的蓬勃發(fā)展提供了有力的工具。二 數(shù)字濾波器2.1什么是數(shù)字濾波器濾波器是指用來(lái)對(duì)輸入信號(hào)進(jìn)行濾波

3、的硬件和軟件。所謂數(shù)字濾波器，是指輸入、輸出均為數(shù)字信號(hào)，通過(guò)一定運(yùn)算關(guān)系改變輸入信號(hào)所含頻率成分的相對(duì)比例或者濾除某些頻率成分的器件。數(shù)字濾波器和模擬濾波器相比，因?yàn)樾盘?hào)的形式和實(shí)現(xiàn)濾波的方法不同，數(shù)字濾波器具有比模擬濾波器精度高、穩(wěn)定、體積小、重量輕、靈活、不要求阻抗匹配等優(yōu)點(diǎn)。一般用兩種方法來(lái)實(shí)現(xiàn)數(shù)字濾波器:一是采用通用計(jì)算機(jī)，把濾波器所要完成的運(yùn)算編成程序通過(guò)計(jì)算機(jī)來(lái)執(zhí)行，也就是采用計(jì)算機(jī)軟件來(lái)實(shí)現(xiàn); 二是設(shè)計(jì)專用的數(shù)字處理硬件。MATLAB 的信號(hào)處理工具箱是專門應(yīng)用于信號(hào)處理領(lǐng)域的專用工具箱，它的兩個(gè)基本組成就是濾波器的設(shè)計(jì)與實(shí)現(xiàn)部分以及譜分析部分。工具箱提供了豐富而簡(jiǎn)便的設(shè)計(jì)，使

4、原來(lái)繁瑣的程序設(shè)計(jì)簡(jiǎn)化成函數(shù)的調(diào)用。只要以正確的指標(biāo)參數(shù)調(diào)用相應(yīng)的濾波器設(shè)計(jì)程序或工具箱函數(shù)，便可以得到正確的設(shè)計(jì)結(jié)果，使用非常方便。2.2數(shù)字濾波器的分類數(shù)字濾波器從功能上分類帶阻濾波器。從濾波器的網(wǎng)絡(luò)結(jié)構(gòu)或者從單位脈沖響應(yīng)分類長(zhǎng)單位沖激響應(yīng)濾波器 和FIR函數(shù)分別為:H (z H (z 第一個(gè)公式中的H (z稱為(N-1階FIR 濾波器函數(shù)。2.3數(shù)字濾波器的設(shè)計(jì)要求濾波器的指標(biāo)常常在頻域給出。數(shù)，所以通常表示為：H (e jw =其中， |稱為幅頻特性函數(shù)，示信號(hào)通過(guò)該濾波器后各頻率成分的衰減情況，器后在時(shí)間上的延時(shí)情況。一般|來(lái)描述設(shè)計(jì)指標(biāo)，相頻特性一般不作要求。IIR 濾波器指標(biāo)參數(shù)

5、如下圖所示。圖中，阻帶邊界頻率；1和2分別為通帶波紋和阻帶波紋；允許的衰減一般用可分為低通濾波器、高通濾波器、帶通濾波器、:可分為IIR 濾波器(即無(wú)限濾波器(即有限長(zhǎng)單位沖激響應(yīng)濾波器 。它們的M b -r r z r =0N 1+a k z -k k =1N -1h (n z -n n =0N 階IIR 濾波器函數(shù)，第二個(gè)公式中的H (z 數(shù)字濾波器的頻響特性函數(shù)一般為復(fù)函H (z |z =e jw =H (e jw e j (w (w稱為相頻特性函數(shù)。幅頻特性表而相頻特性反映各頻率通過(guò)濾波IIR 數(shù)字濾波器，通常只用幅頻響應(yīng)函數(shù) p 和s 分別為通帶邊界頻率和dB :=稱為數(shù)表示，通帶內(nèi)

6、所允許的最大衰減（dB ）和阻帶內(nèi)允許的最小衰減（dB ）分別為p 和s 表示： 一般要求： 低通濾波器的技術(shù)要求2.4數(shù)字濾波器設(shè)計(jì)方法概述IIR 數(shù)字濾波器的設(shè)計(jì)步驟流程圖如下： 步驟流程圖IIR 濾波器設(shè)計(jì)方法有兩類，經(jīng)常用到的一類設(shè)計(jì)方法是借助于模擬濾波器的設(shè)計(jì)方法進(jìn)行的。其設(shè)計(jì)思路是:先設(shè)計(jì)模擬濾波器得到傳輸函數(shù)Ha(s，然后將Ha(s按某種方法轉(zhuǎn)換成數(shù)字濾波器的系統(tǒng)函數(shù)H (Z。這一類方法是基于模擬濾波器的設(shè)計(jì)方法相對(duì)比較成熟，它不僅有完整的設(shè)計(jì)公式，也有完整的圖表供查閱。更可以直接調(diào)用MATLAB 中的對(duì)應(yīng)的函數(shù)進(jìn)行設(shè)計(jì)。另一類是直接在頻域或者時(shí)域中進(jìn)行設(shè)計(jì)的，設(shè)計(jì)時(shí)必須用計(jì)算機(jī)

7、作輔助設(shè)計(jì)，直接調(diào)用MATLAB 中的一些程序或者函數(shù)可以很方便地設(shè)計(jì)出所需要的濾波器。三 IIR 濾波器的設(shè)計(jì)3.1典型的IIR 數(shù)字濾波器的設(shè)計(jì)模擬濾波器的理論和設(shè)計(jì)方法已發(fā)展得相當(dāng)成熟，且有一些典型的模擬濾波器供我們選擇，如巴特沃斯(Butterworth 濾波器、切比雪夫(Chebyshev濾波器、橢圓(Cauer 濾波器、貝塞爾(Bessel 濾波器等，這些典型的濾波器各有特點(diǎn)。用MATLAB 進(jìn)行典型的數(shù)字濾波器的設(shè)計(jì)，一般步驟如下:（1） 將設(shè)計(jì)指標(biāo)歸一化處理。如果采用雙線性變換法，還需進(jìn)行預(yù)畸變。（2） 根據(jù)歸一化頻率，確定最小階數(shù)N 和頻率參數(shù)Wn ?？晒┻x用的階數(shù)擇函數(shù)有:

8、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） 根據(jù)第2步的頻率參數(shù)Wn ，模擬低通濾波原型轉(zhuǎn)換模擬低通、高通、帶通、帶阻濾波器，可用函數(shù)分別是：lp21p,lp

9、2hp,lp2bp,lp2bs 。（5） 運(yùn)用脈沖響應(yīng)不變法或雙線性變法把模濾波器轉(zhuǎn)數(shù)字濾波器，調(diào)用的函數(shù)是impinvar 和bilinear 。脈沖響應(yīng)不變法適用于采樣頻率大于4倍截止頻率的銳截止低通帶通濾波器，而雙線性變換法適合于相位特性要求不高的各型濾波器。（6） 根據(jù)輸出的分子分母系數(shù)，用tf 函數(shù)生成H （z ）的表達(dá)式，再用freqz函數(shù)驗(yàn)證設(shè)計(jì)結(jié)果。設(shè)計(jì)巴特沃思數(shù)字低通濾波器和橢圓數(shù)字低通濾波器，要求通帶邊界頻率fp=2.1kHZ，通帶最大衰減Rp=0.5dB；阻帶邊界頻率fs=8kHZ，阻帶最小衰減Rs=30dB,采樣頻率為Fs=20kHZ。用脈沖響應(yīng)不變法設(shè)計(jì)的巴特沃思數(shù)字

10、低通濾波器的M 程序如下： 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,Rs,'s'%確定butterworth 的最小介數(shù)N 和頻率參數(shù)Wnz,p,k=buttap(N; %設(shè)計(jì)模擬低通原型的零極點(diǎn)增益參數(shù)bp,ap=zp2tf(z,p,k; %將零極點(diǎn)增益轉(zhuǎn)換成分子分母參數(shù)bs,as=lp2lp(bp,ap,Wn*pi*Fs;%將低通原型轉(zhuǎn)換為模擬低通bz,az=impinvar(bs,as,Fs; %用脈沖

11、響應(yīng)不變法進(jìn)行模數(shù)變換 sys=tf(bz,az,T; %給出傳輸函數(shù)H(ZH,W=freqz(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)行后的波形如下： 運(yùn)行結(jié)果：N =4bz = 0.0000 0.0999 0.1914 0.

12、0252az= 1.0000 -1.4336 1.0984 -0.4115 0.0627用雙線性變換法設(shè)計(jì)的橢圓數(shù)字低通濾波器的M 程序如下：fs=20000;wp=2*pi*2100/fs;ws=2*pi*8000/fs;Rp=0.5;Rs=30;Ts=1/fs;Wp=2/Ts*tan(wp/2;Ws=2/Ts*tan(ws/2; %按頻率轉(zhuǎn)換公式進(jìn)行轉(zhuǎn)換N,Wn=ellipord(Wp,Ws,Rp,Rs,'s' %計(jì)算模擬濾波器的最小階數(shù)z,p,k=ellipap(N,Rp,Rs;%設(shè)計(jì)模擬原型濾波器Bap,Aap=zp2tf(z,p,k; %零點(diǎn)極點(diǎn)增益形式轉(zhuǎn)換為傳遞函數(shù)

13、形式b,a=lp2lp(Bap,Aap,Wn; %低通轉(zhuǎn)換為低通濾波器的頻率轉(zhuǎn)化bz,az=bilinear(b,a,fs; %運(yùn)用雙線性變換法得到數(shù)字濾波器傳遞函數(shù)H,f=freqz(bz,az,512,fs;subplot(2,1,1;plot(f,20*log10(abs(H;title('N=2 頻率響應(yīng)'grid on;xlabel('頻率/Hz'ylabel('振幅/dB'subplot(2,1,2;plot(f,abs(H; grid on;xlabel('頻率/Hz'ylabel('振幅/H' 運(yùn)行

14、結(jié)果：N=2bz= 0.1213 0.1662 0.1213az= 1.0000 -0.9889 0.42183.2完全濾波器設(shè)計(jì)除了典型設(shè)計(jì)以外，MATLAB 信號(hào)處理工具箱提供了幾個(gè)直接設(shè)計(jì)IIR 數(shù)字濾波器的函數(shù)，直接調(diào)用就可以設(shè)計(jì)濾波器，這為設(shè)計(jì)通用濾波器提供了方便。 設(shè)計(jì)Butterworth 濾波器用函數(shù)butter(，可以設(shè)計(jì)低通、高通、帶通和帶阻的數(shù)字和模擬濾波器，其特性是通帶內(nèi)的幅度響應(yīng)最大限度的平滑，但損失了截止頻率處的下降斜度。設(shè)計(jì)Chebyshev I型濾波器用函數(shù)chebyl(?？梢栽O(shè)計(jì)低通、高通、帶通和帶阻的數(shù)字和模擬Chebyshev I 型濾波器，其通帶內(nèi)為等波

15、紋，阻帶內(nèi)為單調(diào)。Chebyshev I型濾波器的下降斜度比II 型大，但其代價(jià)目是通帶內(nèi)波紋較大。 設(shè)計(jì)Chebyshev II 型濾波器用函數(shù)cheby2(?？梢栽O(shè)計(jì)低通、高通、帶通和帶阻的數(shù)字和模擬Chebyshev II 型濾波器，其通帶內(nèi)為單調(diào)，阻帶內(nèi)等波紋。Chebyshev II型濾波器的下降斜度比I 型小，但其阻帶內(nèi)波紋較大。設(shè)計(jì)橢圓濾波器用函數(shù)ellip(，與chebyl, cheby2類似，可以設(shè)計(jì)低通、高通、帶通和帶阻的數(shù)字和模擬濾波器。與Butterworth 和chebyshev 濾波器相比，ellip 函數(shù)可以得到下降斜度更大的濾波器，得通帶和阻帶均為等波紋。一般情

16、況下，橢圓濾波器能以最低的階實(shí)現(xiàn)指定的性能指標(biāo)。在使用各類濾波器函數(shù)時(shí)應(yīng)當(dāng)注意以下重點(diǎn):A 、階數(shù)和固有頻率的選擇:N,Wn=buttord(Wp,Ws,Rp,Rs可得到符合要求性質(zhì)的濾波器的最小階數(shù)N 以及數(shù)字Butterworth 濾波器的固有頻率Wn(即3dB 。設(shè)計(jì)的要求是在通帶內(nèi)的衰減不超過(guò)Rp, 在阻帶內(nèi)的衰減不小于Rs ，通帶和阻帶有截止頻率分別是Wp, Ws,它們是歸一化的頻率，范圍是0, 1，對(duì)應(yīng)弧度。B 、有關(guān)濾波器設(shè)計(jì)當(dāng)中的頻率歸一化問(wèn)題:信號(hào)處理工具箱中經(jīng)常使用的頻率是Nyquist 頻率，它被定義為采樣頻率的一半，在濾波器的階數(shù)選擇和設(shè)計(jì)中的截止頻率均使用Nyquis

17、t 頻率進(jìn)行歸一化處理。例如對(duì)于一個(gè)采樣頻率為1000 Hz 的系統(tǒng)，400Hz 的歸一化即為400/500=0.8。歸一化頻率的范圍在0, 1之間。如果要將歸一化頻率轉(zhuǎn)換為角頻率，則將歸一化頻率乘以; 如果要將歸一化頻率轉(zhuǎn)換為Hz, 則將歸一化頻率乘以采樣頻率的一半。C 、設(shè)計(jì)一個(gè)N 階的低通Butterworth 濾波器使用函數(shù)B,A=butter(N, Wn ，返回濾波器系數(shù)矩陣B,A。其中固有頻率Wn 必須是歸一化頻率。它的最大值是采樣頻率的一半。Fs 缺省時(shí)默認(rèn)為2Hz 。如果Wn=Wl,W2是一個(gè)兩元素的向量，則函數(shù)將設(shè)計(jì)出一個(gè)2N 階的帶通濾波器，通帶為W1,W2。設(shè)計(jì)Cheby

18、shev I 型和Chebyshev II 型數(shù)字低通濾波器，要求通帶邊界頻率fp=2.1kHZ，通帶最大衰減Rp=0.5dB；阻帶邊界頻率fs=8kHZ，阻帶最小衰減Rs=30dB,采樣頻率為Fs=20kHZ。Chebyshev I型的M 程序如下：Fs=20000; %抽樣頻率20KHzFlp=2100;Fls=8000;Wp=2*Flp/Fs;Ws=2*Fls/Fs;Rp=0.5; %歸一化的通帶截止頻率 %歸一化的阻帶截止頻率 %通帶最大衰減（單位：dB ）Rs=30; %阻帶最小衰減（單位：dB ）N,Wn=cheb1ord(Wp,Ws,Rp,Rs; %返回最小階數(shù)和截止頻率b,a=

19、cheby1(N,Rp,Wn; %返回H(z的分子分母系數(shù)hw,w=freqz(b,a;subplot(2,1,1;plot(w/pi,20*log10(abs(hw;grid on;xlabel('/'ylabel('幅度（dB ）'title('切比雪夫I 型幅頻響應(yīng)'subplot(2,1,2;plot(w/pi,abs(hw;grid on;xlabel('/'ylabel('幅度（H ）'運(yùn)行后的波形如下： 運(yùn)行結(jié)果：N=2b= 0.1007 0.2014 0.1007a= 1.0000 -0.9872

20、0.4140Chebyshev II型M 程序如下：Fs=20000; %抽樣頻率20KHzFlp=2100;Fls=8000;Wp=2*Flp/Fs; %歸一化的通帶截止頻率Ws=2*Fls/Fs; %歸一化的阻帶截止頻率Rp=0.5;Rs=30; %通帶最大衰減（單位：dB ） %阻帶最小衰減（單位：dB ）N,Wn=cheb2ord(Wp,Ws,Rp,Rs; %返回最小階數(shù)和截止頻率b,a=cheby2(N,Rs,Wn; %返回H(z的分子分母系數(shù)hw,w=freqz(b,a;subplot(2,1,1;plot(w/pi,20*log10(abs(hw;grid on;xlabel(&

21、#39;/'ylabel('幅度（dB ）'title('切比雪夫II 型幅頻響應(yīng)'subplot(2,1,2;plot(w/pi,abs(hw;grid on;xlabel('/'ylabel('幅度（H ）'運(yùn)行后的波形： 運(yùn)行結(jié)果：N=2b= 0.2357 0.4241 0.2357a= 1.0000 -0.2996 0.19503.3 結(jié)果分析從頻率響應(yīng)圖中可以看出:巴特沃斯濾波器具有單調(diào)下降的幅頻特性，通帶內(nèi)平滑; 切比雪I 型濾波器的幅頻特性在通帶內(nèi)有波動(dòng)，阻帶內(nèi)單調(diào);chebyshev II 型濾波器的幅頻特

22、性在阻帶內(nèi)有波動(dòng)，通帶內(nèi)單調(diào); 橢圓濾波器的選擇性相對(duì)前三種是最好的，下降斜度比較大，通帶和阻帶內(nèi)均為等波紋，同樣的性能指標(biāo)，橢圓濾波器可以最低的階數(shù)來(lái)實(shí)現(xiàn)。這樣根據(jù)不同的要求可以選用不同類型的濾波器。四 對(duì)語(yǔ)音信號(hào)的濾波利用函數(shù)wavread 對(duì)語(yǔ)音信號(hào)的采集：i=1;x,fs,bits=wavread('1.wav'%x：語(yǔ)音數(shù)據(jù)；fs ：采樣頻率；bits ：采樣點(diǎn)數(shù)sound(x,fs,bits; %話音回放N=length(x;n=0:N-1;figure(i;subplot(2,1,1;plot(n,x; %畫出原始語(yǔ)音信號(hào)的波形xlabel('n'

23、ylabel('x(n'title('原始語(yǔ)音信號(hào)'subplot(2,1,2;H,f=freqz(x,1,512,fs;plot(f,20*log10(abs(H; %畫出原始語(yǔ)音信號(hào)的頻譜xlabel('n'ylabel('x(n'title('原始語(yǔ)音信號(hào)的頻譜 '運(yùn)行結(jié)果：fs=44100bits=16加入上面設(shè)計(jì)的橢圓濾波器對(duì)其進(jìn)行濾波：fs=20000;wp=2*pi*2100/fs;ws=2*pi*8000/fs;Rp=0.5;Rs=30;Ts=1/fs;Wp=2/Ts*tan(wp/2;Ws=2/T

24、s*tan(ws/2; %按頻率轉(zhuǎn)換公式進(jìn)行轉(zhuǎn)換N,Wn=ellipord(Wp,Ws,Rp,Rs,'s' %計(jì)算模擬濾波器的最小階數(shù)z,p,k=ellipap(N,Rp,Rs;%設(shè)計(jì)模擬原型濾波器Bap,Aap=zp2tf(z,p,k; %零點(diǎn)極點(diǎn)增益形式轉(zhuǎn)換為傳遞函數(shù)形式b,a=lp2lp(Bap,Aap,Wn; %低通轉(zhuǎn)換為低通濾波器的頻率轉(zhuǎn)化bz,az=bilinear(b,a,fs; %運(yùn)用雙線性變換法得到數(shù)字濾波器傳遞函數(shù)H,f=freqz(bz,az,512,fs;plot(f,20*log10(abs(H;title('N=2 頻率響應(yīng)'grid on;xlabel('頻率/Hz'ylabel('振幅 /dB'%-濾波前后波形比較-%i =i+1;figure(i;y=filter(bz,az,x;subplot(2,2,1;plot(n,x; %畫出原語(yǔ)音信號(hào)的波