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<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分别是通带和阻带区的波纹系数。

%******************************************************************* %
%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-->');
%******************************************************************* %
1、ButterWorth模拟和数字滤波
（1）butterd函数：ButterWorth滤波阶数的选择。

调用格式：[n,Wn]=butterd(Wp,Ws,Rp,Rs)，在给定滤波性能的情况下（通带临界频率Wp、阻带临界频率Ws、通带内最大衰减Rp和阻带内最小衰减Rs），计算ButterWorth滤波的阶数n和截止频率Wn。

相同参数条件下的模拟滤波则调用格式为：[n,Wn]=butterd(Wp,Ws,Rp,Rs,’s’)
（2）butter函数：ButterWorth滤波。

调用格式：[b,a]=butter(n,Wn)，根据阶数n和截止频率Wn计算ButterWorth滤波分子分母系数（b为分子系数的矢量形式，a为分母系数的矢量形式）。

相同参数条件下的模拟滤波则调用格式为：[b,a]=butter(n,Wn,’s’)
[例] 采样频率为1Hz，通带临界频率fp =0.2Hz，通带内衰减小于1dB（αp=1）；阻带临界频率fs=0.3Hz，阻带内衰减大于25dB（αs=25）。

一个数字滤波满足以上参数。

[n,Wn]=buttord(0.2,0.3,1,25);
[b,a]=butter(n,Wn);
freqz(b,a,512,1); freqz(b,a,N,fs);%fs为采样频率。