用雙線性變換法設(shè)計原型低通為橢圓型的數(shù)字IIR帶通濾波器

1、 課 程 設(shè) 計課程設(shè)計名稱：用雙線性變換法設(shè)計原型低通為 橢圓型的數(shù)字IIR帶通濾波器 專 業(yè) 班 級 ： xxxxxx 學 生 姓 名 ： xxx 學 號 ： xxxxxxxxxxxx 指 導 教 師 ： xxx 課程設(shè)計時間： 2014.6.16至2014.6.20 電子信息工程 專業(yè)課程設(shè)計任務(wù)書學生姓名XXX專業(yè)班級XXXX 學號XXXXXX題 目用雙線性變換法設(shè)計原型低通為橢圓型的數(shù)字IIR低通濾波器課題性質(zhì)其他課題來源自擬指導教師XX同組姓名主要內(nèi)容用雙線性變換法設(shè)計原型低通為橢圓型的數(shù)字IIR低通濾波器，要求通帶邊界頻率為400Hz，阻帶邊界頻率分別為500Hz，通帶最大衰減1

2、dB，阻帶最小衰減40dB，抽樣頻率為2000Hz，用MATLAB畫出幅頻特性，畫出并分析濾波器傳輸函數(shù)的零極點；信號經(jīng)過該濾波器，其中300Hz，600Hz，濾波器的輸出是什么？用Matlab驗證你的結(jié)論并給出的圖形。任務(wù)要求1、掌握用雙線性變換法設(shè)計原型低通為橢圓型的數(shù)字IIR低通濾波器的原理和設(shè)計方法。2、求出所設(shè)計濾波器的Z變換。3、用MATLAB畫出幅頻特性圖。4、驗證所設(shè)計的濾波器。參考文獻1、程佩青著，數(shù)字信號處理教程，清華大學出版社，20012、Sanjit K. Mitra著，孫洪，余翔宇譯，數(shù)字信號處理實驗指導書（MATLAB版），電子工業(yè)出版社，2005年1月3、郭仕劍等

3、，MATLAB 7.x數(shù)字信號處理，人民郵電出版社，2006年4、胡廣書，數(shù)字信號處理 理論算法與實現(xiàn)，清華大學出版社，2003年審查意見指導教師簽字：教研室主任簽字： 2014 年06 月 12 日 說明：本表由指導教師填寫，由教研室主任審核后下達給選題學生，裝訂在設(shè)計（論文）首頁一需求分析和技術(shù)要求 1、需求分析濾波器從廣義上來說對特定的頻點或頻點以外的頻率進行有效濾波的電路，這種電路保留輸入信號中的有用信息，濾除不需要的信息，從而達到信號的檢測、提取、識別等不同的目的。如果處理的信號是時域離散信號，那么相應的處理系統(tǒng)就稱為數(shù)字濾波器，由于在實際工作中被處理的信號都是幅度量化的數(shù)字信號，因

4、此，數(shù)字濾波器實際上是用有限精度的算法實現(xiàn)一個線性時不變的時域離散系統(tǒng)。數(shù)字濾波器的種類很多，分類方法也不同，可以從功能上分類，也可以從實現(xiàn)方法上及設(shè)計方法上來分類等等。濾波器在功能上總的可分為四類，即低通(LP)、高通(HP)、帶通(BP)、帶阻(BS)濾波器等，從實現(xiàn)方法上，由有限長沖激響應的數(shù)字濾波器被稱為FIR濾波器，具有無限長沖激響應的數(shù)字濾波器稱作IIR濾波器。橢圓濾波器（Elliptic filter）又稱考爾濾波器（Causer filter），是在通帶和阻帶等波紋的一種濾波器。橢圓濾波器相比其他類型的濾波器，在階數(shù)相同的條件下有著最小的通帶和阻帶波動。它在通帶和阻帶的波動相同

5、，這一點區(qū)別于在通帶和阻帶都平坦的巴特沃斯濾波器，以及阻帶平坦、通帶等波紋的切比雪夫濾波器。利用雙線性變換法將模擬傳輸信號Ha(s)變換為數(shù)字傳輸函數(shù)G(z)，從而是z域的數(shù)字傳輸函數(shù)保留s域的模擬傳輸函數(shù)的基本性質(zhì)。設(shè)計成的IIR數(shù)字低通濾波器能夠去掉信號中不必要的高頻成分，降低采樣頻率，避免頻率混淆，去掉高頻干擾。2、技術(shù)要求 1、掌握用雙線性變換法設(shè)計原型低通為橢圓型的數(shù)字IIR低通濾波器的原理和設(shè)計方法。2、求出所設(shè)計濾波器的Z變換。3、用MATLAB畫出幅頻特性圖。4、驗證所設(shè)計的濾波器。二設(shè)計原理及設(shè)計思路1、設(shè)計原理數(shù)字濾波器是一種用來過濾時間離散信號的數(shù)字系統(tǒng)，通過對抽樣數(shù)據(jù)進

6、行數(shù)學處理來達到頻域濾波的目的。數(shù)字濾波器也是具有一定傳輸選擇特性的數(shù)字信號處理裝置，其輸入、輸出均為數(shù)字信號，實質(zhì)上是一個由有限精度算法實現(xiàn)的線性時不變離散系統(tǒng)。IIR數(shù)字濾波器采用遞歸型結(jié)構(gòu)，即結(jié)構(gòu)上帶有反饋環(huán)路。IIR濾波器運算結(jié)構(gòu)通常由延時、乘以系數(shù)和相加等基本運算組成，可以組合成直接型、正準型、級聯(lián)型、并聯(lián)型四種結(jié)構(gòu)形式，都具有反饋回路。數(shù)字濾波器根據(jù)其沖激響應函數(shù)的時域特性，可分為兩種，即無限長沖激響應(IIR)數(shù)字濾波器和有限長沖激響應(FIR)數(shù)字濾波器。IIR 數(shù)字濾波器的特征是，具有無限持續(xù)時間沖激響應，需要用遞歸模型來實現(xiàn)，其差分方程為： （2.1） （2.2）設(shè)計IIR

7、濾波器的任務(wù)就是尋求一個物理上可實現(xiàn)的系統(tǒng)函數(shù)H(z)，使其頻率響應H(z)滿足所希望得到的頻域指標，即符合給定的通帶截止頻率、阻帶截止頻率、通帶衰減系數(shù)和阻帶衰減系數(shù)。本次課程設(shè)計采用雙線性變換法。雙線性變換法設(shè)計IIR數(shù)字濾波器雙線性變換法主要是采用非線性頻率壓縮方法，將整個頻率軸上的頻率范圍壓縮到-/T/T之間，再用z=esT轉(zhuǎn)換到Z平面上。也就是說，第一步先將整個S平面壓縮映射到S1平面的-/T/T一條橫帶里；第二步再通過標準變換關(guān)系z=es1T將此橫帶變換到整個Z平面上去。這樣就使S平面與Z平面建立了一一對應的單值關(guān)系，消除了多值變換性，也就消除了頻譜混疊現(xiàn)象，映射關(guān)系如圖1所示。

8、圖1雙線性變換的映射關(guān)系為了將S平面的整個虛軸j壓縮到S1平面j1軸上的-/T到/T段上，可以通過以下的正切變換實現(xiàn) (2.3)式中,T仍是采樣間隔。當1由-/T經(jīng)過0變化到/T時，由-經(jīng)過0變化到+，也即映射了整個j軸。將式（2.3）寫成(2.4)將此關(guān)系解析延拓到整個S平面和S1平面，令j=s，j1=s1，則得 （2.5）再將S1平面通過以下標準變換關(guān)系映射到Z平面z=es1T (2.6) 從而得到S平面和Z平面的單值映射關(guān)系為： （2.7） （2.8） 式（2.7）與式（2.8）是S平面與Z平面之間的單值映射關(guān)系，這種變換都是兩個線性函數(shù)之比，因此稱為雙線性變換式（2.6）與式（2.7）

9、的雙線性變換符合映射變換應滿足的兩點要求。首先,把z=ej，可得 （2.9）即S平面的虛軸映射到Z平面的單位圓。其次，將s=+j代入式（2.8），得 (2.10)因此(2.11)由此看出，當0時，|z|0時，|z|1。也就是說，S平面的左半平面映射到Z平面的單位圓內(nèi)，S平面的右半平面映射到Z平面的單位圓外，S平面的虛軸映射到Z平面的單位圓上。因此，穩(wěn)定的模擬濾波器經(jīng)雙線性變換后所得的數(shù)字濾波器也一定是穩(wěn)定的。2、設(shè)計思路具體實現(xiàn)方法如下：（1）確定數(shù)字低通濾波器的技術(shù)指標：通帶截止頻率p、通帶衰減Rp、阻帶截止頻率s、阻帶衰減Rs。（2）將數(shù)字低通濾波器的技術(shù)指標轉(zhuǎn)換成模擬低通濾波器的技術(shù)指標

10、。如果采用脈沖響不變法，邊界頻率的轉(zhuǎn)換關(guān)系為: (2.12)如果采用雙線性變換法，邊界頻率的轉(zhuǎn)換關(guān)系為(2.13)（3）按照模擬低通濾波器的技術(shù)指標設(shè)計模擬低通濾波器。（4）利用雙線性變換法將模擬濾波器Ha(s)，從s平面轉(zhuǎn)換到z平面，得到數(shù)字低通濾波器系統(tǒng)函數(shù)H(z)。（5）數(shù)字低通技術(shù)指標為：p=0.4rad， Rp=1dB； s=0.5rad， Rs=40dB（6）模擬低通的技術(shù)指標為：通帶截止頻率p=(2/T)*tan(p/2)=tan(0.392)=0.0068阻帶截止頻率s=(2/T)*tan(s/2)=tan(0.628)=0.011通帶最小衰減Rp=1dB和阻帶最大衰減Rs=4

11、0dB；歸一化截止角頻率p=2pi*Fs/Ft； s=2pi*Fs/Ft（7） 利用模擬橢圓濾波器設(shè)計數(shù)字濾波器。通帶截止頻率為：p=0.4*pi；阻帶截止頻率為：s=0.5*pi；通帶最大衰減為：Rp=1；阻帶最大衰減為：Rs=40；設(shè)定周期為1s；模擬低通濾波器的生成：num1,den1 =ellip(N,Rp,Rs,Wn,s)；最后實現(xiàn)輸入輸出、幅頻特性、相頻特性的圖形。三程序流程圖確定數(shù)字低通濾波器的各項性能指標開始數(shù)字低通濾波器的技術(shù)指標轉(zhuǎn)換成模擬低通濾波器的技術(shù)指標按照模擬低通濾波器的技術(shù)指標設(shè)計模擬低通濾波器利用雙線性變換法將模擬濾波器轉(zhuǎn)變?yōu)閿?shù)字低通濾波器利用模擬橢圓濾波器設(shè)計數(shù)

12、字濾波器結(jié)束圖2 濾波器設(shè)計流程方框圖四程序源代碼MATLAB 程序如下：clear all close all Coeff_LF = 1; % 低頻信號系數(shù) F1(300 Hz)Coeff_HF = 1; % 高頻信號系數(shù) F2(600 Hz) Coeff_NOISE = 0; % 噪聲系數(shù), 如0(沒有噪聲), 1, sqrt(2), 3等等. F1 = 300; % HzF2 = 600; % HzFS = 2000; % Hz 抽樣頻率 % 時間定義EndTime = 1/FS * 1023;Time = 0:1/FS:EndTime;N = length(Time); % 產(chǎn)生正弦信

13、號 x1(t),x2(t)和 x(t)SinSignal_x1= Coeff_LF*sin(2*pi*F1*Time)% x1(t)SinSignal_x2= Coeff_HF*sin(2*pi*F2*Time)% x2(t)SinSignalSum = SinSignal_x1 + SinSignal_x2+ Coeff_NOISE*randn(1, N);SinSignal = SinSignalSum;% x(t) % 繪制信號 x1(t)figuresubplot(3,1,1)plot(Time, SinSignal_x1)xlabel(時間)xlim(0 0.5)title(時域信號

14、 x1(t) subplot(3,1,2)plot(Time, SinSignal_x1)xlabel(時間)xlim(0 0.05)title(時域信號 x1(t) FFTSinSignal_x1 = 2/N*fft(SinSignal_x1)subplot(3,1,3)FnSin_x1 = abs(FFTSinSignal_x1); FW = 2*pi*(Time)/EndTime; plot(FW(2:length(FW)*FS/pi/2, FnSin_x1(2:length(FW)hold onplot(FW(2:length(FW)-2*pi)*FS/pi/2, FnSin_x1(2

15、:length(FW)xlim(-FS/2 FS/2) xlabel(頻率) % 繪制 x2(t)figuresubplot(3,1,1)plot(Time, SinSignal_x2)xlabel(時間)xlim(0 0.5)title(時域信號 x(2) subplot(3,1,2)plot(Time, SinSignal_x2)xlabel(時間)xlim(0 0.05)title(時域信號 x2(t) FFTSinSignal_x2 = 2/N*fft(SinSignal_x2)subplot(3,1,3)FnSin_x2 = abs(FFTSinSignal_x2); plot(FW

16、(2:length(FW)*FS/pi/2, FnSin_x2(2:length(FW)hold onplot(FW(2:length(FW)-2*pi)*FS/pi/2, FnSin_x2(2:length(FW)xlim(-FS/2 FS/2) xlabel(頻率) % 繪制信號 x(t)figuresubplot(3,1,1)plot(Time, SinSignal)xlabel(時間)xlim(0 0.5)title(時域信號) subplot(3,1,2)plot(Time, SinSignal)xlabel(時間)xlim(0 0.05)title(時域信號 x(t) FFTSin

17、Signal = 2/N*fft(SinSignal)subplot(3,1,3)FnSin = abs(FFTSinSignal); plot(FW(2:length(FW)*FS/pi/2, FnSin(2:length(FW)hold onplot(FW(2:length(FW)-2*pi)*FS/pi/2, FnSin(2:length(FW)xlim(-FS/2 FS/2) xlabel(頻率) %-設(shè)計低通濾波器 - Rp = 1; % 通帶衰減 /dBRs = 40; % 阻帶衰減 /dBOmegaP1_1=400; % 通帶邊界頻率OmegaS1_1=500; % 阻帶邊界頻率

18、FS=2000; % 抽樣頻率 Wp1=2*pi*OmegaP1_1/FS; % 將模擬頻率變?yōu)閿?shù)字角頻率Ws1=2*pi*OmegaS1_1/FS; % 將模擬頻率變?yōu)閿?shù)字角頻率OmegaP1=tan(Wp1/2); % 非線性變換OmegaS1=tan(Ws1/2); % 非線性變換Eta_P=OmegaP1/OmegaP1; % 歸一化Eta_S=OmegaS1/OmegaP1; % 歸一化 Lemta_P_EquivalentLowPass=1/Eta_P; % 變?yōu)榈刃У屯▍?shù)Lemta_S_EquivalentLowPass=1/Eta_S; % 變?yōu)榈刃У屯▍?shù) % 估計濾波器階

19、數(shù)N, Wn = Ellipord(Lemta_P_EquivalentLowPass, Lemta_S_EquivalentLowPass,Rp, Rs,s); % 時間濾波器num1,den1 = ellip(N,Rp,Rs,Wn,s);num2,den2=lp2lp(num1,den1,OmegaP1);num,den=bilinear(num2,den2,0.5); % 計算增益響應w = 0:pi/255:pi;h = freqz(num,den,w);g = 20*log10(abs(h); % 繪制增益響應figure;plot(w/pi,g);grid axis(0 1 -60

20、 5);xlabel(omega /pi); ylabel(增益 /dB);title(橢圓低通濾波器的增益響應); % 繪制零點和極點figure;z,p,k=tf2zp(num,den);zplane(z,p);title(零點和極點) figuresubplot(2,1,1)YLP = filter(num,den, SinSignal);plot(Time, YLP)xlim(0 0.05)xlabel(時間)title(低通濾波后的時域信號) FFTSinSignalSum = 2 / N * fft(YLP)subplot(2,1,2)FnSinSum = abs(FFTSinSignalSum); plot(FW(2:length(FW)*FS/pi/2, FnSinSum(2:length(FW)hold onplot(FW(2:length(FW)-2*pi)*FS/pi/2, FnSinSum(2:length(FW)xlim(-FS/2 FS/2) ylim(0 1)xlabel(頻率) %return%-End Program-五仿真結(jié)果圖 六 參考資料1程佩青著.數(shù)字信號處理教程M.清華大學出版社，20012Sanjit K. Mi