【摘要】本设计利用Matlab语言对滤波器设计建立一个软件仿真平台,系统由巴特沃斯滤波器、切比雪夫Ⅰ型滤波器、切比雪夫Ⅱ型滤波器、椭圆滤波器等组成,设计结果由时域冲击响应图形及频谱图形表示,可任意修改各系统参数以分析研究滤波器设计结果,形象直观。本设计开发目的是给滤波器学习分析与研究提供一个软件平台,可作为《信号与系统》、《自动控制原理》、《通信原理》、《数字信号处理》等课程教学辅助。
【关键词】滤波器 MatLab语言 仿真 多媒体教学
【中图分类号】TP31 【文献标识码】A 【文章编号】1009-9646(2008)09(a)-0149-04
滤波器的功能就是允许某一部分频率的信号顺利通过,而另外一部分频率的信号则受到较大的抑制,它实质上是一个选频电路。滤波器中,把信号能够通过的频率范围,称为通频带或通带;反之,信号受到很大衰减或完全被抑制的频率范围称为阻带;通带和阻带之间的分界频率称为截止频率;理想滤波器在通带内的电压增益为常数,在阻带内的电压增益为零;实际滤波器的通带和阻带之间存在一定频率范围称过渡带,在通带内的电压增益不是常数而是一定的波动范围称通带波动,在阻带内的电压增益不是零而是一定的波动范围称阻带波动,滤波器指标如图1所示。由于滤波器设计数学模型复杂,不易理解学习,本文采用Matlab程序,根据给定滤波器参数进行设计,得出设计结果并绘出频域和时域特性图,这种图形化、可视化,形象直观便于理解和学习研究,可提高教学效率。本设计已应用于教学实践中,得到学生普遍好评。
1 系统设计与仿真
本设计开发目的是给滤波器设计的学习分析与研究提供一个软件平台,也是为“信号与系统”、“数字信号处理”、“自动控制原理”、“通信原理”等课程教学提供多媒体教学辅助平台,方便学生对系统的各种特征进行观察,分析比较,提高教学效率。因此,本设计采用菜单形式,把各滤波器分别独立模块,方便模块间进行比较观察分析。
1.1 巴特沃斯低通/高通滤波器
巴特沃斯模拟低通滤波器原型的平方幅度响应函数为:,式中ωc为低通滤波器的截止频率,N为滤波器的阶数。Butterworth低通滤波器特点是通带内具有最大平坦的幅频特性,且随频率增大平滑单调下降;阶数N越高,特性越接近矩形,过渡带越窄;传递函数无零点,极点等距离分布在以|s|=ωc为半径的圆周上。
程序算法:输入滤波器设计参数即通带边界频率fp、阻带边界频率fs、通带最大衰减Rp、阻带最小衰减As;边界频率转化为角频率op=2*pi*fp、os=2*pi*fs;计算最小阶数和截止频率[N,oc]=buttord(op,os,Rp,Rs,"s");计算低通滤波器系统函数分子分母系数,低通[b,a]=butter(N,oc,"s"),高通[b,a]=butter(N,oc,"high","s");计算归一化模拟滤波器零点、极点和增益[z,p,k]=tf2zp(b,a);确定显示频谱特性的角频率范围w;计算频率响应H=freqs(b,a,w);计算单位冲击响应sys=tf(b,a),[h,t]=impulse(sys);计算单位阶跃响应[h1,t1]=step(sys);用plot函数绘出增益系数频谱特性图、分贝系数频谱特性图、单位冲击响应特性图、单位阶跃响应特性图。程序运行仿真结果如图2所示
仿真结果分析:分别绘出增益系数频谱特性图和分贝系数频谱特性图,目的是让两种频谱特性图放在一起进行对比,两种指标的表现能力及优缺点清晰可见。由巴特沃斯低通滤波器的幅频特性曲线图可看出在通带内幅频特性比较平坦,且随着频率增大单调下降;传递函数无零点,极点分布在圆周上;由图2所示给出设计参数,设计结果频谱曲线可看出通带最大衰减Rp=0.63,通带指标有富余;阻带最小衰减,As=30,刚好满足指标要求。任意改变设计参数,仿真结果与滤波器设计理论相符合。
1.2 巴特沃斯带通/带阻滤波器
程序算法:输入设计参数fp1、fp2、fs1、fs2、Rp、As;op1=2*pi*fp1;op2=2*pi*fp2;os1=2*pi*fs1;os2=2*pi*fs2;op=[op1,op2],os=[os1,os2];[N,oc]=buttord(op,os,Rp,As,"s");if 带通 [b,a]=butter(N,oc,"s") else [b,a]=butter(N,oc,"stop","s") end;[z,p,k]=tf2zp(b,a);w=linspace(1,fs2+1000,1000)*2*pi;H=freqs(b,a,w);wx=[op(1),os(1)];Hx=freqs(b,a,wx);dbHx=20*log10(abs(Hx));sys=tf(b,a);[h,t]=impulse(sys);[h1,t1]=step(sys);axes(handles.axes1);plot(w/(2*3.14),abs(H),"-r");axes(handles.axes2);plot(w/(2*3.14),20*log10(abs(H)),"-r");axes(handles.axes3);plot(t,real(h));axes(handles.axes4);plot(t1,real(h1));程序运行结果如图3所示。
运行结果分析:由图上可看出通带位于阻带之间,通带的幅频特性很平稳不波动,阻带明显单调下降。由图3频谱曲线图可知,设计条件为fp1=1000Hz,fp2=2000Hz,fs1=500Hz,fs2=2500Hz,Rp=1dB,As=100dB,则设计结果Rp=0.97<1,As=244>100都满足指标要求。改变设计参数,仿真结果与滤波器设计理论相符合。
1.3 切比雪夫Ⅰ型低通/高通滤波器
ChebyshevⅠ的平方幅度响应函数为
式中为小于1的正数,表示通带内幅值波纹情况;为截止频率,N为切比雪夫多项式的阶数。
程序算法:输入设计参数fp、fs、Rp、As;op=2*pi*fp1;os=2*pi*fs;[N,oc]=cheb1ord(op,os,Rp,As,"s");if 低通 [b,a]=cheby1(N,Rp,oc,"s") else [b,a]=cheby1(N,Rp,oc,"high","s") end;[z,p,k]=tf2zp(b,a);w0=0:0.1:3.14;w=w0*oc;H=freqs(b,a,w);wx=[op,os];Hx=freqs(b,a,wx);dbHx=20*log10(abs(Hx));sys=tf(b,a);[h,t]=impulse(sys) ;[h1,t1]=step(sys);axes(handles.axes1);plot(w/(2*3.14),abs(H),"-r");axes(handles.axes2);plot(w/(2*3.14),20*log10(abs(H)),"-r");axes(handles.axes3);plot(t,real(h));axes(handles.axes4);plot(t1,real(h1));程序运行结果如图4所示。
仿真结果分析:由仿真结果可见切比雪夫Ⅰ低通滤波器通带内具有等波纹起伏特性,阻带内单调下降且衰减更大,随着N的增大接近矩形。传递函数无零点,极点分布于椭圆上。由衰减曲线图可知通带最大衰减和阻带最小衰减均满足指标要求,分贝幅频特性能更清晰显示滤波器频谱特性。改变设计参数,仿真结果与滤波器设计理论相符合。
1.4 切比雪夫Ⅰ型带通/带阻滤波器
程序算法:输入设计参数fp1、fp2、fs1、fs2、Rp、As;op1=2*pi*fp1;op2=2*pi*fp2;os1=2*pi*fs1;os2=2*pi*fs2;op=[op1,op2],os=[os1,os2];[N,oc]=cheb1ord(os,op,Rp,As,"s");if 带通 [b,a]=cheby1(N,Rp,oc,"s") else [b,a]=cheby1(N,Rp,oc,"stop","s")end;[z,p,k]=tf2zp(b,a);w=linspace(1,fs2+1000,1000)*2*pi;H=freqs(b,a,w);wx=[op(1),os(1)];Hx=freqs(b,a,wx);dbHx=20*log10(abs(Hx));sys=tf(b,a);[h,t]=impulse(sys);[h1,t1]=step(sys);axes(handles.axes1);plot(w/(2*3.14),abs(H),"-r");axes(handles.axes2);plot(w/(2*3.14),20*log10(abs(H)),"-r");axes(handles.axes3);plot(t,real(h));axes(handles.axes4);plot(t1,real(h1));程序运行结果如图5所示。
仿真结果分析:由仿真结果可见切比雪夫Ⅰ带通滤波器通带内具有等波纹起伏特性,阻带内单调下降且衰减更大,传递函数无零点,也可推出极点分布于椭圆上。由频谱曲线图可知通带最大衰减和阻带最小衰减均满足指标要求。改变参数,仿真结果与滤波器设计理论相符合。
1.5 切比雪夫Ⅱ型低通/高通滤波器
ChebyshevⅡ型的平方幅度响应函数为
式中为小于1的正数,表示阻带内幅值波纹情况,为截止频率,N为切比雪夫多项式的阶数。
程序算法:输入设计参数fp、fs、Rp、As;op=2*pi*fp1;os=2*pi*fs;[N,oc]=cheb2ord(op,os,Rp,As,"s");if 低通 [b,a]=cheby2(N,As,oc,"s") else [b,a]=cheby2(N,As,oc,"high","s") end;[z,p,k]=tf2zp(b,a);w0=0:0.1:3.14;w=w0*oc;H=freqs(b,a,w);wx=[op,os];Hx=freqs(b,a,wx);dbHx=20*log10(abs(Hx));sys=tf(b,a);[h,t]=impulse(sys) ;[h1,t1]=step(sys);axes(handles.axes1);plot(w/(2*3.14),abs(H),"-r");axes(handles.axes2);plot(w/(2*3.14),20*log10(abs(H)),"-r");axes(handles.axes3);plot(t,real(h));axes(handles.axes4);plot(t1,real(h1));程序运行结果如图6所示。
仿真结果分析:由仿真结果可见切比雪夫Ⅱ滤波器阻带内具有等波纹起伏特性,通带内单调、平滑,随着N的增大接近矩形。传递函数有零点、极点。由频谱特性曲线图可知通带最大衰减和阻带最小衰减均满足指标要求,分贝幅频特性能更清晰显示滤波器频谱特性。改变参数,仿真结果与滤波器设计理论相符合。
1.6 切比雪夫Ⅱ型带通/带阻滤波器
程序算法:输入设计参数fp1、fp2、fs1、fs2、Rp、As;op1=2*pi*fp1;op2=2*pi*fp2;os1=2*pi*fs1;os2=2*pi*fs2;op=[op1,op2],os=[os1,os2];[N,oc]=cheb2ord(os,op,Rp,Rs,"s");if 带通 [b,a]=cheby2(N,As,oc,"s") else [b,a]=cheby2(N,As,oc,"stop","s")end;[z,p,k]=tf2zp(b,a);w=linspace(1,fs2+1000,1000)*2*pi;H=freqs(b,a,w);wx=[op(1),os(1)];Hx=freqs(b,a,wx);dbHx=20*log10(abs(Hx));sys=tf(b,a);[h,t]=impulse(sys);[h1,t1]=step(sys);axes(handles.axes1);plot(w/(2*3.14),abs(H),"-r");axes(handles.axes2);plot(w/(2*3.14),20*log10(abs(H)),"-r");axes(handles.axes3);plot(t,real(h));axes(handles.axes4);plot(t1,real(h1));程序运行结果如图7所示。
仿真结果分析:由仿真结果可见切比雪夫Ⅱ带通滤波器通带内具有单调下降特性,阻带内具有等波纹起伏特性,传递函数有零点和极点。由衰减曲线图可知通带最大衰减和阻带最小衰减均满足指标要求。改变参数,仿真结果与滤波器设计理论相符合。
1.7 椭圆滤波器低通/高通滤波器
椭圆滤波器的平方幅度响应函数为:
式中为小于1的正数,表示波纹情况,为截止频率,N为多项式的阶数。
程序算法:输入设计参数fp、fs、Rp、As;op=2*pi*fp1;os=2*pi*fs;[N,oc]=ellipord(op,os,Rp,As,"s");if 低通 [b,a]=ellip(N,Rp,As,oc,"s") else [b,a]= ellip(N,Rp,As,oc,"high","s") end;[z,p,k]=tf2zp(b,a);w0=0:0.1:3.14;w=w0*oc;H=freqs(b,a,w);wx=[op,os];Hx=freqs(b,a,wx);dbHx=20*log10(abs(Hx));sys=tf(b,a);[h,t]=impulse(sys) ;[h1,t1]=step(sys);axes(handles.axes1);plot(w/(2*3.14),abs(H),"-r");axes(handles.axes2);plot(w/(2*3.14),20*log10(abs(H)),"-r");axes(handles.axes3);plot(t,real(h));axes(handles.axes4);plot(t1,real(h1));程序运行结果如图8所示。
仿真结果分析:由仿真结果可见椭圆滤波器的特性是通带和阻带内都具有等波纹起伏特性,传递函数有零点,也有极点。由频谱曲线图可知通带最大衰减和阻带最小衰减均满足指标要求,分贝幅频特性能更清晰显示滤波器频谱特性。改变设计参数,仿真结果与滤波器设计理论相符合。
1.8 椭圆滤波器带通/带阻滤波器
程序算法:输入设计参数fp1、fp2、fs1、fs2、Rp、As;op1=2*pi*fp1;op2=2*pi*fp2;os1=2*pi*fs1;os2=2*pi*fs2;op=[op1,op2],os=[os1,os2];[N,oc]=ellipord(op,os,Rp,As,"s");if 带通 [b,a]=ellip(N,Rp,As,oc,"s") else [b,a]=ellip(N,Rp,As,oc,"stop","s") end;[z,p,k]=tf2zp(b,a);w=linspace(1,fs2+1000,1000)*2*pi;H=freqs(b,a,w);wx=[op(1),os(1)];Hx=freqs(b,a,wx);dbHx=20*log10(abs(Hx));sys=tf(b,a);[h,t]=impulse(sys);[h1,t1]=step(sys);axes(handles.axes1);plot(w/(2*3.14),abs(H),"-r");axes(handles.axes2);plot(w/(2*3.14),20*log10(abs(H)),"-r");axes(handles.axes3);plot(t,real(h));axes(handles.axes4);plot(t1,real(h1));程序运行结果如图9所示。
仿真结果分析:由仿真结果可看出椭圆滤波器通带内具有等幅波动特性,阻带内也具有等波纹起伏特性,传递函数有零点和极点。由衰减曲线图可知通带最大衰减和阻带最小衰减均满足指标要求。更改参数,仿真结果与滤波器设计理论相符合。
2 软、硬件运行环境及调试运行
硬件要求:PⅡ以上CPU;256M以上内存;20G以上硬盘空间。
软件要求:Windows2000及以上版本操作系统;采用可视化面向对象程序设计的Matlab7.0为系统开发平台。
本系统经过调试修改及实践验证,已能全面满足设计需要,大大提高了对系统学习研究观察的工作效率,显示较高的实用价值,可在“信号与系统”、“通信原理”、“自动控制原理”和“数字信号处理”等课程教学辅助中推广使用。
3 结论
本系统具有如下技术特点:①各滤波器设计的仿真结果满足设计指标要求,频谱特性符合滤波器理论。②界面友好,使用方便,操作简单,运行效率高。③可方便改变滤波器参数,快速绘出与之对应的时域波形,频谱特性波形,方便对滤波器的观察与研究。④采用菜单形式,模块之间的联系清晰,调用方便。各系统采用独立模块,方便不同系统间比较观察研究。⑤对滤波器设计方法进行分类归纳,并进行简单明确的说明,方便“信号与系统”、“通信原理”、“自动控制原理”的教学辅助,也方便利用本软件进行滤波器设计的自学与研究。⑥滤波器系统特性用时域图形,频谱图形表示,形象直观,对学习研究滤波器各特征之间相互关系很有帮助。⑦作者利用此平台作为多媒体教学辅助,把抽象的滤波器设计数学模型用形象的图形来表征,大大提高教学效率,获得学生好评。
参考文献
[1] 郑君里.信号与系统[M].北京:高等教育出版社,2000.
[2] 刘卫国.Matlab程序设计[M].北京:高等教育出版社,2006.
[3] 陈金西,等.离散信号可视化平台的设计与实现[J].中国医疗设备,2008(3):10-14
[4] 陈怀琛.数字信号处理教程—MATLAB释义与实现[M].北京:电子工业出版社,2004.
[5] 陈金西,等.信号可视化平台的设计开发[J].高校教育研究,2008(7):206-208.
注：“本文中所涉及到的图表、注解、公式等内容请以PDF格式阅读原文。”