基于某MATLAB地椭圆滤波器设计

东北大学秦皇岛分校电子信息系综合课程设计基于matlab基于matlab的数字椭滤波器设计专业名称通信工程班级学号4090923学生姓名指导教师设计时间2011.12.19~2012.1.6实用标准文案实用标准文案精彩文档精彩文档课程设计任务书专业：通信工程学号：4090923学生姓名（签名）：设计题目：基于matlab的数字滤波器设计一、设计实验条件笔记本电脑一台；matlab7.0软件。二、设计任务及要求产生一个连续信号，包含低频率，中频，高频分量，对其进行采样，进行频谱分析；设计低通，带通，高通滤波器对信号进行滤波处理，观察滤波后信号的频谱；熟悉基于MATLAB设计椭圆滤波器过程中常用到的函数和命令；了解椭圆型滤波器的基本理论和设计思想；三、设计报告的内容设计题目与设计任务（设计任务书）前言（绪论）（设计的目的、意义等）设计主体（各部分设计内容、分析、结论等）结束语（设计的收获、体会等）参考文献四、设计时间与安排1、设计时间：3周2、设计时间安排：熟悉实验设备、收集资料：6天设计图纸、实验、计算、程序编写调试：6天编写课程设计报告：2天答辩：1天2.前言随着信息时代和数字时代的到来，数字信号信号处理已经成为当前一门极其重要的学科和技术领域。在数字信号处理中起着重要作用并已获得广泛应用的是数字滤波器(DF,DigitalFilter)。数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。MATLAB是英文MATrixLABoratory(矩阵实验室)的缩写。它是美国的MathWorks公司推出的一套用于科学计算和图形处理可视化、高性能语言与软件环境。它的信号处理工具箱包含了各种经典的和现代的数字信号除了技术,是一款非常优秀的算法研究与辅助设计的工具。在现带通信系统中,由于信号中经常混有各种复杂成分,所以很多信号的处理和分析都是基于滤波器而进行的。而滤波器的种类很多,从功能上能将滤波器分为低、带、高、带阻类型。从实现方法上可分为FIR、IIR等。从设计方法上可分为Chebyshev(切比雪夫)，Butterworth(巴特沃兹)。而本次课程设计选题为对采样信号进行频谱分析,利用设计的滤波器对采样信号进行滤波处理,并对仿真结果进行分析和处理。目的在于了解滤波器的基本理论和设计思想,掌握在基于MATLAB设计滤波器过程中常用到的工具和命令。学会设计滤波器的具体步骤，利用MATLAB滤波滤波器设计函数直接实现滤波器。3.设计主体3.1信号发生模块设计3.1.1.DSB信号在条幅信号中，若果将载波抑制，则得到抑制载波双边带信号，简称双边带信号(DSB)。其时域表达式为：S(t)=m(t)-coswtDSBc频域表达式为：S(w)=-[M(w—w)+M(w+w)]DSB2cc从时域表达式中可以看出，DSB信号是调制信号与载波信号直接相乘得到的。而从频域表达式中可以看出,DSB信号的带宽是调制信号带宽的两倍且包含上下两个相互对称的频带，DSB信号的频带宽度是2w。m

在本次设计中，我们使用三种不同频率的单频信号作为调制信号信号，分别为100Hz、50Hz和25Hz。采用1KHz、500Hz和250Hz作为高频、中频和低频的载波信号。在MATLAB中将三路信号的表达式相加，很容易得到所需的混合信号。3.1.2.FFT算法FFT（快速傅里叶变换）使得DFT（离散傅里叶变换）的计算大大简化，在实际中有着广泛的应用oFFT在算法的思想是不断地把长序列的DFT分解成短序列的DFT，并利用旋转因子的周期性和对称性来合并DFT中的某些项，减少DFT的运算次数，大致上可以分为两大类，即按时间抽选法（DIT）和按频率抽选法（DIF）。下面简要介绍按时间抽选法。对于一个长度为N的离散信号序列x[n]，其DFT变换为：X(k)=£TX[n]Wnk，其中Wnk=e-jNnko

NNn=0对于任意0<m<N-1，X(m)=£1x[n]Wnm=x[0]W0m+X[1]W1m+…+x[N一1]W(N-1)moNNNNn=0若序列x[n]的长度为N=2M(如果不是，则只需要添加多个0值进行填充),N^-1A(N^-1A(k)+B(k)=乙a[n]W2n+Nn=0^-1^-1^-1b[n]W[n]W■厶b[n]WknN-n=0n=0N2n=0首先将x（n）按奇偶分成两个子序列，即：1a[n]=x[2n]N[b[n]=x[2n+1]n=0,1'2=一1令A（k）,B（k）分别表示a[n],b[n]的N2各点的傅里叶变换，即:经过分析可得到NX(k)=A(k)+WkB(k)，k=0,1,...,——1N2NNX(一+k)=A(k)—WkB(k)，k=0,1,...,——12N2由于N=2m—1仍然是偶数，故又可将a[n]和b[n]分解成奇偶部分，采用上面同2样的算法，重复进行下去，直到分组信号只剩下一个值时，便得到了信号的傅里叶变换。在MATLAB中，使用函数fft可以方便地计算信号的傅里叶变换。格式：y二fft(x)y=fft(x,n)说明：y二fft(x)当x的长度为2的幕次方时，fft函数采用基2的FFT算法，否则采用稍慢的混合基算法。y=fft(x,n)采用n点FFT，当x的长度小于n时，fft函数在x的尾部补零,当x的长度大于n时，fft函数会截断序列X。3.1.3程序及仿真结果functionst=yzzbN=800;%信号长度N为800；Fs=10000;T=1/Fs;Tp=N*T;%采样频率Fs为10KHz，Tp为采样时间；t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;fc1=Fs/10;%第一路条幅信号载波频率fc仁1KHz；fm1=fc1/10;%第路条幅信号调制信号频率fm1-100Hz；fc2=Fs/20;%第二路条幅信号载波频率fc2=500Hz；fm2=fc2/10;%第二路条幅信号调制信号频率fm2=50Hz；fc3=Fs/40;%第三路条幅信号载波频率fc3=250Hz；fm3=fc3/10;%第三路条幅信号调制信号频率fm3=25Hz；xt仁cos(2*pi*fm1*t).*cos(2*pi*fc1*t);%产生第一路条幅信号;xt2二cos(2*pi*fm2*t).*cos(2*pi*fc2*t);%产生第二路条幅信号；xt3二cos(2*pi*fm3*t).*cos(2*pi*fc3*t);%产生第三路条幅信号;st=xt1+xt2+xt3%三路信号相加，得到复合信号；fxt=fft(st,N);%用FFT算法计算st的频谱；%以下为绘图命令，绘制复合信号的波形和频谱subplot(2,1,1);plot(t,st);grid;xlabel('t/s');ylabel('s(t)');axis([O,Tp,min(st),max(st)]);title('(a)s(t)的波形')subplot(2,1,2);stem(f,abs(fxt)/max(abs(fxt)),'.');grid;title('(b)s(t)的频谱')axis([0,Fs/8,0,1.2]);xlabel('f/Hz');ylabel('幅度');3.1.4分析从仿真图中可以很清楚地看出信号包含6个频率分量，分别为225Hz,275Hz,450Hz,550Hz,900Hz,1100Hz。其中225Hz，275Hz关于载波频率f对称，两频率间隔50Hz，即2fo450Hz,550Hz关于载波频率f，两频率间隔100Hz，即2f，m3c2m2900Hz,1100Hz关于载波频率f对称，间隔频率为200Hz,即2f，可以看出这是clml三路混合的DSB信号，与前面分析一致。3.2数字滤波器设计概述3.2.1.FIR滤波器与IIR滤波器的比较无限长单位冲击响应（IIR）滤波器有以下几个特点：系统的单位冲击响应是无限长的，系统函数在有限z平面（0<|Z<）上有极点存在；结构上存在着输出到输入的反馈。有限长单位冲击响应（IIR）滤波器有以下几个特点：系统的单位冲击响应是有限长的，系统函数在z平面只有零点，而全部极点都在z=0处；结构上主要是非递归结构。在相同的技术指标下，IIR滤波器由于存在着输出对输入的反馈，所以可用比FIR滤波器较少的阶数来满足指标要求，所用的存储单元少，运算次数少，较为经济。FIR滤波器可得到严格的线性相位，而IIR滤波器做不到这一点，IIR滤波器的选择性越好，其相位的非线性愈严重。如果IIR滤波器要求得到线性相位，必须加全通网络进行校正，这样会大大增加滤波器的阶数。FIR滤波器主要采用非递归结构，有限精度运算中它都是稳定的，有限精度运算误差也较小oIIR滤波器极点必须在z平面单位圆内才能稳定，运算中的四舍五入处理有时会引起寄生振荡。从设计上看，IIR滤波器可以利用模拟滤波器设计的现成的闭合公式、数据和表格,因而计算工作量小，对计算工具要求不高。而FIR滤波器只有借助于计算机。在本次设计中我们选择IIR滤波器，除了IIR滤波器的上述优点外，由于IIR滤波器计算量小，我们可以通过手动计算一部分结果，再用计算机仿真，对比结果，以充分理解滤波器设计的过程。同时也能更好地掌握各种模拟滤波器。3.2.2IIR数字滤波器设计中的模拟原型滤波器比较典型的模拟滤波器有巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器和椭圆(Ellipse)滤波器等。其中，巴特沃斯滤波器又叫最平坦响应滤波器,顾名思义，它的响应最为平坦，通带内没有波纹，其频率响应在通带和阻带中都是单调的，且在靠近零频处最平坦，而在趋向阻带时衰减单调增大，缺点是从通带到阻带的过渡带宽，对于带外干扰信号的衰减作用弱。切比雪夫滤波器又分为切比雪夫I型滤波器和切比雪夫II型滤波器。切比雪夫I型滤波器在整个通带内纹波最小，在阻带内随频率单调递增；切比雪夫II型滤波器在通带内随频率光滑且单调递增，零频处最为平坦，在整个阻带内的纹波最小，它们的过渡带较巴特沃斯滤波器陡峭。巴特沃斯滤波器和切比雪夫滤波器的传输函数都是一个常数除以一个多项式，为全极点网络，所有的零点在无穷处，仅在无限大阻带处衰减为无限大，而椭圆函数滤波器在有限频率上既有零点又有极点。极零点在通带内产生等纹波，即它在整个通带和阻带上都具有最小的等纹波，这一点区别于在通带和阻带都平坦的巴特沃斯滤波器，以及通带平坦、阻带等波纹或是阻带平坦、通带等波纹的切比雪夫滤波器。同时，阻带内的有限传输零点减少了过渡区，可获得极为陡峭的衰减曲线。在本次设计中由于载波频率相差较小，所以我们采用过渡带较窄的椭圆滤波器，以达到理想效果3.2.3IIR数字滤波器的设方法及流程在MATLAB中设计IIR数字滤波器的步骤示意图如图二：图二经典法程序设计的具体步骤和函数总结如下：一、设计归一化模拟低通原型滤波器。1•确定所需数字滤波器H(z)的技术指标。低通、高通为w、①、R、Rpsps带通、带阻为wwwwRRulsuslps2•转为相对应类型的模拟滤波器H(s)的指标00，RR保持不变。psps转换时若采用冲击响应不变法则0二w..T.21转换时若采用双线性变换法则0=2tan1wT2注意：冲击响应不变法仅适合低通和带通滤波器，双线性变换法对所有类型滤波器均适合。H(s)转换为归一化模拟低通滤波器H(p)的指标九九LPps如果设计的是高通滤波器则九=00c如果设计的是带通滤波器则九二mt，其中n=0Bn如果设计的是带阻滤波器，则九二n，其中n=0B0转换后的技术指标设计使用滤波器阶数选择函数，确定最小阶数N和固有频率,即边缘频率Wn,根据选用的模拟低通滤波器的类型可以分别用函数buttord、cheb1ord、cheb2ord、ellipord。5•运用最小阶数N产生模拟滤波器原型，模拟低通滤波器的创建函数有：butter、cheby1、cheby2、ellip。二、运用固有频率Wn把模拟低通滤波器原型转换为模拟高通、带通、带阻滤波器,可分别用函数lp2hp、Ip2bp、lp2bs。三、运用冲击响应不变法或双线性变换法把模拟滤波器转换成数字滤波器，分别用函数impinvar和bilinear来实现。四、画出所设计的滤波器幅频图，验证设计结果是否满足指标要求，画图函数为freqz。3.2.4MATLAB中的相关函数Matlab的信号处理软件提供了设计椭圆滤波器的函数：ellipord函数和ellip函数等。.Ellipord函数的功能是求滤波器的最小阶数，其调用格式为：[n,Wn]=ellipord(Wp,Ws,Rp,Rs,'s')其中：n-椭圆滤波器最小阶数；Wp-椭圆滤波器通带截止角频率；Ws-椭圆滤波器阻带起始角频率；Rp-通带波纹(dB)；Rs-阻带最小衰减(dB)；.Ellip函数的功能是用来设计椭圆滤波器，其调用格式：[b,a]=ellip(n,Rp,Rs,Wp)[b,a]=ellip(n,Rp,Rs,Wp,'ftype')其中：‘ftype'=‘high'为高通滤波器，‘ftype'=‘low'为低通滤波器，‘ftype'=‘stop'为带阻滤波器。返回长度为n+1的滤波器系数行向量b和a,进而求得传递函数H(z):H_B(z)_b(l)+b(2)z-i+••"(n+1)z-»zA(z)1+a(2)z-i+・・・a(n+1)z-nellipap函数是求椭圆模拟低通滤波器的原型.格式：[z,p,k]=ellipap(n,Rp,Rs)其中：n-椭圆滤波器最小阶数，Rp-通带波纹(dB)；Rs-阻带最小衰减(dB)。用此函数可以设计出椭圆低通滤波器原型。其传递函数为：H(s)=k[s—z(DMs—z(2)]...[s—z(n)]

s-[s—p(1)][s—p(2)]...[s—p(n)]lp2hp函数是低通到高通模拟滤波器的转换函数。格式：[BT,AT]=lp2hp(a,b,Wp)说明：a，b是低通滤波器原型的分子分母系数，Wp是所需的高通滤波器的通带截止频率。lp2hp可将截止频率为1rad/s的模拟低通滤波器原型转变为截止频率为Wp的高通滤波器。lp2bp函数是低通到带通模拟滤波器的转换函数。格式：[BT,AT]=lp2hp(a,b,Wo,Bw)说明：a，b是低通滤波器原型的分子分母系数，Wo是所需的带通滤波器的中心频率，Bw是其带宽。lp2hp可将截止频率为1rad/s的模拟低通滤波器原型转变为具有指定带宽Bw和中心频率Wo的带通滤波器。.filter函数功能：利用IIR滤波器和FIR滤波器对数据进行滤波。格式：y=filter(b,a,x)说明：filter采用数字滤波器对数据进行滤波，其实现采用移位直接II型结构，因而适用于IIR和FIR滤波器。滤波器的系统函数为b+bz-1+bz—2+...+bz—mH(z)=012m—1+az—i+az—2+...+az—n12n即滤波器系数a=[aaa...a],b=[bb...b],输入序列矢量为x。这里,012n01m标准形式为a=1，如果输入矢量a时，a工1，则MATLAB将自动进行归一化系数的00操作；如果a=0，则给出出错信息。y二filter(b,a,x)利用给定系数矢量a和b对0x中的数据进行滤波，结果放入y矢量中，y的长度取max(N,M)。刀.freqz函数功能：离散时间系统的频率响应。格式：[h,w]=freqz(b,a,n)说明：freqz用于计算数字滤波器H(Z)的频率响应函数H(e」3)。[h,w]=freqz(b,a,n)可得到数字滤波器的n点幅频响应值，这n个点均匀地分布在[0,n]上,并将这n个频点的频率记录在w中，相应的频响值记录在h中。要求n为大于零的整数,最好为2的整数次幕，以便采用FFT计算，提高速度。缺省时n=512。3.3高通滤波器设计3.3.1经典法混合信号的高频段包含900Hz和1100Hz两个频率，所以高通滤波器的通带频率可以为800Hz,阻带频率可以为700Hz。对应的数字滤波器的指标为：①二f•2•兀F二0.50256①二f-2•兀：F二0.43982R二0.1R二60ppSssSPs2o2oQ=—tan(—)=0.51351Q=_tan(—)=0.44705pT2sT2QQQ九二◎这里取Q=Q即九=1，则九二c二p二1.1487pQcppsQQpssFs=10000；T=1；[N,Wn]二ellipord(1,1.1487,0.1,60,'s');%确定高通滤波器的阶数和边缘频率[z0,p0,k0]=ellipap(N,0.1,60);%得到系统的左半平面归一化零极点p=p0*Wn;z=z0*Wn;%将零极点乘以Wn,得到非归一化零极点k=k0*Wn^N;%将"乘以Wn八N，得到非归一化kb=k*real(poly(z));%有零点计算模拟低通滤波器分子系数向量a=real(poly(p));%有极点计算模拟低通滤波器分母系数向量[bb,aa]=lp2hp(b,a,wp);%转换为模拟高通滤波器%由双线性不变法得到数字高通滤波器的分子分母系数向量[bt,at]=bilinear(bb,aa,1/T);[h,w]=freqz(bt,at);%求数字高通的频率响应%以下为绘图命令，绘制高通滤波器的幅度响应plot(w*Fs/(2*pi),20*log10(abs(h)/max(abs(h))));axis([0,1200,-65,0])title('椭圆滤波器的幅度响应');xlabel('f/Hz');ylabel('分贝数');grid3.3.2直接法Fs=10000;wp=800*2/Fs;ws=700*2/Fs;Rp=0.1;Rs=60;%数字高通滤波器的设计指标[n,Wn]=ellipord(wp,ws,Rp,Rs,'s');%确定高通滤波器的阶数和边缘频率[B,A]=ellip(n,Rp,Rs,Wn,'high');%确定滤波器系数[h,w]二freqz(B,A,512);%求数字高通滤波器的频率响应%以下为绘图命令，绘制高通滤波器的幅度响应figure(2);plot(w*Fs/(2*pi),20*log10(abs(h)/max(abs(h))));axis([0,1200,-65,3]);title('数字高通滤波器的幅度响应');xlabel('f(Hz)');ylabel('分贝数');grid同样可以得到如图三所示的滤波器的幅频响应曲线：图三3.3.3仿真结果及分析用设计的滤波器对混合信号进行滤波，程序如下：N=800;Tp二N*T;t=0:T:(N-1)*T;k=0:NT;f二k/Tp;y2=filter(B,A,st);%调用滤波器实现函数，对st进行滤波Y2=fft(y2,N);%用FFT算法计算经过滤波的第一路信号频谱%以下为绘图命令，绘制第一路信号的波形和频谱figure(3),subplot(2,1,1)plot(t,y2);grid;xlabel('t/s');ylabel('y2(t)');axis([O,Tp,min(y2),max(y2)]);title('第1调幅信号的波形')gridon;subplot(2,1,2),stem(f,abs(Y2)/max(abs(Y2)),'.')axis([0,Fs/8,0,1.2]);title('第1路调幅信号的频谱')xlabel('f/Hz');ylabel('幅度');Grid第一路调幅信号的波形及频谱如图四0001a020.01004-0.K0.06&.0?0OBVs第1略谓幅信号的頻谱iiiiiiiir—T—-flttltt+ltt'StTtQ」亠d.q丄―〜■«■■■亠亠亠亠丄-亠亠■■■q亠=■■■亠—―」亠亠上亠」-q亠■■■亠亠d-亠勵4<*母*甲曾寸甲中寸呷t甲寸*甲1■电十中1■丨中甘寸中曹彎・020040D60000010001200f/Hz图四由计算数据可得滤波器阶数N=9,Wn=1；从滤波器的幅频响应曲线可以看出，通带范围为f>800Hz，阻带截止频率为700Hz，在0到700Hz的范围内，最大衰减大于60dB，在通带范围内衰减为0dB,且通带内衰减幅度几乎没有波动。在通过滤波器后的信号波形及频谱可以看出，该频谱只含900Hz和1100Hz两个频率信号。即原信号的高频部分。可以认为该高通滤波器的设计符合要求。3.4带通滤波器设计3.4.1经典法Fs=10000;T=1;%数字高通滤波器的设计指标wp=400*2*pi/Fs;ws=600*2*pi/Fs;wsu=350*2*pi/Fs;wsl=650*2*pi/Fs;Rp=0.1;Rs=60;%用双线性变换法转换为模拟高通滤波器的指标0megap=2/T*tan(wp/2);0megas=2/T*tan(ws/2);0megasu=2/T*tan(wsu/2);Omegasl=2/T*tan(wsl/2);%转化为归一化模拟低通滤波器的指标Omega0=sqrt(Omegap.*Omegas);

B=Omegap-Omegas;yitasu=Omegasu/B;yita0=Omega0/B;Ianmdas=(yitasu.八2-yitaO.八2)/yitasu;lanmdap=1;%确定滤波器的阶数和边缘频率[N,Wn]=ellipord(lanmdap,lanmdas,Rp,Rs,'s')[z0,p0,k0]=ellipap(N,Rp,Rs);%得到系统的左半平面归一化零极点p=p0*Wn;z=z0*Wn;%将零极点乘以Wn,得到非归一化零极点k=k0*Wn^N;%将"乘以Wn八N，得到非归一化kb=k*real(poly(z));%有零点计算模拟低通滤波器分子系数向量a=real(poly(p));w0=(Omegap+Omegas)/2;%有极点计算模拟低通滤波器分母系数向量[bb,aa]=lp2bp(b,a,w0,B);%转换为模拟高通滤波器%由双线性不变法得到数字高通滤波器的分子分母系数向量[bt,at]=bilinear(bb,aa,1/T);[h,w]=freqz(bt,at);%求数字高通的频率响应[h,w]=freqz(bt,at);%求数字高通的频率响应plot(w*Fs/(2*pi),20*log10(abs(h)/max(abs(h))));axis([0,1200,-65,0])title('椭圆滤波器的幅度响应')；xlabel('pirad/s');ylabel('分贝数');grid经典法设计的带通滤波器的幅度响应如图五刑同帯適鴻浜器的嘔侵响应图五3.4.2直接法Fs=10000;T=1/Fs;wp=[400600]*2/Fs;ws二[350650]*2/Fs;Rp=0.1;Rs=60;%数字带通滤波器的设计指标[n,Wn]=ellipord(wp,ws,Rp,Rs,'s');%确定带通滤波器的阶数和边缘频率[B,A]=ellip(n,Rp,Rs,Wn);%确定滤波器系数[h,w]二freqz(B,A);%求数字带通滤波器的频率响应%以下为绘图命令，绘制带通滤波器的幅度响应figure(4);plot(w*Fs/(2*pi),20*log10(abs(h)/max(abs(h))));axis([0,1200,-70,0]);title('数字带通滤波器的幅度响应');xlabel('f(Hz)');ylabel('分贝数');grid得到带通滤波器的幅频响应如图六

数字带通淤渡器的幅度响应图六3.4.3分析用设计的滤波器对混合信号进行滤波，得到的第二路信号的波形和频谱如图从MATLAB的计算结果中可以知道利用经典法设计出的滤波器N=6,而用直接法设计的滤波器N=7。用经典设计法设计出的滤波器通带内平坦，衰减几乎为OdB,

而用直接设计法设计出的滤波器通带内衰减大于O.ldB,但其阻带衰减更大，特性优于经典设计法。但两种滤波器均能正常滤出中频信号。在通过滤波器后的信号波形及频谱可以看出，该频谱只含450Hz和550Hz两个频率信号。即原信号的高频部分。可以认为该高通滤波器的设计符合要求。3.5低通滤波器设计3.5.1经典法Fs=10000;wp=300*2*pi/Fs;ws=400*2*pi/Fs;Rp=0.1;Rs=60;T=1;Omegap=2/T*tan(wp/2);Omegas=2/T*tan(ws/2);[N,Omegac]=ellipord(Omegap,Omegas,Rp,Rs,'s')[z0,p0,k0]=ellipap(N,Rp,Rs);p=p0*Omegac;z=z0*Omegac;k二k0*0megac八N;b=k*real(poly(z));a=real(poly(p));[bt,at]=bilinear(b,a,1/T);[h,w]=fre