數字信號處理--用雙線性變換法設計IIR 數字濾波器

1、 數字信號處理實驗內容用雙線性變換法設計IIR 數字濾波器一、實驗目的：1、熟悉用雙線性變換法設計IIR 數字濾波器的原理與方法。2、掌握數字濾波器的計算機仿真方法。3、熟悉Batterworth濾波器設計方法及特點二、實驗原理：（一）、IIR數字濾波器的設計步驟： 按照一定規(guī)則把給定的濾波器技術指標轉換為模擬低通濾波器的技術指標； 根據模擬濾波器技術指標設計為響應的模擬低通濾波器； 跟據脈沖響應不變法和雙線性不變法把模擬濾波器轉換為數字濾波器； 如果要設計的濾波器是高通、帶通或帶阻濾波器，則首先把它們的技術指標轉化為模擬低通濾波器的技術指標，設計為數字低通濾波器，最后通過頻率轉換的方法來得到

2、所要的濾波器。在MATLAB中，經典法設計數字濾波器主要采用以下步驟：模擬濾波器原型buttap,cheb1ap頻率變換模擬離散化bilinear,impinvarIIR數字濾波器IIR數字濾波器設計步驟（二）、用模擬濾波器設計數字濾波器的方法1、沖激響應不變法：沖激響應不變法是從時域出發(fā)，要求數字濾波器的沖激響應h(n) 對應于模擬濾波器h(t) 的等間隔抽樣。優(yōu)點：時域逼近良好；保持線性關系。缺點：頻域響應混疊。只適用于限帶低通濾波器和帶通濾波器2、雙線性變換法優(yōu)點：克服了頻域混疊缺點：高頻時會引起畸變1）沖激響應不變法impinvar格式：BZ,AZ= impinvar（B,A,Fs）功

3、能：把具有B,A模擬濾波器傳遞函數模型轉換為采樣頻率為Fs的數字濾波器的傳遞函數模型BZ,AZ，Fs默認值為1。例：一個4階的Butterworth模擬低通濾波器的系統(tǒng)函數如下：試用沖激響應不變法求出Butterworth模擬低通數字濾波器的系統(tǒng)函數。num=1;den=1,sqrt(5),2,sqrt(2),1;num1,den1=impinvar(num,den)2）雙線性變換法bilinear格式一：Zd,Pd,Kd= bilinear(Z,P,K,Fs)功能：把模擬濾波器的零極點模型轉換成數字濾波器的零極點模型，Fs是采樣頻率格式二：numd,dend= bilinear(num,de

4、n,Fs)功能：把模擬濾波器的傳遞函數模型轉換為數字濾波器的傳遞函數模型。例：一個三階的模擬Butterworth模擬低通濾波器的系統(tǒng)函數如下：，試用雙線性變換法求出數字Butterworth數字低通濾波器的系統(tǒng)函數。num=1;den=1,sqrt(3),sqrt(2),1;num1,den1=bilinear(num,den,1)3） IIR數字濾波器的頻率變換實現步驟： 按一定的規(guī)則將數字濾波器的技術指標轉換為模擬低通濾波器的技術指標 根據轉換后的技術指標使用濾波器階數函數，確定濾波器的最小階數N和截止頻率Wc 利用最小階數N產生模擬低通原型 利用截止頻率Wc把模擬低通濾波器原型轉換為模

5、擬低通、高通、帶通、帶阻濾波器 利用沖激響應不變法或雙線性變換法把模擬濾波器轉換為數字濾波器表一 IIR濾波器階次估計函數名功能說明buttord計算Butterworth濾波器的階次及截止頻率cheb1ord計算Chebyshev濾波器的階次cheb2ord計算Chebyshev濾波器的階次ellipord計算橢圓濾波器的最小階次表二 模擬低通濾波器原型設計函數名功能說明buttapButterworth模擬低通濾波器原型設計cleb1apChebyshev模擬低通濾波器原型設計cheb2apChebyshev模擬低通濾波器原型設計ellipap橢圓模擬低通濾波器原型設計表三 模擬濾波器變換

6、函數函數名功能說明lp2bp模擬低通轉換為帶通lp2bs模擬低通轉換為帶阻lp2hp模擬低通轉換為高通lp2lp改變模擬低通的截止頻率（三）、數字濾波器的設計1 數字濾波器的設計參數濾波器的4個重要的通帶、阻帶參數為：通帶截止頻率（Hz） ：阻帶起始頻率（Hz）：通帶內波動（dB），即通帶內所允許的最大衰減；：阻帶內最小衰減設采樣速率（即奈奎斯特速率）為，將上述參數中的頻率參數轉化為歸一化角頻率參數：通帶截止角頻率（rad/s） ，； ：阻帶起始角頻率（rad/s） ，通過以上參數就可以進行離散濾波器的設計。2、巴特沃斯濾波器設計1） 巴特沃斯濾波器階數的選擇：在已知設計參數，之后，可利用“b

7、uttord”命令可求出所需要的濾波器的階數和3dB截止頻率，其格式為：n，Wn=buttordWp，Ws，Rp，Rs，其中Wp，Ws，Rp，Rs分別為通帶截止頻率、阻帶起始頻率、通帶內波動、阻帶內最小衰減。返回值n為濾波器的最低階數，Wn為3dB截止頻率。2） 巴特沃斯濾波器系數計算：由巴特沃斯濾波器的階數n以及3dB截止頻率Wn可以計算出對應傳遞函數H(z)的分子分母系數，MATLAB提供的命令如下：(a)巴特沃斯低通濾波器系數計算：b，a=butter(n,Wn)，其中b為H(z)的分子多項式系數，a為H(z)的分母多項式系數(b)巴特沃斯高通濾波器系數計算：b，a=butter(n，W

8、n，High) (c)巴特沃斯帶通濾波器系數計算：b，a=butter(n，W1，W2)，其中W1，W2為截止頻率，是2元向量，需要注意的是該函數返回的是2*n階濾波器系數。(d)巴特沃斯帶阻濾波器系數計算：b，a=butter(ceil(n/2)，W1，W2，stop)，其中W1，W2為截止頻率，是2元向量，需要注意的是該函數返回的也是2*n階濾波器系數。三、巴特沃斯濾波器設計實例：例題1：采樣速率為8000Hz，要求設計一個低通濾波器，=2100Hz，=2500Hz，=3dB，=25dB。用直接設計法程序如下：fn=8000;%采樣頻率 fp=2100; %通帶截止頻率fs=2500; %

9、阻帶起始頻率Rp=3; %通帶最大衰減Rs=25;%阻帶最小衰減Wp=fp/(fn/2);%計算歸一化角頻率Ws=fs/(fn/2);n,Wn=buttord(Wp,Ws,Rp,Rs);%計算階數和截止頻率b,a=butter(n,Wn);%計算H(z)分子、分母多項式系數H,F=freqz(b,a,1000,8000);%計算H(z)的幅頻響應,freqz(b,a,計算點數,采樣速率)subplot(2,1,1)plot(F,20*log10(abs(H) %畫出幅頻特性圖xlabel(Frequency(Hz); ylabel(Magnitude(dB) title(低通濾波器)axis(

10、0 4000 -30 3);grid onpha=angle(H)*180/pi;subplot(2,1,2)plot(F,pha);grid on %畫出相頻特性圖xlabel(Frequency(Hz); ylabel(phase);用雙線性變換法wp=2100*2*pi; %利用 ws=2500*2*pi;Rp=3;Rs=25;Fs=8000;Ts=1/Fs;%選擇濾波器的最小階數N,Wn=buttord(wp,ws,Rp,Rs,s);%創(chuàng)建butterworth模擬濾波器Z,P,K=buttap(N);%把濾波器零極點模型轉化為傳遞函數模型Bap,Aap=zp2tf(Z,P,K); %

11、把模擬濾波器原型轉換成截至頻率為Wn的低通濾波器b,a=lp2lp(Bap,Aap,Wn);%用雙線性變換法實現模擬濾波器到數字濾波器的轉換bz,az=bilinear(b,a,Fs);%繪制頻率響應曲線H,W=freqz(bz,az);plot(W*Fs/(2*pi),abs(H);gridxlabel(頻率/Hz)ylabel(幅度)例題2：模擬原型直接變換法設計數字濾波器：已知四階歸一化低通巴特沃斯模擬濾波器系統(tǒng)函數為，編寫MATLAB程序實現從設計3dB截止頻率為，設采樣周期為T=1，的四階低通巴特沃斯數字濾波器。程序如下：步驟一：將設計內容題所給歸一化巴特沃斯低通濾波器以3dB截止頻

12、率為進行去歸一化。步驟二：用雙線性變化法將低通模擬濾波器變換為低通數字濾波器設計程序如下：clear; T=1; fs=1/T; N=4;wc=pi/2; omegach=2*tan(wc/2)/T;%模擬濾波器的截止頻率M=1; N=1,2.6131,3.4142,2.6131,1;h,w=freqs(M,N,512); %模擬濾波器的幅頻響應subplot(2,1,1);plot(w,20*log10(abs(h); axis(0,10,-90,0),grid on; xlabel(Hz);ylabel(幅度); title(歸一化模擬低通濾波器);Ms,Ns=lp2lp(M,N,omeg

13、ach); %對低通濾波器進行頻率變換hs,ws=freqs(Ms,Ns,512); %模擬濾波器的幅頻響應subplot(2,1,2);plot(ws,20*log10(abs(hs); grid;axis(0,10,-90,0);xlabel(Hz);ylabel(幅度); title(去歸一化模擬低通濾波器);Mz,Nz=bilinear(Ms,Ns,1/T); %對模擬濾波器雙線性變換h1,w1=freqz(Mz,Nz); %數字濾波器的幅頻響應figureplot(w1/pi,20*log10(abs(h1); grid;xlabel(/);ylabel(幅度(dB); title(

14、數字低通濾波器);axis(0,1,-160,0)四、實驗內容：1、采樣速率為10000Hz，要求設計一個巴特沃斯帶阻濾波器，=1000Hz，1500Hz，=1200Hz，1300Hz，=3dB，=30dB。提示：b,a=butter(N,Wc,stop)2、 采樣速率為10000Hz，要求設計一個帶通濾波器，=1000Hz，1500Hz，=600Hz，1900Hz，=3dB，=20dB。3、已知四階歸一化低通巴特沃斯模擬濾波器系統(tǒng)函數為，編寫MATLAB程序實現從設計3dB截止頻率為的四階高通巴特沃斯數字濾波器。4、設計低通濾波器，把輸入信號中的三個信號分離出來。要求：畫出濾波前后信號的波形

15、及頻譜及低通濾波器的幅頻響應。提示：b,a=butter(n,Wn);%計算H(z)分子、分母多項式系數y=filter(b,a,x);%對輸入的信號進行濾波五、思考題：用雙線性變換法設計數字濾波器過程中，變換公式中T 的取值對設計結果有無影響？為什么？程序如下：五、實驗結果：1、程序如下：fn=10000;%采樣頻率 fp=1000,1500; %通帶截止頻率fs=600,1900; %阻帶起始頻率Rp=3; %通帶最大衰減Rs=20;%阻帶最小衰減Wp=fp/(fn/2);%計算歸一化角頻率Ws=fs/(fn/2);n,Wn=buttord(Wp,Ws,Rp,Rs);%計算階數和截止頻率b

16、,a=butter(n,Wn,stop);%計算H(z)分子、分母多項式系數H,F=freqz(b,a,1000,8000);%計算H(z)的幅頻響應,freqz(b,a,計算點數,采樣速率)subplot(2,1,1)plot(F,20*log10(abs(H) %畫出幅頻特性圖xlabel(Frequency(Hz); ylabel(Magnitude(dB) title(帶阻濾波器)axis(0 4000 -30 3);grid onpha=angle(H)*180/pi;subplot(2,1,2)plot(F,pha);grid on %畫出相頻特性圖xlabel(Frequency

17、(Hz); ylabel(phase);2、程序如下：fn=10000;%采樣頻率 fp=1000,1500; %通帶截止頻率fs=600,1900; %阻帶起始頻率Rp=3; %通帶最大衰減Rs=20;%阻帶最小衰減Wp=fp/(fn/2);%計算歸一化角頻率Ws=fs/(fn/2);n,Wn=buttord(Wp,Ws,Rp,Rs);%計算階數和截止頻率b,a=butter(n,Wn);%計算H(z)分子、分母多項式系數H,F=freqz(b,a,1000,8000);%計算H(z)的幅頻響應,freqz(b,a,計算點數,采樣速率)subplot(2,1,1)plot(F,20*log1

18、0(abs(H) %畫出幅頻特性圖xlabel(Frequency(Hz); ylabel(Magnitude(dB) title(帶通濾波器)axis(0 4000 -30 3);grid onpha=angle(H)*180/pi;subplot(2,1,2)plot(F,pha);grid on %畫出相頻特性圖xlabel(Frequency(Hz); ylabel(phase);3、程序如下：clear; T=1; fs=1/T; N=4;wc=pi/4; omegach=2*tan(wc/2)/T;%模擬濾波器的截止頻率M=1; N=1,sqrt(5),2,sqrt(2),1;h,w=freqs(M,N,512); %模