图5-9 椭圆模拟原型滤波器平方幅频图

 

程序运行结果见图5-9。可见阶数为4的椭圆滤波器的过渡带已相当窄（陡），但这种特性的获得是以牺牲通带和阻带的单调平滑特性为代价的。可以看到滤波器的阶数越高平方幅频响应越接近于矩形。

5.3.5  Bessel滤波器

前面讲过的各类原型滤波器均没有绘出其相位随频率的变化特性（相频特性）。在后面的数字信号处理学习中将会看到它们的相位特性是非线性的。本节所介绍的Bessel滤波器就能最大限度地减少相频特性的非线性，使得通带内通过的信号形状不变(拷贝不走样)。

Bessel模拟低通滤波器的特点是在零频时具有最平坦的群延迟，并在整个通带内群延迟几乎不变。在零频时的群延迟为。由于这一特点，Bessel模拟滤波器通带内保持信号形状不变。但数字Bessel滤波器没有平坦特性，因此MATLAB信号处理工具箱只有模拟Bessel滤波器设计函数。

函数besselap用于设计Bessel模拟低通滤波器原型，调用格式为：

[z,p,k]=besselap(N)

式中，N为滤波器的阶数，应小于25。z,p,k为滤波器的零点、极点和增益。

滤波器的传递函数无零点，具有与（5-10）式相同的形式。下面用实例观看Bessel滤波器的幅频和相频特性。

【例5-5】绘制5阶和10阶Bessel低通滤波器原型的平方幅频和相频图。

%Samp5_5

clf

n=0:0.01:2;   %设置频率点

for ii=1:2

   switch ii

   case 1,pos=1;N=5;

   case 2,pos=3;N=10;

   end

   [z,p,k]=besselap(N);   %设计Bessel模拟滤波器

   [b,a]=zp2tf(z,p,k);   %将零点极点增益形式转化为传递函数形式

   [H,w]=freqs(b,a,n);   %求得传递函数的复数频率响应

   magH2=(abs(H)).^2;

   phaH=unwrap(angle(H));    %求得函数的相位角并进行解缠绕运算

   phaH=phaH*180/pi;       %将相位角由弧度转化为度

   posplot=['2,2,' num2str(pos)];   %设置绘图位置字符串

   subplot(posplot);

   plot(w,magH2); grid on       %绘出平方幅频响应

   title(['N=' num2str(N)]);

   xlabel('w/wc');   ylabel('Bessel |H(jw)|^2');

   grid on

   subplot(['2,2,' num2str(pos+1)]);

   plot(w,phaH);

   xlabel('w/wc');ylabel('Bessel 相位/^o');

 title(['N=' num2str(N)]);

   grid on

end

图5-10 Bessel模拟原型滤波器相频图

 

可见，Bessel滤波器具有最优线性相频的特点，但这个特点的获得是以牺牲窄过渡带为代价的,即滤波器的幅频平方特性与矩形特性相差甚远。

对所有的模拟原型滤波器做一总结可知：Butterworth滤波器在通带和阻带内均具有平滑单调的特点，但在相同过渡带宽的条件下，该滤波器所需的阶数最多。Chebyshev I和II型滤波器在通带或阻带内具有波纹，但在相同过渡带宽的条件下，该滤波器所需的阶数比Butterworth滤波器要少。椭圆滤波器在通带和阻带内均有波纹出现，但在相同过渡带宽的条件下，该滤波器所需的阶数最少。Bessel滤波器具有最宽的过渡带，但具有最优的线性相频特性。因此没有绝对“好”的滤波器，要根据解决问题的不同选择不同的滤波器，因此，每一种滤波器的设计方法我们都要熟练掌握。

5.4  频 率 变 换

前面所讲的模拟原型滤波器均是截止频率为1的滤波器，在实际设计中是很难遇到的，然而它是设计其他各类滤波器的基础。我们通常遇到的是截止频率任意的低通滤波器、高通滤波器、带通滤波器和带阻滤波器。如何以由低通原型模拟滤波器为基础设计这些滤波器呢？这就要用到我们今天要讲的频率变换。所谓频率变换是指各类滤波器（低通、高通、带通、带阻）和低通滤波器原型的传递函数中频率自变量之间的变换关系。通过频率变换，我们可以从模拟低通滤波器原型获得模拟的低通滤波器、高通滤波器、带通滤波器和带阻滤波器，再借助于s域至z域的变换关系又可以设计各类后面所讲的无限冲激响应数字滤波器，这是滤波器设计的重要方法之一。

MATLAB信号处理工具箱有lp2lp,lp2hp,lp2bp,lp2bs四个频率变换函数。下面分别叙述其使用方法及各参量的意义。

（1）       函数lp2lp用于实现由低通模拟原型滤波器至低通滤波器的频率变换，调用格式为

[bt,at]=lp2lp(b,a,)

其中，a,b为模拟原型滤波器的分母和分子多项式的系数，为低通滤波器所期望的截止频率（rad/s），若给定的单位为Hz，应乘以2_{http://www.fzjpk.com/dizhen/word/chap5.files/image103.gif}。bt,at为返回的低通滤波器的分母和分子多项式的系数。该函数将模拟原型滤波器传递函数执行下面变换：

                                                    (5-15)

式中，H(p)为低通原型滤波器传递函数，H(s)为低通滤波器传递函数。该项操作可以执行模拟原型滤波器由截止频率为1到指定截止频率的变换，其原理讨论已超出本课程的范围，可参看其他信号处理参考书。下面的例子说明如何进行模拟原型低通滤波器变换为截止频率不为1的模拟低通滤波器。

【例5-6】将4阶椭圆模拟原型滤波器变换为截止频率为0.5的椭圆模拟低通滤波器，其中通带波纹Rp=2dB,阻带衰减Rs=30dB。

%Samp5_6

Rp=2;Rs=30;   %模拟原型滤波器的通带波纹为2dB,阻带衰减为30dB。

[z,p,k]=ellipap(4,Rp,Rs);   %设计椭圆滤波器

[b,a]=zp2tf(z,p,k);      %由零点极点增益形式转换为传递函数形式

[H,w]=freqs(b,a,0:0.01:2);   %给出复数频率响应

subplot(2,1,1),plot(w,abs(H).^2);   %绘出平方幅频函数

xlabel('w/wc');ylabel('椭圆 |H(jw)|^2');

title('原型低通椭圆滤波器(wc=1)')

[bt,at]=lp2lp(b,a,0.5);  %将模拟原型低通滤波器的截止频率变换为0.5

[Ht,wt]=freqs(bt,at,0:0.01:2);  %给出复数频率响应

subplot(2,1,2),plot(wt,abs(Ht).^2); %绘出平方幅频函数

xlabel('w/wc');ylabel('椭圆 |H(jw)|^2');

title('低通椭圆滤波器(wc=0.5)')

         图5-11 将4阶椭圆模拟原型滤波器变换为截止频率为0.5的椭圆模拟低通滤波器的结果

 

程序的运行结果见图5-11。清楚地表明将截止频率由1变换到0.5。

（2）       函数lp2hp用于实现由低通模拟滤波器至高通滤波器的频率变换。调用格式：

[bt,at]=lp2hp(b,a,)

式中，_{http://www.fzjpk.com/dizhen/word/chap5.files/image103.gif}。该函数将模拟原型滤波器传递函数执行下面变换：

                                                   (5-16)

【例5-7】将6阶Chebyshev I型原型滤波器变换为截止频率为0.8的模拟高通滤波器，其中通带波纹Rp=0.5dB。

%Samp5_7

Rp=0.5;   %设置滤波器的通带波纹为0.5dB

[z,p,k]=cheb1ap(6,0.5);   %设计Chebyshev I型模拟原型滤波器

[b,a]=zp2tf(z,p,k);    %由零点极点增益形式转化为传递函数形式

[H,w]=freqs(b,a,0:0.01:2);   %计算传递函数的复数频率响应

subplot(2,1,1),plot(w,abs(H).^2); %绘制传递函数的平方幅频响应

xlabel('w/wc');ylabel('Chebyshev I |H(jw)|^2');

title('Chebyshev I 低通原型滤波器(wc=1)')

[bt,at]=lp2hp(b,a,0.8);   %由低通原型滤波器转换为截止频率为0.8的高通滤波器

[Ht,wt]=freqs(bt,at,0:0.01:2);   %计算滤波器的复数频率响应

subplot(2,1,2),plot(wt,abs(Ht).^2);  %绘制传递函数的平方幅频响应

xlabel('w/wc');ylabel('Chebyshev I |H(jw)|^2');

title('Chebyshev I 高通滤波器(wc=0.8)')

图5-12 将6阶Chebyshev I型原型滤波器变换为截止频率为0.8的模拟高通滤波器的结果图

 

程序运行结果见图5-12。可见实现了从低通滤波器到高通滤波器的转换，并且截止频率也符合例题要求。

（3）       函数lp2bp用于实现由低通模拟原型滤波器至带通滤波器的频率变换。调用格式：

[bt,at]=lp2bp(b,a,_{http://www.fzjpk.com/dizhen/word/chap5.files/image101.gif},Bw)

式中，_{http://www.fzjpk.com/dizhen/word/chap5.files/image101.gif}为带通滤波器的中心频率(rad/s)，Bw为带通滤波器带宽(rad/s)。而

                                          (5-17)

式中, 为带通滤波器的下边界频率，为带通滤波器上边界频率。若给定的边界频率为Hz需乘以2_{http://www.fzjpk.com/dizhen/word/chap5.files/image103.gif}。

该函数将模拟原型滤波器传递函数执行下面变换运算：

                                              (5-18)

下面用例题说明截止频率和滤波器类型的转换。这里要注意，输出的带通滤波器阶数为模拟原型滤波器阶数的2倍。

【例5-8】将6阶Chebyshev II型原型滤波器变换为模拟带通滤波器，其中上边界截止频率为0.8 rad/s,下边界截止频率为1.4 rad/s,阻带衰减Rs=20dB。

%Samp5_8

Rs=20;    %滤波器的阻带衰减为20dB

[z,p,k]=cheb2ap(6,Rs);   %设计6阶阻带衰减为20dB的Chebyshev II型原型滤波器

[b,a]=zp2tf(z,p,k);    %将零点极点增益形式转换为传递函数形式

[H,w]=freqs(b,a,0:0.01:2);   %计算模拟原型滤波器的复数频率响应

subplot(2,1,1),plot(w,abs(H).^2);  %绘制平方幅频响应

xlabel('w/wc');ylabel('Chebyshev II |H(jw)|^2');

title('Chebyshev II 型低通原型滤波器(wc=1)')

w1=0.8;w2=1.4;   %给定欲设计滤波器通带的下限和上限频率

w0=sqrt(w1*w2);bw=w2-w1; %计算中心点频率和频带宽度

[bt,at]=lp2bp(b,a,w0,bw);  %频率转换

[Ht,wt]=freqs(bt,at,0:0.01:2);  %计算滤波器的复数频率响应

subplot(2,1,2),plot(wt,abs(Ht).^2); %绘制平方幅频响应

xlabel('w/wc');ylabel('Chebyshev II |H(jw)|^2');

title('Chebyshev II型带通滤波器(wc=0.8~1.4)')

图5-13 将6阶Chebyshev II型原型滤波器变换为模拟带通滤波器的结果图

 

程序运行结果为图5-13。可见该程序实现了由低通原型滤波器到带通滤波器的转换。

（4）       函数lp2bs用于实现由低通模拟原型滤波器至带阻滤波器的频率变换。调用格式

[bt,at]=lp2bs(b,a,_{http://www.fzjpk.com/dizhen/word/chap5.files/image101.gif},Bw)

式中，_{http://www.fzjpk.com/dizhen/word/chap5.files/image101.gif}为带阻滤波器的中心频率(rad/s)，Bw为带阻滤波器带宽(rad/s)。而

                                          (5-19)

式中, 为带阻滤波器的下边界频率，为带阻滤波器上边界频率。若给定的边界频率为Hz需乘以2_{http://www.fzjpk.com/dizhen/word/chap5.files/image103.gif}。

该函数将模拟原型滤波器传递函数执行下面变换运算：

                     (5-20)

 

注意：输出的带阻滤波器和带通滤波器是滤波器原型阶数的2倍。

【例5-9】 将6阶Butterworth原型滤波器变换为模拟带阻滤波器，其中上边界频率为0.7 rad/s,下边界频率为1.5 rad/s。

%Samp5_9

[z,p,k]=buttap(6);   %设计阶数为6的Butterworth模拟原型滤波器

[b,a]=zp2tf(z,p,k);  %将零点极点增益形式转换为传递函数形式

[H,w]=freqs(b,a,0:0.01:2); %计算指定频率点的复数频率响应

subplot(2,1,1),plot(w,abs(H).^2); %绘制平方幅频响应

xlabel('w/wc');ylabel('Butterworth |H(jw)|^2');

title('Btterworth低通原型滤波器(wc=1)')

w1=0.7;w2=1.5;   %带阻滤波器的下限及上限频率

w0=sqrt(w1*w2);bw=w2-w1;   %计算中心点频率和阻带宽度

[bt,at]=lp2bs(b,a,w0,bw);   %频率转换

[Ht,wt]=freqs(bt,at,0:0.01:2);  %计算带阻滤波器指定点的复数频率响应

subplot(2,1,2),plot(wt,abs(Ht).^2);%绘制平方幅频响应

xlabel('w/wc');ylabel('Butterworth |H(jw)|^2');

title('Butterworth带阻滤波器(wc=0.7~1.5) ')

图5-14 将6阶Butterworth原型滤波器变换为模拟带阻滤波器的结果图

 

程序运行结果为图5-10。可见，该程序实现了由低通原型滤波器到带阻滤波器之间的转换，并且频带参数与所要求的完全一致。

5.5  滤波器最小阶数选择

前面所述的模拟滤波器设计中，滤波器阶数是 我们在编程时任意指定的。其实它是决定滤波器品质的主要参数之一。通常在满足性能指标的前提下，阶数应该尽可能小，以满足易于实现、提高运算速度的要求。 而在滤波器阶数和滤波器性能之间存在一定的函数关系，我们通过这一函数关系可以求出满足滤波性能指标的最低阶数。下面介绍常用低通滤波器最小阶数的确定原 理及MATLAB实现；接着介绍MATLAB信号处理工具箱中用来计算最小阶数和截止频率的工具函数。

5.5.1  滤波器最小阶数选择原理：以Butterworth低通模拟滤波器为例

模拟低通滤波器的设计指标为：通带边界频率_{http://www.fzjpk.com/dizhen/word/chap5.files/image046.gif}，通带波纹Rp(dB)、阻带衰减Rs(dB)。

当_{http://www.fzjpk.com/dizhen/word/chap5.files/image132.gif}时

                                 (5-21)

以截止频率_{http://www.fzjpk.com/dizhen/word/chap5.files/image139.gif}，则上式可写成：

                      (5-22)

同理，当_{http://www.fzjpk.com/dizhen/word/chap5.files/image143.gif}时，

                                            (5-23)

综合上面两式可得：

                     (5-24)

式中，N应向上取整，则N为该滤波器的最小阶数。_{http://www.fzjpk.com/dizhen/word/chap5.files/image066.gif}表示为：

                         (5-25)

或

                                                   (5-26)

下面用MATLAB编程可以计算滤波器的最小阶数N和截止频率_{http://www.fzjpk.com/dizhen/word/chap5.files/image066.gif}。

【例5-10】设计一个模拟Butterworth滤波器。设计指标为通带边界频率_{http://www.fzjpk.com/dizhen/word/chap5.files/image159.gif}，通带波纹1dB，阻带衰减16dB,试确定最小阶数N和截止频率。

%Samp5_10

wp=200*pi;ws=300*pi;Rp=1;Rs=16;%滤波器的通带阻带边界频率、通带波纹和阻带衰减

N=ceil(log10((10^(Rp/10)-1)/(10^(Rs/10)-1))/(2*log10(wp/ws)))

% （5-22）式，这里ceil函数表示向上取整

Wcp=wp/((10^(Rp/10)-1)^(1/(2*N)))   %（5-23）式

Wcs=ws/((10^(Rs/10)-1)^(1/(2*N)))   %（5-24）式

给程序的输出为：

N =   7

Wcp =  691.9839

Wcs =  725.7292

程序输出的两个截止频率分别为运用通带波纹和阻带衰减得到的。为使信号保留较多的有用信息，对于低通滤波器可选用较大的截止频率，即725.7292。

其他几种滤波器的最小阶数选择思路与此相同，由于篇幅所限，我们略去。

5.5.2  滤波器最小阶数选择函数

上面给出了Butterworth滤波器的最小阶数和截止频率的选择公式及程序。其实MATLAB工具箱中运用滤波器的最小阶数选择公式给出了滤波器最小阶数选择函数。几种滤波器最小阶数的选择函数如下：

[n,wc]=buttord(wp,ws,Rp,Rs, 's');     Butterworth 滤波器

[n,wc]=cheb1ord(wp,ws,Rp,Rs, 's');     Chebyshev I 滤波器

[n,wc]=cheb2ord(wp,ws,Rp,Rs, 's');     Chebyshev II 滤波器

[n,wc]=ellipord(wp,ws,Rp,Rs, 's');     Elliptic 滤波器

式中，wp为通带边界频率，ws为阻带边界频率，单位为rad/s。Rp,Rs分别为通带波纹和阻带衰减，单位为dB。二者分别表示通带内的最大允许幅值损失和阻带下降的分贝数。's'表示模拟滤波器（缺省时，该函数适用于数字滤波器）；函数返回值n为模拟滤波器的最小阶数；wc为模拟滤波器的截止频率，单位为rad/s。这四个函数适用于低通、高通、带通、带阻滤波器。

若wp<ws，对应于低通模拟滤波器，当wp>ws时对应于高通模拟滤波器，对于带通和带阻滤波器存在两个过渡带，wp和ws均应该为含有两个元素的向量，分别表示两个过渡带的边界频率。即：wp=[通带下界频率, 通带上界频率]，ws=[阻带下界频率, 阻带上界频率]。对于带通滤波器，这四个频带界线的大小排列为:阻带下界频率<通带下界频率<通带上界频率<阻带上界频率；对于带阻滤波器，这四个频带界线的大小排列为:通带下界频率<阻带下界频率<阻带上界频率<通带上界频率。这时返回值wc包括两个元素（第一个元素小于第二个元素），分别为通带和阻带之间的界线频率。函数自动判断是带通还是带阻滤波器

【例5-11】设计一个模拟低通Butterworth滤波器和Chebyshev II型滤波器，满足通带边界频率_{http://www.fzjpk.com/dizhen/word/chap5.files/image159.gif}，通带波纹Rp=1dB，阻带衰减Rs=16dB,试确定最小阶数N和通带下降3dB的截止频率。

%Samp5_11

Wp=200*pi;Ws=300*pi;Rp=1;Rs=16;%滤波器的通带阻带截止频率、通带波纹和阻带衰减

disp('Butterworth滤波器:'); %在命令窗口显示括号内的文本

[N,Wc]=buttord(Wp,Ws,Rp,Rs,'s') %计算Butterworth滤波器的最小阶数和截止频率

disp('Chebyshev II型滤波器:');

[N,Wc]=cheb2ord(Wp,Ws,Rp,Rs,'s') %计算ChebyshevII滤波器的最小阶数和截止频率

程序的输出为：

Butterworth滤波器:

N =7

Wc =725.7292

Chebyshev II型滤波器:

N =4

Wc =  839.8376

程序输出的Butterworth滤波器的最小阶数和截止频率与例5-10的程序输出相同。证明调用求取最小阶数和截止频率的程序就是采用了（5-24）~（5-26）的公式。运用Chebyshev II型滤波器所需的滤波器的最小阶数远小于Butterworth滤波器，说明运用Chebyshev II型滤波器比Butterworth滤波器更容易实现、计算速度更快。大家可以自己设计一些例子运用本节给出的求取滤波器最小阶数的函数求解滤波器最小阶数和截止频率。

5.6  模拟滤波器的性能测试

滤波器设计好之后，一般要进行各方面的测试。在正式设计滤波器之前，我们先介绍如何测试滤波器的性能。对于模拟滤波器，在本章第二节我们已采用函数freqs来求模拟滤波器的频率响应，这里我们详细介绍该函数，若其调用格式为：

[h[,w]]=freqs(b,a[,n])

式中，b,a分别为模拟滤波器传递函数分子和分母多项式系数； h为对应频率点的传递函数值。上面的表示中，[,w]和[,n]表示可有可无的参数，本书中的该类表示均为可有可无的参数, 由用户根据需要确定。w为频率点的值,n为频率点数。若n=128，则用128个频率点来给出此模拟滤波器的频率特性（给定频率点的传递函数值），缺省时为200。若该函数不写输出变量，则执行后绘出该滤波器的幅频响应和相位响应图。此函数模拟滤波器的传递函数形式为：

                       (5-27)

MATLAB工具箱还提供了两个函数abs和angle，由频率响应求幅频响应和相频响应。其中angle的输出单位为rad。可采用rad2deg函数转化为度。另外注意函数的幅频响应经常用分贝(dB)来表示。

求出的幅频响应后，通过下式转换为分贝:(dB)。

【例5-12】已知模拟滤波器的传递函数为：，绘制滤波器的幅频响应和相频响应。

%Samp5_12

b=[0.2 0.3 1];   %滤波器传递函数分子多项式系数

a=[1 0.4 1]; %滤波器传递函数分母多项式系数

figure(1),freqs(b,a);  %第一种输出方法

[h,w]=freqs(b,a); %计算滤波器的复数频率响应

mag=abs(h);pha=angle(h);  %得到滤波器的幅频和相频响应

figure(2),subplot(2,1,1),loglog(w,mag);  %运用双对数坐标绘制幅频响应

grid on;xlabel('角频率/rad\cdots^-^1');ylabel('振幅');

subplot(2,1,2),semilogx(w,pha*180/pi) %运用半对数坐标绘相频响应

grid on;xlabel('角频率/rad\cdots^-^1');

ylabel('相位/^o');

 

图5-15 例5-12中滤波器的幅频响应和相频响应

 

程序运行结果的两幅图完全一致（图5-15），说明freqs的绘图程序与本例绘图程序完全一致。

我们知道，除了用传递函数描述滤波器特性外，还可用脉冲(冲激)响应来描述滤波器，因为在模拟滤波器中，脉冲响应与传递函数是Laplace变换对。此外还可以用阶跃响应（输入一个阶跃时系统的输出）来描述滤波器特性。下面我们介绍在MATLAB中如何得到模拟滤波器的脉冲响应和阶跃响应。

将滤波器的传递函数表示成分子和分母多项式系数的形式，如分子和分母多项式的系数为[b,a]，则在MATLAB中用H=[tf(b,a)]来表示此模拟滤波器，采用[[y,t]=]impulse(H)给出该系统的模拟脉冲响应。采用[[y,t]=]step(H)来得到该系统的阶跃响应。这两个函数与freqs一样，若没有输出则程序自动绘图模拟该滤波器的脉冲响应或阶跃响应。输出值y,t分别为该滤波器的脉冲响应或阶跃响应及其对应的时间序列。

另外，我们还可以运用一个输入信号模拟该滤波器的输出。若给定滤波器的输入值序列和对应的时间序列为u,t，则系统的输出可用y=lsim(H,u,t)来模拟, y为对应t的输出。若该函数没有输出变量则程序自动绘图显示。

【例5-13】设计一个5阶的Chebyshev I型带通滤波器，通带波纹3dB，下边界频率100Hz，上边界频率500Hz，绘制幅频响应图。给出该滤波器的脉冲响应、阶跃响应。假定输入sin(2*pi*30*t)+0.5*cos(2*pi*300*t)+2*sin(2*pi*800*t)的信号，求模拟滤波器的输出并给出模拟输入信号和模拟输出信号的Fourier振幅谱。

%Samp5_13

N=5;Rp=3;f1=100;f2=500;  %滤波器阶数、边界频率（Hz）

w1=2*pi*f1;w2=2*pi*f2;    %边界频率（rad/s）

[z,p,k]=cheb1ap(N,Rp);   %设计Chebyshev I型原型低通滤波器

[b,a]=zp2tf(z,p,k);  %转换为传递函数形式

Wo=sqrt(w1*w2);    %中心频率

Bw=w2-w1;    %频带宽度

[bt,at]=lp2bp(b,a,Wo,Bw);   %频率转换

[h,w]=freqs(bt,at);    %计算复数频率响应

figure(1)

subplot(2,2,1),semilogy(w/2/pi,abs(h));   %绘制幅频响应

xlabel('频率/Hz');grid on;title('幅频图');

subplot(2,2,2),plot(w/2/pi,angle(h)*180/pi);%绘制相频响应

xlabel('频率/Hz');

ylabel('相位/^o');grid on;title('相频图')

H=[tf(bt,at)];  %在MATLAB中表示此滤波器

[h1,t1]=impulse(H);  %输出绘出系统的脉冲响应图

subplot(2,2,3),plot(t1,h1);xlabel('时间/s');title('脉冲响应')

[h2,t2]=step(H);     %输绘出系统的阶跃响应图

subplot(2,2,4),plot(t2,h2);xlabel('时间/s');title('阶跃响应')

figure(2)

dt=1/2000;

t=0:dt:0.1;    %给出模拟滤波器输出的时间范围

u=sin(2*pi*30*t)+0.5*cos(2*pi*300*t)+2*sin(2*pi*800*t);  %模拟输入信号

subplot(2,2,1),plot(t,u)  %绘制模拟输入信号

xlabel('时间/s');title('输入信号')  

[ys,ts]=lsim(H,u,t);   %模拟系统的输入u时的输出

subplot(2,2,2),plot(ts,ys)  %绘制模拟输出信号

xlabel('时间/s'),title('输出信号'); 

subplot(2,2,3),plot((0:length(u)-1)/(length(u)*dt),abs(fft(u))*2/ length(u));%绘输入信号振幅谱

title('输入信号振幅谱'),xlabel('频率/Hz')

Subplot(2,2,4),

Y=fft(ys);

plot((0:length(Y)-1)/(length(Y)*dt),abs(Y)*2/length(Y)); %绘制输出信号振幅谱

title('输出信号振幅谱')

xlabel('频率/Hz')

图5-16 例5-13设计的Chebyshev I型滤波器的幅频响应、相频响应、脉冲响应和阶跃响应

 

图5-17 例5-13所设计滤波器模拟输入和输出信号的时间域图形和振幅谱

 

程序输出结果见图5-16和图5-17。图5-16给出了该程序得到的幅频图、相频图、脉冲响应和阶跃响应。幅频图清楚地给出了的通带范围和阻带范围。图5-17给出了模拟输入、输出信号的时间域波形及振幅谱。输入信号中含有30Hz,300Hz和500Hz的波，由滤波器的幅频特性（图5-16左上图）可知，30Hz和500Hz全在阻带内，只有300Hz的波可通过滤波器。从输入信号和模拟输出信号的时间域波形和振幅谱，可以看到达到了预期的效果。注意，由该滤波器的相频特性（图5-16的相频图）可知，该滤波器并不是线性相位。

5.7  模拟滤波器的设计

5.7.1  模拟滤波器设计步骤

用户对设计的滤波器提出设计要求，我们可以针对滤波器的设计要求设计滤波器。通常用户对模拟滤波器提出的要求有：（1）滤波器的性能指标，包括截止频率_{http://www.fzjpk.com/dizhen/word/chap5.files/image101.gif}(对于低通和高通)或上下边界频率、，通带波纹、阻带衰减等；（2）滤波器的类型，通常为Butterworth、Chebyshev I、 Chebyshev II、 Elliptic或Bessel滤波器。我们根据滤波器的类型通常按下列步骤设计滤波器。

（1）       给定模拟滤波器的性能指标，如截止频率_{http://www.fzjpk.com/dizhen/word/chap5.files/image101.gif}(对于低通和高通)或上下边界频率、；通带波纹、阻带衰减以及滤波器类型等（用户给定）。

（2）       确定滤波器阶数（根据本章第5节介绍的方法和函数）。

（3）       设计模拟低通原型滤波器。MATLAB信号处理工具箱的滤波器原型设计函数有butterap、cheb1ap、cheb2ap、ellipap、besselap。（本章第5.3节介绍的方法和函数）

（4）       按频率变换设计模拟滤波器（低通、高通、带通、带阻）。MATLAB信号处理工具箱的频率变换函数有lp2lp、lp2hp、lp2bp、lp2bs。（本章第5.4节介绍的方法和函数）。

下面我们用例子说明如何根据用户的要求设计滤波器。

【例5-14】设计带通Chebyshev I型模拟滤波器Ws1=0.2*pi（rad/s）; Ws2=0.8*pi（rad/s）; Rs=60dB; Wp1=0.35*pi（rad/s）; Wp2=0.65*pi（rad/s）; Rp=1dB。

%Samp5_14

wp=[0.35 0.65]*pi;   %给出滤波器的通带边界频率，用弧度表示

ws=[0.2 0.8]*pi;   %给出滤波器的阻带边界频率，，用弧度表示

Rp=1;Rs=60;   %给出滤波器的通带波纹和阻带衰减

[n,wn]=cheb1ord(wp,ws,Rp,Rs,'s');  %根据滤波器的要求求滤波器的最小阶数

[z,p,k]=cheb1ap(n,Rp);  %用求得的滤波器的最小阶数设计Chebyshev I型滤波器

[b,a]=zp2tf(z,p,k); %将零点极点增益形式转换为传递函数形式

Wo=sqrt(wn(1)*wn(2)); %计算中心点频率

Bw=wn(2)-wn(1);  %计算频带宽度

[bt,at]=lp2bp(b,a,Wo,Bw);  %模拟原型低通滤波器转换为带通滤波器

[h,w]=freqs(bt,at,128);  %计算滤波器的复数频率响应

plot(w/pi,20*log10(abs(h)));   %绘制滤波器的振幅响应

xlabel('角频率/pi');

ylabel('|H(jw)|/dB')

grid on;hold on;

plot([0.2 0.2],ylim, 'r'); plot([0.8 0.8],ylim, 'r');%绘出阻带界限以显示达到的阻带衰减

%plot([0.35 0.35],ylim); plot([0.65 0.65],ylim);%可绘出通带界限以显示达到的通带衰减

 

图5-18  例5-14设计的Chebyshev I型带通滤波器的幅频响应

 

程序的运行结果见图5-18。可以看到，设计的滤波器的通带范围确实为0.35*pi~0.65*pi,并且在阻带边界处下降分贝数大于60dB，满足了用户的设计要求。程序中ylim函数是提取现行坐标轴y轴的下限和上限。程序用它绘出了阻带边界。

5.7.2  模拟滤波器设计函数

上面滤波器的设计步骤比较麻烦，根据设计要求求解滤波器的最小阶数和边界频率之后需要设计模拟原型滤波器并进行频率转换。其实MATLAB将这一系列的过程组合成了更为方便的设计函数：butter,cheby1,cheby2,ellip,besself。 这些函数称为模拟滤波器完全设计函数。用户在求得滤波器的最小阶数和截止频率之后只需调用一次完全设计函数就可以自动完成所有设计过程，编程十分简单。这 些工具函数适用于模拟滤波器的设计，但同样也适用于数字滤波器。本节只讨论这些函数在模拟滤波器设计中的应用。但要注意，MATLAB是将上述一系列的步骤复合而已，并不是一种新的设计模拟滤波器的方法。

[b,a]=butter(n,wn[,'ftype'], 's')

[z,p,k]=butter(n,wn[,'ftype'], 's')

[b,a]=cheby1(n,Rp,wn[,'ftype'], 's')

[z,p,k]=cheby1(n,Rp,wn[,'ftype'], 's')

[b,a]=cheby2(n,Rs,wn[,'ftype'], 's')

[z,p,k]=cheby2(n,Rs,wn[,'ftype'], 's')

[b,a]=ellip(n,Rp,Rs,wn[,'ftype'], 's')

[z,p,k]=cheby2(n,Rp,Rs,wn[,'ftype'], 's')

[b,a]=besself(n,wn[,'ftype'], 's')

[z,p,k]=besself(n,wn[,'ftype'], 's')

在上面的调用方式中，n为滤波器的阶数，wn为滤波器的截止频率，单位rad/s（wn>0）；’s’为模拟滤波器，缺省时为数字滤波。

   ‘ftype’滤波器的类型可取为：

    ‘high’高通滤波器，截止频率为wn。

    ‘stop’带阻滤波器，截止频率为wn=[w1,w2] (w<1>w2)。

‘ftype’缺省时为低通或带通滤波器。

当设计带通滤波器时，截止频率也为wn=[w1,w2] (w1<>w2)。a,b分别为滤波器的传递函数分子和分母多项式系数向量；z,p,k分别为滤波器的零点、极点和增益。Rp,Rs分别为所设计滤波器的通带波纹和阻带衰减, 单位为dB。

滤波器的传递函数具有下面的形式：

                        (5-26)

若滤波器为带通或带阻型，则滤波器的阶数为2n，否则阶数为n。下面用一些例子说明滤波器的设计。

【例5-15】设计一个Butterworth模拟带通滤波器，设计指标为：通带频率：1000~2000Hz，两侧过渡带宽500Hz，通带波纹1dB，阻带衰减100dB。假设一个信号，其中f1=100Hz,f2=1500Hz,f3=2900Hz。信号的采样频率为10000Hz。试将原信号与通过该滤波器的模拟信号进行比较。

%Samp5_15

wp=[1000 2000]*2*pi;ws=[500 2500]*2*pi;Rp=1;Rs=100; %滤波器设计参数，对于给定Hz应乘以2_{http://www.fzjpk.com/dizhen/word/chap5.files/image103.gif}

[N,Wn]=buttord(wp,ws,Rp,Rs,'s'); %求得滤波器的最小阶数和截止频率 

w=linspace(1,3000,1000)*2*pi;   %设置绘制频率响应的频率点

[b,a]=butter(N,Wn,'s');    %设计模拟Butterworth滤波器

H=freqs(b,a,w);    %计算给定频率点的复数频率响应

magH=abs(H);phaH=unwrap(angle(H)); %计算幅频响应和相频响应

plot(w/(2*pi),20*log10(magH));  %以频率为横坐标绘制幅频响应

xlabel('频率/Hz');ylabel('振幅/dB');

title('Butterworth模拟带通滤波器');

hold on;plot([1000 1000],ylim,'r');plot([2000 2000],ylim,'r');%绘通阻带边界grid on

figure(2)

dt=1/10000;  %模拟信号采样间隔

f1=100;f2=1500;f3=2900;%输入信号的三个频率成分

t=0:dt:0.04;  %给定模拟时间段

x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t); %输入信号

H=[tf(b,a)];     %滤波器在MATLAB系统中的表示

[y,t1]=lsim(H,x,t);   %模拟输出

subplot(2,1,1),plot(t,x),title('输入信号') %绘出输入信号

subplot(2,1,2),plot(t1,y)   %绘制输出信号

title('输出信号'),xlabel('时间/s')

图5-19   例5-15设计的模拟滤波器幅频响应，两条竖线表示阻带边界

 

图5-20  例5-15中滤波器的输入信号和输出信号

 

程序的运行结果为图5-19和图5-20。可见，滤波器的幅频响应完全符合要求。在阻带边界确实下降到100dB。大家修改程序可以观看通带边界处下降的分贝数不超过1个dB。模拟的输入信号中含有三个频率成分，根据图5-19可知，2900Hz和100Hz的频率均在阻带内，只有1000Hz的频率在通带内，通过模拟系统后，将2900Hz的高频成分和100Hz的低频成分滤除，输出只含有1500Hz频率成分的振动。

【例5-16】 设计一个Chebyshev I型模拟带阻滤波器，设计指标为：阻通带频率：1000～2000Hz，两侧过渡带宽500Hz，通带波纹1dB，阻带衰减大于50dB。假设一个信号，其中f1=100Hz,f2=1500Hz,f3=2900Hz。信号的采样频率为20000Hz。试将原信号与通过该滤波器的模拟信号进行比较。

%Samp5_16

close all;figure(1)

ws=[1000 2000]*2*pi;wp=[500 2500]*2*pi;Rp=1;Rs=50;  %滤波器设计参数，乘以2_{http://www.fzjpk.com/dizhen/word/chap5.files/image103.gif}转换为弧度/s

[N,Wn]=cheb1ord(wp,ws,Rp,Rs,'s');   %求模拟滤波器的最小阶数和截止频率

w=linspace(1,3000,1000)*2*pi;   %设置计算复数频率响应的频率点

[b,a]=cheby1(N,Rp,Wn,'stop','s');  %设计带阻滤波器

H=freqs(b,a,w);    %计算给定频率点的复数频率响应

magH=abs(H);phaH=unwrap(angle(H));  %求滤波器的幅频响应和相频响应

plot(w/(2*pi),20*log10(magH));  %以频率为横轴绘制幅频响应

xlabel('频率/Hz');ylabel('振幅/dB');

title('Butterworth模拟带阻滤波器');

hold on;plot([1000 1000],ylim, 'r');plot([2000 2000],ylim, 'r');%绘阻带边界

grid on

figure(2)

dt=1/20000;

f1=100;f2=1500;f3=2900;   %输入信号含有的频率

t=0:dt:0.04;%时间序列

x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t); %输入信号

H=[tf(b,a)];   %滤波器在MATLAB系统中的表示

[y,t1]=lsim(H,x,t);   %模拟输出信号

subplot(2,1,1),plot(t,x),title('输入信号')  %绘输入信号

subplot(2,1,2),plot(t1,y)   %绘制输出信号

title('输出信号'),xlabel('时间/s')

图5-21 例5-16设计的带阻滤波器的幅频响应

图中，两条竖线表示阻带边界

 

图5-22 例5-16检测滤波器的输入输出信号

 

程序的输出结果见图5-21和图5-22。可见设计的带阻滤波器完全符合要求，阻带边界达到50dB的衰减。模拟的输出信号完全滤除了1500Hz的信号，模拟输出信号中只含用100Hz和2900Hz的信号。

【例5-17】设计一个Chebyshev II型模拟带阻滤波器，设计指标与上例相同。

%Samp5_17

ws=[1000 2000]*2*pi;wp=[500 2500]*2*pi;Rp=1;Rs=50;  %滤波器设计参数

[N,Wn]=cheb2ord(wp,ws,Rp,Rs,'s');   %求滤波器的最小阶数和截止频率

w=linspace(1,3000,1000)*2*pi;   %设置计算复数频率响应的频率点

[b,a]=cheby2(N,Rs,Wn,'stop','s');  %设计带阻滤波器

H=freqs(b,a,w);   %计算给定频率点的复数频率响应

magH=abs(H);phaH=unwrap(angle(H));  %计算幅频响应和相频响应

plot(w/(2*pi),20*log10(magH));  %以频率为横轴绘制幅频响应

xlabel('频率/Hz');ylabel('振幅/dB');

title('Chebyshev II 模拟带阻滤波器')

hold on; plot([1000 1000],ylim, 'r');plot([2000 2000],ylim, 'r');%阻带边界

grid on

图5-23  例5-17设计滤波器的幅频响应图

 

程序运行结果为图5-23。可以看出滤波器在阻带边界处达到了下降50dB的要求，此滤波器在阻带内有振荡，这正是Chebyshev II型滤波器的特点。

【例5-18】设计一个高通椭圆滤波器，设计性能指标为：通带边界频率wp=1500Hz,阻带边界频率ws=1000Hz,通带波纹Rp=1dB,阻带衰减Rs=30dB。假设一个信号，其中f1=400Hz,f2=1600Hz。信号的采样频率为10000Hz。试将原信号与通过该滤波器的模拟信号进行比较。

%Samp5_18

wp=1500*2*pi;ws=1000*2*pi;Rp=1;Rs=30;   %滤波器设计参数

[N,Wn]=ellipord(wp,ws,Rp,Rs,'s');   %求滤波器的最小阶数和截止频率

w=linspace(1,3000,1000)*2*pi;   %给出计算复数频率响应的频率点

[b,a]=ellip(N,Rp,Rs,Wn,'high','s');  %设计高通椭圆滤波器

H=freqs(b,a,w);   %计算给定频率的复数频率响应

magH=abs(H);phaH=unwrap(angle(H));  %求幅频响应和相频响应

figure(1)

plot(w/(2*pi),20*log10(magH));  %以频率为横轴绘幅频响应

xlabel('频率/Hz');ylabel('振幅/dB');

title('椭圆模拟高通滤波器');

hold on;plot([1000 1000],ylim, 'r'); %阻带边界

grid on

figure(2)   %给出另一个图形窗口

dt=1/10000;   %采样间隔

f1=400;f2=1600;%信号中所含频率成分

t=0:dt:0.04;%时间序列

x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t);  %输入信号

H=[tf(b,a)];   %滤波器在MATLAB系统中的表示

[y,t1]=lsim(H,x,t);  %模拟输出

subplot(2,1,1),plot(t,x),title('输入信号') %绘输入信号

subplot(2,1,2),plot(t1,y)  %绘制输出信号

title('输出信号'),xlabel('时间/s')

图5-24  例5-18设计滤波器的幅频响应图

 

图5-25例5-18滤波器的输入和输出信号

 

程序运行结果见图5-24和图5-25。可见所设计的高通滤波器完全符合要求，在阻带边界下降的分贝数达30dB。并且通带和阻带均有波纹，这正是椭圆滤波器的特点。将该滤波器模拟输入400Hz和1600Hz的信号后，输出结果完全滤除了400Hz低频成分，只含有1600Hz的信号。

【例5-19】设计一个五阶Bessel低通模拟滤波器，截止频率为200Hz1000rad/s,绘制滤波器的频率特性图。假设一个信号，其中f1=100Hz,f2=1000Hz。信号的采样频率为10000Hz。试将原信号与通过该滤波器的模拟信号进行比较。

%Samp5_19

N=5;  %滤波器阶数

Wn=200*2*pi1000;  %边界频率

[b,a]=besself(N,Wn); %设计Bessel滤波器

figure(1)

[H,w]=freqs(b,a,512);   %绘出滤波器的复数频率特性

magH=abs(H);phaH=unwrap(angle(H));  %求幅频响应和相频响应

subplot(2,1,1),plot(w/(2*pi),20*log10(magH));  %以频率为横轴绘幅频响应

grid on;xlabel('频率/Hz');ylabel('振幅/dB');

subplot(2,1,2),plot(w/2/pi,angle(H)*180/pi);%绘制相频响应

grid on;xlabel('频率/Hz');ylabel('相位/^o')

figure(2)    %新的绘图框

dt=1/10000;   %信号采样间隔

f1=100;f2=1000;  %输入两个频率成分

t=0:dt:0.1;

x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号

H=[tf(b,a)];   %滤波器在MATLAB系统中的表示

[y,t1]=lsim(H,x,t);  %模拟滤波器的输出

subplot(2,1,1),plot(t,x),title('输入信号')  %绘出输入信号

subplot(2,1,2),plot(t1,y)  %绘出输出信号

title('输出信号'),xlabel('时间/s')

http://www.fzjpk.com/dizhen/word/chap5.files/image205.gif

图5-26 例5-19设计Bessel滤波器的频率特性

 

图5-27 例5-19中Bessel滤波器的输入和输出

 

 

程序运行结果为图5-26和图5-27。虽然模拟Bessel滤波器的幅频特性不如其他滤波器好，但其相频特性为线性的。因此图5-27中的输出信号与输入的低频信号形状不变，只是有延迟。

第五章  习题

1.     模拟滤波器可分为哪几种类型？设计时要求的主要性能指标有哪些？

2.     设计程序研究本章第二节设计Butterworth、Chebyshev I,II、椭圆滤波器例题中滤波器的非线性相位特性，并与Bessel滤波器的线性相位特性进行比较。

3.     设计一个Butterworth模拟低通滤波器，要求满足下面指标：截止频率wc=1000Hz,通带波纹小于3dB,阻带始点频率ws=1500Hz，阻带衰减大于35dB。（1）确定滤波器的最小阶数。（2）滤波器的传递函数；（3）绘制滤波器的平方幅频和相频特性图。（4）当输入信号为x(t)=sin(2*pi*f1*t)+cos(2*pi*f2*t),f1=500Hz,f2=2000Hz时，绘制滤波器输出信号曲线。（5）绘制滤波器的单位脉冲响应曲线。

4.     确定三阶Chebyshev I型低通模拟滤波器的传递函数，性能指标为：通带波纹为1dB，截止频率为wc=2rad/s。

5.     低通模拟滤波器满足下列条件：边界频率为500Hz和750Hz,通带波纹为2dB,阻带衰减为30dB。试求满足此要求的Butterworth,Chebyshev I, ChebyshevII和椭圆Elliptic滤波器的最小阶数。

6.     设计一个高通滤波器，边界频率为100Hz和150Hz，通带波纹不大于2dB,阻带衰减不小于30dB。试设计满足此要求的Butterworth, Chebyshev I, ChebyshevII和Elliptic椭圆滤波器并绘制它们的平方幅频响应图。绘制此系统的脉冲响应和阶跃响应。

7.     设计满足下列条件的带通Butterworth滤波器：中心频率1kHz,带宽为200Hz,通带波纹在1dB以下，阻带衰减在40dB以上。确定滤波器的阶数、系统函数并绘制幅频特性、相频特性和零极点图。

8.     FSS-3地震计零极点为：零点：[0.0,0.0,1.127E3, -6.25E+01+j*6.307089E+02, -6.25E+01-j*6.307089E+02, -1.827E+02]。极点：[-4.442212-j*4.443554, -4.442212+j*4.443554,-1.3645E+02-j*1.137163E+02,-1.3645E+02+j*1.137163E+02,-5.51E+01-j*4.409807E+02, -5.51E+01+j*4.409807E+02]。归一化因子：-7.342E-2。增益：4.405E+2。将增益与归一化因子相乘作为计算滤波器的增益。试绘出该地震计的频率特性。

9.     FSS-3型地震仪系统传递函数为：零点：[1.355083E-4+j*1.043789E-2, 1.355083E-4-j*1.043789E-2, 1.684441E2+j*7.147554E1, 1.684441E2-j*7.147554E1, 6.572567+j*2.38467E2, 6.572567-j*2.38467E2, 8.87029E1+j*1.476857E+2, 8.87029E1-j*1.476857E+2, 3.297725E1+j*1.584205E2, 3.297725E1-j*1.584205E2, 1.018952E1+j*1.653736E2, 1.018952E1-j*1.653736E2, 2.924976+j*1.530131E2, 2.924976-j*1.530131E2, 9.947939+j*1.550421E2, 9.947939-j*1.550421E2]; 极点：[-4.442212+j*4.443554, -4.442212-j*4.443554, -8.49856E3, -7.043488E1+j*3.404871E2, -7.043488E1-j*3.404871E2, -4.773638+j*1.452003E2, -4.773638-j*1.452003E2, -6.337485+j*1.39173E2, -6.337485-j*1.39173E2, -7.902612+j*1.330283E2, -7.902612-j*1.330283E2, -2.083218E1+j*1.247569E2, -2.083218E1-j*1.247569E2, -2.792995E1+j*1.203601E2, -2.792995E1-j*1.203601E2, -3.405637E01+j*9.91384E1, -3.405637E01-j*9.91384E1, -4.302115E1+j*7.238421e1, -4.302115E1-j*7.238421e1, -4.992166e1+j*3.912072e1, -4.992166e1-j*3.912072e1, -5.206292e1];归一化因子：为1.264053E12; 增益为：4.405E2, 试绘出该地震系统的频率特性。

10.  SLJ-100型强震加速度计传递函数及其零极点：零点：[-3.601E02+j*2.532257488882203E3, -3.601E02+j*2.532257488882203E3]；极点：[-2.299999999999999E2+ j*4.516635916254485E02, -2.299999999999999E2-j*4.516635916254485E02,  -5.636E02]增益：k=21.955。试绘制其频率特性。