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

1、第一章 緒論數(shù)字濾波是數(shù)字信號(hào)處理的重要內(nèi)容，是由乘法器、加法器和單位延時(shí)器組成的一種運(yùn)算過(guò)程，其功能是對(duì)輸入離散信號(hào)進(jìn)行運(yùn)算處理，以達(dá)到改變信號(hào)頻譜的目的。數(shù)字濾波器根據(jù)頻域特性可分為低、高通、帶通和帶阻四個(gè)基本類(lèi)型；根據(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)展和成就。鑒于此，數(shù)字濾波器的設(shè)計(jì)就顯得尤為重要

2、。 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)行濾波的硬件和軟件。所謂數(shù)字濾波器，是指

3、輸入、輸出均為數(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ì)專(zhuān)用的數(shù)字處理硬件。 MATLAB的信號(hào)處理工具箱是專(zhuān)門(mén)應(yīng)用于信號(hào)處理領(lǐng)域的專(zhuān)用工具箱，它的兩個(gè)基本組成就是濾波器的設(shè)計(jì)與實(shí)現(xiàn)部分以及譜分析部分。工具箱提供了豐富而簡(jiǎn)便的設(shè)計(jì)，使原來(lái)繁瑣的程序設(shè)計(jì)簡(jiǎn)化成函數(shù)的調(diào)用

4、。只要以正確的指標(biāo)參數(shù)調(diào)用相應(yīng)的濾波器設(shè)計(jì)程序或工具箱函數(shù)，便可以得到正確的設(shè)計(jì)結(jié)果，使用非常方便。2.2數(shù)字濾波器的分類(lèi) 數(shù)字濾波器從功能上分類(lèi):可分為低通濾波器、高通濾波器、帶通濾波器、帶阻濾波器。 從濾波器的網(wǎng)絡(luò)結(jié)構(gòu)或者從單位脈沖響應(yīng)分類(lèi):可分為IIR濾波器(即無(wú)限長(zhǎng)單位沖激響應(yīng)濾波器)和FIR濾波器(即有限長(zhǎng)單位沖激響應(yīng)濾波器)。它們的函數(shù)分別為:第一個(gè)公式中的H (z)稱(chēng)為N階IIR濾波器函數(shù)，第二個(gè)公式中的H (z)稱(chēng)為(N-1)階FIR濾波器函數(shù)。2.3數(shù)字濾波器的設(shè)計(jì)要求 濾波器的指標(biāo)常常在頻域給出。數(shù)字濾波器的頻響特性函數(shù)一般為復(fù)函數(shù)，所以通常表示為：其中，|稱(chēng)為幅頻特性函數(shù)

5、，(w)稱(chēng)為相頻特性函數(shù)。幅頻特性表示信號(hào)通過(guò)該濾波器后各頻率成分的衰減情況，而相頻特性反映各頻率通過(guò)濾波器后在時(shí)間上的延時(shí)情況。一般IIR數(shù)字濾波器，通常只用幅頻響應(yīng)函數(shù)|來(lái)描述設(shè)計(jì)指標(biāo)，相頻特性一般不作要求。IIR濾波器指標(biāo)參數(shù)如下圖所示。圖中，p和s分別為通帶邊界頻率和阻帶邊界頻率；1和2分別為通帶波紋和阻帶波紋；允許的衰減一般用dB數(shù)表示，通帶內(nèi)所允許的最大衰減（dB）和阻帶內(nèi)允許的最小衰減（dB）分別為p和s表示： 一般要求：低通濾波器的技術(shù)要求2.4數(shù)字濾波器設(shè)計(jì)方法概述IIR數(shù)字濾波器的設(shè)計(jì)步驟流程圖如下：步驟流程圖 IIR濾波器設(shè)計(jì)方法有兩類(lèi)，經(jīng)常用到的一類(lèi)設(shè)計(jì)方法是借助于模擬

6、濾波器的設(shè)計(jì)方法進(jìn)行的。其設(shè)計(jì)思路是:先設(shè)計(jì)模擬濾波器得到傳輸函數(shù)Ha(s)，然后將Ha(s)按某種方法轉(zhuǎn)換成數(shù)字濾波器的系統(tǒng)函數(shù)H (Z)。這一類(lèi)方法是基于模擬濾波器的設(shè)計(jì)方法相對(duì)比較成熟，它不僅有完整的設(shè)計(jì)公式，也有完整的圖表供查閱。更可以直接調(diào)用MATLAB中的對(duì)應(yīng)的函數(shù)進(jìn)行設(shè)計(jì)。另一類(lèi)是直接在頻域或者時(shí)域中進(jìn)行設(shè)計(jì)的，設(shè)計(jì)時(shí)必須用計(jì)算機(jī)作輔助設(shè)計(jì)，直接調(diào)用MATLAB中的一些程序或者函數(shù)可以很方便地設(shè)計(jì)出所需要的濾波器。第三章IIR濾波器的設(shè)計(jì)3.1典型的IIR數(shù)字濾波器的設(shè)計(jì)模擬濾波器的理論和設(shè)計(jì)方法已發(fā)展得相當(dāng)成熟，且有一些典型的模擬濾波器供我們選擇，如巴特沃斯(Butterwor

7、th )濾波器、切比雪夫(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ù)有:buttord，cheblord，cheb2ord，ellipord等。（3） 運(yùn)用最小階數(shù)N設(shè)計(jì)模擬低通濾波器原型。模擬低通濾波器的創(chuàng)建函數(shù)有:buttap,cheblap, cheb2ap,ellipap和besselap，這些函數(shù)輸出的是零極點(diǎn)式形式

8、，還要用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,lp2hp,lp2bp,lp2bs。（5） 運(yùn)用脈沖響應(yīng)不變法或雙線性變法把模濾波器轉(zhuǎn)數(shù)字濾波器，調(diào)用的函數(shù)是impinvar和bilinear。脈沖響應(yīng)不變法適用于采樣頻率大于4倍截止頻率的銳截止低通帶通濾波器，而雙線性變換法適合于相位特性要求不高的各型濾波器。（6） 根據(jù)輸出

9、的分子分母系數(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ù)字低通濾波器的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');%確定bu

10、tterworth的最小介數(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); %用脈沖響應(yīng)不變法進(jìn)行模數(shù)變換 sys=tf(bz,az,T); %給出傳輸函數(shù)H(Z)H,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)格xl

11、abel('頻率/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.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=

12、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ù)形式 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

13、,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)行結(jié)果：N=2bz= 0.1213 0.1662 0.1213az= 1.0000 -0.9889 0.42183.2完全濾波器設(shè)計(jì)除了典型設(shè)計(jì)以外，

14、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)為等波紋，阻帶內(nèi)為單調(diào)。Chebyshev I型濾波器的下降斜度比II型大，但其代價(jià)目是通帶內(nèi)波紋較大。 設(shè)計(jì)Chebyshev II型濾波器用函數(shù)cheb

15、y2()。可以設(shè)計(jì)低通、高通、帶通和帶阻的數(shù)字和模擬Chebyshev II型濾波器，其通帶內(nèi)為單調(diào)，阻帶內(nèi)等波紋。Chebyshev II型濾波器的下降斜度比I型小，但其阻帶內(nèi)波紋較大。 設(shè)計(jì)橢圓濾波器用函數(shù)ellip()，與chebyl, cheby2類(lèi)似，可以設(shè)計(jì)低通、高通、帶通和帶阻的數(shù)字和模擬濾波器。與Butterworth和chebyshev濾波器相比，ellip函數(shù)可以得到下降斜度更大的濾波器，得通帶和阻帶均為等波紋。一般情況下，橢圓濾波器能以最低的階實(shí)現(xiàn)指定的性能指標(biāo)。在使用各類(lèi)濾波器函數(shù)時(shí)應(yīng)當(dāng)注意以下重點(diǎn):A、階數(shù)和固有頻率的選擇:N,Wn=buttord(Wp,Ws,Rp,

16、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ì)中的截止頻率均使用Nyquist頻率進(jìn)行歸一化處理。例如對(duì)于一個(gè)采樣頻率為1000 Hz的系統(tǒng)，400Hz的歸一化即為400/500=0.8。歸一化頻率的范圍在0, 1之間。如果要將歸一化頻率轉(zhuǎn)換為

17、角頻率，則將歸一化頻率乘以;如果要將歸一化頻率轉(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ì)Chebyshev I型和Chebyshev II型數(shù)字低通濾波器，要求通帶邊界頻率fp=2.1kHZ，通帶最大衰減Rp=0.5dB；阻帶邊界頻率fs=8kHZ，阻帶最小衰減Rs=30dB,采樣頻率為

18、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=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;x

19、label('/');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 0.4140Chebyshev II型M程序如下：Fs=20000; %抽樣頻率20KHzFlp=2100; Fls=8000;Wp=2*Flp/Fs; %歸一化的通帶截止

20、頻率Ws=2*Fls/Fs; %歸一化的阻帶截止頻率Rp=0.5; %通帶最大衰減（單位：dB）Rs=30; %阻帶最小衰減（單位：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('/');ylabel('幅度（dB）')title('切比雪夫II型幅頻響應(yīng)');subplot(2,1,2);pl

21、ot(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型濾波器的幅頻特性在阻帶內(nèi)有波動(dòng)，通帶內(nèi)單調(diào);橢圓濾波器的選擇性相對(duì)前三種是最好的，下降斜度比較大，通帶和阻帶內(nèi)均為等波紋，同樣的性能指標(biāo)，橢圓濾波器可以最低的階數(shù)來(lái)實(shí)現(xiàn)。

22、這樣根據(jù)不同的要求可以選用不同類(lèi)型的濾波器。第四章 對(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); %畫(huà)出原始語(yǔ)音信號(hào)的波形 xlabel('n'); ylabel('x(n)'); title('原始語(yǔ)音信號(hào)'); subplot(2,1,

23、2); H,f=freqz(x,1,512,fs); plot(f,20*log10(abs(H); %畫(huà)出原始語(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/Ts*tan(ws/2); %按頻率轉(zhuǎn)換公式進(jìn)行轉(zhuǎn)換 N,Wn=el

24、lipord(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(&

25、#39;頻率/Hz'); ylabel('振幅/dB'); %-濾波前后波形比較-% i =i+1; figure(i); y=filter(bz,az,x); subplot(2,2,1); plot(n,x); %畫(huà)出原語(yǔ)音信號(hào)的波形 xlabel('n'); ylabel('x(n)'); title('原始語(yǔ)音信號(hào)'); subplot(2,2,3); plot(n,y); %IIRLowPass后的波形 title('IIRLowPass后的波形'); %-濾波前后語(yǔ)音信號(hào)頻譜-%subplot(2,2,2); H,f=freqz(x,1,512,fs); plot(f,20*log10(abs(H); %濾波前語(yǔ)音信號(hào)頻譜 xlabel('頻率/Hz'); ylabel('振幅/dB'); title('IIR低通濾波前語(yǔ)音信號(hào)頻譜'); grid on; subplot(2,2,4); H,f=freqz(y,1,512,fs); plot(f,20*log10(abs(H); %濾波后語(yǔ)音信號(hào)頻譜 xlabel('頻率/Hz'); ylabel('振幅/dB'); title(&#