傅里叶分析：用圆画羊驼¶

傅里叶变换是一种将输入信号从时域变换到频域的方法。它的应用范围极其广泛，比如音频工程师可以用傅里叶变换挑出歌曲中某个不需要的频率，再用逆傅里叶变换得到去掉该频率后的声音。在数字图像处理中，我们也能借助傅里叶变换做图像滤波，比如高斯模糊，以及用 JPEG 格式压缩图像。最后，正如本文将要展开的：画画！我说的画画是什么意思？思路是这样的：我们可以取一条路径，它代表某个东西，比如一条鱼、π 符号，或者几乎任何我们能画出来的图形——把铅笔落在纸上，一笔画完，中途不抬笔。有了这条路径，我们就能把一堆以不同频率绕圆旋转的向量首尾相接，最后一个向量的末端会描出我们最初的草图，或者至少是与之非常接近的东西。下面举几个例子说明我在讲什么：

图 1 用圆画一条直线¶

所以你可以用一种与直线尽可能相反的图形画出一条直线。修改这些圆就能改变画出的形状。比如，减小外圆的半径、提高它的旋转速度，我们就能画出一个正方形。

图 2 一个正方形¶

那么，如果再加上第三个圆会怎样？它能让我们的曲线更精细。用三个圆，我们可以画出一条看起来像鱼的曲线。

图 3 一条鱼。¶

想看仅用彼此相连的圆能做到多极端的事，这里有一头羊驼，用了 $1024$ 个圆来绘制，每个圆都有不同的频率、半径和起始角度。一般来说，（几乎）每条闭合曲线都能用圆画出来，这些圆称为本轮。

图 4 一只羊驼？！¶

一般来说，我们添加的圆越多，能画出的图形就越复杂。同样，对于同一幅画，添加的圆越多，它就越接近原始的“输入”图形。要理解这段动画是怎么来的，我们需要深入傅里叶变换（以及傅里叶级数）背后的一些数学原理和直觉。如前所述，我们可以在一个圆上调整三个条件来画出不同的曲线：它旋转的速率（频率）、圆的大小（半径），以及它的起始角度（相位）。这三个条件都可以用单个复数来表示。

下面这个例子你可以自己动手试试，直观感受一下增加圆对最终图形的影响。

傅里叶变换¶

但首先，我大概应该先简单介绍一下傅里叶变换是什么。它是一种从时域到频域的变换。这是什么意思？如果我们有一个函数 $\sin (2 π \times 3 t)$ ，那么对这个 $\sin (2 π \times 3 t)$ 函数做傅里叶变换，结果就是在频率 $\pm 3 hz$ 处的一个尖峰。

图 5 $\sin (2 π \times 3 t)$ 的图像¶

图 6 $\sin (2 π \times 3 t)$ 在频域中的幅值。¶

不出所料，在频率 $3 hz$ 和 $- 3 hz$ 处出现了尖峰。本质上，我们要做的是找出我们想要处理的对象中存在哪些频率，然后把这些正弦波组合成对原始对象的一个近似。在接下来的几节里，我会尝试展示傅里叶变换是如何工作的。

复数¶

这一节会快速复习一下什么是复数，如果你已经熟悉，可以跳过。一个复数由两部分组成：“实”部和“虚”部。它把我们早已熟悉的数拿过来，本质上又加了一个轴来表示它们。复数写作 $3 + 4 i$ ，其中 $3$ 是实部， $4 i$ 是虚部。这两部分相加，构成一个复数。复数在这里还有一个有用的性质：它们可以很方便地看作二维向量，其中实部是 $x$ ，虚部是 $y$ ，对应二维欧几里得空间中的标准点 $(x, y)$ 。由于三角学与复指数之间的联系，这种表示在这里非常方便。把两者联系起来的公式称为欧拉恒等式。本文后面的一节会进一步讨论它。

线性代数¶

线性代数是理解如何利用向量在圆上旋转来绘制图形的重要概念，具体来说，就是如何用其他向量表示某个向量。在标准的二维 $(x, y)$ 欧几里得空间中，所有向量都可以用 $𝐱$ 和 $𝐲$ 单位向量表示。这两个向量通常记作 $\hat{𝐢}$ 和 $\hat{𝐣}$ （读作“i-hat”和“j-hat”），其中 $\hat{𝐢} = (1, 0)$ （ $x$ 轴）， $\hat{𝐣} = (0, 1)$ （ $y$ 轴）。要做到这一点，我们有两种可用的运算：标量乘法和加法。例如，向量 $\vec{a} = (5, - 8)$ 可以表示为

(1)¶

\begin{matrix} \begin{matrix} \begin{matrix} \vec{a} & = 5 \cdot \hat{𝐢} - 8 \cdot \hat{𝐣} \\ = 5 \cdot (1, 0) + - 8 \cdot (0, 1) \\ = (5, 0) + (0, - 8) \\ = (5, - 8) \end{matrix} \end{matrix} \end{matrix}

这种用标量乘法与向量加法组合出另一个向量的方式，称为一组向量的线性组合。在这个例子里，这组向量就是单位向量 $\hat{𝐢}$ 和 $\hat{𝐣}$ 。

图 7 两个向量 $\vec{u}$ 和 $\vec{v}$ ，夹角为 $90^{\circ}$ ，长度均为 $1$ ¶

但 $\hat{𝐢}$ 和 $\hat{𝐣}$ 并不是唯一能用来表示其他向量的向量。事实上，任意两个向量，只要彼此不成倍数关系，都可以通过缩放再相加得到任意其他向量。例如 $(1, - 1)$ 和 $(- 2, 2)$ 就不行，因为 $(- 2,2) = - 2 \times (1, - 1)$ 。在上图中，向量 $\vec{u} = (\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}})$ 和 $\vec{v} = (- \frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}})$ 长度为 1，且互相正交（夹角为 $90^{\circ}$ ）。这两个向量谁都无法通过缩放变成对方，因此它们线性无关，可以作为整个二维空间的基向量。

点积（有时也叫“内积”）用来确定：用选定的向量表示某个向量时，各向量需要缩放的倍数。要用向量 $\vec{𝐰} = (0.7, 1.2)$ （见上图）表示成向量 $\vec{u}$ 和 $\vec{v}$ 的线性组合 $\vec{w} = c_{1} \cdot \vec{u} + c_{2} \cdot \vec{v}$ ，我们需要求出这两个向量各自要乘上的数。具体来说，点积给出的正是这两个常数系数 $c_{1}$ 和 $c_{2}$ 。换个角度理解，点积回答的是“ $\vec{u}$ 里有多少在 $\vec{w}$ 中？”这个问题。等后面讲到 Fourier 变换背后的数学时，这个直觉会很重要。总之，把 $\vec{w}$ 与 $\vec{u}$ 做点积就能算出 $c_{1}$ ，而且顺序无关紧要（例如 $⟨ \vec{u}, \vec{w} ⟩ = ⟨ \vec{w}, \vec{u} ⟩$ ）。

(2)¶

\begin{matrix} \begin{matrix} \begin{matrix} c_{1} & = ⟨ \vec{u}, \vec{w} ⟩ \\ = (u_{x} \times w_{x}) + (u_{y} \times w_{y}) \\ = (\frac{1}{\sqrt{2}} \times 0.7) + (\frac{1}{\sqrt{2}} \times 1.2) \\ = 1.3435 \end{matrix} \end{matrix} \end{matrix}

也就是说， $\vec{u}$ 要乘以 $1.3435$ ，即被拉伸到原长度的 $1.3435$ 倍。于是 $\vec{w}$ 中属于 $\vec{u}$ 的分量（ $c_{1}$ ）就是 $1.3435$ 。

(3)¶

\begin{matrix} \begin{matrix} \begin{matrix} c_{2} & = ⟨ \vec{v}, \vec{w} ⟩ \\ = (v_{x} \times w_{x}) + (v_{y} \times w_{y}) \\ = (- \frac{1}{\sqrt{2}} \times 0.7) + (\frac{1}{\sqrt{2}} \times 1.2) \\ = 0.353553 \end{matrix} \end{matrix} \end{matrix}

同样地，我们需要将 $\vec{v}$ 缩放 $0.353553$ 倍，也就是说它会被压缩到原尺寸的 $0.353553$ 倍。因此 $\vec{w}$ 中 $\vec{v}$ 的分量（ $c_{2}$ ）为 $0.353553$ 。

最终得到的线性组合为

(4)¶

\begin{matrix} \begin{matrix} \begin{matrix} \vec{w} & = 1.3435 \cdot \vec{u} + 0.353553 \cdot \vec{v} \\ = (0.95, 0.95) + (- 0.25, 0.25) \\ = (0.7, 1.2) \end{matrix} \end{matrix} \end{matrix}

这确实等于我们最初的向量 $\vec{w}$ 。

图 8 $\vec{u}$ 与 $\vec{v}$ 的线性组合得到 $\vec{w}$ ¶

任意两个向量的点积定义为

(5)¶

\begin{matrix} \begin{matrix} \begin{matrix} ⟨ \vec{v_{1}}, \vec{v_{2}} ⟩ & = {v_{1}}_{x} \times {v_{2}}_{x} + {v_{1}}_{y} \times {v_{2}}_{y} \\ = ‖ \vec{v_{1}} ‖ \times ‖ \vec{v_{2}} ‖ \times \cos (θ) \end{matrix} \end{matrix} \end{matrix}

其中 $θ$ 是 $v_{1}$ 与 $v_{2}$ 之间的夹角， $‖ \vec{v} ‖$ 是向量的模。模也称为范数，定义为 $‖ \vec{v} ‖ = \sqrt{v_{x}^{2} + v_{y}^{2}}$ 。这个式子大概看着眼熟，因为它就是我们中小学数学课上都学过的勾股定理。

图 9 直角三角形¶

点积还有一些重要的性质，正好对我们有用。首先，任意向量与自身的点积等于该向量模的平方（ $\vec{v_{1}} \cdot \vec{v_{1}} = ‖ \vec{v_{1}} ‖^{2}$ ）。例如，

(6)¶

\begin{matrix} ⟨ [5, 7], [5, 7] ⟩ = {\sqrt{5^{2} + 7^{2}}}^{2} = {\sqrt{74}}^{2} = 74 \end{matrix}

其次，两个互相垂直（相交成直角）的向量的点积为 $0$ 。这是因为 $\cos (\frac{π}{2}) = 0$ 。因此，当 $\vec{v_{1}}, \vec{v_{2}}$ 垂直时，可得

(7)¶

\begin{matrix} ⟨ \vec{v_{1}}, \vec{v_{2}} ⟩ = ‖ \vec{v_{1}} ‖ \times ‖ \vec{v_{2}} ‖ \times 0 = 0 \end{matrix}

按照我们之前对“ $\vec{v_{1}}$ 中有多少个 $\vec{v_{2}}$ ”的直观理解，这个结果为 $0$ 是合理的，因为这两个向量正交。

连续函数的内积¶

上一节讲的是两个向量的点积，但这个定义可以推广到连续函数。关键区别在于求和的方式不同。如前所述，二维点积的公式是

(8)¶

\begin{matrix} ⟨ \vec{a}, \vec{b} ⟩ = a_{x} \times b_{x} + a_{y} \times b_{y} \end{matrix}

更准确地说，是把两个向量中对应的项相加；此时向量长度有限（这里是 $2$ ），所以求和是离散的。这个公式在任意维度下还可以写成

(9)¶

\begin{matrix} ⟨ \vec{a}, \vec{b} ⟩ = \sum_{n = 1}^{N} a_{n} \times b_{n} \end{matrix}

其中 $N$ 是维度。这里的 sigma 求和是在离散区间上进行的。例如在三维空间中，求和区间是 $[1, 3]$ 。不过，这个求和只适用于 $ℝ^{n}$ 中的向量。长度为 $n$ 的复向量 $ℂ^{n}$ 的公式与之非常相似，只需对第二个向量取共轭。

(10)¶

\begin{matrix} ⟨ \vec{a}, \vec{b} ⟩ = \sum_{n = 1}^{N} a_{n} \times b_{n} \end{matrix}

需要注意的是，这里的复共轭使得该公式“对第二个参数反线性”，也就是说 $\vec{a} \cdot (\vec{b} + \vec{c}) = \vec{a} \cdot \vec{b} + \vec{a} \cdot \vec{c}$ 。但它在标量乘法下不是线性的，我们必须把标量的共轭提出来， $\vec{a} \cdot (b \times \vec{c}) = b \times (\vec{a} \cdot \vec{c})$ 。不过它对第一个参数是线性的，也就是说，从内积的第一个参数中提出标量倍数时，不需要取该标量的共轭。一般来说，它对取复共轭的那个参数反线性，对另一个参数线性。

我们需要一种能在连续区间上进行的求和，且求和区间可以无限小。这听起来正是积分的用武之地。

(11)¶

\begin{matrix} ⟨ f (t), g (t) ⟩ = \frac{1}{T} \int_{0}^{T} f (t) g (t) d t \end{matrix}

唯一的限制是函数 $f (t)$ 和 $g (t)$ 必须在区间 $[0, T]$ 上平方可积，也就是说，该函数绝对值平方的积分在 $[0, T]$ 上有限。

(12)¶

\begin{matrix} f (t) \in L^{2} (0, T) if \int_{0}^{T} | f (t) |^{2} d t is finite. \end{matrix}

注意 $L^{2}$ 是平方可积函数的集合， $L^{2} (0, T)$ 是区间 $[0, T]$ 上平方可积函数的集合。

为了让它更接近我们之前对 2D 空间中点积的定义，看起来更相似，我们可以把这个积分写成黎曼和（sigma 记号），如下所示

(13)¶

\begin{matrix} \frac{1}{T} \int_{0}^{T} f (t) g (t) d t = \lim_{N \to \infty} \frac{1}{T} \sum_{k = 0}^{N} f (\frac{T \times k}{N}) \times g (\frac{T \times k}{N}) \times \frac{T}{N} \end{matrix}

其中 $\frac{T}{N}$ 相当于这里的 $d t$ ，也就是每一步所乘的那个无穷小量。函数 $f$ 和 $g$ 的自变量，是积分上界 $T$ 乘以当前求和步数与极限趋近值之比，这样就能保证 $f$ 和 $g$ 的参数始终落在积分区间内。当步数较小（ $k = 0, 1, 2$ 等）时， $\frac{T \times k}{N}$ 会非常接近 $0$ ，因此自变量也非常接近 $0$ 。而当 $k$ 趋近于 $N$ 时， $\frac{T \times k}{N}$ 会接近 $1$ ，从而接近 $T$ 。

用复指数表示正弦与余弦¶

复指数让我们能方便地处理正弦和余弦。三角学与复平面之间存在着一个出人意料的联系，描述这一联系的公式叫做欧拉恒等式。下面是欧拉恒等式的公式，正如在复数一节中所承诺的那样。

(14)¶

\begin{matrix} e^{i x} = \cos (x) + i \times \sin (x) \end{matrix}

其中 $x$ 表示绕圆走多远。例如， $e^{\frac{π}{4} i}$ 绕圆走了 $\frac{π}{4}$ 弧度，最终落在 $\frac{\sqrt{2}}{2} + \frac{\sqrt{2}}{2} i$ 。此外， $e^{π i}$ 绕圆走了 $π$ 弧度，最终落在 $- 1 + 0 i$ ，这正是欧拉另一个著名的恒等式：

(15)¶

\begin{matrix} e^{π i} = - 1 \end{matrix}

图 10 $e^{i x}$ 与三角函数的联系¶

为什么会这样？把 $e$ 提升到虚数次幂，为什么会绕单位圆旋转？答案在物理里。若时刻 $t$ 的位置 $p$ 为

(16)¶

\begin{matrix} p (t) = e^{i t} \end{matrix}

那么该点的速度就是 $p (t)$ 的一阶导数。

(17)¶

\begin{matrix} v (t) = \frac{d p}{d t} = i \times e^{i t} \end{matrix}

这就够了！速度与位置形式相同，只是从实数变成了虚数。乘以 $i$ 相当于逆时针旋转 $90^{\circ}$ 。例如，若 $a \in ℝ$ ，那么 $i \times a$ 就是原来那个 $a$ 落到虚轴上的结果。但如果 $a$ 本来就是虚数呢？此时 $a = i \times k$ ，其中 $k$ 是某个实数。再乘以 $i$ ，得到 $i \times a = i^{2} \times k$ 。由 $i$ 的定义可知 $i^{2}$ 等于 $- 1$ ，所以 $i \times a = - k$ 。这只是一个实数的相反数，也就是旋转 $180^{\circ}$ 。既然乘以 $i$ 相当于旋转 $90^{\circ}$ ，那么速度函数始终垂直于当前位置的方向。这样画出的曲线是一个圆，这正是 $f (t) = e^{i t} where t \in [0, 2 π]$ 画出圆的原因。

另外注意

(18)¶

\begin{matrix} \frac{1}{2 π} \int_{0}^{2 π} \cos (t) \sin (t) d t = 0 \end{matrix}

因为 $\cos t$ 与 $\sin t$ 是正交函数，这使得 $\cos$ 和 $\sin$ 构成傅里叶变换的一组基。

图 11 演示 $\cos$ 与 $\sin$ 是正交函数。¶

本质上，我们做的是从二维复数空间 $(1,0), (0, 1) i$ 到以函数为基 $(. . ., e^{i \frac{- 4 π}{T} t}, e^{i \frac{- 2 π}{T} t}, 1, e^{i \frac{2 π}{T} t}, e^{i \frac{4 π}{T} t}, . . .)$ 的无穷维空间的基变换。这意味着我们从一个产生复数的函数，映射到一个产生无穷维复向量的新函数。利用欧拉公式把基转换为 $(. . ., \cos (\frac{- 2 π}{T}), \sin (\frac{- 2 π}{T}), 1, \cos (\frac{2 π}{T}), \sin (\frac{2 π}{T}), . . .)$ ，这也可以表示为 $ℝ$ 上的无穷维基。由于这是一个线性运算，并且是从一个函数空间映射到另一个函数空间，我们称之为线性算子。

逼近方波¶

图 12 方波。¶

方波差不多算是信号处理界的“hello world”。方波定义为 $[- π, π]$ 上的周期波，其中半个周期为“开”，即 $1$ ，半个周期为“关”，即 $0$ 。

(19)¶

\begin{matrix} square (t) = sgn (\sin (t)) \end{matrix}

其中 sgn 是符号函数，定义如下

(20)¶

\begin{split} \begin{equation}\text{sgn}(x)=\begin{cases}1 &\text{if}\; x>0 \\ 0 &\text{if}\; x=0 \\ -1 &\text{if}\; x ## 画羊驼¶ 只要能把想画的东西表示成一条相当平滑的闭合曲线，我们就能用上一节展示的同样原理画出几乎任何东西。至于如何得到表示羊驼的函数，有两种做法。我们可以在矢量图形程序中画一只羊驼，导出为 SVG，然后直接用这个 SVG，因为 SVG 本来就是数学定义的曲线（直线、贝塞尔曲线等）。或者，我们可以对这个 SVG 按等间距采样，得到一组采样点，再对这些采样点做离散傅里叶变换（DFT），从而得到傅里叶系数。本文采用第二种做法，因为它更直接。 ``` const numSamples = 400; let pts = []; function entry() { const svgUrl = './images/llama.svg'; const response = await fetch(svgUrl); const text = await response.text(); pts = await sampleSvgPoints(text, numSamples); } entry(); // Might not sample exactly numPts number of sample points, but it will be // close. Not a concern using the O(N^2) DFT. Using FFT, that can be // concerning since the ideal number of point samples is an exponent of 2. async function sampleSvgPoints(xml: string, numPts: number): Promise { const parser = new DOMParser(); const doc = parser.parseFromString(xml, 'application/xml'); let samplePts = []; const paths = Array.from(doc.querySelectorAll('path')); const totalLength = paths.reduce((acc, v) => acc + v.getTotalLength(), 0); const step = (totalLength / numPts); paths.forEach(path => { const sample = []; const n = M.floor(path.getTotalLength()); for (let i = 0; i numPts) samplePts.pop(); return samplePts; } function toPt(pt: SVGPoint): [number, number] { return [pt.x, pt.y]; } ``` 对 SVG 采样取点接下来，我们对这些采样点做离散傅里叶变换（DFT），求出所需的傅里叶系数，就像上面处理方波那样。离散傅里叶变换的公式如下 (33)¶\[ \begin{equation}c_k=\frac{1}{N}\sum_{j=0}^{N}{p(j)e^{-i\frac{2\pi}{N}kj}}\end{equation} ParseError: Expected & or \\ or \cr or \end at position 121: …\text{if}\; x #̲# 画羊驼¶ 只要能把想画的…

其中 $N$ 是待变换的采样点列表中的采样点总数。在我们的场景里， $p (k)$ 本质上就是按索引取数组元素。例如，若 $p$ 是采样点列表，在 Javascript 中就是 p[k]。DFT 与傅里叶级数其实是一回事，区别在于前者是对等间距的有限个采样点求和，而不是让步长趋于 $0$ 取极限，等价地说，是让求和的项数趋于 $\infty$ 。DFT 也是把时域上 $N$ 个等间距的复数采样点转换到频域表示的方法，也就是求出每个整数倍正弦频率的强度。这正是把基变换直接应用到 $N$ 维基 $(. . ., e^{i \frac{- 4 π}{T} t}, e^{i \frac{- 2 π}{T} t}, 1, e^{i \frac{2 π}{T} t}, e^{i \frac{4 π}{T} t}, . . .)$ 上，几节前在 Complex Exponentials as Sines and Cosines 一节中已经解释过。

用两层循环就能相当容易地实现 DFT：外层循环计算每个系数 $c_{k}$ ，内层循环做数值积分求和，再除以 $N$ ，得到每个频率 $k$ 的平均强度。这种做法非常慢，但在点数不太多时还能用。其复杂度为 $𝒪 (k j)$ ，最坏情况下是 $𝒪 (N^{2})$ 。用更高效的算法可以做得好得多，比如 Radix-2 快速傅里叶变换（FFT），平均复杂度为 $𝒪 (N \log (N))$ ，不过这篇文章已经太长了，那可以留作以后的话题。

import * as M from 'mathjs';

export type DFTData = M.Complex & {
 freq: number;
 radius: number;
 phase: number;
};

export async function dft(cps: M.Complex[]): Promise {
 const X = [];
 const N = cps.length;
 for (let i = 0; i

$𝒪 (N^{2})$ DFT implementation in TypeScript.

既然拿到了系数，就可以把曲线参数化，从而把它画出来。

function parameterize(cn: M.Complex[]) {
 return async (t: number): Promise =>
 await sigma(k =>
 M.multiply(M.complex(cn[k - 1].re, cn[k - 1].im), 
 M.exp(M.multiply(M.i, (k - 1) * t))), [1, numSamples]) as M.Complex;
}

那么，用这组系数去逼近一只羊驼，应该看到的是：正弦项越多，结果越像羊驼。下面这段动画能看出来。这里的「模式」（modes）指的是用到的正频率个数，另外还用了 $0$ 频率以及对应的负频率。

Fig. 14 $1 \to 500$ modes of a llama approximation.¶

从动画里能看出，大部分活是在低频段干完的，因为它很快就有了羊驼的样子；再加上越来越多的高频项，小的细节才会慢慢补上。

In order to recreate the image at the beginning of the article, with the circles moving around in organized chaos drawing the outline of the llama, we need some code to go through and draw each frequency’s circle. To do that we need the radius and phase which we have included in the DFTData type in the first code block. To do the actual graphics representation, we’ll be using the p5.js library. Note that the code depends on variables defined in the above code blocks. I’ll also provide the full source at the end.

function draw() {
 for (let i = 0; i M.pi * 2)
 points.pop();

time += ((2 * M.pi) / pts.length);
&#125;

这样就能画出文章开头展示的那段动画。

Fourier 动画的完整源码在这里：

adekau/fourier —— 用傅里叶变换画图。

](https://github.com/adekau/fourier)

Index Notes on Classical & Intuistionistic Logic