數(shù)字信號處理--用雙線性變換法設(shè)計(jì)IIR 數(shù)字濾波器

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

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

3、能：把具有B,A模擬濾波器傳遞函數(shù)模型轉(zhuǎn)換為采樣頻率為Fs的數(shù)字濾波器的傳遞函數(shù)模型BZ,AZ，F(xiàn)s默認(rèn)值為1。例：一個4階的Butterworth模擬低通濾波器的系統(tǒng)函數(shù)如下：試用沖激響應(yīng)不變法求出Butterworth模擬低通數(shù)字濾波器的系統(tǒng)函數(shù)。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)功能：把模擬濾波器的零極點(diǎn)模型轉(zhuǎn)換成數(shù)字濾波器的零極點(diǎn)模型，F(xiàn)s是采樣頻率格式二：numd,dend= bilinear(num,de

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

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

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

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

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

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

10、39;) 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');用雙線性變換法wp=2100*2*pi; %利用 ws=2500*2*pi;Rp=3;Rs=25;Fs=8000;Ts=1/Fs;%選擇濾波器的最小階數(shù)N,Wn=buttord(wp,ws,Rp,Rs,'s');%創(chuàng)建butterworth

11、模擬濾波器Z,P,K=buttap(N);%把濾波器零極點(diǎn)模型轉(zhuǎn)化為傳遞函數(shù)模型Bap,Aap=zp2tf(Z,P,K); %把模擬濾波器原型轉(zhuǎn)換成截至頻率為Wn的低通濾波器b,a=lp2lp(Bap,Aap,Wn);%用雙線性變換法實(shí)現(xiàn)模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換bz,az=bilinear(b,a,Fs);%繪制頻率響應(yīng)曲線H,W=freqz(bz,az);plot(W*Fs/(2*pi),abs(H);gridxlabel('頻率/Hz')ylabel('幅度')例題2：模擬原型直接變換法設(shè)計(jì)數(shù)字濾波器：已知四階歸一化低通巴特沃斯模擬濾波器系統(tǒng)函數(shù)為，編寫

12、MATLAB程序?qū)崿F(xiàn)從設(shè)計(jì)3dB截止頻率為，設(shè)采樣周期為T=1，的四階低通巴特沃斯數(shù)字濾波器。程序如下：步驟一：將設(shè)計(jì)內(nèi)容題所給歸一化巴特沃斯低通濾波器以3dB截止頻率為進(jìn)行去歸一化。步驟二：用雙線性變化法將低通模擬濾波器變換為低通數(shù)字濾波器設(shè)計(jì)程序如下：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); %模擬濾波器的幅頻響應(yīng)subplot(2,1,1);plot(w,20*log10(abs(h); axis

13、(0,10,-90,0),grid on; xlabel('Hz');ylabel('幅度'); title('歸一化模擬低通濾波器');Ms,Ns=lp2lp(M,N,omegach); %對低通濾波器進(jìn)行頻率變換hs,ws=freqs(Ms,Ns,512); %模擬濾波器的幅頻響應(yīng)subplot(2,1,2);plot(ws,20*log10(abs(hs); grid;axis(0,10,-90,0);xlabel('Hz');ylabel('幅度'); title('去歸一化模擬低通濾波器'

14、);Mz,Nz=bilinear(Ms,Ns,1/T); %對模擬濾波器雙線性變換h1,w1=freqz(Mz,Nz); %數(shù)字濾波器的幅頻響應(yīng)figureplot(w1/pi,20*log10(abs(h1); grid;xlabel('/');ylabel('幅度(dB)'); title('數(shù)字低通濾波器');axis(0,1,-160,0)四、實(shí)驗(yàn)內(nèi)容：1、采樣速率為10000Hz，要求設(shè)計(jì)一個巴特沃斯帶阻濾波器，=1000Hz，1500Hz，=1200Hz，1300Hz，=3dB，=30dB。提示：b,a=butter(N,Wc,

15、9;stop')2、 采樣速率為10000Hz，要求設(shè)計(jì)一個帶通濾波器，=1000Hz，1500Hz，=600Hz，1900Hz，=3dB，=20dB。3、已知四階歸一化低通巴特沃斯模擬濾波器系統(tǒng)函數(shù)為，編寫MATLAB程序?qū)崿F(xiàn)從設(shè)計(jì)3dB截止頻率為的四階高通巴特沃斯數(shù)字濾波器。4、設(shè)計(jì)低通濾波器，把輸入信號中的三個信號分離出來。要求：畫出濾波前后信號的波形及頻譜及低通濾波器的幅頻響應(yīng)。提示：b,a=butter(n,Wn);%計(jì)算H(z)分子、分母多項(xiàng)式系數(shù)y=filter(b,a,x);%對輸入的信號進(jìn)行濾波五、思考題：用雙線性變換法設(shè)計(jì)數(shù)字濾波器過程中，變換公式中T 的取值對設(shè)計(jì)

16、結(jié)果有無影響？為什么？程序如下：五、實(shí)驗(yàn)結(jié)果：1、程序如下：fn=10000;%采樣頻率 fp=1000,1500; %通帶截止頻率fs=600,1900; %阻帶起始頻率Rp=3; %通帶最大衰減Rs=20;%阻帶最小衰減Wp=fp/(fn/2);%計(jì)算歸一化角頻率Ws=fs/(fn/2);n,Wn=buttord(Wp,Ws,Rp,Rs);%計(jì)算階數(shù)和截止頻率b,a=butter(n,Wn,'stop');%計(jì)算H(z)分子、分母多項(xiàng)式系數(shù)H,F=freqz(b,a,1000,8000);%計(jì)算H(z)的幅頻響應(yīng),freqz(b,a,計(jì)算點(diǎn)數(shù),采樣速率)subplot(2,

17、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(Hz)'); ylabel('phase');2、程序如下：fn=10000;%采樣頻率 fp=1000,150

18、0; %通帶截止頻率fs=600,1900; %阻帶起始頻率Rp=3; %通帶最大衰減Rs=20;%阻帶最小衰減Wp=fp/(fn/2);%計(jì)算歸一化角頻率Ws=fs/(fn/2);n,Wn=buttord(Wp,Ws,Rp,Rs);%計(jì)算階數(shù)和截止頻率b,a=butter(n,Wn);%計(jì)算H(z)分子、分母多項(xiàng)式系數(shù)H,F=freqz(b,a,1000,8000);%計(jì)算H(z)的幅頻響應(yīng),freqz(b,a,計(jì)算點(diǎn)數(shù),采樣速率)subplot(2,1,1)plot(F,20*log10(abs(H) %畫出幅頻特性圖xlabel('Frequency(Hz)'); yla

19、bel('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); %模擬濾波器的幅頻響應(yīng)subplot(2,1,1);plot(w,20*log10(abs(h); axis