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

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

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

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

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

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

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

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

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

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

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

12、將模擬濾波器轉變為數字低通濾波器利用模擬橢圓濾波器設計數字濾波器結束圖2 濾波器設計流程方框圖四程序源代碼MATLAB 程序如下：clear all close all Coeff_LF = 1; % 低頻信號系數 F1(300 Hz)Coeff_HF = 1; % 高頻信號系數 F2(600 Hz) Coeff_NOISE = 0; % 噪聲系數, 如0(沒有噪聲), 1, sqrt(2), 3等等. F1 = 300; % HzF2 = 600; % HzFS = 2000; % Hz 抽樣頻率 % 時間定義EndTime = 1/FS * 1023;Time = 0:1/FS:EndTi

13、me;N = length(Time); % 產生正弦信號 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)xla

14、bel('時間')xlim(0 0.5)title('時域信號 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, FnSi

15、n_x1(2:length(FW)hold onplot(FW(2:length(FW)-2*pi)*FS/pi/2, FnSin_x1(2: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

16、)title('時域信號 x2(t)') FFTSinSignal_x2 = 2/N*fft(SinSignal_x2)subplot(3,1,3)FnSin_x2 = abs(FFTSinSignal_x2); plot(FW(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)p

17、lot(Time, SinSignal)xlabel('時間')xlim(0 0.5)title('時域信號') subplot(3,1,2)plot(Time, SinSignal)xlabel('時間')xlim(0 0.05)title('時域信號 x(t)') FFTSinSignal = 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

18、(FW(2:length(FW)-2*pi)*FS/pi/2, FnSin(2:length(FW)xlim(-FS/2 FS/2) xlabel('頻率') %-設計低通濾波器 - Rp = 1; % 通帶衰減 /dBRs = 40; % 阻帶衰減 /dBOmegaP1_1=400; % 通帶邊界頻率OmegaS1_1=500; % 阻帶邊界頻率FS=2000; % 抽樣頻率 Wp1=2*pi*OmegaP1_1/FS; % 將模擬頻率變為數字角頻率Ws1=2*pi*OmegaS1_1/FS; % 將模擬頻率變為數字角頻率OmegaP1=tan(Wp1/2); % 非線性變換

19、OmegaS1=tan(Ws1/2); % 非線性變換Eta_P=OmegaP1/OmegaP1; % 歸一化Eta_S=OmegaS1/OmegaP1; % 歸一化 Lemta_P_EquivalentLowPass=1/Eta_P; % 變為等效低通參數Lemta_S_EquivalentLowPass=1/Eta_S; % 變為等效低通參數 % 估計濾波器階數N, Wn = Ellipord(Lemta_P_EquivalentLowPass, Lemta_S_EquivalentLowPass,Rp, Rs,'s'); % 時間濾波器num1,den1 = ellip(

20、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 5);xlabel('omega /pi'); ylabel('增益 /dB');title('橢圓低通濾波器的增益響應'); % 繪制零點和極

21、點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-五仿真結果圖 六 參考資料1程佩青著.數字信號處理教程