Matlab的滤波器设计

基于Matlab的滤波

一、巴特沃斯IIR滤波的

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

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

butter函数的用法为：

[b,a]=butter(n,Wn,/ftype/) % ’ftype’must be one of 'high', 'stop', 'low',

'bandpass', 'z' or 's'

其中n代表滤波阶数，W n代表滤波的截止频率，这两个参数可使用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

3．带通滤波：Wp和Ws为二元矢量且WpWs，如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分别是通带和阻带区的波纹系数。

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

%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;