椭圆滤波器的设计精品资料.doc

燕山大学课程设计说明书第一章 摘 要本文简单、直观地介绍了椭圆低通滤波器的基本理论和设计思想，阐述了设计椭圆低通滤波器的具体步骤，利用MATLAB产生一个包含低频、中频、高频分量的连续信号，并实现对连续信号进行的采样。文中还对采样信号进行频谱分析，利用设计的椭圆滤波器对采样信号进行滤波处理，并对仿真结果进行分析和处理。详细介绍了在用MATLAB设计椭圆滤波器用到的工具和命令。 第二章 引 言 信号处理是科学研究和工程技术许多领域都需要进行的一个重要环节，传统上对信号的处理大都采用模拟系统实现。随着人们对信号处理要求的日益提高，以及模拟信号处理中一些不可克服的缺点，对信号的许多处理而采用数字的方法进行。近年来由于大规模集成电路和计算机技术的进步，信号的数字处理技术得到了飞速发展。数字信号处理系统无论在性能、可靠性、体积、耗电量、成本等诸多方面都比模拟信号处理系统优越的多，使得许多以往采用模拟信号处理的系统越来越多地被数字处理系统所代替,数字信号处理技术在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。在数字信号处理中,数字滤波器十分重要并已获得广泛应用,数字滤波器与模拟滤波器比较,具有精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配以及实现模拟滤波器无法实现的特殊滤波功能等优点。在各种滤波器中,椭圆滤波器具有其独特的优点。本次设计中所用到数学软件为MATLAB。MATLAB和Mathematica、Maple并称为三大数学软件，它是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。第三章 椭圆滤波器的基本理论3.1椭圆滤波器的概述常用数字滤波器的类型有巴特沃斯(Butterworth) ,切比雪夫(Chebyshev) 及椭圆型 滤波器，其中椭圆滤波器（Elliptic filter）又称考尔滤波器（Cauer filter），是一种性能优越的滤波器。从传递函数来看，巴特沃斯和切比雪夫滤波器的传输函数都是一个常数除以一个多项式, 为全极点网络, 仅在无限大阻带处衰减为无限大, 而椭圆函数滤波器在有限频率上既有零点又有极点。极零点在通带内产生等波纹, 阻带内的有限传输零点减少了过渡区, 可获得极为陡峭的衰减曲线。也就是说在阶数相同的条件下,椭圆滤波器相比于其他类型的滤波器，能获得更窄的过渡带宽和较小的阻带波动, 就这点而言, 椭圆滤波器是最优的。它陡峭的过渡带特性是用通带和阻带的起伏为代价来换取的，并且在通带和阻带的波动相同，这一点区别于在通带和阻带都平坦的巴特沃斯滤波器，以及通带平坦、阻带等波纹或是阻带平坦、通带等波纹的切比雪夫滤波器。 总结起来，椭圆滤波器具有以下特点：1、 椭圆低通滤波器是一种零、极点型滤波器，它在有限频率范围内存在传输零点和极点。2、 椭圆低通滤波器的通带和阻带都具有等波纹特性，因此通带，阻带逼近特性良好。 3、 对于同样的性能要求，它比前两种滤波器所需用的阶数都低，而且它的过渡带比较窄。但是椭圆滤波器传输函数是一种较复杂的逼近函数, 利用传统的设计方法进行电路网络综合要进行繁琐的计算, 还要根据计算结果进行查表, 整个设计, 调整都十分困难和繁琐。而用MATLAB设计椭圆滤波器可以大大简化设计过程。3.2椭圆滤波器设计的数学推导椭圆滤波器的振幅平方函数为 ： （1）其中是雅可比(Jacobi) 椭圆函数，雅可比椭圆函数是阶数N的有理函数，N=5时的特性曲线如图1所示。图1 N=5时雅可比椭圆函数的特性曲线 由图1 可见,在归一化通带内( - 1 1) ,() 在(0 ,1) 间振荡,而超过 后,() 在( , ) 间振荡。L 越大，也变大。这一特点使滤波器同时在通带和阻带具有任意衰减量。L 是一个表示波纹性质的参量。雅可比椭圆函数还具有以下性质： （2）阶数N等于通带和阻带内最大点和最小点的总和，为与通带衰减有关的参数。系统函数和阶数N是由系统下面的性能指标来确定的，主要有：截止频率,通带内最大衰减和阻带截止频率以及阻带内最小衰减。假定是频率归一化的基准频率，即 （3）定义频率的选择性因数k为 （4）则截止频率分别归一化为 （5）再次假定 （6） （7） （8） （9）则得到椭圆滤波器的阶数N为 （10）这时，令归一化的基准频率为，则得到归一化后的椭圆低通滤波器的系统函数为 （11）式中， 所以，实际的椭圆低通滤波器就可以由归一化的系统函数来得到 （12） 图2 为典型N 为奇数的椭圆滤波器的幅度特性, 当，和A 确定后,阶次N 即可确定,进而可以设计出椭圆滤波器。 图2 椭圆滤波器的幅度特性3.3关于归一化的讨论归一化是一种简化计算的方式，主要是为了数据处理方便提出来的，即将有量纲的表达式，经过变换，化为无量纲的表达式，成为纯量。比如，复数阻抗可以归一化书写：Z = R + jL = R(1 + jL/R) ，复数部分变成了纯数量了，没有量纲。归一化方法（Normalization Method）把数据映射到01范围之内处理，更加便捷快速，应该归到数字信号处理范畴之内。其具体作用是归纳统一样本的统计分布性。归一化在01之间是统计的概率分布，归一化在-1+1之间是统计的坐标分布。归一化有同一、统一和合一的意思。无论是为了建模还是为了计算，其基本度量单位要同一。但是归一化处理并不总是合适的，根据输出值的分布情况，标准化等其它统计变换方法有时可能更好。具体情况还需具体分析。第四章 设计过程4.1椭圆滤波器设计结构图椭圆滤波器设计结构图如图3所示：图3 结构框图4.2设计椭圆数字滤波器的步骤由于模拟滤波器的设计方法非常成熟, 许多典型系统有成熟的公式、图表可以查阅,便于设计;因此设计数字滤波器的主要方法是:首先设计一个合适的模拟滤波器, 然后将他“ 变换”成满足给定指标的数字滤波器。设计椭圆数字滤波器的步骤:(1) 确定数字滤波器性能指标、 ;(2) 将数字滤波器性能指标转换成相应的模拟滤波器性能指标;(3) 设计满足指标要求的模拟滤波器;(4) 通过变换将模拟滤波器转换成数字滤波器。4.3数字椭圆低通滤波器的MATLAB实现Matlab 是MathWorks 公司于1984 年正式推出的一套集数值计算、符号运算及图形处理等强大功能于一体的科学计算语言。作为强大的科学计算平台,他几乎能够满足所有的计算需求。他的应用范围涵盖了当今几乎所有的领域,如电子、半导体制造、医学研究、航空航天、汽车制造、分子模型、影视制作、建筑等行业。Matlab 具有以下优势和特点:友好的工作平台和编程环境,简单易用的程序语言,强大的科学计算及数据处理能力,出色的图形处理功能,应用广泛的模块集和工具箱,实用的程序接口和发布平台,模块化的设计和系统级的仿真。随着Matlab的不断完善, 尤其是Matlab 的信号处理工具箱( SignalProcessing Toolbox) 的推出,如今Matlab 已经成为数字信号处理DSP(Digital Signal Processing) 应用中分析和仿真设计的主要工具。4.3.1设计椭圆滤波器所用函数Matlab 的信号处理工具箱提供了设计椭圆滤波器的函数:ellipord 函数、ellip 函数和ellipap函数。1. ellipord 函数的功能是求滤波器的最小阶数和截止频率,其调用格式: N, = ellipord( , , , ) 可以得到数字椭圆型滤波器的最小阶数N和截止频率 ,并使滤波器在通带内(0 ,) 的波纹系数小于通带最大衰减 ,阻带内( ,1) 的波纹系数大于阻带最小衰减。其中是椭圆滤波器通带截止角频率，是椭圆滤波器阻带起始角频率。 根据本次任务书的设计要求，需要产生一个连续信号，包含低频5Hz，中频15Hz，高频30Hz的三个分量，并对其进行采样，采样频率为100Hz，采样点数为100。设计低通滤波器对信号进行滤波处理，滤除中频和高频信号。由于已知参数有限，对于设计中所用到的参数可取= 0.1，=40，通带截止频率Wp=5Hz，阻带截止频率 Ws=10Hz，归一化处理wp=2*Wp/Fs; ws=2*Ws/Fs 。根据程序：Wp=5;Ws=10;Fs=100;rp=0.1;rs=40;wp=2*Wp/Fs;ws=2*Ws/Fs;n,wn=ellipord(wp,ws,rp,rs)可得出：n =4wn =0.1000 2. ellip 函数的功能是设计滤波器,其调用格式:b,a = ellip ( N, , ,) , 利用ellipord 函数得到的最小阶数N和截止频率,可以设计低通滤波器。其中，b、a分别为椭圆滤波器传输函数的分子、分母多项式。3. ellipap函数的功能是直接返回椭圆滤波器的零点z、极点p和增益k，其调用格式：z , p , k=ellipap(N , )4.3.2频谱分析所用函数Matlab的信号处理工具箱提供了频谱分析函数:fft函数、filter函数和freqz函数。 1. fft函数的功能是对信号进行快速傅里叶变换，其调用格式： Y = fft(X) Y = fft(X,n)Y = fft(X,dim) Y = fft(X,n,dim)matlab的fft序号是从1到n，大多数采用从0到n-1，Y=fft（x）之后，这个Y是一个复数，它的模值应该除以（length(x)2），才能得到各个频率信号实际幅值。fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析，只需考察0Nyquist频率范为内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率01进行。另外，振幅的大小与所用采样点数有关，例如：采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz与15Hz振动幅值之比均为4：1，与真实振幅0.5：2是一致的。为了与真实振幅对应，需要将变换后结果乘以2除以N。 2. Freqz函数功能是用来求幅频响应，其调用格式：h,w=freqz(b,a,n)h,f=freqz(b,a,n,Fs)h=freqz(b,a,w)h=freqz(b,a,f,Fs)freqz(b,a,n)说明: freqz 用于计算数字滤波器H(Z)的频率响应函数H(ej)。h,w=freqz(b,a,n)可得到数字滤波器的n点复频响应值，这n个点均匀地分布在0,上,并将这n个频点的频率记录在w中，相应的频响值记录在h中。要求n为大于零的整数,最好为2的整数次幂,以便采用FFT计算，提高速度。缺省时n =512。h,f=freqz(b,a,n,)用于对在0,/2上等间隔采样n点，采样点频率及相应频响值分别记录在f 和h中。由用户指定（以HZ为单位）值。h=freqz(b,a,w)用于对在0,上进行采样，采样频率点由矢量w指定。h=freqz(b,a,f,Fs) 用于对在0,上采样，采样频率点由矢量f指定。freqz(b,a,n) 用于在当前图形窗口中绘制幅频和相频特性曲线。3. filter函数功能是利用IIR滤波器和FIR滤波器对数据进行滤波，其调用格式： y,zf=filter(b,a,x) y=filter(b,a,x,zi) y=filter(b,a,x)说明：filter采用数字滤波器对数据进行滤波，其实现采用移位直接型结构，因而适用于IIR和FIR滤波器。滤波器的系统函数为 即滤波器系数a=.,b=.,输入序列矢量为x。这里，标准形式为=1，如果输入矢量a时，1，则MATLAB将自动进行归一化系数的操作；如果=0，则给出出错信息。y=filter(b,a,x)利用给定系数矢量a和b对x中的数据进行滤波，结果放入y矢量中，y的长度取max(N,M)。y=filter(b,a,x,zi)可在zi中指定x的初始状态。y,zf=filter(b,a,x)除得到矢量y外，还得到x的最终状态矢量zf。4.3.3 其他命令函数在设计过程中出以上功能函数外，还用到了很多其他的函数。例如： clc 清除命令窗口中的内容 Clear 清除内存中的变量和函数 Plot 绘制线性图形 Abs 取模 等。4.4 椭圆低通滤波器的设计程序- 28 - %原始混合信号的产生及对其采样Fs=100;t=(1:100)/Fs;s1=sin(2*pi*t*5);s2=sin(2*pi*t*15);s3=sin(2*pi*t*30);s=s1+s2+s3;figure(1);subplot(2,1,1)plot(t,s)xlabel('时间/t')ylabel('幅值')title('原始输入模拟信号')subplot(2,1,2)stem(t,s)xlabel('时间/t')ylabel('幅值')title('采样后的输入信号')%椭圆低通滤波器的设计%求取最小阶数和截止频率Wp=5;Ws=10;Fs=100;rp=0.1;rs=40;wp=2*Wp/Fs;ws=2*Ws/Fs;n,wn=ellipord(wp,ws,rp,rs)b,a=ellip(4,0.1,40,5*2/Fs);H,w=freqz(b,a,512);figure(2);plot(w*Fs/(2*pi),abs(H);xlabel('频率/Hz');ylabel('幅度');title('椭圆低通滤波器的幅度响应')grid;%对滤波后的信号进行分析和变换sf=filter(b,a,s);figure(3);plot(t,sf);xlabel('时间/t');ylabel('幅值');title('滤波后的信号')axis(0 1 -1 1);S=fft(s,512);SF=fft(sf,512);w=(0:255)/256*(Fs/2);figure(4);subplot(2,1,1)plot(w,abs(S(1:256);xlabel('频率/Hz');ylabel('幅度');title('滤波前的傅立叶变换图')grid;subplot(2,1,2)plot(w,abs(SF(1:256);xlabel('频率/Hz');ylabel('幅度');title('滤波后的傅立叶变换图')grid;第五章 仿真图5.1原始信号及其采样仿真图，如图3所示图3 原始输入信号及其采样图5.2信号通过椭圆低通滤波器的仿真图，如图4所示图4 信号通过椭圆低通滤波器的仿真图5.3 椭圆低通滤波器的幅度特性，如图5所示图5 椭圆低通滤波器的幅度特性5.4 对信号进行傅里叶变换的仿真，如图6所示图6 信号的傅里叶变换第六章 结论 在最常用的巴特沃斯、切比雪夫、椭圆函数3 种滤波器中,实现同样性能指标所需的阶数及阻带衰减，椭圆滤波器给出的设计阶数比前两种低,而且频率特性较好,过渡带较窄,但是椭圆滤波器在通带上的非线性相位响应最明显。数字滤波器的应用十分广泛，利用Matlab，很容易地设计地设计出数字椭圆滤波器，数字椭圆滤波器可以用较少的阶数获得很高的选择特性，在相位要求不敏感的场合，如语音信号处理等，适合用数字椭圆滤波器。在设计过程中可以对比滤波器的特性，随时更改参数，已达到滤波器设计的最优。心得体会整个设计过程，让我感触最深的就是功能的强大性与掌握编程各种函数和语句的重要性。当然首先要了解所要编程运行的对象的原理。在课程设计的过程中，我深深的感受到我所学的东西太少了，需要学习的东西太多了。在一周的课程设计时间里，我每天都过的很充实，查资料、读程序，反复琢磨，学习的过程是艰辛的，但是同时也是快乐的。通过实实在在的课程设计，发现自己在课堂上所学的知识对于解决实际问题来说，是远远不够的，要想掌握技术，还需要加倍的努力。不过，在这短短一周的时间里，还是让我对数字信号处理有了更深了解，也学到了新的知识让我掌握了visio、MATLAB等软件的简单使用，明白了软件仿真对设计的重要性。再次，这次课程设计让我充分认识到团队合作的重要性，只有分工协作才能保证整个项目的有条不紊。在整个课程设计过程中我再次认识到坚持、耐心、细心等品质的重要，这对今后的学习和工作是有很大帮助的。总之，这次课程设计让我获益良多。参考文献【1】matlab信号处理详解 陈亚勇等编著 人民邮电出版社，2001年【2】电子滤波器设计 宁彦卿等译 科学出版社. 2008年【3】 数字信号处理原理及应用 谢平 王娜 李洪彬编著 机械工业出版社 2008年【4】 数字椭圆滤波器的Matlab设计与实现 王靖 李永全 现代电子技术2007年第6 期总第245附录资料：matlab画二次曲面一、螺旋线1.静态螺旋线a=0:0.1:20*pi;h=plot3(a.*cos(a),a.*sin(a),2.*a,'b','linewidth',2);axis(-50,50,-50,50,0,150);grid onset(h,'erasemode','none','markersize',22);xlabel('x轴');ylabel('y轴');zlabel('z轴');title('静态螺旋线'); 2.动态螺旋线t=0:0.1:10*pi;i=1;h=plot3(sin(t(i),cos(t(i),t(i),'*','erasemode','none');grid onaxis(-2 2 -2 2 0 35)for i=2:length(t) set(h,'xdata',sin(t(i),'ydata',cos(t(i),'zdata',t(i); drawnow pause(0.01)endtitle('动态螺旋线');(图略) 3.圆柱螺旋线t=0:0.1:10*pi;x=r.*cos(t);y=r.*sin(t);z=t;plot3(x,y,z,'h','linewidth',2);grid onaxis('square')xlabel('x轴');ylabel('y轴');zlabel('z轴');title('圆柱螺旋线') 二、旋转抛物面b=0:0.2:2*pi;X,Y=meshgrid(-6:0.1:6);Z=(X.2+Y.2)./4;meshc(X,Y,Z);axis('square')xlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('旋转抛物面')或直接用：ezsurfc('(X.2+Y.2)./4') 三、椭圆柱面load clownezsurf('(2*cos(u)','4*sin(u)','v',0,2*pi,0,2*pi)view(-105,40) %视角处理shading interp %灯光处理colormap(map) %颜色处理grid on %添加网格线axis equal %使x,y轴比例一致xlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('椭圆柱面') %添加标题四、椭圆抛物面b=0:0.2:2*pi;X,Y=meshgrid(-6:0.1:6);Z=X.2./9+Y.2./4;meshc(X,Y,Z);axis('square')xlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('椭圆抛物面')或直接用：ezsurfc('X.2./9+Y.2./4')五、'双叶双曲面ezsurf('8*tan(u)*cos(v)','8.*tan(u)*sin(v)','2.*sec(u)',-pi./2,3*pi./2,0,2*pi)axis equalgrid onaxis squarexlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('双叶双曲面')六、双曲柱面load clownezsurf('2*sec(u)','2*tan(u)','v',-pi/2,pi/2,-3*pi,3*pi)hold on %在原来的图上继续作图ezsurf('2*sec(u)','2*tan(u)','v',pi/2,3*pi/2,-3*pi,3*pi)colormap(map)shading interpview(-15,30)axis equalgrid onaxis equalxlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('双曲柱面')七、双曲抛物面（马鞍面）X,Y=meshgrid(-7:0.1:7);Z=X.2./8-Y.2./6;meshc(X,Y,Z);view(85,20)axis('square')xlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('双曲抛物面')或直接用：ezsurfc('X.2./8-Y.2./6') 八、抛物柱面X,Y=meshgrid(-7:0.1:7);Z=Y.2./8;h=mesh(Z);rotate(h,1 0 1,180) %旋转处理%axis(-8,8,-8,8,-2,6);axis('square')xlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('抛物柱面')或直接用：ezsurfc('Y.2./8') 九、环面ezmesh('(5+2*cos(u)*cos(v)','(5+2*cos(u)*sin(v)','2*sin(u)',0,2*pi,0,2*pi)axis equalgrid onxlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('环面')十、椭球ezsurfc('(5*cos(u)*sin(v)','(3*sin(u)*sin(v)','4*cos(v)',0,2*pi,0,2*pi)axis equalgrid onxlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('椭球')十一、单叶双曲面ezsurf('4*sec(u)*cos(v)','2.*sec(u)*sin(v)','3.*tan(u)',-pi./2,pi./2,0,2*pi)axis equalgrid onxlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('单叶双曲面')十二、旋转单叶双曲面load clownezsurf('8*sec(u)*cos(v)','8.*sec(u)*sin(v)','2.*tan(u)',-pi./2,pi./2,0,2*pi)colormap(map)view(-175,30)%alpha(.2) %透明处理axis equalgrid onaxis squarexlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('旋转单叶双曲面')十三、圆柱面subplot(1,2,1)ezsurf('(2*cos(u)','2*sin(u)','v',0,2*pi,0,2*pi)grid onshading interpaxis equalxlabel('x轴');ylabel('y轴');zlabel('z轴');title('圆柱面')subplot(1,2,2)cylinder(30)shading interpaxis squaretitle('调用cylinder函数所得圆柱面')十四、二次锥面clc,clear;P=1,0,0; 0,cos(45*pi/180),sin(45*pi/180); 0,-sin(45*pi/180),cos(45*pi/180);for k2 = 1:31 for k1 = 1:31 x(k1,k2) = (k2-1)*cos ( (k1-1)*12*pi/180); y(k1,k2) = (k2-1)*sin ( (k1-1)*12*pi/180); z(k1,k2) = sqrt(x(k1,k2)2+y(k1,k2)2); Vxyz = P*x(k1,k2),y(k1,k2),z(k1,k2)' x1(k1,k2)=Vxyz(1); y1(k1,k2)=Vxyz(2); z1(k1,k2)=Vxyz(3); endendsurf(x,y,z)hold on;surf(x1,y1,z1);shading flat;