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

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

3、的硬件和軟件。所謂數字濾波器，是指輸入、輸出均為數字信號，通過一定運算關系改變輸入信號所含頻率成分的相對比例或者濾除某些頻率成分的器件。數字濾波器和模擬濾波器相比，因為信號的形式和實現濾波的方法不同，數字濾波器具有比模擬濾波器精度高、穩定、體積小、重量輕、靈活、不要求阻抗匹配等優點。一般用兩種方法來實現數字濾波器:一是采用通用計算機，把濾波器所要完成的運算編成程序通過計算機來執行，也就是采用計算機軟件來實現; 二是設計專用的數字處理硬件。MATLAB 的信號處理工具箱是專門應用于信號處理領域的專用工具箱，它的兩個基本組成就是濾波器的設計與實現部分以及譜分析部分。工具箱提供了豐富而簡便的設計，使

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

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

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

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

8、buttord，cheblord ，cheb2ord ，ellipord 等。（3） 運用最小階數N 設計模擬低通濾波器原型。模擬低通濾波器的創建函數有:buttap,cheblap, cheb2ap,ellipap 和besselap ，這些函數輸出的是零極點式形式，還要用zp2tf 函數轉換成分子分母多項式形式。如果想根據最小階數直接設計模擬低通濾波器原型，可用butter,chebyl,cheby2,ellip,bessel 等函數，只是注意要將函數中的Wn 設為1。（4） 根據第2步的頻率參數Wn ，模擬低通濾波原型轉換模擬低通、高通、帶通、帶阻濾波器，可用函數分別是：lp21p,lp

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

10、低通濾波器的M 程序如下： fp=2100;fs=8000;Fs=20000;Rp=0.5;Rs=30;T=1/Fs; %設計指標W1p=fp/Fs*2;W1s=fs/Fs*2;%求歸一化頻率N,Wn=buttord(W1p,W1s,Rp,Rs,'s'%確定butterworth 的最小介數N 和頻率參數Wnz,p,k=buttap(N; %設計模擬低通原型的零極點增益參數bp,ap=zp2tf(z,p,k; %將零極點增益轉換成分子分母參數bs,as=lp2lp(bp,ap,Wn*pi*Fs;%將低通原型轉換為模擬低通bz,az=impinvar(bs,as,Fs; %用脈沖

11、響應不變法進行模數變換 sys=tf(bz,az,T; %給出傳輸函數H(ZH,W=freqz(bz,az,512,Fs; %生成頻率響應參數subplot(2,1,1;plot(W,20*log10(abs(H; %繪制幅頻響應grid on; %加坐標網格xlabel('頻率/Hz'ylabel('振幅/dB'subplot(2,1,2;plot(W,abs(H; grid on;xlabel('頻率/Hz'ylabel('振幅/H'運行后的波形如下： 運行結果：N =4bz = 0.0000 0.0999 0.1914 0.

12、0252az= 1.0000 -1.4336 1.0984 -0.4115 0.0627用雙線性變換法設計的橢圓數字低通濾波器的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; %按頻率轉換公式進行轉換N,Wn=ellipord(Wp,Ws,Rp,Rs,'s' %計算模擬濾波器的最小階數z,p,k=ellipap(N,Rp,Rs;%設計模擬原型濾波器Bap,Aap=zp2tf(z,p,k; %零點極點增益形式轉換為傳遞函數

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

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

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

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

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

18、shev I 型和Chebyshev II 型數字低通濾波器，要求通帶邊界頻率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; %返回最小階數和截止頻率b,a=

19、cheby1(N,Rp,Wn; %返回H(z的分子分母系數hw,w=freqz(b,a;subplot(2,1,1;plot(w/pi,20*log10(abs(hw;grid on;xlabel('/'ylabel('幅度（dB ）'title('切比雪夫I 型幅頻響應'subplot(2,1,2;plot(w/pi,abs(hw;grid on;xlabel('/'ylabel('幅度（H ）'運行后的波形如下： 運行結果：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; %返回最小階數和截止頻率b,a=cheby2(N,Rs,Wn; %返回H(z的分子分母系數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 型幅頻響應'subplot(2,1,2;plot(w/pi,abs(hw;grid on;xlabel('/'ylabel('幅度（H ）'運行后的波形： 運行結果：N=2b= 0.2357 0.4241 0.2357a= 1.0000 -0.2996 0.19503.3 結果分析從頻率響應圖中可以看出:巴特沃斯濾波器具有單調下降的幅頻特性，通帶內平滑; 切比雪I 型濾波器的幅頻特性在通帶內有波動，阻帶內單調;chebyshev II 型濾波器的幅頻特

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

23、ylabel('x(n'title('原始語音信號'subplot(2,1,2;H,f=freqz(x,1,512,fs;plot(f,20*log10(abs(H; %畫出原始語音信號的頻譜xlabel('n'ylabel('x(n'title('原始語音信號的頻譜 '運行結果：fs=44100bits=16加入上面設計的橢圓濾波器對其進行濾波：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; %按頻率轉換公式進行轉換N,Wn=ellipord(Wp,Ws,Rp,Rs,'s' %計算模擬濾波器的最小階數z,p,k=ellipap(N,Rp,Rs;%設計模擬原型濾波器Bap,Aap=zp2tf(z,p,k; %零點極點增益形式轉換為傳遞函數形式b,a=lp2lp(Bap,Aap,Wn; %低通轉換為低通濾波器的頻率轉化bz,az=bilinear(b,a,fs; %運用雙線性變換法得到數字濾波器傳遞函數H,f=freqz(bz,az,512,fs;plot(f,20*log10(abs(H;title('N=2 頻率響應'grid on;xlabel('頻率/Hz'ylabel('振幅 /dB'%-濾波前后波形比較-%i =i+1;figure(i;y=filter(bz,az,x;subplot(2,2,1;plot(n,x; %畫出原語音信號的波