手把手教你编写傅里叶动画

嵌入式ARM 2022-12-30 12:00

测量工具选型不当可能会毁掉你的PCB调试成果 工程师调试电路的秘密武器

先来看几个比较有艺术性的动画：

上面三个绘制动画（苹果logo、爱心、中国结）是我用一千个圆，让他们各自按照自己的方向与速度进行滚动，然后把他们首尾进行相连，同时把最后一个圆上某一点的路径记录下来，这便是图中绿色的轮廓。说到这里，你有没有感觉：卧槽，怎么做的？我也想做！

这篇文章我会花大量的篇章来介绍这个方法（读者看到动态图可能会有轻微的卡顿，是因为微信公众号限制动态图的帧率，但实际渲染远比这平滑）。

事情先回溯到两百多年前，法国数学家傅里叶提出一个观点：“任何一个连续周期性函数都可以用正弦函数与余弦函数的和来表示”。虽然这句话在当时没有引起重视，但时至今日，我们已经没有任何权利去怀疑这一说辞的分量，因为它的分量已经远远超过当年傅里叶所表达问题的本身。

从傅里叶级数衍生出来的傅里叶变换，它的重要性在任何一门工程领域尤其是电子信息领域不言而喻，毫不夸张地说，我们每天能够用诸如微信这种即时通讯工具跟好友聊天，除了感谢腾讯之外，还应该感谢一下傅里叶同学，因为没有他的那句话，绝对不可能有电子通讯的今天。

但是，我相信大部分学过《高等数学》的人，对傅里叶变换的概念肯定还停留在只会做题的阶段，并不了解傅里叶变换的本质，而傅里叶变换却是我做图像处理中经常用到的一种算法。所以，我写这篇文章的原因就在于想跟大家聊一聊我自己所了解的傅里叶变换，但要说傅里叶变换却绕不开傅里叶级数，所以本篇文章我会主要讲解傅里叶级数，而傅里叶变换在实际工程中的应用例如图像降噪、图像分类、图像压缩等技术，我会在另外一篇篇章中介绍。

如果读完本文您觉得对你有用，我希望您能够把它分享出去，让更多的人看到。同时，如果你发现有错误，我也希望的到您的指正。

提到「变换」这个词，想必大家一定不会陌生，没错，我们在大学时候高数里面的变换实在是太多了，诸如拉普拉斯变换、泰勒展开等。可是大家有没有想过这个问题：我们们为什么要弄这些千奇百怪的「变换」呢，难道是为了好玩吗？

答案显然不是，之所以我们需要这么多变换，最主要的原因还是我们懒啊，这些变换能使我们把非常复杂的问题简单化，甚至在大多数情况下如果我们不通过这些变换我们的问题是不能解的。

仔细回想一下，其实我们在生活中已经不知不觉地应用了各种变换，例如电子秤并非直接测量实际的体重，而是根据内部压力传感器电阻的变化通过一定的转换来输出实际的数值的。类似的事情，其实在古代已经有了，例如阿基米德通过测量溢出水质量来测量皇冠的密度，曹冲通过石头的重量来测量大象的体重等等，这些都是一种变换。

小曹冲通过测量吃同样水深的石头来测量大象的体重，这也是一种变换。

那么，我今天要说的傅里叶变换到底有什么好处呢？换句话说，我们为什么要进行傅里叶变换呢？

原因很简单：傅里叶变换能把一个周期性函数分解成多个正余弦函数的叠加，这种变换对我们分析信号有诸多好处，尤其在数字信号处理或者电工学中有重大应用，例如脉冲信号一般都是矩形波，如果我们把它进行傅里叶变换，我们就能看到组成这个矩形波的各谐波的频率、相位、振幅等等，由于正余弦函数的正交性，不仅方便我们计算，而且对我们分析问题带来很大方便。

那么，如何理解傅里叶级数呢？如果你你记得公式的话，它是长这样子的：

$f(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty}(a_n\cos \frac{n\pi t}{l}+b_n\sin \frac{n\pi t}{l})$

其中 $f(x)$ 是以 $2l$ 为周期的函数，

余弦系数 $a_n:$

$a_n=\frac{1}{l}\int_{0}^{2l}f(t)\cos \frac{n \pi t}{l}dt(n=0,1,2...)$

正弦系数 $b_n:$

$b_n=\frac{1}{l}\int_{0}^{2l}f(t)\sin \frac{n \pi t}{l}dt(n=1,2,3...)$

可这个公式不太完美，过于复杂，而且不能从直观上看出傅里叶公式的奥妙之处，因此我们换种方式，用欧拉公式去化简它，而傅里叶级数的精髓也在这里。但是，要引入欧拉公式之前，我们要从虚数说起。

一、虚数的几何意义

虚数 $i$ 究竟代表什么意思呢，我在高中时候曾经对这个公式 $i^2=-1$ 大惑不解，究竟一个什么样的数平方居然会等于负数呢，不知你是否跟我有同样的疑问，而这个问题直到我大学才算彻底明白。

其实虚数的真正含义，代表着旋转：

假如一根数轴上有一个点 $B(+1)$ ，如果我们对这根数轴正半部分做两次逆时针旋转，那么 $B$ 就会变成 $B'(-1)$ ：

这样我们可以说：

$(+1)*(旋转90^\circ)*(旋转90^\circ)=-1$

如果我们把逆时针旋转 $90^\circ$ 用 $i$ 来表示，那么上式就是：

$1 * i *i=-1$

也就是： $i^2=-1$

明白了吧，虚数代表一个旋转量，也就是说，只要你看到虚数，就应该想到旋转。

二、复数的几何意义

我们进一步拓宽到复数领域，如果把上图中的纵轴表示为虚数轴，横轴代表实数轴，那么任何一个复数 $a+bi$ 都能在上述坐标轴中找到唯一的坐标，这样能够极大方便我们计算向量的旋转问题。

例如如果一个复数 $3+2i$ ,如果我们想将它的方向逆时针旋转 $45^\circ$ 那么我们只需要将其乘以一个 $45^\circ$ 的复数 $1+i$ 即可：

$(3+2i)(1+i)=1+5i$

那么新的方向就是 $1+5i$ ：

也就是说，对于任意一个复数向量，我们如果想对其进行旋转，那么只需要乘以对应角度的复数即可。这里给出一个简单粗暴的证明方法：

如上图，如果两个向量的长度分别为 $||v||=r_1,||u||=r_2$ ,我们令：

$a+bi=r_1(\cos \alpha+i\sin\alpha)$
$c+di=r_2(\cos \beta+i\sin\beta)$

咦，这个公式右边看起来怎么那么像欧？没错，我们就是为了引入欧拉公式，下面会讲到。上式左边相乘，得到：

$(a+bi)(b+ci)=r_1r_2(\cos \alpha+i\sin\alpha)(\cos \beta+i\sin\beta)$

利用和差公式化简右边，可得到：

$(a+bi)(b+ci)=r_1r_2(\cos(\alpha+\beta)+i\sin(\alpha+\beta))$

这就说明，两个复数相乘，结果就等于旋转半径相乘、旋转角度相加。

三、欧拉公式的几何意义：一种旋转运动

我们再引入欧拉公式，根据泰勒公式（具体推导过程请参考我另外一篇文章来看看比尔盖茨当年写的BASIC解释器源代码吧，你就知道泰勒级数有什么用了）：

$e^x=1+x+\frac{x^2}{2!}+...+\frac{x^n}{n!}$

用 $ix$ 代替 $x$ 代入上式，得到：

$e^{ix}=1+ix-\frac{x^2}{2!}-\frac{ix^3}{3!}+\frac{x^4}{4!}+\frac{ix^5}{5!}...$

我们把带有虚数项与不带虚数项分开：

$e^{ix}=1-\frac{x^2}{2!}+\frac{x^4}{4!}+...+(-1)^{n}\frac{x^{2n}}{(2n)!}+$

$+i(x-\frac{x^3}{3!}+\frac{x^5}{5!}+...+(-1)^{n-1}\frac{x^{2n-1}}{(2n-1)!})$

由于正弦的泰勒展开为：

$\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}+...+(-1)^{n-1}\frac{x^{2n-1}}{(2n-1)!}$

余弦泰勒展开为：

$\cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}+...+(-1)^{n}\frac{x^{2n}}{(2n)!}$

将上面两个式子带入 $e^{ix}$ ，最终，我们得到：

$e^{ix}=\cos x+i\sin x$

这便是著名的欧拉公式。特别地，当 $x=\pi$ 时，上式可化为：

$e^{i\pi}+1=0$

这就是著名的欧拉恒等式，它被称为上帝最美公式，因为仅仅一个公式就包含了自然底数 $e$ ，圆周率 $\pi$ ，实数 $0$ ，虚数 $i$ 这么多重要的数学元素。额，扯远了。。。

但是，不管是欧拉公式还是欧拉恒等式，但凡你第一次看到诸如这种 $e^{ix}$ 虚数在指数上面，一定会感觉如鲠在咽，至少不会比指数为无理数那样更舒服，因为你至少不能把它理解成 $i$ 个 $e$ 相乘吧，这在虚数范围内毫无意义。那么 $e^{ix}$ 究竟是什么意思呢？答案是它代表着一种运动，而且是旋转的运动，为什么呢，请往下看：

我们知道，自然底数 $e$ 的定义是：

$e^x=\lim_{n\to\infty} \left( 1+\frac{x}{n} \right)^n$

当 $x=1$ 时，我们就得到了 $e$ .我们把这个式子扩展到复数领域，用 $i\theta$ 代替 $x$ ，得到：

$e^{i\theta}=\lim_{n\to\infty} \left( 1+\frac{\theta}{n}i \right)^n$

也就是说， $e^{i\theta}$ 相当于 $n$ 个复数向量 $(1+\frac{\theta}{n}i)$ 相乘，而我们上面已经证明了，复数相乘的几何意义就是对这个复数平面向量旋转一定的角度，同时长度进行了一定的伸缩。所以 $n$ 次相乘就代表着将复向量 $(1+\frac{\theta}{n}i)$ 进行了 $n$ 次旋转，每次旋转的角度是 $\arctan(\frac{\theta}{n})$ 。

为了更加直观理解这一操作，我们再引入复平面单位圆，这里我们先令 $\theta=\pi,n=10$ ，这样，复平面内这个复向量的原始位置就是 $(1+\frac{\pi}{10}i)$ 。

由于连续进行了10次相乘，每次的相乘，都意味着对这个向量进行旋转操作，同进行伸缩。所以10次相乘就是这个向量从原始位置不断旋转，连续十次，每次的角度都为 $\arctan(\pi/10)$ 弧度：

可见，每次经过一次旋转，长度都为上一次的 $\sqrt{1+(\pi/10)^2}$ 倍，同时角度增加 $\arctan(\pi/10)$ 。当 $n=10$ 时候，旋转的角度基本上接近了 $\pi$ 。

我们再增大 $n$ 试试，当 $n=50$ 或者 $n=100$ 时候，看旋转情况：

$n=50$ ：

$n=100:$

可见，随着 $n$ 不断增大，复向量 $z_1$ 的长度不断逼近1，同时最终角度趋向于 $\pi$ 。不难看出，当 $n\to\infty$ 时候，复向量最终会旋转 $\pi$ 弧度，长度等于1。

因此，我们得出一个重要结论： $e^{i\theta}$ 的实际几何意义就是副向量绕复平面单位圆做 $\theta$ 角度的旋转。如果我们用时间 $t$ 来替换 $\theta$ ，那么 $e^{it}$ 就意味着在一个单位时间内复向量做了 $t$ 弧度的旋转。

这意味着，如果 $e$ 的指数不一样的话，那么他们在单位时间内旋转的角度是不一样的，也就是说指数的大小决定了他们的旋转速度。并且，指数的正负还能决定它们的旋转方向，为什么呢？

我们假设一个复向量的方向为 $e^{it}$ ,对其求导：

$\frac{d(e^{it})}{dt}=ite^{it}$

可见，求导后上式右边乘了一个虚数 $i$ ，前面我们说了，虚数的本质就是做90度的旋转，那么意味着这个瞬时旋转速度肯定垂直于复向量的方向。而且指数上的虚数 $i$ 的系数正负，能够决定它的旋转方向，例如 $e^{-it}$ 则意味着初始向量朝相反的方向做圆周运动。

$e^{it}$ 逆时针旋转：

$e^{-2it}$ 顺时针旋转，速度为原来的两倍：

如果初始向量不在1的位置，而是另外一个位置角度，例如初始向量的位置在 $\pi/4$ 的位置，而且长度为2，根据前面我们证明的复数相乘的结论，我们只需要在 $e^{it}$ 前面乘以相应长度为2，角度为 $\pi/4$ 的复向量即可： $(2+2i)e^{it}$

初始角度： $\pi/4$ ，长度：2。

值得注意的是，一个复向量绕单位元旋转一周后，最终位置与初始位置一样，这意味着在任何整个周期区间内，对 $e^{ix}$ 积分恒为0：

$\bbox[yellow]{\int_{0}^{2\pi}e^{it}dt=0\qquad (1)}$

这个结论既有趣又重要，它在后面能为我们大大简化计算参数系数的计算量。

而且，根据上面推导出来的欧拉公式：

$e^{it}=\cos t+i\sin t$

用 $-t$ 替换上面的 $t$ ，两者联立，进一步推导出： $\cos t=\frac{e^{it} + e^{-it}}{2},\sin t=\frac{e^{it} - e^{-it}}{2i}$

四、用欧拉公式化简傅里叶级数，窥其本质：

好了，我们利用上面两个正余弦函数的虚数形式，再化简傅里叶级数：

$f(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty}(a_n\cos \frac{n\pi t}{l}+b_n\sin \frac{n\pi t}{l})$

$=\frac{a_0}{2}+\sum_{n=1}^{\infty}[\frac{a_n}{2}(e^{i\frac{n\pi t}{l}}+e^{-i\frac{n\pi t}{l}})$

$-\frac{ib_n}{2}(e^{i\frac{n\pi t}{l}}-e^{-i\frac{n\pi t}{l}})]$

$=\frac{a_0}{2}+\sum_{n=1}^{\infty}[\frac{a_n-ib_n}{2}e^{i\frac{n\pi t}{l}}+\frac{a_n+ib_n}{2}e^{-i\frac{n\pi t}{l}}]$

$=c_0+\sum_{n=1}^{\infty}(c_ne^{i\frac{n\pi t}{l}}+c_{-n}e^{-i\frac{n\pi t}{l}})$

$=\sum_{n=-\infty}^{+\infty}c_ne^{i\frac{n\pi t}{l}}$

最终，我们得到傅里叶级数的复数形式为： $\bbox[yellow] { f(t)=\sum_{n=-\infty}^{+\infty}c_ne^{i\frac{n\pi t}{l}} }$ 嗯，你没有看错，刚才那么一大坨的公式居然被我们用欧拉公式给简化到了如此简洁的程度。根据前面的结论： $e^{it}$ 的实际几何意义就是复向量绕圆周做 $t$ 速度的旋转运动，前面的 $c_n$ 系数只是决定了其旋转半径以及初始旋转角度。

因此不难看出，傅里叶级数的本质几何意义就是无数个不同半径的复向量以不同速度、不同方向进行旋转然后进行向量叠加的结果。本文刚开始的三个动画我就是通过一千个复向量做圆周运动进行累加，只需要把最后一个向量的坐标记录下来即可。

$\bbox[yellow] { c_ne^{i\frac{n\pi t}{l}} \begin{cases} c_n, & \text{决定旋转半径及初始角度} \\[2ex] e^{i\frac{n\pi t}{l}}, & \text{决定旋转速度} \end{cases} }$

这种叠加的结果让人觉得不可思议，甚至说恐怖也丝毫不为过，因为只要规定合适的参数，它几乎可以拟合出来我们任何想要的曲线，更神奇的是，有些函数明明带有明显的跳跃不连续点，傅里叶级数居然也能拟合出来，前面我说了方波函数可通过傅里叶级数来获得，即便它带有跳跃点也无所谓：

下面是一个方波函数：我们先用一个圆旋，并记录其某一点纵坐标路径，没错，正弦函数就是这样来的：

咋一看，有那么几分相似。我们再增加到三个圆：

再看，貌似有几分相像了，但还是有差距。

再将圆数量增加到50个：

可见，此时方波函数几乎已经与我们的路径拟合在一块了，如果不是你亲眼看到，我相信你很难相信方波函数居然可以通过正余弦函数叠加来近似得到，而更恐怖的是方波的间断点竟然也能拟合出来。并且，随着圆数量的增加，我们拟合的结果会越来越准确。实在是太恐怖了，而且通过合适的圆数量它几乎能拟合出任何你想要的的图像路径。

不过，要说明的是，上面每个圆的旋转方向与速度都是计算出来的，具体是如何计算的呢，请继续往下看。

五、傅里叶动画的制作

有了上面的知识，我们做文章开头那几个动画就有思路了。我只是通过计算，把上面一连串的 $c_n$ 给计算出来，计算的结果就是他们的方向与速度，还有长度，并让他们各自按照自己的圆心旋转即可。但是，这里有两个问题：

（1）函数 $f(t)$ 究竟是如何来的？

（2）每个复向量前面的系数 $c_n$ 又是如何计算出来的？

下面，我们仔细剖析这两个问题。

5.1，获取 $f(t)$ 方法：

首先来说 $f(t)$ 是如何来的。注意这里的 $f(t)$ 并不是一组连续的实数，而是一组二维向量的坐标值。为了画出前面那三幅优美的动画，我们需要对时间 $t$ 与像素的坐标值 $(a,b)$ 做一一映射，为了方便，我们用复数的形式表示每一个像素点的坐标，也即任意给定一个时间 $t$ ，都能找到唯一的一个复数坐标与其对应： $t\to (a+bi)=re^{i\theta}$ 那么究竟如何找出这些像素点坐标的集合呢？在这里，我使用了opencv来读取一张图像轮廓的坐标数据集。不过在这里需要实现配置一下opencv环境，具体方法请参考我另外一篇文章：手把手教你安装OpenCV与配置环境。对图像边缘轮廓的提取大致可分为下面五个步骤：

1）对图像进行灰度转化；

2）将转化后的灰度图像进行高斯模糊处理，目的是让轮廓更加平滑；

3）再将处理后的图像进行二值化处理，实现黑白分明的轮廓；

4）利用Canny()算子提取轮廓；

5）使用findContours()接口提取图像轮廓坐标数据集合。

根据这个方法，我们几乎可以提取任何图片的轮廓：

c++实现代码如下：

//函数功能：对图像进行轮廓坐标提取//返回值：int//作者：@刘亚曦#include "stdafx.h"#include #includeusing namespace cv;using namespace std;int main(){  Mat src = imread("自己图片的路径");  Mat grayImage;  cvtColor(src, grayImage, CV_BGR2GRAY);  //灰度处理  GaussianBlur(grayImage, grayImage, Size(3, 3), 0, 0);  //高斯模糊处理  threshold(grayImage, grayImage, 128, 255, CV_THRESH_BINARY);  //二值化处理  Mat cannyImage;  Canny(grayImage, cannyImage, 128, 255, 3);    //提取边缘算子  vector<vector> contours;  vector hierarchy;  //contours为输出的轮廓数据集合  findContours(cannyImage, contours, hierarchy, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE, Point(0, 0));  return 0;}

这行代码要根据你图片的实际情况来稍微调整几个参数，例如灰度阈值、高斯核半径、canny算子阈值等等。这几个参数在一定程度上决定了你提取图片轮廓的质量。不过在这里这几个参数暂时不做重点讲解，等我有时间会在另外一篇文章中专门讲解一下几种边缘检测算法与参数的意义。大家按照我上面代码里面的参数设置即可，基本上能适应大多数图像了。

这样，我们就能得到我想所要的图片轮廓坐标集合，上面代码中，contours参数就是储存输出值。而且这里有趣的是，它提取的图片轮廓像素坐标是有一定顺序的，一般是以逆时针排列的，这对我们后面处理积分运算极大方便。这就是 $f(t)$ 的来历。不过我们再使用过程中，为了方便，还是要把所有的实数坐标转化为复数坐标。

5.2，计算 $c_n$ 系数：

那么究竟如何计算 $c_n$ 的实际值呢？我们往下看：

$f(t)=...+c_{-2}e^{i\frac{-2 \pi t}{l}}+c_{-1}e^{i\frac{-1 \pi t}{l}}+$

$c_0+c_{1}e^{i\frac{1 \pi t}{l}}+c_{2}e^{i\frac{2 \pi t}{l}}...\qquad (2)$ 首先我们对其两边取区间为单个 $f(t)$ 周期也即 $[0,2l]$ 的积分：

$\int_{0}^{2l}f(t)dt=...\int_{0}^{2l}c_{-2}e^{i\frac{-2 \pi t}{l}}+\int_{0}^{2l}c_{-1}e^{i\frac{-1 \pi t}{l}}$

$+\int_{0}^{2l}c_0+c_{1}e^{i\frac{1 \pi t}{l}}+\int_{0}^{2l}c_{2}e^{i\frac{2 \pi t}{l}}...\qquad (3)$ 根据我们上面(1)式的结论，在一个整周期内，复向量的积分恒为零，也就意味着上式右边出了 $c_0$ 那一选项，其余都为零：

$...\int_{0}^{2l}c_{-2}e^{i\frac{-2 \pi t}{l}}=0;\int_{0}^{2l}c_{-1}e^{i\frac{-1 \pi t}{l}} = 0$

$\int_{0}^{2l}c_{1}e^{i\frac{1 \pi t}{l}}=0;\int_{0}^{2l}c_{2}e^{i\frac{2 \pi t}{l}}=0...$ 所以式(2又可化简为：

$\int_{0}^{2l}f(t)dt=\int_{0}^{2l}c_0$ 这样我们便能得出： $c_0=\frac{1}{2l}\int_{0}^{2l}f(t)dt$ 最后，我们求除了 $c_0$ 以外的其他参数系数。其实方法也一样，为了求某一项的系数，我们可以把本项乘以对应指数的负数来消掉本项的指数，然后两边再次积分，例如为了求 $c_{2}e^{i\frac{2 \pi t}{l}}$ 本项的 $c_2$ ，我们首先对(2)式两边分别乘以 $e^{-i\frac{2 \pi t}{l}}$ ：

$f(t)e^{-i\frac{2 \pi t}{l}}=...+c_{-2}e^{i\frac{-2 \pi t}{l}}e^{-i\frac{2 \pi t}{l}}+...$

$+c_0e^{-i\frac{2 \pi t}{l}}+c_{2}e^{i\frac{2 \pi t}{l}}e^{-i\frac{2 \pi t}{l}}...$

两边再次积分，注意，再根据上面的结论，上式除了 $c_2$ 项外，其余全部为0，因此我们得到： $\int_{0}^{2l}f(t)e^{-i\frac{2 \pi t}{l}}=\int_{0}^{2l}c_{2}$ 求得： $c_2=\frac{1}{2l}\int_{0}^{2l}f(t)e^{-i\frac{2 \pi t}{l}}dt$ 进一步，我们得到： $\bbox[yellow]{ c_n=\frac{1}{2l}\int_{0}^{2l}f(t)e^{-i\frac{n \pi t}{l}}dt }$ 至此，我们已经将所有的系数求出来了。剩下的就是通过代码讲这些计算方法实现，并通过动画模拟出来。下面我们开始编写代码，为了动画能够在浏览器上直接运行，这次我采用JS代码实现。

5.3，动画引擎代码的编写与实现：

其实这个框架非常简单，我们首先根据上面推导出来的公式，由于 $c_n$ 决定了每个复向量的半径大小与初始角度大小，所以我们必须先计算出来它。而剩下的幂指数部分只是决定了其旋转方向与速度而已，并且指数部分只有一个变量 $n$ ，剩下的都是已知变量，很好计算。等全部计算出来后，将所有圆绘制出来，并根据实际时间 $t$ 来计算他们的实际位置，重新计算并刷新屏幕，就形成了连续的动画。

代码部分：

在计算过程中，我们要不断计算复数的相加，复数的相乘等方法，所以我们先封装几个常用的工具函数方便我们后面调用，我们返回值都为一个长度为2的数组，第一个代表实部值，第二个值代表虚部。实现代码如下：

//返回类型：[]//作者：@刘亚曦//功能：指数转化为复数function exp(theta) {    return [Math.cos(theta), Math.sin(theta)];}//功能：两个复数相乘function mul(za, zb) {    return [za[0] * zb[0] - za[1] * zb[1], za[0] * zb[1] + za[1] * zb[0]];}//功能：两个复数相加function add(za, zb) {    return [za[0] + zb[0], za[1] + zb[1]];}

我们再定义一个数组 $K$ ,用它来储存 $c_n$ 下标 $n$ 的值。方便我们后面做积分运算。在这里要特别注意，数组内的值严格按照[0,-1,1,-2,2...]进行排列，并且不可交换。代码如下：

//功能：储存下标值//作者：@刘亚曦//返回值：[]var circleCount = 1000; //圆的数量var K = [];for (var i = 0; i < circleCount; i++) {    K[i] = (1 + i >> 1) * (i & 1 ? -1 : 1);}

然后我们开始计算 $c_n$ 的值，根据我们前面推导出来的运算公式： $c_n=\frac{1}{2l}\int_{0}^{2l}f(t)e^{-i\frac{n \pi t}{l}}dt$ 可实际上我们一张图像上的像素坐标并不真的“连续”，所以说在这里并不存在实际意义的积分运算。但是我们可以通过穷举法，每次让 $\Delta t$ 偏移一个最小的单位，不断与 $f(t)$ 进行相乘，最后全部累加做和，就能模拟出积分运算。这里我们取最小间隔时间刷新单位，让 $\Delta t=1$ 。由此，我们可以用下面公式来模拟积分运算： $\int_{0}^{2l}f(t)e^{-it}dt=[f(0)e^0\Delta t+f(1)e^1\Delta t+...+f(2l)e^{2l}\Delta t]$ 这里的 $2l$ 就是我们路径数组的长度，也就是在上面我们计算 $f(t)$ 时输出的数组长度。

我们再定义一个getCn()函数，来计算 $c_n$ 的值，注意在这里，上面我们给出的函数周期为 $2l$ ，显然我们这里整个轮廓坐标要形成一个闭环，如果我们每次让 $\Delta t$ 偏移一个单位，这样我们实际的周期应该为我们轮廓数据集数组的长度。轮廓坐标保存在path二维数组里面，因此这里的 $2l=$ path.length。

不过要注意，在计算 $e^{-i\frac{n \pi t}{l}}$ 的时候，我们要调用前面封装的exp函数，将指数形式转化为复数形式，方便我们运算，代码如下：

//函数功能：穷举法计算复向量系数c_n//返回值：[]//作者：@刘亚曦function getCn() {    var z = [0, 0];    var N = path.length;        //路径坐标数组    for (var j = 0; j < K.length; j++) {        for (var i = 0; i < N; i++) {            var za = [path[i][0], path[i][1]];          //f(t)            var zb = exp(K[j] * i * 2 * -Math.PI / N);  //调用指数转复数函数            z = add(z, mul(za, zb));                    //调用相乘函数        }        z[0] = z[0] / N;        z[1] = z[1] / N;        Cn[j] = [z[0], z[1]];    }}

至此，我们就将全部的所需要的 $c_n$ 参数计算出来了。

最后，我们定义三个绘画函数，用来绘制我们的圆、路径、以及连接圆心的直线，在主函数里面循环刷新，即可形成动画。其实这三个函数非常简单，前面我们的 $c_n$ 数组参数已经计算出来了，这个参数决定了每个圆的半径。

我们还需要每个圆的另外两个参数才能实现绘制，一个是其旋转速度，另一个是其中心坐标。我上面说了，其速度是有公式的指数部分 $e^{-i\frac{n \pi t}{l}}$ 决定的，而这里的 $n$ 我们已经知道了([0,-1,1,-2,2...])， $l$ 也知道了，就剩下一个 $t$ 了，这个 $t$ 就是我们的时间“尺度”，上面我们说了，让 $\Delta t$ 每次偏移一个单位，所以这里的 $t$ 我们每次计算的时候都要自增1。因此，我们早在主循环函数再定义一个time的参数，每次循环后实现time=time+1操作，不过如果你感觉整体绘画速度偏慢的话，可以适当增加time参数的自增速度，例如让其自增2，可以提高整体速度。这样随着time增加，我们的所有圆的瞬时圆心就可以计算出来，每次循环的时候无非重新刷新屏幕再重新把所有的圆绘制出来而已，形成肉眼看到的动画。实现代码如下：

//功能：画圆函数//作者：@刘亚曦//返回值：voidfunction DrawCircles() {    let p = [center_x, center_y];    var a = 2 * Math.PI * time / path.length;    for (var i = 0; i < Cn.length; i++) {        context.beginPath();        var r = Math.hypot(Cn[i][0], Cn[i][1]);        context.arc(p[0], p[1], r, 0, 2 * Math.PI);        context.lineWidth = 1.0;        context.strokeStyle = "rgba(255,128,32,0.7)";        if (i > 0) {            context.stroke();   //第一个圆不画        }        p = add(p, mul(Cn[i], exp(a * K[i])));    }}

再定义绘制路径函数，这个几乎跟上面画圆函数没有区别，因为我们已经计算出来了圆心的实时坐标，直接将他们连起来即可：

//绘制连接圆心函数//作者：@刘亚曦//返回值：voidfunction DrawLines() {    context.beginPath();    let p = [center_x, center_y];    var a = 2 * Math.PI * time / path.length;    for (var i = 0; i < Cn.length; i++) {        if (i == 1) {            context.moveTo(p[0], p[1]);     //第一个线不画        }        p = add(p, mul(Cn[i], exp(a * K[i])))        context.lineTo(p[0], p[1]);    }    context.lineWidth = 1;    context.strokeStyle = "rgba(255,255,255,0.6)";    context.stroke();}

最后，来绘制路径函数。这个函数比较特殊，我重点说明一下。这里的路径指的是最后一个圆上面一个点的实时坐标，具体是哪一点呢，这是由最后一个圆对应的 $c_n$ 系数决定的，因为 $c_n$ 决定了每个圆的半径与初始角度。这样，我们在计算路径的时候，只需根据当前的时间time参数，遍历每个圆上面计算的点，将其的首尾坐标进行相加即可。比如第一个圆圆心计算出来了，再根据它上面一点作为圆心再计算下一个圆，以此类推，即可算出最后一个圆的路径。

这里还要说明一点，为了形成完整的轮廓，我们要把最后一点的坐标记录下来保存在数组里面，但是为了防止内存过度占用，我们要限制数组的长度，也就是只绘制以当前时间为基准的后面一定数量的坐标。我使用了&运算来实现这样的操作，这样我们数组当前保存的坐标集合就是当前最新的数组集合。实现代码如下：

//功能：绘制路径函数//作者：@刘亚曦//返回值：voidfunction DrawPath() {    let p = [center_x, center_y];    var a = 2 * Math.PI * time / path.length;    for (var i = 0; i < Cn.length; i++) {        p = add(p, mul(Cn[i], exp(a * K[i])));    }    var x = p[0];    var y = p[1];    valuePointer++;    values_x[valuePointer & pointCount] = x;    values_y[valuePointer & pointCount] = y;    context.beginPath();    context.strokeStyle = "rgba(0,255,0,1)";    context.moveTo(x, y);    for (var i = 1; i <= pointCount; ++i) {        context.lineTo(values_x[(valuePointer - i) & pointCount], values_y[(valuePointer - i) & pointCount]);    }    context.stroke();}

注意，上面的代码我都跳过了第一个圆的绘制，读者看到的圆是加上是少了一个的。为什么要这样呢？这是因为 $c_0$ 系数的原因：在我们计算的 $c_n$ 数组中，只有 $c_0$ 比较特殊，它的幂函数的指数部分为0，意味着虽然它也是一个圆，但它却不做任何滚动。为了实现更好的观赏性，我就把它取消了。但读者可以把我那行跳过第一个圆的代码注释掉，下去可以自己试试完整圆的效果。

上面几个函数都封装好了，剩下的就更简单了，我们在主循环函数里面不断去调用它实时刷新屏幕即可，代码如下：

//循环刷新函数//作者@刘亚曦(function frame() {    context.clearRect(0, 0, canvas.width, canvas.height);    context.fillStyle = "#000000";    context.fillRect(0, 0, canvas.width, canvas.height);    DrawCircles();    DrawPath();    DrawLines();    time = time + 2;    window.requestAnimationFrame(frame);})();