實(shí)驗(yàn)四 IIR數(shù)字濾波器的設(shè)計(jì)與濾波

1、實(shí)驗(yàn)四 IIR數(shù)字濾波器的設(shè)計(jì)與濾波一、巴特沃斯模擬濾波器的設(shè)計(jì)1 模擬濾波器的設(shè)計(jì)參數(shù)模擬濾波器的4個(gè)重要的通帶、阻帶參數(shù)為：或Omegap：通帶截止頻率 或Omegas：阻帶截至頻率：通帶內(nèi)波動(dòng)(dB)，即通帶內(nèi)所允許的最大衰減；：阻帶內(nèi)最小衰減通過以上參數(shù)就可以進(jìn)行模擬濾波器的設(shè)計(jì)。2 巴特沃斯模擬濾波器設(shè)計(jì)1） 巴特沃斯濾波器階數(shù)的選擇：在已知設(shè)計(jì)參數(shù)，之后，可利用“buttord”命令可求出所需要的濾波器的階數(shù)和3dB截止頻率，其格式為：N，Omegac=buttordfp，fs，Rp，Rs，s，其中fp，fs，Rp，Rs分別為通帶截止頻率、阻帶起始頻率、通帶內(nèi)波動(dòng)、阻帶內(nèi)最小衰減。

2、返回值N為濾波器的最低階數(shù)，Wc為3dB截止頻率。2） 巴特沃斯濾波器系數(shù)計(jì)算：由巴特沃斯濾波器的階數(shù)N以及3dB截止頻率Omegac可以計(jì)算出對(duì)應(yīng)傳遞函數(shù)H(z)的分子分母系數(shù)，MATLAB提供的命令如下：l 巴特沃斯低通濾波器系數(shù)計(jì)算：b，a=butter(N,Omegac)，其中b為H(z)的分子多項(xiàng)式系數(shù)，a為H(z)的分母多項(xiàng)式系數(shù)l 巴特沃斯高通濾波器系數(shù)計(jì)算：b，a=butter(N,Omegac,High) l 巴特沃斯帶通濾波器系數(shù)計(jì)算：b，a=butter(N，Omega1，Omega2)，其中Omega1，Omega2為通帶截止頻率，是2元向量，需要注意的是該函數(shù)返回的是

3、2N階濾波器系數(shù)。l 巴特沃斯帶阻濾波器系數(shù)計(jì)算：b，a=butter(N，Omega1，Omega2，stop)，其中Omega1，Omega2為通帶截止頻率，是2元向量，需要注意的是該函數(shù)返回的也是2N階濾波器系數(shù)。二、巴特沃斯數(shù)字濾波器的設(shè)計(jì)1. 數(shù)字濾波器的設(shè)計(jì)參數(shù)濾波器的4個(gè)重要的通帶、阻帶參數(shù)為：通帶截止頻率（Hz） ：阻帶起始頻率（Hz）：通帶內(nèi)波動(dòng)（dB），即通帶內(nèi)所允許的最大衰減；：阻帶內(nèi)最小衰減設(shè)采樣速率（即奈奎斯特速率）為，將上述參數(shù)中的頻率參數(shù)轉(zhuǎn)化為歸一化頻率參數(shù)：歸一化通帶截止頻率，； ：歸一化阻帶截至頻率，通過以上參數(shù)就可以進(jìn)行數(shù)字濾波器的設(shè)計(jì)。l 低通濾波器情況：

4、采樣頻率為8000Hz，要求通帶截止頻率為1500Hz，阻帶起始頻率為2000Hz，通帶內(nèi)波動(dòng)3dB，阻帶內(nèi)最小衰減為50dB，則=1500/4000，=2000/4000，=3dB，=50dB。l 高通濾波器情況：采樣頻率為8000Hz，要求通帶截止頻率為1500Hz，阻帶起始頻率為1000Hz，通帶內(nèi)波動(dòng)3dB，阻帶內(nèi)最小衰減為65dB，則=1500/4000，=1000/4000，=3dB，=65dB。l 帶通濾波器情況：采樣頻率為8000Hz，要求通帶截止頻率為800Hz，1500Hz，阻帶起始頻率為500Hz，1800Hz，通帶內(nèi)波動(dòng)3dB，阻帶內(nèi)最小衰減為45dB，則=800/40

5、00，1500/4000，=500/4000，1800/4000，=3dB，=45dB。l 帶阻濾波器情況：采樣頻率為8000Hz，要求通帶截止頻率為800Hz，1500Hz，阻帶起始頻率為1000Hz，1300Hz，通帶內(nèi)波動(dòng)3dB，阻帶內(nèi)最小衰減為55dB，則=800/4000，1500/4000，=1000/4000，1300/4000，=3dB，=45dB。2. 巴特沃斯數(shù)字濾波器設(shè)計(jì)1) 巴特沃斯濾波器階數(shù)的選擇：在已知設(shè)計(jì)參數(shù)，之后，可利用“buttord”命令可求出所需要的濾波器的階數(shù)和3dB截止頻率，其格式為：N，Wc=buttordWp，Ws，Rp，Rs，其中Wp，Ws，Rp

6、，Rs分別為通帶截止頻率、阻帶起始頻率、通帶內(nèi)波動(dòng)、阻帶內(nèi)最小衰減。返回值N為濾波器的最低階數(shù)，Wc為3dB截止頻率。2) 巴特沃斯濾波器系數(shù)計(jì)算：由巴特沃斯濾波器的階數(shù)N以及3dB截止頻率Wc可以計(jì)算出對(duì)應(yīng)傳遞函數(shù)H(z)的分子分母系數(shù)，MATLAB提供的命令如下：l 巴特沃斯低通濾波器系數(shù)計(jì)算：b，a=butter(N,Wc)，其中b為H(z)的分子多項(xiàng)式系數(shù)，a為H(z)的分母多項(xiàng)式系數(shù)l 巴特沃斯高通濾波器系數(shù)計(jì)算：b，a=butter(N,Wc,High) l 巴特沃斯帶通濾波器系數(shù)計(jì)算：b，a=butter(N，W1，W2)，其中W1，W2為通帶截止頻率，是2元向量，需要注意的是該

7、函數(shù)返回的是2N階濾波器系數(shù)。l 巴特沃斯帶阻濾波器系數(shù)計(jì)算：b，a=butter(N，W1，W2，stop)，其中W1，W2為通帶截止頻率，是2元向量，需要注意的是該函數(shù)返回的也是2N階濾波器系數(shù)。3. 巴特沃斯濾波器設(shè)計(jì)實(shí)例1) 采樣速率為8000Hz，要求設(shè)計(jì)一個(gè)低通濾波器，=2100Hz，=2500Hz，=3dB，=25dB。程序如下：fn=8000; fp=2100; fs=2500; Rp=3; Rs=25;Wp=fp/(fn/2);%計(jì)算歸一化頻率Ws=fs/(fn/2);N,Wc=buttord(Wp,Ws,Rp,Rs);%計(jì)算階數(shù)和截止頻率b,a=butter(N,Wc);%

8、計(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)') title('低通濾波器')axis(0 4000 -30 3);grid onsubplot(2,1,2)pha=angle(H)*180/pi;plot(F,pha);grid on2) 采樣速率為10000Hz，要求設(shè)計(jì)一個(gè)帶

9、阻濾波器，=1000Hz，1500Hz，=1200Hz，1300Hz，=3dB，=30dB。程序如下：fn=10000; fp=1000,1500; fs=1200,1300; Rp=3; Rs=30;Wp=fp/(fn/2);%計(jì)算歸一化角頻率Ws=fs/(fn/2);N,Wc=buttord(Wp,Ws,Rp,Rs);%計(jì)算階數(shù)和截止頻率b,a=butter(N,Wc,'stop');%計(jì)算H(z)分子、分母多項(xiàng)式系數(shù)H,F=freqz(b,a,1000,10000);%計(jì)算H(z)的幅頻響應(yīng),freqz(b,a,計(jì)算點(diǎn)數(shù),采樣速率)subplot(2,1,1)plot(F

10、,20*log10(abs(H)axis(0 5000 -35 3)xlabel('Frequency(Hz)') ylabel('Magnitude(dB)') title('帶阻濾波器')grid onsubplot(2,1,2)pha=angle(H)*180/pi;plot(F,pha)grid on三、IIR數(shù)字濾波器的設(shè)計(jì)與信號(hào)濾波(一)、實(shí)驗(yàn)?zāi)康脑O(shè)計(jì)IIR濾波器，實(shí)現(xiàn)對(duì)存在加性干擾的時(shí)域離散信號(hào)進(jìn)行濾波。(二)、實(shí)驗(yàn)內(nèi)容已知帶加性干擾的信號(hào)用x(n)表示，x(n)=xs(n)+(n)，式中xs(n)是有用的信號(hào)，是一個(gè)00.2rad

11、的帶限信號(hào)。(n)是一個(gè)干擾信號(hào)，其頻譜分布在0.3rad以上。要求設(shè)計(jì)一個(gè)巴特沃斯IIR數(shù)字濾波器對(duì)信號(hào)x(n)進(jìn)行濾波，將干擾(n)濾除。要求在xs(n) 所在的通帶內(nèi)濾波器幅度平坦，在0.2rad處幅度衰減不大于1dB，在噪聲所在的0.3rad 以上的頻帶內(nèi)濾波器幅度衰減大于等于40 dB。(三)、實(shí)驗(yàn)步驟1. 根據(jù)題目要求確定要設(shè)計(jì)的數(shù)字濾波器的技術(shù)指標(biāo)(低通濾波器指標(biāo)：wp=0.2rad ,ws=0.3rad,p=1dB,s=40dB)；2. 用雙線性變換法頻率轉(zhuǎn)換公式，將DF技術(shù)指標(biāo)轉(zhuǎn)換為巴特沃斯AF的技術(shù)指標(biāo)；3. 調(diào)用MATLAB函數(shù)buttord和butter，設(shè)計(jì)該模擬濾波

12、器；4. 用雙線性變換法將設(shè)計(jì)的模擬濾波器轉(zhuǎn)換為數(shù)字濾波器，求出該濾波器的系統(tǒng)函數(shù)并畫出它的直接型網(wǎng)絡(luò)結(jié)構(gòu)圖；5. 畫出該濾波器的幅度特性曲線；6. 根據(jù)1所確定的技術(shù)指標(biāo)，調(diào)用MATLAB函數(shù)buttord和butter，直接設(shè)計(jì)數(shù)字濾波器，觀察設(shè)計(jì)結(jié)果與上面用雙線性變換法的設(shè)計(jì)結(jié)果是否相同。8. 濾波仿真：調(diào)用MATLAB工具箱函數(shù)filter對(duì)下面給出的帶加性干擾的信號(hào)x(n)進(jìn)行濾波，觀察濾波效果(繪制濾波前后信號(hào)的時(shí)域和頻域波形)。實(shí)驗(yàn)信號(hào)為心電圖信號(hào)采樣序列，人體心電圖信號(hào)在測(cè)量過程中往往受到工業(yè)高頻干擾，所以必須經(jīng)過低通濾波處理后，才能作為判斷心臟功能的有用信息。下面給出一實(shí)際心

13、電圖信號(hào)采樣序列樣本x(n)，其中存在高頻干擾。在實(shí)驗(yàn)中，以x(n)作為輸入序列，濾除其中的干擾成分。x(n)：=-4, -2, 0, -4, -6, -4, -2, -4, -6, -6, -4, -4, -6, -6, -2, 6, 12, 8, 0, -16, -38, -60, -84, -90, -66, -32, -4, -2, -4, 8, 12, 12, 10, 6, 6, 6, 4, 0, 0, 0, 0, 0, -2, -4, 0, 0, 0, -2, -2, 0, 0, -2, -2, -2, -2, 0 (四)、實(shí)驗(yàn)用MATLAB函數(shù)簡(jiǎn)介filter功能：一維數(shù)字濾波器

14、直接型實(shí)現(xiàn)函數(shù)。調(diào)用格式：yn=filter(B, A, xn)：按直接型實(shí)現(xiàn)結(jié)構(gòu)對(duì)輸入信號(hào)向量xn濾波，返回濾波器輸出信號(hào)向量yn，調(diào)用參數(shù)B和A分別為濾波器系統(tǒng)函數(shù)的分子和分母多項(xiàng)式系數(shù)向量。其實(shí)質(zhì)是求解差分方程：A(1)y(n)=B(1)x(n)+B(2)x(n-1)+B(M+1)x(n-M)-A(2)y(n-1)-A(N+1)y(n-N如果A(1)不等于1時(shí)，則對(duì)系數(shù)關(guān)于A(1)歸一化后計(jì)算輸出信號(hào)y(n)。當(dāng)A=1時(shí)，對(duì)應(yīng)FIR-DF的直接型實(shí)現(xiàn)。五、實(shí)驗(yàn)報(bào)告要求1. 寫出所設(shè)計(jì)的數(shù)字濾波器的主要技術(shù)指標(biāo)、系統(tǒng)函數(shù)、網(wǎng)絡(luò)結(jié)構(gòu)；2. 寫出兩種IIR數(shù)字濾波器的設(shè)計(jì)原理、步驟；3. 寫出

15、仿真濾波試驗(yàn)的結(jié)果。附錄1：用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器的程序：%輸入信號(hào)及其波形x=-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,. 0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,. 4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0;n=0:55;subplot(2,2,1);stem(n,x,'.');axis(0 60 -100 50);hold on;n=0:60;m=zeros(61);plot(

16、n,m);xlabel('n');ylabel('x(n)');title('心電圖信號(hào)采樣序列x(n)');X=fft(x,256);subplot(2,2,2)k=0:255;plot(2*k/256,abs(X)xlabel('k');ylabel('X(k)');title('心電圖信號(hào)采樣序列的頻譜');%用雙線性變換法濾波器設(shè)計(jì)IIR-DFT=1;Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs=40;Omegap=2*tan(Wp/2)/T; %模擬波器通帶截止頻率Omegas=

17、2*tan(Ws/2)/T; %模擬波器通帶截止頻率Fs=1/T; %采樣頻率N,Wc=buttord(Omegap,Omegas,Rp,Rs,'s'); %計(jì)算巴特沃斯AF階數(shù)N及3dB截止頻率WcZ,P,K=buttap(N); %求butterworth歸一化原型Ha(p)，零極點(diǎn)增益形式b,a=zp2tf(Z,P,K); % 轉(zhuǎn)換成歸化系統(tǒng)函數(shù)形式，Ha(p)為分子分母多項(xiàng)式形式bt,at=lp2lp(b,a,Wc); %去歸一化轉(zhuǎn)換為實(shí)際低通Ha(s)% b1,a1=butter(N,Omegac,'s');%可以代替以上三個(gè)步驟 bb,ab=bilinear(bt,at,Fs); %雙線性變換得H(z)，F(xiàn)s為采樣頻率，改為脈沖響應(yīng)不變法，觀%察實(shí)驗(yàn)結(jié)果有何區(qū)別%求數(shù)字濾波器的幅頻特性H,w=freqz(bb,ab,100); %對(duì)應(yīng)采樣間隔為0.02*pimag=abs(H);db=20*log10(mag+eps)/max(mag);subplot(2,2,3);plot(w/pi,db);axis(0,0.6,-50,10);grid ontitle('濾波器的幅頻響應(yīng)曲線');%用設(shè)計(jì)濾波器對(duì)輸入信號(hào)進(jìn)行濾波y=filter(bb,ab,x);