濾波器設(shè)計―fdatool工具及代碼設(shè)計

1、黎美琪201300800610 通信工程2班實驗一基于FDATool的FIR及IIR濾波器設(shè)計一、實驗?zāi)康?1.學(xué)會使用fdatool設(shè)計濾波器2.分析比較不同濾波器的特性二、實驗條件Matlab 2013b pc機三、實驗過程知識儲備:使用FDATool設(shè)計FIR濾波器的具體步驟1.1 濾波器指標若需要設(shè)計一個16階的FIR濾波器(h(0=0,給定的參數(shù)如下:(1 低通濾波器注:階數(shù)問題這個選項直接影響濾波器的性能,階數(shù)越高,性能越好,但是相應(yīng)在FPGA實現(xiàn)耗用的資源需要增多。在這個設(shè)置中提供2個選項:Specify order和Minimum order,Specify order是工程師

2、自己確定濾波器的階數(shù),Minimum order是讓工具自動確定達到期望的頻率相應(yīng)所需要的最小階數(shù)。(2 采樣頻率Fs為8kHz,濾波器Fpass為3kHz,Fstop為3050hz(3 輸入序列位寬為9位(最高位為符號位在此利用MATLAB來完成FIR濾波器系數(shù)的確定。1.2 打開MATLAB的FDAToolMATLAB集成了一套功能強大的濾波器設(shè)計工具FDATool(Filter Design & Analysis Tool,可以完成多種濾波器的設(shè)計、分析和性能評估。在MATLAB主命令窗口內(nèi)鍵入“fdatool”,打開FDATool程序界面,如圖B.2所示。 1.3 濾波器設(shè)計F

3、DATool界面左下側(cè)排列了一組工具按鈕,其功能分別如下所述:濾波器轉(zhuǎn)換(TransForm Filter設(shè)置量化參數(shù)(Set Quantization Parameters實現(xiàn)模型(Realize Model導(dǎo)入濾波器(Import Filter多速率濾波器(Multirate Filter零極點編輯器(Pole-zero Editor設(shè)計濾波器(Design Filter選擇其中的選擇Design Filter按鈕,進入設(shè)計濾波器界面,進行下列選擇,如圖B.3所示。濾波器類型(Filer Type為低通(Low Pass設(shè)計方法(Design Method為FIR/IIR,分別采用Equi

4、ripple、Least-squares、Window、Constr.LeastPth-norm 、Constrained Equiripple、Constr.Band Equiripple(FIR濾波器設(shè)計和 Butterworth、Chebyshev Type I、Chebyshev Type I、Elliptic、Maximally flat、Least Pth-norm、Constr.LeastPth-norm(IIR濾波器設(shè)計。濾波器階數(shù)(Filter order定制為15Fs =8kHz,Fpass=3kHz ,Fstop=3050hz最后單擊Design Filter 圖標,讓M

5、ATLAB 計算FIR 濾波器系數(shù)并作相關(guān)分析。其系統(tǒng)函數(shù)H(z可用下式來表示: 顯然上式可以寫成: 即可以看成是一個15 階的FIR 濾波器的輸出結(jié)果經(jīng)過了一個單位延時單元,所以在FDATool 中,把它看成15 階FIR 濾波器來計算參數(shù)。1.4濾波器分析計算完FIR 濾波器系數(shù)以后,往往需要對設(shè)計好的FIR 濾波器進行相關(guān)的性能分析,以便了解該濾波器是否滿足設(shè)計要求。分析操作步驟如下:選擇選擇工具欄中“Magnitude and Phase Response”,啟動幅、相頻比較分析;選擇工具欄中的圓圖圖標,啟動零極點位置分析(主要分析濾波器的穩(wěn)定性,但是通過濾波器基本信息欄也可以直接看濾

6、波器的穩(wěn)定性。除此之外,工具欄 的不同圖標還可以對濾波器進行群延時、相位延時、單位脈沖響應(yīng)、單位階躍響應(yīng)、濾波器的系數(shù)的分析(此實驗中沒有對各方面都展開分析。在圖形的左側(cè),顯示了當(dāng)下濾波器的基本信息,如: 1.5量化可以看到,FDATool計算出的值是一個有符號的小數(shù),如果建立的FIR濾波器模型需要一個整數(shù)作為濾波器系數(shù),就必須進行量化,并對得到的系數(shù)進行歸一化。為此,單擊FDATool 左下側(cè)的工具按鈕進行量化參數(shù)設(shè)置。量化參數(shù)有三種方式:雙精度、單精度和定點。(注意:當(dāng)使用DSP處理器時,FTR濾波器的系數(shù)要量化,量化誤差將導(dǎo)致所設(shè)計的濾波器形狀發(fā)生很大的變化;對于IIR濾波器,會改變?yōu)V波

7、器的特性,甚至導(dǎo)致濾波器不穩(wěn)定。1.6 導(dǎo)出濾波器系數(shù)為導(dǎo)出設(shè)計好的濾波器系數(shù),選擇FDATool菜單的“File”“Export”命令,打開Export(導(dǎo)出對話框,選擇導(dǎo)出到工作區(qū)(Workplace。這時FIR濾波器系數(shù)就存入到一個一維變量Num(名稱按照順序取為a-g中了,不過這時Num中的元素是以小數(shù)形式出現(xiàn)的。對于IIR濾波器,導(dǎo)出得到矩陣SOS(名稱按照順序取為h、k-p和G(名稱按照順序取為h1、k1-p1。導(dǎo)出后的數(shù)據(jù)在workplace中的顯示如下: 其中FIR濾波器系數(shù) 對于IIR濾波器:b,a=sos2tf(SOS;得到濾波器系統(tǒng)函數(shù)的分子分母多項式系數(shù),只是差了一個比

8、例因子k。K=cumprod(G; k=K(end;濾波輸出結(jié)果filteredpWave可以這樣求:filteredpWave=filter(b,a,pWave*k;導(dǎo)出的數(shù)與傳輸函數(shù)的系數(shù)的關(guān)系a,b含義,在MATLAB中函數(shù)如下公式: 這種用法主要出現(xiàn)在信號處理方向,控制方向上主要使用z,而不是z-1,這是習(xí)慣問題,都一樣。其次,了解了其具體函數(shù),則它的傳遞函數(shù)就好寫了,如下:h = tf(b,a,0.1,'variable','z-1' %0.1表示采樣時間,你可以根據(jù)自己設(shè)計的濾波器來設(shè)定。h表示的傳遞函數(shù)就是上面公式中的H(z。當(dāng)然,你也可以用“Di

9、screte Transfer Fcn”模塊來建立傳遞函數(shù),由于這個模塊是直接以Z的形式,所以a,b系統(tǒng)應(yīng)該倒過來寫,相信這個你應(yīng)該理解。 按照上述操作,得到的濾波器設(shè)計相關(guān)圖形如下:FIR濾波器設(shè)計濾波器幅度、相位響應(yīng)濾波器零極點圖形濾波器單位脈沖響應(yīng)Equiripple Least-squares Window(Kaiser Beta=0.5 濾波器幅度、相位響應(yīng)濾波器零極點圖形濾波器單位脈沖響應(yīng)Constr.LeastPth-norm Constrained Equiripple Constr.Band Equiripple II R濾波器設(shè)計濾波器幅度、相位響應(yīng)濾波器零極點圖形濾波器單

10、位脈沖響應(yīng)Butterworth(min order175 Chebyshev Type I(min order 32 Chebyshev Type II(min order 32 Elliptic(min order 12 Maximally flat LeastxPth-norm、 Constr.Least Pth-norm 雖然以上的IIR濾波器的極點都是在單位圓上,但是上述濾波器都是穩(wěn)定的,從其基本信息可以看出 圖形數(shù)據(jù)分析:1.通過觀察FIR濾波器的“幅度相位圖”及“脈沖響應(yīng)圖”可以發(fā)現(xiàn),脈沖響應(yīng)關(guān)于中點對稱的FIR濾波器無相位失真。對稱的脈沖響應(yīng)產(chǎn)生了通帶內(nèi)的線性相位關(guān)系。這樣,線

11、性相位關(guān)系也保證無相位失真;但是IIR濾波器的相位均是非線性的,即其不能保證無相位失真。2.對于給定的濾波器形狀,IIR濾波器較FIR濾波器所需要的系數(shù)要少。3.從“零極點圖”可以發(fā)現(xiàn),濾波器的系數(shù)為N個,對應(yīng)的有N-1個零點和N-1個極點,FIR濾波器的極點都位于單位圓內(nèi)(原點,即所有的FIR濾波器都是穩(wěn)定的。而IIR(它是遞歸濾波器的極點由分母多項式所確定,極點分布在單位圓內(nèi)及單位圓上,這意味著不能保證其穩(wěn)定性,但是從基本信息處可以看到上面所討論的幾種IIR濾波器是穩(wěn)定的。4.通常的幾種模擬濾波器的特性:Butterworth濾波器在通帶和阻帶都是單調(diào)的,意味著它們在一個方向上平滑變化;C

12、hebyshev Type I濾波器在阻帶內(nèi)是單調(diào)的,但是在通帶內(nèi)是有波紋的;Chebyshev Type II濾波器在通帶內(nèi)是單調(diào)的,但是在阻帶內(nèi)是有波紋的;Elliptic濾波器在通帶和阻帶內(nèi)都有波紋。5.對于同樣的設(shè)計要求,切比雪夫濾波器(min order 32比巴特沃斯濾波器(min order175的階數(shù)要低得多。補充說明: 實驗二利用函數(shù)法設(shè)計FIR及IIR濾波器一、實驗?zāi)康?.熟悉掌握設(shè)計濾波器的操作步驟2.學(xué)會用頻率響應(yīng)采樣法和窗函數(shù)法設(shè)計FIR濾波器,并比較其性能3.學(xué)會用單位脈沖響應(yīng)法和雙線性法設(shè)計IIR濾波器,并比較其性能二、實驗條件Matlab 2013b pc機三、

13、實驗過程濾波器指標:Fs=8000hz;fpass=3000hz;fstop=3050;s=0.1;p=0.1(1窗函數(shù)法設(shè)計FIR低通濾波器(主要可以參見實驗四設(shè)計步驟1.在過渡帶寬度的中間,選擇通帶邊緣頻率(Hz:f1=fpass+(fstop-fpass/2=30252.計算1=2*pi*f1/Fs,并將此值代入理想低通濾波器的脈沖響應(yīng)h1n中:h1n=sin(n*1/(n*pi= sin(n*0.75625*pi/(n*pi3.計算濾波器阻帶衰減為Dbs=|20lg(0.1|=20通帶邊緣增益為Dbp=20lg(1-0.1=-0.91從表1中選擇滿足阻帶衰減及其他濾波器要求的窗函數(shù)矩形

14、窗,用表中的公式計算所需窗的非零項目數(shù)目N=0.91*Fs/(fstop-fpass=145.6選擇奇數(shù)項145,這樣脈沖響應(yīng)完全可以對稱,避免了濾波器產(chǎn)生相位失真,對于|n|<=(N-1/2,設(shè)計窗函數(shù)wn。 表14.對于|n|<=(N-1/2,從式hn=h1nwn計算(有限脈沖響應(yīng),對于其他n值hn=0,此脈沖響應(yīng)是非因果的。5.將脈沖響應(yīng)右移(N-1/2,確保第一個非零值在n=0處,使此低通濾波器為因果的。編寫代碼如下:%LPF design 窗函數(shù)法clc;clear all;close all;omga_d=3025*pi/4000; %理想低通濾波器的頻率omga=0:

15、pi/100:pi;N=145;w1=window(rectwin,N;%矩形窗subplot(411;plot(1:N,w1;axis(1 N 0 1.1;%取值范圍不同,矩形窗的中心位置不同title('窗函數(shù)時域形狀'legend('rectwin'%標注窗函數(shù)名稱M=floor(N/2;%M=72n=1:M;%hd=sin(n*omga_d./(n*omga_d*omga_d/pi;hd=sin(n*omga_d./(n*pi;%理想低通濾波器對應(yīng)的脈沖響應(yīng)的右邊序列hdd=fliplr(hd,1/omga_d,hd; %繪制原脈沖函數(shù)圖像,關(guān)于(N+1

16、/2對稱subplot(412;plot(1:N,hdd;axis(1 N -0.2 0.45;title('脈沖函數(shù)hd(n'h_d1=hdd.*w1'%被窗函數(shù)截取后的脈沖函數(shù)m=1:M;H_d1=2*cos(omga'*m*h_d1(M+2:N'+h_d1(M+1;%對截取后的脈沖函數(shù)做newdtft運算得到Hsubplot(413;plot(omga,H_d1;grid on;axis(0 pi 0 1;title('窗函數(shù)法得濾波器形狀H'subplot(414;plot(abs(fft(h_d1'title('

17、fft求得濾波器形狀'圖形顯示: 雖然矩形窗得到的低通濾波器有較好的特性,但是由于其邊緣陡峭,旁瓣出現(xiàn)振鈴,且變化較為明顯。(2頻率響應(yīng)采樣法設(shè)計FIR低通濾波器。設(shè)計原理及濾波器性能分析:.頻率采樣法是從頻域出發(fā),對給定的理想濾波器的頻響進行N點等間隔采樣 然后以此Hd(k作為實際FIR濾波器的頻率特性采樣值H(k,即令: 由DFT定義,可以用這N個頻域的采樣值H(k來惟一確定FIR的單位脈沖響應(yīng)h(n,即: 圖示如下: 對設(shè)計出的濾波器頻率響應(yīng)特性進行分析。由頻域采樣定理中的內(nèi)插公式可以知道,利用這N個頻域采樣值H(k同樣可以求得FIR濾波器的頻率響應(yīng),這個將逼近理想濾波器的頻響。

18、的內(nèi)插公式為: (1式中(w是內(nèi)插函數(shù): (2把(2式帶入(1式,化簡后可得: (3從式(3可以看到,在各頻率采樣點上,設(shè)計的濾波器,實際的頻率響應(yīng)嚴格地與理想濾波器的頻率響應(yīng)數(shù)值相等,即 但是在采樣點之間的頻率響應(yīng)是由各采樣點的加權(quán)內(nèi)插函數(shù)疊加而形成的,因而有一定的逼近誤差。該誤差大小取決于理想頻率響應(yīng)的形狀,理想頻響特性變化越平緩,內(nèi)插值越接近理想值,逼近誤差越小;反之,如果采樣點之間的理想頻響特性變化越陡,則內(nèi)插值與理想值之間的誤差越大,因而在理想濾波器不連續(xù)點的兩邊,就會產(chǎn)生尖峰,而在通帶和阻帶就會產(chǎn)生波紋。線性相位條件:FIR濾波器的最大優(yōu)點是嚴格的線性相位特性。下面討論為實現(xiàn)線性相

19、位,在頻域內(nèi)采樣得到的Hd(k應(yīng)滿足什么條件。FIR濾波器具有線性相位的條件是h(n是實序列,且滿足h(n=±h(N-1-n,即h(n關(guān)于對稱,其中N為濾波器的長度。以第一類線性相位條件h(n=h(N-1-n(偶對稱為基礎(chǔ)來推導(dǎo)頻域采樣Hd(k滿足的條件。理想濾波器的頻率響應(yīng)可以表示為: (4為實現(xiàn)第一類線性相位條件,相位函數(shù)(和幅度函數(shù)Hg(應(yīng)分別滿足: (5 (6在w=02*pi之間對進行N點等間隔采樣,即,將w 代入式(4、(5、(6并表示為k的函數(shù),可以得到: (7 (8 (9式(8和式(9就是頻率采樣值Hd(k滿足第一類線性相位的條件。設(shè)計步驟:1. 確定希望逼近的理想濾波

20、器的頻率響應(yīng)2. 在頻域內(nèi)對進行k點等間隔采樣,利用頻率采樣設(shè)計公式求頻率采樣值Hd(k,采樣間隔=/k=1/512,這樣在通帶內(nèi)共有c個采樣點,分別是k=0,1,c-1。3. 用離散傅里葉逆變換求得要設(shè)計的實際濾波器的單位脈沖響應(yīng)h(n: 4. 根據(jù)傅里葉變換的定義求得實際濾波器的頻率響應(yīng),驗證是否滿足濾波器技術(shù)指標的要求,主要驗證濾波器的阻帶衰減是否能夠滿足阻帶的要求。Matlab代碼中需要用到的函數(shù)知識>> help linspacelinspace Linearly spaced vector.linspace(X1, X2 generates a row vector o

21、f 100 linearlyequally spaced points between X1 and X2.linspace(X1, X2, N generates N points between X1 and X2.For N = 1, linspace returns X2.編寫代碼如下:%LPF design 頻率響應(yīng)采樣法%待設(shè)計濾波器的基本參數(shù)f_p=3000;%所設(shè)計濾波器的通帶邊緣頻率f_s=3050;%阻帶邊緣頻率Fs=8000;%采樣頻率%All in Hzdelta_p=0.1;%通帶波紋delta_s=0.1;%阻帶波紋omg_p=2*pi*f_p/Fs;%通帶邊緣模擬

22、頻率轉(zhuǎn)化為數(shù)字頻率omg_s=2*pi*f_s/Fs;%阻帶邊緣模擬頻率轉(zhuǎn)化為數(shù)字頻率K=512;%理想濾波器長度M=56;%單位脈沖響應(yīng)單邊長度C=floor(omg_p+omg_s/2/pi*K;%理想濾波器中有值部分長度?Mag=ones(1,C,zeros(1,K-C;%幅度phi=linspace(0,pi,K*M/2;%相位,范圍為0-pi,間隔為pi/K,線性相位特性omg*M/2 H=Mag.*exp(-1j*phi;%|H|*ej*omg*M/2,待設(shè)計的理想濾波器H=H,conj(fliplr(H(1:K-1;%將H轉(zhuǎn)換為因果的h0=fftshift(real(ifft(H

23、;h=h0(K+M/2-M:K+M/2+M;%相移ej*omg*M/2,時域采樣平移M/2個點subplot(211;plot(h;%脈沖響應(yīng)axis(1 120 -0.2 0.8;title('脈沖函數(shù)hd(n'subplot(212;plot(abs(fft(h;%所設(shè)計的濾波器的形狀title('所設(shè)計的濾波器形狀'圖像顯示: 增加過渡帶抽樣點,可加大阻帶衰減,但是導(dǎo)致過渡帶變寬增加N,使抽樣點變密,減小過渡帶帶寬,但是增加了計算量優(yōu)點:頻域直接計算缺點:抽樣頻率只能是2/N或者/N的整數(shù)倍,截止頻率不能任意取值(3單位脈沖響應(yīng)法設(shè)計IIR濾波器核心原理:

24、通過對連續(xù)函數(shù)ha(t等間隔采樣得到離散序列ha(nT。令h(n=ha(nT, T為采樣間隔。它是一種時域上的轉(zhuǎn)換方法。設(shè)計步驟:1. 由模擬濾波器的傳輸函數(shù)Ha(s求出其對應(yīng)的模擬沖激響應(yīng)ha(t2. 對ha(t每T秒進行一次采樣,采樣時間為t=nT,得到數(shù)字脈沖響應(yīng)h(n=ha(nT3. 由h(n進行z變換可以求出數(shù)字濾波器的傳輸函數(shù)H(z4. 將z用ej*Omga替換得到數(shù)字濾波器的頻率響應(yīng)H(Omga編寫代碼如下:(用要求的參量得不到需要設(shè)計的濾波器,這里更換了參數(shù)%LPF design 脈沖響應(yīng)不變法(Butterworthfp=100;fs=300;Rp=0.91;%通帶衰減系數(shù)R

25、s=20;%阻帶衰減系數(shù)Fs=1000;wp=2*pi*fp;ws=2*pi*fs;N,Wn=buttord(wp,ws,Rp,Rs,'s'%求濾波器所需階數(shù)及截止頻率B,A=butter(N,Wn,'s'%計算濾波器傳輸函數(shù)系數(shù)W=0:pi/100:pi;H,W=freqs(B,A,W;%濾波器頻率響應(yīng)%H=20*log10(abs(H;plot(W,H,grid on;xlabel('rad'ylabel('幅度/db'title('所設(shè)計的濾波器形狀'axis(0 pi 0 1.1;運行結(jié)果N = 3 Wn

26、= 876.3856B =1.0e+08 *0 0 0 6.7311A = 1.0e+08 *0.0000 0.0000 0.0154 6.7311圖形顯示: (4雙線性法設(shè)計IIR濾波器設(shè)計步驟1.確定待求通帶邊緣頻率fpass=3000hz,待求阻帶邊緣頻率fstop=3050hz,待求阻帶衰減20db,通帶邊緣增益-0.91db。2.用式子=2*f/Fs把由hz表示的待求邊緣頻率轉(zhuǎn)化為由弧度表示的數(shù)字頻率,得到pass、stop。3.計算預(yù)扭曲模擬頻率以避免雙線性變換帶來的失真。由w=2Fs*tan(/2求得wpass和wstop單位是弧度/秒。4.計算所需濾波器的階數(shù)5.把wpass帶

27、入n階模擬巴特沃斯濾波器傳輸函數(shù)H(s中,并且對其進行雙線性變換s=2Fs*(z-1/(z+1得到n階數(shù)字傳輸函數(shù)H(z。6.由H(z求出濾波器實現(xiàn)所需的差分方程。7.把w=2Fs*tan(/2代入 即可得濾波器的形狀|H(|。編寫代碼如下:%LPF design 雙線性變換法(Butterworth%待設(shè)計濾波器的基本參數(shù)fpass=3000;fstop=3050;Fs=8000;%單位 hz Rs=0.1; Rp=0.1; omg_pass=2*pi*fpass/Fs omg_stop=2*pi*fstop/Fs wpass=2*Fs*tan(omg_pass/2 %阻帶波紋 %通帶波紋

28、%通帶截止頻率數(shù)字頻率 %阻帶截止頻率數(shù)字頻率，單位 rad wstop=2*Fs*tan(omg_stop/2 %預(yù)扭曲模擬頻率，單位 rad/s n=ceil(log10(1/(Rs2-1/(2*log10(wstop/wpass %所需濾波器階數(shù),ceil 向上取整 omg=0:pi/100:pi; H=1./sqrt(2*Fs*tan(omg./2./wpass.(2*n+1; plot(omg,abs(H;grid on; title('所設(shè)計的濾波器形狀' axis(0 pi 0 1.1; xlabel('數(shù)字頻率/（rad）'ylabel('幅度' 計算結(jié)果： omg_pass =2.3562 omg_stop =2.3955 wpass =3.8627e+04 wstop =4.0880e+04 n =41 圖形顯示： 濾波特性較好，通帶內(nèi)沒有波紋，邊緣較陡峭 分析比較 1. 對比實驗