用Excel手搓傅里叶级数（2）

3. 傅里叶级数

上一节，我们验证了 不同周期的正弦函数之间、不同周期的余弦函数之间、正弦函数和余弦函数之间都是正交的，因此可以在傅里叶级数展开时作为 基底。不同幅度的这些函数相加，就可以收敛于某个周期函数。哪个周期的正弦函数、哪个周期的余弦函数 所对应的幅度，也即乘以这些函数的系数，是我们要求取的傅里叶系数。

公式大概长得像下图这样。

其中a0是常数项，也可以视为cos(0)，（根据不同定义）可能乘以一个系数；

与余弦对应的系数是 a_i ;

与正弦对应的系数是 b_i ;

3.1 傅里叶级数，展开，验证

书中的例子，需要展开为傅里叶级数的函数周期为T=1。

（1）原理，以及准备各个函数对应的数据

根据T=1

以及 角速度ω=2π/T

ω=2π/T=2π

以上备用。

我们可以跳过 傅里叶系数是待展开函数在各个基底上的映射 这一原理，查到公式如下。

上式中，n与i同义，指不同角速度/周期的正弦或余弦所对应的一个系数。我们能够注意到 a_n与余弦cos对应，b_n与正弦sin对应，其余形式完全相同。

上式中的方框，在cos和sin之后的那两个方框，内容是 n*ω*t。

把上面备用的

ω=2π

T=1

代入到上面两个式子中，以b_n为例得到下式。

a_n的形式非常类似，仅把sin改为cos。

得到上述公式，距离在Excel中实施只差一步。

这个定积分形式与上一篇中我们所见到的 黎曼和 形式求不同周期的、正弦函数、余弦函数间的相关性的公式非常相似。相似到不过是换了函数而已。验证基底相关性为0时，我们所用的两个函数分别是sin,cos之类的，现在用的是 f(x)和sin。再乘以一个系数，这里是2。

我们回到上一篇中的下面这组数据，对它做傅里叶展开。

这样，我们得到待展开的函数对应的数据。

在上一篇中，我们已经求出 角度、弧度、sin(t)、cos(t)。

按同样的方法，保留角度、弧度，我们求出 sin(2π*n*x)和cos(2π*n*x)，其中n取1，2，3，4，5。

如下图所示，是当n=3时，以弧度作为x，绿框所在的这一列是sin(2π*n*x)。

这样，我们得到共10列不同周期的正弦和余弦函数对应的数据。

常数项，即 a₀，书中的作法是 常数项也是个函数，其数据为值为 y=SQRT(2)/2 ，即0.707左右，的直线。与其他教材上 a₀的定义不同，结果是等价的。

这样，我们有了常数项的函数所对应的数据。

现在，我们有了

待展开的函数 所对应的数据 （波）；
10列不同周期的正弦和余弦函数 所对应的数据；
常数项的函数 所对应的数据。

以上，在Excel中共12列。

（2）傅里叶系数

我在Excel中隐藏无关的行和列，仅展示 b₃这一系数的求取过程。

第一步，求得每一行，即黎曼和中的每个横坐标变化所对应的y。如下图所示。

其中G6是此前做好备用的 sin(2π*3*x)的第6个数据，
N6是待展开函数 的第6个数据，

这一行的结果暂存在b3列的S6单元格。

乘以1/1000，是因为我把整个一周期内的波形等分为1000份，因此每一行对黎曼和的贡献是1/1000。

如果忘了黎曼和，可以回顾下图，在上一篇中出现过。当时用于求内积，这里用于求傅里叶系数。

把b3所有行累加，再根据下图中的公式乘以2。

得到b3这个系数的值，如下。

值为0.424400615。

用同样的方法求出11个系数。

（3）验证

我们只使用当n=1时sin(2π*x)和cos(2π*x)这两个三角函数，因此傅里叶只涉到 a1,b1，以及常数项。这样得到的结果，在下图中表示为 “傅里叶 I”这一列。

“傅里叶 I”这一列的第一行，那个公式的含义从左向右为

蓝色 b1参数* 红色C3 正弦函数 +
紫色a1参数*D3余弦函数+
粉色a0*M3常数项函数。

验证一下对不对。如下图所示，蓝色是待展开的函数，橙色是“傅里叶 I”展开。

看起来好像对，又好像不太对。两个函数看起来有点像，但是区别又很大。

我们把更精细的傅里叶级数也加进来，得到下图。

其中绿色的，是n=1至n=5 傅里叶 V 级数展开的结果，有点像了。而且与刚刚的傅里叶 I级相比，肉眼可见像得多了。

但是我们可能会注意到一个问题，从I级到V级，一共应该有5条曲线，为什么这里只画了3条呢。有帖子提到方波（可以扩大范围至奇函数）的所有sin参数（和cos中的偶数周期倍数）非常小，小到接近0。所以，在这个例子中，每两级傅里叶级数波形极其相似，I被II覆盖，III被IV覆盖了。

不仅理论如此，在实践的数据中我们也可以看到，除了b1,b3,b5这三组，其余的都贡献微小。

肉眼看到的结果不足够有说服力，我们可以再求 均方差。方差越大，展开得到的波形与原函数就越不相似。

先求出每一行的差。我隐藏无关的行，在Excel中看起来如下图所示。

然后得到 差的平方。

最后对 差的平方 求均值，得到均方差。

对I,III,V三个级别的傅里叶展开求均方差。如下图所示，我们能够看到均方差越来越小。这表明，傅里叶级数越高，得么的曲线与原函数越接近。与我们肉眼的感觉是相同的。

（4）换周期，换采样间隔，换原函数

调程序的一大乐事，是在程序写完以后，修改代码，测试不同效果，或者反复程序，居然还不崩溃。在Excel里手搓傅里叶级数也一样，折腾各种变化很有意思。

书里的周期是1，我们换一个。

根据周期T=1求得ω，并把T和ω代入公式，得到下式。

制备出各个周期的正弦函数和余弦函数，备用。

制备出各个傅里叶系数，备用。

每一行。

上图中，每行的贡献为 1/1000*T = 1/1000*2π.

求黎曼和。

验证一下。

周期不同，看不出来什么区别，不够炫酷。

换采样间隔吧。1000份真是很多，翻页需要好久，如果在Excel中不隐藏数据行的话。改为20等份。

以下，周期使用 2 pi，未展开说明。

穷举弧度时，要说明等分为20而不是1000。

求傅里叶系数时，每行贡献1/20而不是1/1000。

其余不变。绘图验证，只有20个点，瘦骨嶙峋，果然看起来漂亮很多。

换个原函数。修改“波”这一列数据，改变原函数。瞬间Excel就给出新的傅里叶展开。

再换个原函数。

再换。

这个过程非常过瘾。改一下波形数据，对应的傅里叶级数就同时跟着变化。夫子步亦步，夫子趋变趋。

4. 待续

还可以有些展开，也非常有趣。不过，那些已经不应该在 用excel手搓傅里叶级数 这个故事中了。我们在以后的故事里再继续。

（1）傅里叶变换，模和幅角

又

谐波的强度为 sqrtp(a^2+b^2)
幅角为 cos (arctan (-b/a))

因此，可以把每对 cos 和 sin 转换为 只用cos。

这样就把原函数变换到了频域上，完成了傅里叶变换。

仍然可以用Excel求出（不是无穷级数的）结果，与其他工具计算的结果对比。

（2）用Geogebra绘制动图，圆周与傅里叶级数展开，本轮和均轮

我们可以使用Geogebra绘制出 红色大圆之上的绿色中圆之上的蓝色小圆上某个点的轨迹。红、绿、蓝三个圆的周期和半径不同，对应了不同的傅里叶级数。

改变横坐标，可以看到红、绿、蓝三个圆转动导致纵坐标描绘出方波。在红圆的圆心，我们可以看到蓝紫色的“行星”轨迹，像火星一样忽快忽慢，有时倒退。

怎么把这个效果做出来呢？

（3）快速傅里叶变换FFT

因为FFT求系数数的时间复杂度更低一些，所以有可能在Excel中手搓出相当多的级数而工作量可以容忍。

不过，如上所说，用Excel手搓傅里叶级数这一场的大幕落下，那些都应该在另一些故事中再继续。