MATLAB大作业

MATLAB⼤作业班级：姓名：学号：2014年4⽉24⽇MATLAB⼤作业

——实现傅⾥叶级数的Gibbs现象⼀、要解决的问题及意义

⽤MATLAB实现傅⾥叶级数的Gibbs现象，更为直观、具体、形象地认识并更为深⼊地了解Gibbs现象。⼆、Gibbs现象的简单介绍Gibbs现象的来由

⾃从法国数学家傅⾥叶(Fourier)发明了傅⾥叶级数，关于级数理论的研究随即⾛向了⼀个新时代的⾥程碑，傅⾥叶级数在纯粹数学、应⽤数学、物理，乃⾄⼯业应⽤上都取得了辉煌的成就，其优点是光滑性好，收敛条件要求低（尤其是⼤范围⼀致收敛⽐幂级数要求低得多），但是它也并不是完美⽆缺的。19世纪末，美国数学家

J.W.Gibbs在研究⼀个函数的三⾓级数展开时，发现了这种级数展开的缺点,即不只在被展开函数的不连续点处，甚⾄在该点附近，其傅⾥叶级数也不收敛到该函数，这就是著名的Gibbs（吉布斯）现象。Gibbs 现象（⼜叫吉布斯效应）

⼀般来说，将周期信号分解得到的三⾓函数形式的傅⾥叶级数的项数是⽆限的。也就是说，通常只有⽆穷项的傅⾥叶级数才能与原函数精确相等。但在实际应⽤中，显然⽆法取⾄⽆穷多项，⽽只能采⽤有限项级数来逼近⽆穷项级数。⽽且，所取项数越多，有限项级数就越逼近原函数，原函数与有限项级数间的⽅均误差就越⼩，⽽且低次谐波分量的系数不会因为所取项数的增加⽽变化。

将具有不连续点的周期函数（如矩形脉冲)进⾏傅⽴叶级数展开后，选取有限项进⾏合成，当选取的项数越多，在所合成的波形中出现的峰起越靠近原信号的不连续点。当选取的项数很⼤时，该峰起值趋于⼀个常数，⼤约等于总跳变值的9%。这种现象称为吉布斯现象。三、数学⽅法求解

连续时间周期信号的傅⾥叶级数

频域分析法即傅⾥叶分析法，它是变换域分析法的基⽯。其中，傅⾥叶级数是变换域分析法的理论基础，傅⾥叶变换作为频域分析法的重要数学⼯具，具有明确的物理意义，在不同的领域得到⼴泛的应⽤。1.连续时间周期信号的分解

根据⾼等数学的知识，任何周期为T的周期函数)(t

f，在满⾜狄⾥赫利条件时，则该

周期信号可以展开成傅⾥叶级数。傅⾥叶级数有三⾓形式和指数形式两种。1）三⾓形式的傅⾥叶级数三⾓形式的傅⾥叶级数为：)12(3,2,1)sin()cos(2

)sin()sin()cos()cos()cos(2)(11213210b -??=Ω+Ω+=?+Ω+Ω+??+Ω+Ω+Ω+=∑∑∞=∞

=n t n t n t t t t t t f n n n nb a a b a a a a式中系数na 、n b

称为傅⾥叶系数，可由下式求得。()()()()()2

22222

012cos 2sin T T T T T T n n a f t dtT a f t n t dt

T b f t n t dt T ---==Ω=Ω()22- 其中，T π2=

Ω为基波频率,Ωn 为n 次谐波频率。如果将()21-式中同频率的正弦和余弦分量合并，则三⾓形式的傅⾥叶级数可表⽰为：()()0

1t cos 2n n n A f A n t ?∞==+Ω+∑ 1,2,3,n =

()23- 上式中001,2,arctann n

n n A a A n a b ??====- ()24- 00

cos ,1,2,sin n n n n n n a A a A n b A ??=??===-

可以看出，傅⾥叶系数n a 和n b

都是n 或()n Ω的函数，其中n a 和n A 是n 或()n Ω的偶函数，即有

n a -=n a ；⽽n b 和n ?是n 或()n Ω的奇函数，即有n b -=n b -。2） 指数形式的傅⾥叶级数 根据欧拉公式：][21)cos()()(n n t n j t n j n e e t n +Ω-+Ω+=+Ω()25-并考虑n A和

n ?奇偶性可将()23-改写为指数形式的傅⾥叶级数：(),0,1,2,3,jn t n n f t F e n ∞Ω=-∞==±±±∑)62(-即周期信号可分解为⼀系列不同频率的虚指数信号之和，式中n F称为傅⾥叶复系数，可由下式求得：()22

1T T jn t n F f t e dt

T-Ω-=()27-2. 连续时间周期信号的傅⾥叶综合

任何满⾜狄⾥赫⾥条件的周期信号，可以表⽰成式()21-或()26-的和式形式，

()21-或()26-式常称为连续周期信号的傅⾥叶级数综合公式。

⼀般来说，傅⾥叶级数系数有⽆限个⾮零值，即任何具有有限个间断点的周期信号都⼀定有⼀个⽆限项⾮零系数的傅⾥叶级数表⽰。但对数值计算来说，这是⽆法实现的。在实际的应⽤中，但我们可以⽤有限项的傅⾥叶级数求和来逼近。

为了⽐较有限项谐波的逼近情况，同时也让⼤家有⼀个更为直观的感受，下⾯给出了周期矩形波信号各次谐波的合成波形，如图1所⽰。

-505

-0.200.20.40.6基波-505-0.500.51

基波+2次谐波

-5

05-0.500.51

1.5基波+2次谐波+3次谐波-505-0.500.51

1.5基波+2次谐波+3次谐波+6次谐波图1 周期矩形脉冲信号的合成

由图1可见，当它所包含的谐波分量越多时，合成波形愈接近于原来的矩形波脉冲。由图1还可以看到，合成波形所包含的谐波分量愈多时，除间断点附近外，它越接近于原

矩形波脉冲。在间断点附近，随着所含谐波次数的增加，合成波形的尖峰愈接近间断点，但尖峰幅度并未明显减少。可以证明，即使合成波形所含谐波次数n →∞时，在间断点处仍有约9%的偏差，这种现象称为吉布斯（Gibbs ）现象。在傅⾥叶级数的项数取得很⼤时，间断点处尖峰下的⾯积⾮常⼩以致趋近于零，因⽽在均⽅的意义上合成波形同原波形的真值之间没有区别。3.吉布斯现象

我们知道满⾜狄⾥赫利条件的周期函数表⽰成的傅⽴叶级数都收敛。狄⾥赫利条件如下：1)在任何周期内，x(t)必须绝对可积；

2)在任⼀有限区间中，x(t)只能取有限个最⼤值或最⼩值； 3)在任何有限区间上，x(t)只能有有限个第⼀类间断点。

所谓的吉布斯现象就是：在x(t)的不可导点上，如果我们只取x(t)等式右边的⽆穷级数中的有限项作和X(t)，那么X(t)在这些点上会有起伏。具体现象如下图所⽰，以下分别为谐波次数为N=50，N=100，N=500合成波的情况。

-2-1.5-1-0.50.511.5200.5

1谐波次数N=50时的合成波形

-2-1.5-1-0.50.511.5200.5

1谐波次数N=100时的合成波形

-2-1.5-1-0.50.511.5200.5

1谐波次数N=500时的合成波形图2 不同时N 值时的合成波

从上⾯的图像中可以看出，当N=500的时候，合成波与原来的⽅波拟合得⾮常好，但是在不可导的点上，即为x=-1.5，x=-0.5，x=0.5，x=1.5这样的点的时候，合成波会有较⼤的波动，这就是⾮常明显的吉布斯现象。四、MATLAB 求解的流程图

五、⼩结

本次MATLAB课程⾄此已经接近尾声，学习MA TLAB的时间虽然很短暂，但在这⼏个星期的学习过程中收获颇多。本次作业的核⼼内容就是利⽤MA TLAB强⼤的图形处理功能，符号运算功能以及数值计算功能，实现傅⾥叶级数的Gibbs现象。在整个作业完成过程中，⾸先，我对所学的信号与系统与数字信号处理有了更深的了解，⽐如傅⽴叶级数、信号频谱等；其次，实现过程是通过MA TLAB软件完成的，MA TLAB 的图形功能强⼤，具有良好的⼈机界⾯，此次设计过程中熟练了MATLAB的编程，掌握了很多函数的作⽤及使⽤⽅法；最后，通过此次实际操作，我对设计所⽤到的软件MATLAB有了更加深刻地了解，MA TLAB不管在数值计算⽅⾯的功能很强⼤，⽽且其图形仿真功能更能满⾜各个领域的需要，因此我们以后更要经常运⽤MATLAB软件，使其成为⾃⼰不可或缺的辅助⼯具。

在编写本次作业的相关源程序的时候，我还浏览了⼤量的⽹站，在⽹上搜索了很多关于MATLAB的资料，在这个过程中，我才惊觉⽹络上也存在着⼤量的有⽤的信息，在以后的学习⾥还应该多注意充分合理有效地利⽤⽹络，通过⽹络来了解学习更多的东西。在收集资料的阶段，我还温习了《数字信号处理》这门课程⾥的相关知识，对已经学过的理论知识有了更进⼀步的认识和理解。

虽然我基本顺利完成了本次⼤作业的要求，但是我感觉到我对MA TLAB的理解我掌握还停留在⽐较浅的层次，要想真正掌握它还需要继续努⼒学习它。

这次⼤作业也使我明⽩了在知识的领域⾥我还存在着许多的不⾜，并且再⼀次深深地体会到理论和实践之间还存在着很⼤的差别。在以后的学习中，还应该更多地注意到将所学知识运⽤到实践中去，不断地训练⾃⼰的实际应⽤能⼒，同时也在实践中更好地理解⾃⼰所学的理论知识。附：源程序:t=-2:0.001:2;

N=1; % N为输⼊要达到的最⾼次谐波的次数c0=0.5;

fN=c0*ones(1,length(t));for n=1:2:N

fN=fN+cos(pi*n*t)*sinc(n/2);endfigureplot(t,fN)

axis([-2 2 -0.2 1.2])

title('(1) N=1','FontSize',12,'FontWeight','bold');grid ont=-2:0.001:2;

N=3; % N为输⼊要达到的最⾼次谐波的次数c0=0.5;

fN=c0*ones(1,length(t));for n=1:2:N

fN=fN+cos(pi*n*t)*sinc(n/2);end

figureplot(t,fN)

axis([-2 2 -0.2 1.2])

title('(2) N=3','FontSize',12,'FontWeight','bold');grid ont=-2:0.001:2;

N=10; % N为输⼊要达到的最⾼次谐波的次数c0=0.5;

fN=c0*ones(1,length(t));for n=1:2:N

fN=fN+cos(pi*n*t)*sinc(n/2);endfigureplot(t,fN)

axis([-2 2 -0.2 1.2])

title('(3) N=10','FontSize',12,'FontWeight','bold'); grid ont=-2:0.001:2;

N=20; % N为输⼊要达到的最⾼次谐波的次数c0=0.5;fN=c0*ones(1,length(t));for n=1:2:N

fN=fN+cos(pi*n*t)*sinc(n/2);endfigureplot(t,fN)

axis([-2 2 -0.2 1.2])

title('(4) N=20','FontSize',12,'FontWeight','bold'); grid ont=-2:0.001:2;

N=50; % N为输⼊要达到的最⾼次谐波的次数c0=0.5;fN=c0*ones(1,length(t));for n=1:2:N

fN=fN+cos(pi*n*t)*sinc(n/2);endfigureplot(t,fN)

axis([-2 2 -0.2 1.2])

title('(5) N=50','FontSize',12,'FontWeight','bold'); grid ont=-2:0.001:2;

N=200; % N为输⼊要达到的最⾼次谐波的次数c0=0.5;fN=c0*ones(1,length(t));for n=1:2:N

fN=fN+cos(pi*n*t)*sinc(n/2);endfigureplot(t,fN)

axis([-2 2 -0.2 1.2])

title('(6) N=200','FontSize',12,'FontWeight','bold'); 结果图：