d．带阻滤波器：
Wp和Ws为二元矢量且Wp>Ws，如Wp=[0.1,0.8],Ws=[0.2,0.7]。
例如，采样频率10Hz，通带截止频率fp=3Hz，阻带截止频率fs=4Hz，通带衰减小于1dB，阻带衰减大于20dB，使用双线性变换法由模拟滤波器原型设计数字滤波器。
MATLAB程序见附录，如图1所示。
图1ButterWorth低通滤波器
3.2契比雪夫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分别是通带和阻带区的波纹系数。
譬如，采样频率为10Hz，设计一个数字低通滤波器，要求其通带临界频率fp=3Hz，通带αp=1dB），阻带临界频率fs=4Hz，阻带内衰减大于15dB（αs=15）内衰减小于1dB（αs=1dB）。
其程序见附录，设计出的滤波器如图2
图2Chebyshev低通（I型）滤波器
3.3数字滤波器的设计
数字滤波器（digitalfilter）是由数字乘法器、加法器和延时单元组成的一种装置。
其功能是对输入离散信号的数字代码进行运算处理，以达到改变信号频谱的目的。
数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。
数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。
IIR数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配。
所以IIR滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。
FIR数字滤波器的单位脉冲响应是有限长序列。
它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题，设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。
在对滤波器实际设计时，整个过程的运算量是很大的。
例如利用窗函数法设计M阶FIR低通滤波器时，首先要根据
（1）式计算出理想低通滤波器的单位冲激响应序列
，然后根据
（2）式计算出M个滤波器系数
。
当滤波器阶数比较高时，计算量比较大，设计过程中改变参数或滤波器类型时都要重新计算。
（1）
（2）
设计完成后对已设计的滤波器的频率响应要进行校核，要得到幅频相频响应特性，运算量也是很大的。
我们平时所要设计的数字滤波器，阶数和类型并不一定是完全给定的，很多时候都是要根据设计要求和滤波效果不断的调整，以达到设计的最优化。
在这种情况下，滤波器的设计就要进行大量复杂的运算，单纯的靠公式计算和编制简单的程序很难在短时间内完成设计。
利用MATLAB强大的计算功能进行计算机辅助设计，可以快速有效的设计数字滤波器，大大的简化了计算量，直观简便。
FIR数字滤波器的窗函数设计法
设理想带阻滤波器频率响应为
例如，利用Kaiser窗函数，设计长度为55的阻带滤波器，使阻带衰减为60dB。
参数β可由式α
β=0.1102（α-8.7）其中α=60，确定
MATLAB程序见附录，如图3所示。
与其他高级语言的程序设计相比,MATLAB环境下可以更方便、快捷地设计出具有严格线性相位的FIR滤波器,节省大量的编程时间,提高编程效率,且参数的修改也十分方便.还可以进一步进行优化设计,相信随着版本的不断更新,MATLAB在数字滤波器技术中必将发挥更大的作用.
IIR数字滤波器的设计方法
IIR数字滤波器最大的优点是给定一组指标时，他的阶数要比相同组的FIR滤波器的低得多。
基于Matlab的IIR滤波器的设计方法主要有完全设计法、模拟原型设计法和直接设计法等。
a.模拟原型法
采用经典低通滤波器作为连续域上的设计模型，通过频域变换得到IIR数字滤波器，最后还要进行离散化处理。
Matlab提供的低通模拟滤波器原型函数包括：
besselap，buttap，cheb1lp，cheb2ap，ellipap；频域变换函数包括：
lp2bp，lp2bs，lp2hp，lp2lp；离散化处理函数有bilinear和impinvar。
图3带阻滤波器幅频特性
b.完全设计法
Matlab信号处理工具箱提供了几个直接设计IIR数字滤波器的函数，直接调用就可以设计滤波器，这为设计通用滤波器提供了方便。
设计Butterworth滤波器用函数butter（），可以设计低通、高通、带通和带阻的数字和模拟滤波器，其特性是通带内的幅度响应最大限度的平滑，但损失了截止频率处的下降斜度。
设计ChebyshevI型滤波器用函数chebyl（）。
可以设计低通、高通、带通和带阻的数字和模拟ChebyshevI型滤披器，其通带内为等波纹，阻带内为单调。
ChebyshevI型滤波器的下降斜度比II型大，但其代价是通带内波纹较大。
设计ChebyshevII型滤波器用函数cheby2（）。
可以设计低通、高通、带通和带阻的数字和模拟ChebyshevII型滤波器，其通带内为单调，阻带内等波纹。
ChebyshevII型滤波器的下降斜度比I型小，但其阻带内波纹较大。
设计椭圆滤波器用函数ellip（），与cheby1，cheby2类似，可以设计低通、高通、带通和带阻的数字和模拟滤波器。
与Butterworth和chebyshev滤波器相比，ellip函数可以得到下降斜度更大的滤波器，得到通带和阻带均为等波纹。
一般情况下，椭圆滤波器能以最低的阶实现指定的性能指标。
c.直接设计法
直接设计方法的思想是基于给定的滤波器参数直接在离散域上寻找合适的数字滤波器，他不限于常规的滤波器类型，如低通、高通、带通和带阻等。
这种方法甚至可以设计多带的频率响应，Matlab提供yulewalk函数用于辅助设计。
设计一个数字低通滤波器满足以下参数，采样频率为1Hz，通带临界频率fp=0.2Hz，通带内衰减小于1dB（αp=1）；阻带临界频率fs=0.3Hz，阻带内衰减大于25dB（αs=25）。
Matlab使用归一化的频率参数（临界频率除以采样频率的1/2）。
这样临界频率参数的取值范围在0和1之间，1代表Fs/2（用角频率表示的时候对应π）。
MATLAB源程序见附录，如图4所示。
图4直接设计的数字高通滤波器
d.通用Butterworth设计方法
使用这种方法设计的Butterworth数字滤波器可以有不同数目的零点和极点，Matlab提供的maxflat函数实现了这一功能。
这个函数与butter函数很相似，但他可以指定两个阶参数，其中归一化和非归一化各一个。
如果这两个参数的值相同，那么他与butter函数的结果就是相同的。
例如butter（）函数设计的滤波器，其MATLAB程序见附录，增益图如图5所示。
e.参数建模法
寻找接近于所需要设计的滤波器的通用模型，时域上的建模函数为lpc，prony，Stmcb；频域上的建模函数有invfreqs和invfreqz。
典型IIR数字滤波器的设双线性变换法为了克服冲激响应不变法的频率混叠现象，需要使s平面与z平
面建立一一对应的单值映射关系，即求出s=f（z），然后将它带入H（s），就可以求得H（z），即
H（z）=H（s）|s=f（z）
为了克服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里，其次再通过上面讨论过的标准变换关系z=es1T将此横带变换到整个z平面上去，这样就使s平面与z平面是一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。
例如，双线性变换法设计一个ChebyshevⅡ高通滤波器，使其幅频特性逼近一个具有以下技术指标的模拟ChebyshevⅡ高通滤波器：
Ws=2*pi*1kHz，Wp=2*pi*1.4kHz，在Ws处的最小衰减为15dB，在Wp处的最大衰减不超过0.3dB，抽样频率为20kHz。
其MATLAB程序见附录，频率响应如图6
图5Butterworth带阻滤波器增益图
图6ChebyshevⅡ高通滤波器的频率响应
可以看到经过离散采样、数字滤波后分离出了83.3hz的频率分量。
之所以选取上面的叠加信号
作为原始信号，是由于在实际工作中是要对已经经过差分滤波的信号进一步做带通滤波，信号的各分量基本同
一致，可以反映实际的情况。
本例设计的滤波器已在实际工作中应用，取得了不错的效果。
结论
采用MATLAB设计滤波器，使原来非常繁琐复杂的程序设计变成了简单的函数调用，为滤波器的设和实现开辟了广阔的天地，尤其是Matlab工具箱使各个领域的研究人员可以直观方便地进行科学研究与工程应用。
其中的信号处理工具箱、图像处理工具箱、小波工具箱等更是为数字滤波研究的蓬勃发展提供了可能。
MATLAB信号处理工具箱为滤波器设计及分析提供了非常优秀的辅助设计工具,在设计数字滤波器时,善于应用MATLAB进行辅助设计,能够大大提高设计效率。
参考文献
[1]陈德树.计算机继电保护原理与技术.北京：
水利电力出版社，1992.
[2]蒋志凯.数字滤波与卡尔曼滤波.北京：
中国科学技术出版社，1993
[3]楼顺天、李博菡.基于MATLAB的系统分析与设计－信号处理.西安：
西安电子
[4]董长虹等.MATLAB信号处理与应用.北京：
国防工业出版社，2005
[5][美]M.H.海因斯著，张建华等译.数字信号处理.北京：
科学出版社，2002
[6]张葛祥，李娜.MATLAB仿真技术与应用.北京：
清华大学出版社，2002
附录
1、ChebyshevⅡ高通滤波器的频率响应
ws=2*pi*1000;
ws1=ws*2*pi;
wp=2*pi*1400;
wp1=wp*2*pi;
rp=0.3;
rs=15;
fs=20000;
[N,Wn]=cheb2ord（wp1,ws1,rp,rs,’s’）;
[z,p,k]=cheb2ap（N,rs）;
[A,B,C,D]=zp2ss（z,p,k）;
[At,Bt,Ct,Dt]=lp2hp（A,B,C,D,Wn）;
[At1,Bt1,Ct1,Dt1]=bilinear（At,Bt,Ct,Dt,fs）;
[num,den]=ss2tf（At1,Bt1,Ct1,Dt1）;
freqz（num,den）;
[H,W]=freqz（num,den）;
plot（W*fs/（2*pi）,abs（H））;grid;
xlabel（‘频率/Hz’）;
ylabel（‘幅值’）;
2、带阻滤波器幅频特性实现
n=55-1;
w=[0.40.6];
beta=0.1102*（60-8.7）;
kaiw=Kaiser（n+1,beta）;
b=fir1（n,w,’stop’,kaiw）;
[h,w]=freqz（b,1,512,2）;
plot（w,20*log10（abs（h）））;grid;
xlabel（‘频率（归一化）’）;
ylabel（（‘幅度dB’））;
3、Chebyshev低通（I型）滤波器
T=0.1;FS=1/T;
fp=3;fs=4;
Rp=1;
As=15;
wp=fp/FS*2*pi;
ws=fs/FS*2*pi;
OmegaP=（2/T）*tan（wp/2）;
OmegaS=（2/T）*tan（ws/2）;
ep=sqrt（10^（Rp/10）-1）;
A=10^（As/20）;
OmegaC=OmegaP;
OmegaR=OmegaS/OmegaP;
g=sqrt（A*A-1）/ep;
N=ceil（log10（g+sqrt（g*g-1））/log10（OmegaR+sqrt（OmegaR*OmegaR-1）））;
[z,p,k]=cheb1ap（N,Rp）;
a=real（poly（p））;
aNn=a（N+1）;
p=p*OmegaC;
a=real（poly（p））;
aNu=a（N+1）;
k=k*aNu/aNn;
b0=k;
B=real（poly（z））;
b=k*B;
[bz,az]=bilinear（b,a,FS）;
%freqz（bz,az,200,FS,'whole'）;
H=freqz（bz,az,200,'whole'）;
plot（abs（H））;
xlabel（‘频率/Hz’）;
ylabel（‘幅值’）;
4、直接设计数字滤波器
FS=1;
[n,Wn]=buttord（0.3/（FS/2）,0.2/（FS/2）,1,25）;
[b,a]=butter（n,Wn,'high'）;
freqz（b,a,512,FS）;
5、ButterWorth低通滤波器
T=0.1;FS=1/T;
fp=3;fs=4;
wp=fp/FS*2*pi;
ws=fs/FS*2*pi;
Rp=1;
As=20;
OmegaP=（2/T）*tan（wp/2）;
OmegaS=（2/T）*tan（ws/2）;
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（OmegaP/OmegaS）））;
OmegaC=OmegaP/（（10^（Rp/10）-1）^（1/（2*N）））;
[z,p,k]=buttap（N）;
p=p*OmegaC;
k=k*OmegaC^N;
B=real（poly（z））;
b0=k;
cs=k*B;
ds=real（poly（p））;
[b,a]=bilinear（cs,ds,FS）;
freqz（b,a,512,FS）;
6、Butterworth带阻滤波器增益图
function[g,w]=gain（num,den）
%ComputesthegainfunctionindBofa
%transferfunctionat256equallyspacedpoints
%ontheunitcircle
%Numeratorcoefficientsareinvectornum
%Denominatorcoefficientsareinvectorden
%Frequencyvaluesarereturnedinvectorw
%Gainvaluesarereturnedinvectorg
w=0:
pi/255:
pi;
h=freqz（num,den,w）;
g=20*log10（abs（h））;
%DesignofaButterworthBandstopDigitalFilter
Ws=[0.40.6];Wp=[0.20.8];Rp=0.4;Rs=50;
%EstimatetheFilterOrder
[N1,Wn1]=buttord（Wp,Ws,Rp,Rs）;
%DesigntheFilter
[num,den]=butter（N1,Wn1,'stop'）;
%Displaythetransferfunction
disp（'NumeratorCoefficientsare'）;disp（num）;
disp（'DenominatorCoefficientsare'）;disp（den）;
%Computethe
升级会员 