Matlab中的滤波器设计

Matlab中的滤波器设计

http://hi.baidu.com/cfin_nifc/item/bf737648734006e11381daba

IIR滤波器有以下几个特点： 

1．IIR数字滤波器的系统函数可以写成封闭函数的形式。 

2．IIR数字滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。 

3．IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。 

4．IIR数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络。 

在MATLAB下设计IIR滤波器可使用Butterworth函数设计出巴特沃斯滤波器，使用Cheby1函数设计出契比雪夫I型滤波器，使用Cheby2设计出契比雪夫II型滤波器，使用ellipord函数设计出椭圆滤波器。下面主要介绍前两个函数的使用。 

与FIR滤波器的设计不同，IIR滤波器设计时的阶数不是由设计者指定，而是根据设计者输入的各个滤波器参数（截止频率、通带滤纹、阻带衰减等），由软件设计出满足这些参数的最低滤波器阶数。在MATLAB下设计不同类型IIR滤波器均有与之对应的函数用于阶数的选择。 

一、巴特沃斯IIR滤波器的设计 

在MATLAB下，设计巴特沃斯IIR滤波器可使用butter函数。 

Butter函数可设计低通、高通、带通和带阻的数字和模拟IIR滤波器，其特性为使通带内的幅度响应最大限度地平坦，但同时损失截止频率处的下降斜度。在期望通带平滑的情况下，可使用butter函数。 

butter函数的用法为： 

[b,a]=butter(n,Wn,/ftype/) 

其中n代表滤波器阶数，Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。buttord函数可在给定滤波器性能的情况下，求出巴特沃斯滤波器的最小阶数n，同时给出对应的截止频率Wn。buttord函数的用法为： 

[n,Wn]= buttord(Wp,Ws,Rp,Rs) 

其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。 

不同类型（高通、低通、带通和带阻）滤波器对应的Wp和Ws值遵循以下规则： 

1．高通滤波器：Wp和Ws为一元矢量且Wp>Ws； 

2．低通滤波器：Wp和Ws为一元矢量且Wp<Ws； 

3．带通滤波器：Wp和Ws为二元矢量且Wp<Ws，如Wp=[0.2,0.7],Ws=[0.1,0.8]; 

4．带阻滤波器：Wp和Ws为二元矢量且Wp>Ws，如Wp=[0.1,0.8],Ws=[0.2,0.7]。 

二、契比雪夫I型IIR滤波器的设计 

在期望通带下降斜率大的场合，应使用椭圆滤波器或契比雪夫滤波器。在MATLAB下可使用cheby1函数设计出契比雪夫I型IIR滤波器。 

cheby1函数可设计低通、高通、带通和带阻契比雪夫I型滤IIR波器，其通带内为等波纹，阻带内为单调。契比雪夫I型的下降斜度比II型大，但其代价是通带内波纹较大。 

cheby1函数的用法为： 

[b,a]=cheby1(n,Rp,Wn,/ftype/) 

在使用cheby1函数设计IIR滤波器之前，可使用cheblord函数求出滤波器阶数n和截止频率Wn。cheblord函数可在给定滤波器性能的情况下，选择契比雪夫I型滤波器的最小阶和截止频率Wn。 

cheblord函数的用法为： 

[n,Wn]=cheblord(Wp,Ws,Rp,Rs) 

其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。

-----------------------------------------------------------------------

MATLAB参考程序和仿真内容 

%*******************************************************************% 

%mode: 1--巴特沃斯低通；2--巴特沃斯高通；3--巴特沃斯带通；4--巴特沃斯带阻 

% 5--契比雪夫低通；6--契比雪夫高通；7--契比雪夫带通；8--契比雪夫带阻 

%fp1,fp2： 通带截止频率，当高通或低通时只有fp1有效 

%fs1, fs2： 阻带截止频率，当高通或低通时只有fs1有效 

%rp: 通带波纹系数 

%as: 阻带衰减系数 

%sample: 采样率 

%h: 返回设计好的滤波器系数 

%*******************************************************************% 

function[b,a]=iirfilt(mode,fp1,fp2,fs1,fs2,rp,as,sample) 

wp1=2*fp1/sample;wp2=2*fp2/sample; 

ws1=2*fs1/sample;ws2=2*fs2/sample; 

%得到巴特沃斯滤波器的最小阶数N和3bd频率wn 

if mode<3[N,wn]=buttord(wp1,ws1,rp,as); 

elseif mode<5[N,wn]=buttord([wp1 wp2],[ws1 ws2],rp,as); 

%得到契比雪夫滤波器的最小阶数N和3bd频率wn 

elseif mode<7[N,wn]=cheb1ord(wp1,ws1,rp,as); 

else[N,wn]=cheblord([wp1 wp2],[ws1 ws2],rp,as); 

end 

%得到滤波器系数的分子b和分母a 

if mode= =1[b,a]=butter(N,wn);end 

if mode= =2[b,a]=butter(N,wn,/high/);end 

if mode= =3[b,a]=butter(N,wn);end 

if mode= =4[b,a]=butter(N,wn,/stop/);end 

if mode= =5[b,a]=cheby1(N,rp,wn);end 

if mode= =6[b,a]=cheby1(N,rp,wn,/high/);end 

if mode= =7[b,a]=cheby1(N,rp,wn);end 

if mode= =8[b,a]=cheby1(N,rp,wn,/stop/);end 

set(gcf,/menubar/,menubar); 

freq_response=freqz(b,a); 

magnitude=20*log10(abs(freq_response)); 

m=0:511; 

f=m*sample/(2*511); 

subplot(3,1,1);plot(f,magnitude);grid; %幅频特性 

axis([0 sample/2 1.1*min(magnitude) 1.1*max(magnitude)]); 

ylabel('Magnitude');xlabel('Frequency-->'); 

phase=angle(freq_response); 

subplot(3,1,2);plot(f,phase);grid; %相频特性 

axis([0 sample/2 1.1*min(phase) 1.1*max(phase)]); 

ylabel('Phase');xlabel('Frequency-->'); 

h=impz(b,a,32); %32点的单位函数响应 

t=1:32; 

subplot(3,1,3);stem(t,h);grid; 

axis([0 32 1.2*min(h) 1.1*max(h)]); 

ylabel('h(n)');xlabel('n-->'); 

%*******************************************************************% 

假设需设计一个巴特沃斯低通IIR滤波器，通带截止频率为2KHz，阻带截止频率为3KHz，通带波纹系数为1，阻带衰减系数为20，采样频率为10KHz，则只需在MATLAB的命令窗口下键入： 

[b,a]=iirfilt(1,2000,3000,2400,2600,1,20,10000) 

程序进行模拟，并且按照如下顺序输出数字滤波器系统函数的系数 

b= b0 b1 ……bn 

a= a0 a1 ……an 

关于软件仿真实验内容，建议在完成大量仿真例子的基础上，选择能够体现实验要求的4个例子进行记录，系统函数只要记录系统的阶数。