IIR数字滤波器的设计及应用——MATLAB
一、實(shí)驗(yàn)?zāi)康?/h2>
(1)熟悉雙線性變換法和雙重映射法設(shè)計(jì)IIR數(shù)字濾波器的原理與方法。
(2)掌握IIR數(shù)字濾波器的MATLAB實(shí)現(xiàn)方法設(shè)計(jì)各種濾波器。
(3)觀察分析濾波器輸入輸出數(shù)據(jù)波形,理解數(shù)字濾波的概念。
二、實(shí)驗(yàn)原理
1、IIR數(shù)字濾波器設(shè)計(jì)的一般方法
IIR數(shù)字濾波器設(shè)計(jì)最常采用的方式是先按技術(shù)指標(biāo)設(shè)計(jì)模擬濾波器,然后采用合適的方法將模擬濾波器離散化為所需要的數(shù)字濾波器。按照這一方式,IR數(shù)字濾波器的設(shè)計(jì)過程為:
(1)確定數(shù)字濾波器實(shí)際需滿足的技術(shù)指標(biāo)參數(shù);
(2)將數(shù)字濾波器設(shè)計(jì)的技術(shù)指標(biāo)轉(zhuǎn)換為模擬濾波器設(shè)計(jì)的技術(shù)指標(biāo);
(3)設(shè)計(jì)模擬濾波器,得到模擬濾波器的系統(tǒng)函數(shù)Ha(s);
(4)采用合適的離散化方法,將所設(shè)計(jì)的模擬濾波器離散化,為所需要的數(shù)字離散化,得到數(shù)字濾波器系統(tǒng)函數(shù)H(z)
2、常用模擬原型濾波器及其仿真計(jì)算方法
在設(shè)計(jì)數(shù)字濾波器之前所設(shè)計(jì)的模擬濾波器稱為原型濾波器,常用的模擬原型濾波器有巴特沃斯模擬濾波器、切比雪夫模擬濾波器和橢圓模擬濾波器三種。Matlab環(huán)境中,三種模擬濾波器的仿真設(shè)計(jì)方法介紹如下:
(1)巴特沃斯模擬濾波器設(shè)計(jì)。MATIAB信號(hào)處理工具箱提供的用于巴特沃什模擬濾波器設(shè)計(jì)的函數(shù)主要是 buttord()和butter()。這兩個(gè)函數(shù)的調(diào)用格式分別如下:
Buttord()函數(shù)的調(diào)用用于計(jì)算巴特沃什模擬濾波器的階次N和3dB截止頻率。輸入?yún)?shù)為四個(gè),含義分別為:
Wp:通帶截止頻率;ws:阻帶截止頻率;Rp:通帶最大衰減,單位為dB;Rs:阻帶最小衰減,單位為dB。
函數(shù)調(diào)用后,可返回兩個(gè)參數(shù),分別是:
N:模擬濾波器的階次;wc:模擬濾波器的3dB截止頻率。在函數(shù)的使用中,當(dāng)wp<ws時(shí),設(shè)計(jì)低通濾波器;當(dāng)wp>ws時(shí),設(shè)計(jì)高通濾波器;當(dāng)wp和ws為二元矢量時(shí),設(shè)計(jì)帶通或帶阻濾波器,這時(shí)wc也是二元矢量。
Buttord()函數(shù)用于計(jì)算巴特沃什模擬濾波器系統(tǒng)函數(shù)中分子和分母多項(xiàng)式的系數(shù)向量B和A。輸入?yún)?shù)為四個(gè),含義分別為:
N:模擬濾波器的階次;
wc:模擬濾波器的3dB截止頻率;當(dāng)設(shè)計(jì)的是低通或高通濾波器時(shí),為一元向量;設(shè)計(jì)的是帶通或帶阻濾波器時(shí),是二元向量,向量的兩個(gè)元素分別為3dB上限截止頻率和下限截止頻率。
ftype:濾波器類型,缺省時(shí)默認(rèn)設(shè)計(jì)低通(wc為一元向量)或帶通(wc為二元向量)濾波器;ftype=high,設(shè)計(jì)高通濾波器,ftype=stop,設(shè)計(jì)帶阻濾波器(wc為二元向量)。
s:固定,代表是模擬濾波器的設(shè)計(jì)。
函數(shù)調(diào)用后,可返回兩個(gè)參數(shù),分別是:
B:系統(tǒng)函數(shù)分子多項(xiàng)式系數(shù)向量;A:系統(tǒng)函數(shù)分母多項(xiàng)式系數(shù)向量。
由系數(shù)向量B和A可以寫出模擬濾波器的系統(tǒng)函數(shù)為:
(2)切比雪夫模擬濾波器設(shè)計(jì)。MATIAB信號(hào)處理工具箱提供的用于第一類切比雪夫模擬濾波器設(shè)計(jì)的函數(shù)主要是 cheblord()和cheby1()。這兩個(gè)函數(shù)的調(diào)用格式如
[N,wc]=cheblord(Wwp, ws,Rp,Rs,'s[B,A]=cheby1(N, we,'ftype','s')這兩個(gè)函數(shù)中的所有參數(shù)基本與前述的巴特沃什模擬濾波器設(shè)計(jì)函數(shù)調(diào)用中的參數(shù)含義一致,唯一的區(qū)別在于 wc指的是第一類切比雪夫模擬濾波器的通帶截止頻率,而不是3dB截止頻率。
(3)橢圓模擬濾波器設(shè)計(jì)MATLAB信號(hào)處理工具箱提供的用于橢圓模擬濾波器設(shè)計(jì)的函數(shù)主要是ellipord()和 ellip()。這兩個(gè)函數(shù)的調(diào)用格式和參數(shù)含義與前述幾個(gè)函數(shù)基本類似,此處不再贅述。
3、模擬濾波器到數(shù)字濾波器的離散化方法
模擬濾波器設(shè)計(jì)完成后,由模擬濾波器到數(shù)字濾波器的離散化,有兩種方式。即沖激響應(yīng)不變法和雙線性映射法。
沖激響應(yīng)不變法即基本思想是:數(shù)字濾波器的單位取樣響應(yīng)是模擬濾波器單位沖激響應(yīng)的采樣值即:
h(n)= h(t)|t=nT
這種離散化方法能夠保證數(shù)字角頻率和模擬角頻率之間嚴(yán)格的線性關(guān)系,從而不改變?yōu)V波器頻率響應(yīng)的形狀,但會(huì)存在頻率寧限濾波器的設(shè)計(jì)(低通和響應(yīng)的混疊失真,因而只適合于進(jìn)行帶限濾波器的設(shè)計(jì)(低通和帶通)。
雙線性映射法主要為克服沖激響應(yīng)不變法的混疊失真現(xiàn)象,其基本思想如下:
(1)先把模擬濾波器的s平面映射為s平面,即限制Ω在-Π/T~Π/T之間取值。
(2)然后采用沖激不變法的方式將s1平面映射到數(shù)字濾波器的z平面。
雙線性映射法克服了混疊現(xiàn)象,但是模擬角頻率和數(shù)字角頻率之間不再滿足嚴(yán)格的線性關(guān)系,故存在著頻率響應(yīng)的畸變?yōu)榭朔l率響應(yīng)畸變的影響,在模擬濾波技術(shù)指標(biāo)確定時(shí)需要進(jìn)行反畸變校正,即:
MATLAB的數(shù)字信號(hào)處理工具箱提供impinvar()函數(shù)和
雙線性映射法。這兩個(gè)函bilinear 函數(shù)用以實(shí)現(xiàn)沖激響應(yīng)不變法和數(shù)的調(diào)用格式如下:
[Bz,Az]= bilinear(B,A,fs)其中,B,A分別是H。(s)的分子、分母多項(xiàng)式的系數(shù)向量,Bz,Az分別是H(Z)的分子、分母多項(xiàng)式的系數(shù)向量,fs為采樣頻率。
4、IIR數(shù)字濾波器的實(shí)現(xiàn)方法
IIR數(shù)字濾波器的設(shè)計(jì)計(jì)算最終是獲取濾波器的系統(tǒng)函數(shù)H(Z),但是應(yīng)用所設(shè)計(jì)的數(shù)字濾波器對(duì)實(shí)際信號(hào)進(jìn)行處理,一般是借助差分方程來實(shí)現(xiàn)的。即通過系統(tǒng)函數(shù)H(Z)獲得濾波器系統(tǒng)的差分方程,然后根據(jù)差分方程利用順序迭代的方式進(jìn)行濾波器實(shí)現(xiàn)程序的編寫。但需注意的是,順序迭代程序之前需要賦初始條件,初始條件的個(gè)數(shù)取決于濾波器的階次N。以二階IIR數(shù)字濾波器為例,其實(shí)現(xiàn)程序如下:
三、實(shí)驗(yàn)步驟、數(shù)據(jù)記錄及處理
代碼示例:
%輸入并顯示濾波前心電圖采樣序列; clc;clear all;close all;%清空工作空間; xn=[-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;Fs=100;%設(shè)置時(shí)域序列與采樣數(shù)頻率; figure('name','心電圖時(shí)域與頻域采樣圖')%定義繪圖框標(biāo)題名; subplot(211);stem(n,xn);grid on;xlabel('n');ylabel('x(n)');title('心電圖采樣序列時(shí)域');%繪制時(shí)域采樣圖; ndft=length(n);%采樣數(shù)據(jù)點(diǎn)數(shù); Xk=fft(xn,ndft);%快速傅里葉變換; mag=abs(Xk);%獲得傅里葉變換后的振幅; k=0:ndft/2-1;%離散時(shí)間序列; f=Fs*k/ndft;%采樣頻率序列; subplot(212);stem(f,mag(1:(ndft/2)),'.');grid on;xlabel('f');title('心電信號(hào)頻譜圖');%繪制采樣頻譜圖; clc;clear all;close all;%清空工作空間; Fs=100;fp=15;fs=23;k1=-2;k2=-20;Ts=1/Fs;%定義信號(hào)采樣頻率;通帶截止頻率;通帶起始頻率; wp=2*pi*fp;ws=2*pi*fs;%分別求出模擬通帶截止頻率wp和模擬阻帶截止頻率ws; %wp=(2/Ts)*tan(wp/2);Ws=(2/Ts)*tan(ws/2);%分別求出數(shù)字通帶截止頻率Wp和數(shù)字阻帶截止頻率Ws; [N,wc]=buttord(wp,ws,k1,k2,'s');%計(jì)算巴特沃什模擬濾波器的階次; N%顯示巴特沃什模擬濾波器的階次 [b1,a1]=butter(N,wc,'s');%計(jì)算H(s)系數(shù)向量b1和a1; [H,m]=freqz(b1,a1);%計(jì)算巴特沃什模擬濾波器的頻率響應(yīng); mag=abs(H);%求巴特沃什模擬濾波器的幅頻響應(yīng)振幅; pha=angle(H);%求巴特沃什模擬濾波器的相頻響應(yīng)相位; db=20*log10((mag+eps)/max(mag));%將幅頻特性轉(zhuǎn)換為對(duì)數(shù)形式; figure('name','巴特沃什模擬濾波器幅頻相頻特性')%定義巴特沃什模擬濾波器幅頻相特性圖形標(biāo)題名 subplot(211);plot(m*100/(2*pi),db);grid on;xlabel('w/2pi');ylabel('db');title('巴特沃什模擬濾波器幅頻特性');%繪制巴特沃什模擬濾波器幅頻特性 subplot(212);plot(m*100/(2*pi),pha);grid on;xlabel('w/2pi');ylabel('rad');title('巴特沃什模擬濾波器相頻特性');%巴特沃什模擬濾波器相頻特性 clear all;clc;close all;%清空工作空間; Fs=100;fp=15;fs=23;Rp=-2;Rs=-20;Ts=1/Fs;%定義信號(hào)采樣頻率;通帶截止頻率;通帶起始頻率; wp=2*pi*fp;ws=2*pi*fs;%分別求出模擬通帶截止頻率wp和模擬阻帶截止頻率ws [N,wc]=buttord(wp,ws,Rp,Rs,'s');%計(jì)算巴特沃什模擬濾波器的階次 N%顯示雙線性數(shù)字濾波器的階次 [b1,a1]=butter(N,wc,'s');%計(jì)算H(s)系數(shù)向量b1和a1 [b,a]=bilinear(b1,a1,Fs);%利用雙線性映射法將巴特沃什模擬濾波器轉(zhuǎn)換為數(shù)字濾波器 [H,m]=freqz(b,a);%計(jì)算巴特沃什數(shù)字濾波器的頻率響應(yīng) mag=abs(H);%求巴特沃什數(shù)字濾波器的幅頻響應(yīng) pha=angle(H);%求巴特沃什數(shù)字濾波器的相頻響應(yīng) db=20*log10((mag+eps)/max(mag));%將幅頻特性轉(zhuǎn)換為對(duì)數(shù)形式 figure('name','雙線性數(shù)字濾波器幅頻相頻特性') subplot(211);plot(m*100/(2*pi),db);grid on;axis([0,56,-100,50]);xlabel('w/2pi');ylabel('db');title('雙線性數(shù)字濾波器幅頻特性');%繪制雙線性數(shù)字濾波器幅頻特性; subplot(313);plot(m*100/(2*pi),pha);grid on;xlabel('w/2pi');ylabel('rad');title('雙線性數(shù)字濾波器相頻特性');%繪制雙線性數(shù)字濾波器相頻特性; clear all;clc;close all;%清空工作空間; %畫出濾波前心電圖采樣序列 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;Fs=100; figure('name','雙線性映射法濾波后心電序列波形') subplot(221);stem(n,x);grid on;axis([0,56,-100,50]);xlabel('n');ylabel('x(n)');title('濾波前心電圖采樣序列'); %畫出雙線性映射法將巴特沃什模擬濾波器轉(zhuǎn)換為數(shù)字濾波器的幅頻特性曲線 Fs=100;fp=15;fs=23;Rp=-2;Rs=-20;Ts=1/Fs; wp=2*pi*fp;ws=2*pi*fs;%分別求出模擬通帶截止頻率wp和模擬阻帶截止頻率ws [N,wc]=buttord(wp,ws,Rp,Rs,'s');%計(jì)算巴特沃什模擬濾波器的階次 N%顯示巴特沃什濾波器的階次 [b1,a1]=butter(N,wc,'s');%計(jì)算H(s)系數(shù)向量b1和a1 [b,a]=bilinear(b1,a1,Fs);%利用雙線性映射法將巴特沃什模擬濾波器轉(zhuǎn)換為數(shù)字濾波器 [H,m]=freqz(b,a);%計(jì)算巴特沃什數(shù)字濾波器的頻率響應(yīng) mag=abs(H);%求巴特沃什數(shù)字濾波器的幅頻響應(yīng) db=20*log10((mag+eps)/max(mag));%將幅頻特性轉(zhuǎn)換為對(duì)數(shù)形式 subplot(222);plot(m*100/(2*pi),db);grid on;axis([0,56,-100,50]);xlabel('w/2pi');ylabel('db');title('雙線性數(shù)字濾波器幅頻特性');%繪制雙線性數(shù)字濾波器幅頻特性 %利用所設(shè)計(jì)出的濾波器差分方程 y(1)=x(1); y(2)=(x(1)+x(2))/2; y(3)=(x(1)+x(2)+x(3))/3; y(4)=(x(1)+x(2)+x(3)+x(4))/4; y(5)=(x(1)+x(2)+x(3)+x(4)+x(5))/5; y(6)=(x(1)+x(2)+x(3)+x(4)+x(5)+x(6))/6; for n=7:length(x); y(n)=b(1)*x(n)+b(2)*x(n-1)+b(3)*x(n-2)+b(4)*x(n-3)+b(5)*x(n-4)+b(6)*x(n-5)+b(7)*x(n-6)-a(2)*y(n-1)-a(3)*y(n-2)-a(4)*y(n-3)-a(5)*y(n-4)-a(6)*y(n-5)-a(7)*y(n-6); end %繪制濾波后心電序列波形 subplot(223); n=0:55;%設(shè)置時(shí)域序列 ndft=length(x);%采樣數(shù)據(jù)點(diǎn)數(shù); stem(n,y,'.');grid on;axis([0,56,-100,50]);hold on; n=0:60;%設(shè)置時(shí)域序列 m=zeros(61); subplot(223);plot(n,m);grid on;xlabel('n');ylabel('y(n)');title('IIR數(shù)字濾波后心電采樣序列'); Xk=fft(y,ndft);%快速傅里葉變換; mag=abs(Xk);%獲得傅里葉變換后的振幅; k=0:ndft/2-1;%離散時(shí)間序列; f=Fs*k/ndft;%采樣頻率序列; subplot(224);plot(f,mag(1:(ndft/2)));grid on;xlabel('f');title('IIR數(shù)字濾波后心電信號(hào)頻譜圖'); clear all;clc;close all;%清空工作空間; Fs=100;fp=15;fs=23;Rp=-2;Rs=-20;Ts=1/Fs;%定義信號(hào)采樣頻率;通帶截止頻率;通帶起始頻率; wp=2*pi*fp;ws=2*pi*fs;%分別求出模擬通帶截止頻率wp和模擬阻帶截止頻率ws [N,wc]=buttord(wp,ws,Rp,Rs,'s');%計(jì)算巴特沃什模擬濾波器的階次 N%顯示巴特沃什濾波器的階次 [b1,a1]=butter(N,wc,'s');%計(jì)算H(s)系數(shù)向量b1和a1 [b,a]=impinvar(b1,a1,Fs);%利用沖激響應(yīng)不變法將巴特沃什模擬濾波器轉(zhuǎn)換為數(shù)字濾波器 [H,m]=freqz(b,a);%計(jì)算巴特沃什數(shù)字濾波器的頻率響應(yīng) mag=abs(H);%求巴特沃什數(shù)字濾波器的幅頻響應(yīng) pha=angle(H);%求巴特沃什數(shù)字濾波器的相頻響應(yīng) db=20*log10((mag+eps)/max(mag));%將幅頻特性轉(zhuǎn)換為對(duì)數(shù)形式 figure('name','沖激響應(yīng)不變法數(shù)字濾波器幅頻相頻特性') subplot(211);plot(m*100/(2*pi),db);grid on;axis([0,56,-100,50]);xlabel('w/2pi');ylabel('db');title('沖激響應(yīng)不變法數(shù)字濾波器幅頻特性'); subplot(313);plot(m*100/(2*pi),pha);grid on;xlabel('w/2pi');ylabel('rad');title('沖激響應(yīng)不變法數(shù)字濾波器相頻特性'); clear all;clc;close all; %畫出濾波前心電圖采樣序列 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;Fs=100; figure('name','沖激響應(yīng)不變法濾波后心電序列波形') subplot(221);stem(n,x);grid on;axis([0,56,-100,50]);xlabel('n');ylabel('x(n)');title('濾波前心電圖采樣序列'); %畫出雙線性映射法將巴特沃什模擬濾波器轉(zhuǎn)換為數(shù)字濾波器的幅頻特性曲線 Fs=100;fp=15;fs=23;Rp=-2;Rs=-20;Ts=1/Fs; wp=2*pi*fp;ws=2*pi*fs;%分別求出模擬通帶截止頻率wp和模擬阻帶截止頻率ws [N,wc]=buttord(wp,ws,Rp,Rs,'s');%計(jì)算巴特沃什模擬濾波器的階次 N%顯示巴特沃什濾波器的階次 [b1,a1]=butter(N,wc,'s');%計(jì)算H(s)系數(shù)向量b1和a1 [b,a]=impinvar(b1,a1,Fs);%利用沖激響應(yīng)不變法將巴特沃什模擬濾波器轉(zhuǎn)換為數(shù)字濾波器 [H,m]=freqz(b,a);%計(jì)算巴特沃什數(shù)字濾波器的頻率響應(yīng) mag=abs(H);%求巴特沃什數(shù)字濾波器的幅頻響應(yīng) db=20*log10((mag+eps)/max(mag));%將幅頻特性轉(zhuǎn)換為對(duì)數(shù)形式 subplot(222);plot(m*100/(2*pi),db);grid on;axis([0,56,-100,50]);xlabel('w/2pi');ylabel('db');title('沖激響應(yīng)不變數(shù)字濾波器幅頻特性'); %利用所設(shè)計(jì)出的濾波器差分方程 y(1)=x(1); y(2)=(x(1)+x(2))/2; y(3)=(x(1)+x(2)+x(3))/3; y(4)=(x(1)+x(2)+x(3)+x(4))/4; y(5)=(x(1)+x(2)+x(3)+x(4)+x(5))/5; y(6)=(x(1)+x(2)+x(3)+x(4)+x(5)+x(6))/6; for n=7:length(x); y(n)=b(1)*x(n)+b(2)*x(n-1)+b(3)*x(n-2)+b(4)*x(n-3)+b(5)*x(n-4)+b(6)*x(n-5)+b(7)*x(n-6)-a(2)*y(n-1)-a(3)*y(n-2)-a(4)*y(n-3)-a(5)*y(n-4)-a(6)*y(n-5)-a(7)*y(n-6); end %繪制濾波后心電序列波形 subplot(223); n=0:55; ndft=length(x); stem(n,y,'.');grid on;axis([0,56,-100,50]); hold on; n=0:60; m=zeros(61); subplot(223);plot(n,m);grid on; xlabel('n');ylabel('y(n)');title('IIR數(shù)字濾波后心電采樣序列'); Xk=fft(y,ndft); mag=abs(Xk); k=0:ndft/2-1; f=Fs*k/ndft; subplot(224); plot(f,mag(1:(ndft/2)));grid on; xlabel('f');title('IIR數(shù)字濾波后心電信號(hào)頻譜圖');四、分析
分析略
五、思考題
思考題略
六、總結(jié)
總結(jié)略
總結(jié)
以上是生活随笔為你收集整理的IIR数字滤波器的设计及应用——MATLAB的全部?jī)?nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: ps图像压缩插件:TinyPNG and
- 下一篇: Evernote是什么软件?印象笔记fo