傅里叶级数实现本轮画图¶
大一下，高数拓展作业。
将一张 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')