吉布斯现象的matlab实现

信号与系统实验
实验者信息：
一：实验题目
运用matlab 验证吉布斯现象 二：实验原理
对于具有不连续点（跳变点）的波形，所取级数项数越多，近似波形的方均误差虽可减小，但在跳变点处的峰起（上冲）值不能减小，此峰起随项数增多向跳变点靠近。

（详见课本279页） 三：实验内容
1.计算N 取不同值时的方均误差N E 的值
2.用matlab 画出N 取不同值时
N
0n 1n 1n 1S(t)=a [a cos(n t)b sin(n t)]ωω=++∑的图。

四：实验过程
我们取课本99页的函数f (t)作为实验的目标函数
-E/2 (T /2t T /4)f(t)=E/2(T /4t T /4)
-E/2(T /4t T /2)-<<-⎧⎪
-<<⎨⎪<<⎩
matlab 有专门产生矩形波的函数square(t)，为了显示对称波形，
将它移动四分之一个周期，即square(t+T/2)，然后再乘上我们此题中的振幅E
2
,因此得到
f_t=(E/2)*square(t+T/2);
已知周期函数f (t)的傅里叶级数为
0n 1n 1n 1f(t)=a [a cos(n t)b sin(n t)]
ωω∞
=++∑
因为此题的函数既是偶函数，又是奇谐函数。

因此在它的傅里叶级数中只可能含有奇次谐波的余弦项。

所以n 2E n a sin()n 2ππ=
于是
1112E
11
f (t)[cos(t)cos(3t)cos(5t)]35ωωωπ=
-+-
若取傅里叶级数的前(2N 1)+项来逼近周期函数f (t)，则有限项傅里叶级数为
N
0n 1n 1n 1S(t)=a [a cos(n t)b sin(n t)]
ωω=++∑
进一步简化为
2N 1N n 12E sin(nt)
S T n -==∑ 为了得到对称波形，将它移动四分之一个周期，即
2N 1
N n 1
n T
sin(nt )2E
2S T
n
-=⨯+
=
∑
N S 就是代码中的S
这样用S(t)逼近f (t)引起的误差函数为N (t)f (t)S(t)ε=- 方均误差等于
010
t T 2
2
N N N t 11E (t)(t)dt T εε+==
⎰
在程序中我们取t=-T:0.0001:T;将t 等分成多份，所以不能直接用积分函数int （即不连续），所以依据定义的方法来求积分（极限法）
1T
2
N
1t 0
t T lim((t )t)ε∆→=-∆∑因为均分所以积分变为1T
2N
t T
t ε
=-∆⨯
∑其中
1T
2N
t T
ε
=-∑为
所有分量1t 对应的2
N ε值的总和
求N E 时应除以总周期即2T ，t
2T
∆即是将区间分为多少等份的倒数，matlab 中为
1
length(t)
，其中length(t)为将区间分成的等份数。

综上 ，我们得1T
2N
t T
N sum((f_t-S).^2)
E length(t)
length(t)
ε
=-=
=
∑ （后面是用在
matlab 中的计算方法）
第一问：
n 2N 1=-，n 为n E 的计数，
运行matlab 分别令N 1;N 3;N 5=== 得出123E 0.7577;E 0.3975;E 0.2678===…………如下表1
表1 方均误差
第二问：
-4
-3-2-1
01234
-3-2
-1
1
2
3
时间t
相应的函数值f (t )和S (t )
有限级数 S1
-4
-3-2-1
01234
-3-2
-1
1
2
3
时间t
相应的函数值f (t )和S (t )
有限级数 S2
-4
-3-2-1
01234
-3-2
-1
1
2
3
时间t
相应的函数值f (t )和S (t )
有限级数 S3
-4
-3-2-1
01234
时间t
相应的函数值f (t )和S (t )
有限级数 S20
放大跳变点附近的图像得到下图
五：实验结论与分析
由第一问的表1和第二问的图像可知，随着N 的增大，即随着所取级数项数的增多，进似方均误差减小，且峰起随着项数的增多向跳变点靠近，并且从有限级数S20跳变点附近的放大图可知峰起值趋于相同（峰起值趋近于跳变值的9%，这里不再予以验证）。

代码如下：
clear; E=4; T=pi;
t=-T:0.0001:T;%周期 N=input('N'); i=1;
for i=1:N%设出i ，让它变化，从而绘出S1，S2，S3……SN a=0;
1.3
1.35
1.4
1.45
1.5 1.55 1.6
1.65
1.7
1.75
2.052.12.152.22.252.32.35时间t
相应的函数值f (t )和S (t )
有限级数 S20
for n=1:2:i
b=(E/2)*sin(n*t+n*T/2)/n;%移动四分之一周期以显示对称波形a=a+b;
end
S(i,:)=4*a/T;
f_t=(E/2)*square(t+T/2);
plot(t,S,t,f_t)
xlabel('时间t')
ylabel('相应的函数值f(t)和S(t)')
title('有限级数S20')
end
E_N=sum((f_t-S(N,:)).^2)/length(t)。