一、實驗目的

（1）熟悉雙線性變換法和雙重映射法設計IIR數字濾波器的原理與方法。
（2）掌握IIR數字濾波器的MATLAB實現方法設計各種濾波器。
（3）觀察分析濾波器輸入輸出數據波形，理解數字濾波的概念。

二、實驗原理

1、IIR數字濾波器設計的一般方法

IIR數字濾波器設計最常采用的方式是先按技術指標設計模擬濾波器，然后采用合適的方法將模擬濾波器離散化為所需要的數字濾波器。按照這一方式,IR數字濾波器的設計過程為:
(1）確定數字濾波器實際需滿足的技術指標參數;
(2）將數字濾波器設計的技術指標轉換為模擬濾波器設計的技術指標;
(3）設計模擬濾波器,得到模擬濾波器的系統函數Ha(s);
(4）采用合適的離散化方法，將所設計的模擬濾波器離散化，為所需要的數字離散化，得到數字濾波器系統函數H(z)

2、常用模擬原型濾波器及其仿真計算方法

在設計數字濾波器之前所設計的模擬濾波器稱為原型濾波器，常用的模擬原型濾波器有巴特沃斯模擬濾波器、切比雪夫模擬濾波器和橢圓模擬濾波器三種。Matlab環境中，三種模擬濾波器的仿真設計方法介紹如下:
(1）巴特沃斯模擬濾波器設計。MATIAB信號處理工具箱提供的用于巴特沃什模擬濾波器設計的函數主要是 buttord()和butter()。這兩個函數的調用格式分別如下:

[N,wc]=buttord(wp，ws,Rp，Rs，'s')

Buttord()函數的調用用于計算巴特沃什模擬濾波器的階次N和3dB截止頻率。輸入參數為四個，含義分別為:
Wp:通帶截止頻率;ws:阻帶截止頻率;Rp:通帶最大衰減,單位為dB;Rs:阻帶最小衰減，單位為dB。
函數調用后,可返回兩個參數，分別是:
N:模擬濾波器的階次;wc:模擬濾波器的3dB截止頻率。在函數的使用中,當wp<ws時，設計低通濾波器;當wp>ws時,設計高通濾波器;當wp和ws為二元矢量時,設計帶通或帶阻濾波器，這時wc也是二元矢量。

[B,A]=butter(N, wc,'ftype','s’

Buttord()函數用于計算巴特沃什模擬濾波器系統函數中分子和分母多項式的系數向量B和A。輸入參數為四個,含義分別為:
N:模擬濾波器的階次;
wc:模擬濾波器的3dB截止頻率;當設計的是低通或高通濾波器時,為一元向量;設計的是帶通或帶阻濾波器時，是二元向量，向量的兩個元素分別為3dB上限截止頻率和下限截止頻率。
ftype:濾波器類型，缺省時默認設計低通（wc為一元向量)或帶通（wc為二元向量）濾波器;ftype=high，設計高通濾波器，ftype=stop，設計帶阻濾波器（wc為二元向量）。
s:固定，代表是模擬濾波器的設計。
函數調用后，可返回兩個參數，分別是:
B:系統函數分子多項式系數向量;A:系統函數分母多項式系數向量。
由系數向量B和A可以寫出模擬濾波器的系統函數為:

（2）切比雪夫模擬濾波器設計。MATIAB信號處理工具箱提供的用于第一類切比雪夫模擬濾波器設計的函數主要是 cheblord()和cheby1()。這兩個函數的調用格式如

[N,wc]=cheblord(Wwp, ws，Rp,Rs，'s[B,A]=cheby1(N, we,'ftype'，'s')

這兩個函數中的所有參數基本與前述的巴特沃什模擬濾波器設計函數調用中的參數含義一致，唯一的區別在于 wc指的是第一類切比雪夫模擬濾波器的通帶截止頻率，而不是3dB截止頻率。
(3）橢圓模擬濾波器設計MATLAB信號處理工具箱提供的用于橢圓模擬濾波器設計的函數主要是ellipord()和 ellip()。這兩個函數的調用格式和參數含義與前述幾個函數基本類似,此處不再贅述。
3、模擬濾波器到數字濾波器的離散化方法
模擬濾波器設計完成后，由模擬濾波器到數字濾波器的離散化，有兩種方式。即沖激響應不變法和雙線性映射法。
沖激響應不變法即基本思想是:數字濾波器的單位取樣響應是模擬濾波器單位沖激響應的采樣值即:

h(n)= h(t)|t=nT

這種離散化方法能夠保證數字角頻率和模擬角頻率之間嚴格的線性關系，從而不改變濾波器頻率響應的形狀,但會存在頻率寧限濾波器的設計(低通和響應的混疊失真，因而只適合于進行帶限濾波器的設計（低通和帶通）。
雙線性映射法主要為克服沖激響應不變法的混疊失真現象,其基本思想如下:
(1）先把模擬濾波器的s平面映射為s平面,即限制Ω在-Π/T~Π/T之間取值。
(2）然后采用沖激不變法的方式將s1平面映射到數字濾波器的z平面。
雙線性映射法克服了混疊現象,但是模擬角頻率和數字角頻率之間不再滿足嚴格的線性關系,故存在著頻率響應的畸變為克服頻率響應畸變的影響,在模擬濾波技術指標確定時需要進行反畸變校正，即:

MATLAB的數字信號處理工具箱提供impinvar()函數和
雙線性映射法。這兩個函bilinear 函數用以實現沖激響應不變法和數的調用格式如下:

[Bz,Az]= impinvar(B,A,fs)

[Bz,Az]= bilinear(B,A,fs)其中，B,A分別是H。(s)的分子、分母多項式的系數向量,Bz，Az分別是H(Z)的分子、分母多項式的系數向量,fs為采樣頻率。
4、IIR數字濾波器的實現方法
IIR數字濾波器的設計計算最終是獲取濾波器的系統函數H(Z)，但是應用所設計的數字濾波器對實際信號進行處理,一般是借助差分方程來實現的。即通過系統函數H(Z)獲得濾波器系統的差分方程，然后根據差分方程利用順序迭代的方式進行濾波器實現程序的編寫。但需注意的是，順序迭代程序之前需要賦初始條件,初始條件的個數取決于濾波器的階次N。以二階IIR數字濾波器為例，其實現程序如下:

(1)=x(1); y(2)=(x(1)+x(2))/2;for n=3:length(x) y(n)=b(1)*x(n)+b(2)*x(n-1)+b(3)*x(n-2)-a(2)*y(n-1)-a(3)*y(n-2); end;

三、實驗步驟、數據記錄及處理

代碼示例：

%輸入并顯示濾波前心電圖采樣序列； 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;%設置時域序列與采樣數頻率； figure('name','心電圖時域與頻域采樣圖')%定義繪圖框標題名； subplot(211);stem(n,xn);grid on;xlabel('n');ylabel('x(n)');title('心電圖采樣序列時域');%繪制時域采樣圖； ndft=length(n);%采樣數據點數； Xk=fft(xn,ndft);%快速傅里葉變換； mag=abs(Xk);%獲得傅里葉變換后的振幅； k=0:ndft/2-1;%離散時間序列； f=Fs*k/ndft;%采樣頻率序列； subplot(212);stem(f,mag(1:(ndft/2)),'.');grid on;xlabel('f');title('心電信號頻譜圖');%繪制采樣頻譜圖；

clc;clear all;close all;%清空工作空間； Fs=100;fp=15;fs=23;k1=-2;k2=-20;Ts=1/Fs;%定義信號采樣頻率；通帶截止頻率；通帶起始頻率； wp=2*pi*fp;ws=2*pi*fs;%分別求出模擬通帶截止頻率wp和模擬阻帶截止頻率ws； %wp=(2/Ts)*tan(wp/2);Ws=(2/Ts)*tan(ws/2);%分別求出數字通帶截止頻率Wp和數字阻帶截止頻率Ws； [N,wc]=buttord(wp,ws,k1,k2,'s');%計算巴特沃什模擬濾波器的階次； N%顯示巴特沃什模擬濾波器的階次 [b1,a1]=butter(N,wc,'s');%計算H（s）系數向量b1和a1； [H,m]=freqz(b1,a1);%計算巴特沃什模擬濾波器的頻率響應； mag=abs(H);%求巴特沃什模擬濾波器的幅頻響應振幅； pha=angle(H);%求巴特沃什模擬濾波器的相頻響應相位； db=20*log10((mag+eps)/max(mag));%將幅頻特性轉換為對數形式； figure('name','巴特沃什模擬濾波器幅頻相頻特性')%定義巴特沃什模擬濾波器幅頻相特性圖形標題名 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;%定義信號采樣頻率；通帶截止頻率；通帶起始頻率； wp=2*pi*fp;ws=2*pi*fs;%分別求出模擬通帶截止頻率wp和模擬阻帶截止頻率ws [N,wc]=buttord(wp,ws,Rp,Rs,'s');%計算巴特沃什模擬濾波器的階次 N%顯示雙線性數字濾波器的階次 [b1,a1]=butter(N,wc,'s');%計算H（s）系數向量b1和a1 [b,a]=bilinear(b1,a1,Fs);%利用雙線性映射法將巴特沃什模擬濾波器轉換為數字濾波器 [H,m]=freqz(b,a);%計算巴特沃什數字濾波器的頻率響應 mag=abs(H);%求巴特沃什數字濾波器的幅頻響應 pha=angle(H);%求巴特沃什數字濾波器的相頻響應 db=20*log10((mag+eps)/max(mag));%將幅頻特性轉換為對數形式 figure('name','雙線性數字濾波器幅頻相頻特性') subplot(211);plot(m*100/(2*pi),db);grid on;axis([0,56,-100,50]);xlabel('w/2pi');ylabel('db');title('雙線性數字濾波器幅頻特性');%繪制雙線性數字濾波器幅頻特性; subplot(313);plot(m*100/(2*pi),pha);grid on;xlabel('w/2pi');ylabel('rad');title('雙線性數字濾波器相頻特性');%繪制雙線性數字濾波器相頻特性;

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('濾波前心電圖采樣序列'); %畫出雙線性映射法將巴特沃什模擬濾波器轉換為數字濾波器的幅頻特性曲線 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');%計算巴特沃什模擬濾波器的階次 N%顯示巴特沃什濾波器的階次 [b1,a1]=butter(N,wc,'s');%計算H（s）系數向量b1和a1 [b,a]=bilinear(b1,a1,Fs);%利用雙線性映射法將巴特沃什模擬濾波器轉換為數字濾波器 [H,m]=freqz(b,a);%計算巴特沃什數字濾波器的頻率響應 mag=abs(H);%求巴特沃什數字濾波器的幅頻響應 db=20*log10((mag+eps)/max(mag));%將幅頻特性轉換為對數形式 subplot(222);plot(m*100/(2*pi),db);grid on;axis([0,56,-100,50]);xlabel('w/2pi');ylabel('db');title('雙線性數字濾波器幅頻特性');%繪制雙線性數字濾波器幅頻特性 %利用所設計出的濾波器差分方程 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數字濾波后心電采樣序列'); 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數字濾波后心電信號頻譜圖');

clear all;clc;close all;%清空工作空間； 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');%計算巴特沃什模擬濾波器的階次 N%顯示巴特沃什濾波器的階次 [b1,a1]=butter(N,wc,'s');%計算H（s）系數向量b1和a1 [b,a]=impinvar(b1,a1,Fs);%利用沖激響應不變法將巴特沃什模擬濾波器轉換為數字濾波器 [H,m]=freqz(b,a);%計算巴特沃什數字濾波器的頻率響應 mag=abs(H);%求巴特沃什數字濾波器的幅頻響應 pha=angle(H);%求巴特沃什數字濾波器的相頻響應 db=20*log10((mag+eps)/max(mag));%將幅頻特性轉換為對數形式 figure('name','沖激響應不變法數字濾波器幅頻相頻特性') subplot(211);plot(m*100/(2*pi),db);grid on;axis([0,56,-100,50]);xlabel('w/2pi');ylabel('db');title('沖激響應不變法數字濾波器幅頻特性'); subplot(313);plot(m*100/(2*pi),pha);grid on;xlabel('w/2pi');ylabel('rad');title('沖激響應不變法數字濾波器相頻特性');

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('濾波前心電圖采樣序列'); %畫出雙線性映射法將巴特沃什模擬濾波器轉換為數字濾波器的幅頻特性曲線 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');%計算巴特沃什模擬濾波器的階次 N%顯示巴特沃什濾波器的階次 [b1,a1]=butter(N,wc,'s');%計算H（s）系數向量b1和a1 [b,a]=impinvar(b1,a1,Fs);%利用沖激響應不變法將巴特沃什模擬濾波器轉換為數字濾波器 [H,m]=freqz(b,a);%計算巴特沃什數字濾波器的頻率響應 mag=abs(H);%求巴特沃什數字濾波器的幅頻響應 db=20*log10((mag+eps)/max(mag));%將幅頻特性轉換為對數形式 subplot(222);plot(m*100/(2*pi),db);grid on;axis([0,56,-100,50]);xlabel('w/2pi');ylabel('db');title('沖激響應不變數字濾波器幅頻特性'); %利用所設計出的濾波器差分方程 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數字濾波后心電采樣序列'); 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數字濾波后心電信號頻譜圖');

四、分析

分析略

五、思考題

思考題略

六、總結

總結略

總結

以上是生活随笔為你收集整理的IIR数字滤波器的设计及应用——MATLAB的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。