傅里叶级数实现本轮画图¶

将一张 SVG 图转化为关于时间的周期函数 \(f: \mathbb{R} \to \mathbb{C}\)，再展开成复数形式的傅里叶级数。级数的每一项都是一个本轮（Epicycle），让它们随时间运动，绘制出的轨迹和原 SVG 图一样。

将 SVG 转化为函数¶

对于一条二维曲线，其上每一点都可以用一个复数表示。

考虑到一张 SVG 图里可能有多条曲线，为了写代码方便，可以定义这样一个函数，对于图中的第 \(i\) 条曲线（\(i=1,2,\cdots\)）：

\(f(i-1)\) 为该曲线的起点。
\(f(x)\)，\(x \in (i-1,i)\) 表示

\[ \dfrac{\text{从} f(i-1) \text{沿该曲线到} f(x) \text{的距离}}{\text{该曲线的长度}} = x-(i-1) \]

于是，假设 SVG 图中有 \(n\) 条曲线，\(f(x)\) 的定义域为 \([0,n)\)。然后对 \(f(x)\) 做周期延拓，就得到了一个周期为 \(n\) 的函数 \(f: \mathbb{R} \to \mathbb{C}\)。

计算傅里叶级数的系数¶

原始的计算公式是一个积分

\[ c_n = \frac{1}{T} \int_{0}^{T} f(t) \exp \left(-\frac{2 \pi nt}{T}i \right) dt \]

这里用最简单的 Riemann Sum 来近似。假设把区间分成 \(N\) 份，改写后化简得

\[ c_n = \frac{1}{N} \sum_{j=1}^{N} f(t_j) \exp \left(-\frac{2 \pi nt_j}{T}i \right) \]

注意到，这个形式可以改写成两个向量的内积。令

\[ g(n, t) = \exp \left(-\frac{2 \pi nt}{T}i \right) \]

则有

\[ c_n = \frac{1}{N} \begin{bmatrix} g(n,t_1) & g(n,t_2) & \cdots & g(n,t_N) \end{bmatrix} \begin{bmatrix} f(t_1)\\ f(t_2)\\ \vdots\\ f(t_N) \end{bmatrix} \]

那么，所有系数构成的列向量就表示为

\[ \begin{bmatrix} \vdots\\ c_0\\ c_1\\ \vdots\\ c_n\\ \vdots\\ \end{bmatrix} = \dfrac{1}{N} \begin{bmatrix} \vdots & \vdots & & \vdots\\ g(0,t_1) & g(0,t_2) & \cdots & g(0,t_N)\\ g(1,t_1) & g(1,t_2) & \cdots & g(1,t_N)\\ \vdots & \vdots & \ddots & \vdots\\ g(n,t_1) & g(n,t_2) & \cdots & g(n,t_N)\\ \vdots & \vdots & & \vdots\\ \end{bmatrix} \begin{bmatrix} f(t_1)\\ f(t_2)\\ \vdots\\ f(t_N) \end{bmatrix} \]

写成这种形式后，可以很方便地用 numpy 写代码。

还有一个好处是减少了重复计算。列向量 \((f(t_1),f(t_2),\cdots,f(t_N))^T\) 是每次计算系数时的公共部分，只要算一次就行。如果只看一开始的式子，不一定能看出来。

绘图原理¶

根据傅里叶级数的公式

\[ f(x)=\sum_{n=-\infty}^{+\infty} c_n \exp(\frac{2\pi nx}{T} i) \]

让 \(x\) 随时间变化，级数每一项 \(c_n \exp(\dfrac{2\pi nx}{T} i)\) 的轨迹都是一个圆。\(\dfrac{2\pi n}{T}\) 是角速度。\(c_n\) 是一个定值，决定了初相和圆的半径。

把所有项加起来，整体的轨迹就是原 SVG 图。实际画图时只取一部分项做近似。

代码¶

绘图使用 3b1b 的 Manim 实现。配置字段写在了 TestScene2 的 CONFIG 中。提示：

如果觉得画得不准，可以增加级数的项数，或者提高积分的采样数。
如果画出的图不光滑，可以增加动画的时长。
实时渲染可能比较卡，保存成视频就流畅了。

from manimlib import *
import numpy as np


class FourierTerm(object):
    """ 表示傅里叶级数中的一项，即 c * exp(2 * PI * 1j * n * x / t) """

    def __init__(self, c: complex, n: int, t: float) -> None:
        self.c = c
        self.n = n
        self.t = t

    def evalAsRealVector(self, x: float) -> np.ndarray:
        """ 给定 x ，计算该项的值，并将结果用一个三维的实向量表示 """
        v = self.c * np.exp(2 * np.pi * 1j * self.n * x / self.t)
        return np.array([v.real, v.imag, 0]) # 把复数转成对应的实向量

    @property
    def magnitude(self) -> float:
        """ 该项的模长 """

        # 该项等于 c * exp(2 * PI * 1j * n * x / t)
        # exp(2 * PI * 1j * n * x / t) 的模长恒为 1
        # 所以整体的模长就等于 c 的模长
        return abs(self.c)


class SVGFunction(object):
    def __init__(self, svg: SVGMobject) -> None:
        # 获取 svg 图中的所有曲线
        self.curves = svg.family_members_with_points()

        # 函数的周期，规定为曲线的数量
        self.period = len(self.curves)

    def sample(self, x: float) -> complex:
        # 根据 x 找到对应的曲线，然后按比例采样
        p, i = np.modf(x % self.period)
        point = self.curves[int(i)].point_from_proportion(p)
        return point[0] + point[1] * 1j # 用复数表示这个点

    def expand(self, sampleCount: int, termRange) -> list[FourierTerm]:
        g = lambda n, t: np.exp(-2 * np.pi * n * t * 1j / self.period)

        samples = np.linspace(0, self.period, sampleCount + 1)[:-1] # 取指定数量的等间隔点
        matF = np.array([self.sample(x) for x in samples])
        matG = np.array([g(n, samples) for n in termRange])
        matC = np.matmul(matG, matF) / sampleCount

        results = [FourierTerm(c, n, self.period) for c, n in zip(matC, termRange)]
        results.sort(key=lambda t: abs(t.n))
        return results


class TestScene2(Scene):
    CONFIG = {
        'svgPath': 'test.svg',         # 图片路径
        'svgScale': 2.5,               # 图片缩放系数
        'termRange': range(-150, 151), # 选取的项
        'integralSampleCount': 10000,  # 积分采样数

        'animateTime': 80,             # 动画时长
        'repeatCount': 2,              # 重复次数
        'lineStrokeWidth': 3,
        'colorGradient': [RED, ORANGE, YELLOW, GREEN, BLUE, PURPLE]
    }

    def construct(self) -> None:
        self.camera.background_rgba = [0, 0, 0, 1] # bg

        debugLabels = VGroup(
            Tex('Epicycle Count: {}'.format(len(self.termRange)), font_size=22),
            Tex('Integral Sample Count: {}'.format(self.integralSampleCount), font_size=22)
        ).arrange(RIGHT).move_to(TOP + RIGHT_SIDE - np.array([0.15, 0.15, 0]), aligned_edge=UR)
        self.add(debugLabels)

        svg = SVGMobject(self.svgPath, stroke_color=WHITE)
        svg.scale(self.svgScale)
        self.add(svg)
        self.doAnimation(svg)

    def doAnimation(self, svg: SVGMobject):
        func = SVGFunction(svg)
        fourierTerms = func.expand(int(self.integralSampleCount), self.termRange)

        samplePoints = []
        arrows = VGroup()
        epicycles = VGroup()

        # 添加箭头和本轮
        fouriers = VGroup(arrows, epicycles)
        for term in fourierTerms:
            arrows.add(Arrow(buff=0, stroke_width=2, color=self.getTermColor(term)))
            epicycles.add(Circle(radius=term.magnitude, buff=0, stroke_width=1, stroke_opacity=0.4, color=self.getTermColor(term)))
        self.add(fouriers)

        self.playTime = 0
        def updateArrowsAndEpicycles(dt):
            parentPos = np.zeros(3)

            self.playTime += dt
            progress = min(1, self.playTime / self.animateTime)
            x = interpolate(0, func.period * self.repeatCount, progress)

            for arr, epic, term in zip(arrows, epicycles, fourierTerms):
                start = parentPos
                end = start + term.evalAsRealVector(x)
                arr.put_start_and_end_on(start, end)
                epic.move_to(start)
                parentPos = end

            samplePoints.append(parentPos)

            if len(samplePoints) > 1:
                self.add(Line(
                    start=samplePoints[-2],
                    end=samplePoints[-1],
                    stroke_width=self.lineStrokeWidth,
                    buff=0,
                    color=self.getTermColor(fourierTerms[-1])
                ))

        # 随时间更新
        fouriers.add_updater(lambda m, dt: updateArrowsAndEpicycles(dt))
        self.wait(self.animateTime + 1)

    def getTermColor(self, term: FourierTerm):
        return self.colorGradient[abs(term.n) % len(self.colorGradient)]

if __name__ == '__main__':
    moduleName = os.path.basename(__file__)
    os.system(f'manimgl {moduleName} TestScene2')