使用 Bird-Meertens Formalism 推导 FFT

BMF 是一种非常高妙的程序演算。我们用它来简单推导一下 FFT。

$\gdef\concat{\mathbin{+\mkern-8mu+} } \gdef\map{\operatorname{\ast} } \gdef\red{\operatorname{/} } \gdef\single{\operatorname{[\cdot]} } \gdef\comp{\operatorname{\circ} } \gdef\tri{\operatorname{\triangle} } \gdef\app{\operatorname{\$} } \gdef\repeat{\operatorname{@} } \gdef\fk{\operatorname{\zeta} } \gdef\sq{\operatorname{\square} } \gdef\len{\operatorname{\ell} } \gdef\key#1{ {\color{red}\, #1 \,} }$

基本运算符号

$f \map$ 表示将一个列表中所有元素变为应用了 $f$ 之后的结果（map）。
$\oplus \red$ 表示将列表中的元素按照具有结合律的运算 $\oplus$ 合并（reduce）。
$\single$ 表示将某个元素转化为一个长度为 1 的仅包含它自己的列表（singleton）。
$\concat$ 表示连接两个序列。
$\Tau$ 表示二维列表的转置（transpose）。
$\Gamma_{\oplus}$ 表示将两个等长的列表对应项通过二元运算 $\oplus$ 计算出结果形成的新的列表（zip with）。
$\len x$ 表示 $x$ 的长度。
$f \tri x$ 的变换规则如下： $(f\tri x)_i = (f\comp)^{i}(x_i)$ ，其中 $f$ 是一个映射， $(f \comp)^i$ 表示迭代 $i$ 次。
$f \sq x$ 的变换规则如下： $(f \sq x)_i = (f \comp)^{\len x}(x_i)$ 。
$f \app x$ 表示 $f(x)$ 。
$n \repeat x$ 表示一个长度为 $n$ 全是 $x$ 的列表。
$n \fk z = \omega_n\times z$ 。

基本定理

map 对 $\circ$ 有分配律： $(f \circ g)\map = (f\map)\circ (g\map)$ 。

关于 $\Tau$ 的一系列性质：

$\Tau$ 的一个等价形式： $\Tau = \Gamma_{\concat} \red\circ \single \map \map$ 。
$\Tau$ 对 map 有分配律： $f \ast \ast \circ \Tau = \Tau \circ f \ast \ast$ 。
$\Tau$ 对 reduce 的影响： $\Gamma_{\oplus}\red = \oplus \red \map \comp \Tau$ 。因此可以把前者理解为竖着的 reduce。

我们称 $\map$ ， $\tri$ ， $\sq$ 为 pointwise 运算，下面的 $\oplus$ ， $\otimes$ 均为 pointwise 运算：

关于 $\map$ 交换律：若 $f \circ g = g \circ h$ ，那么 $f \oplus\comp g \map = g \map \comp h \operatorname\oplus$ 。
若 $f\circ g = g\circ f$ ，那么 $f \oplus \comp g \operatorname\otimes = g \otimes \comp f \operatorname\oplus$ ，并且有 $g\oplus \comp f \operatorname\oplus = (g\comp f)\oplus = (f\comp g)\oplus = f\oplus\comp g \operatorname\oplus$ 。而且这也可以推导出 $f\oplus\comp f\operatorname\otimes = f\otimes \comp f \operatorname\oplus$ 。
$f\oplus\otimes\comp \Tau = \Tau\comp f\otimes\oplus$ 。

有关 $n\fk$ 的一个性质是 $(nm\fk)^m = n\fk$ 。

有关 $n\repeat$ 的性质：

分解： $n\repeat\concat\red = \concat\red \comp q\repeat \comp p\repeat \circ\concat\red$ ，其中 $p\times q = n$ 。
与 $\concat\red$ 交换： $n\repeat\comp\concat\red = \concat\red\map\comp n\repeat$ 。
假设操作对象是一个矩阵（二维列表），那么 $k\repeat = \Tau\comp k\repeat\map$ 。 $\Tau$ 对于高维列表来说相当于交换第一二维的坐标，因此 $y = (k\repeat\map) x$ 得到的结果是 $y_{i,t, j} = x_{i, j}$ （ $0\le t < k$ ），再转置一下就能得到 $y_{t, i, j} = x_{i, j}$ 。

描述 DFT

对于序列 $x_0, \ldots, x_{n - 1}$ 的离散傅里叶变换可以描述为

$\mathcal{F}_n = +\red\map \comp n\fk\tri\tri \comp n\repeat$

理解方式如下：

$(n\fk)\tri$ 表示将序列 $a_j$ 变成 $a_jw_n^j$ 。
$(n \fk)\tri\tri$ 表示将第 $i$ 个序列变成 $a_jw_n^{ij}$ 。
$+\red\map$ 表示分别求和。

分治

观察上面的 specification，我们首先要处理的是形如 $f\tri\tri$ 的变换。考虑分治。

先考虑一个简单的情况：对一维列表做 $f\tri$ 变换。将这个过程分治，相当于将其切成若干段，不妨考虑切出来的段等长。那么把原列表做 $f \tri$ 变换就可以转换为，把矩阵每行分别做变换，再连起来。写成 BMF 就是

$f\tri\comp\concat\red = \concat\red\comp f\sq\tri\comp f \tri \map$

同理， $f\tri\tri$ 是针对二维列表的变换：

$f\tri\tri$ 相当于对第 $i$ 行第 $j$ 列的元素迭代 $i\times j$ 次。
$f \map\tri$ 相当于对第 $i$ 行第 $j$ 列的元素迭代 $i$ 次。
$f \tri\map$ 相当于对第 $i$ 行第 $j$ 列的元素迭代 $j$ 次。

这个过程也可以分治。这里我们的分治方法是：对于一行，按照一维的情况来做。此外，我们也可以将所有行分成若干组，每一组行数相等。分割后我们拿到的其实是一个四维列表。把四维列表拼成二维列表的变换是 $\concat/\circ\concat \red \map\map$ 。

接下来做一些 BMF 推导：

$\begin{aligned} & \key{f\tri\tri\circ \concat\red}\circ\concat\red\map\map \\ = & \concat\red\comp f\tri\sq\tri\comp f\tri\tri\key{\map}\comp \concat\red\map\key{\map}\\ = & \concat\red\comp f\tri\sq\tri\comp (\key{f\tri}\tri\circ \key{\concat\red}\map)\map\\ = & \concat\red\comp f\tri\sq\tri\comp (\concat\red\map \comp (f\sq\tri\comp f \tri \map) \tri)\key{\map} \\ = & \concat\red \comp f\tri\sq\tri \comp \concat\red\map\map \comp (f\sq\tri\comp f \tri \map) \tri\map \\ \end{aligned}$

这里我们单独把 $f\tri \sq \tri \comp \concat \red \map \map$ 拿出来给大家变换一下：

$f\tri\comp\concat\red = \concat\red\comp f\sq\tri\comp f \tri \map$ 。
$\key{f\tri} \sq \comp \key{\concat\red} \map = \concat\red\map \comp (f\sq\tri\comp f \tri \map)\sq$ 。
$\key{f\tri\sq}\tri\comp \key{\concat\red\map}\map = \concat\red\map\map \comp (f\sq\tri\comp f \tri \map)\sq\tri$ 。

根据基本定理， $f\sq\tri\comp f\tri\map = f\tri\map\comp f\sq\tri$ （可交换）因此

$\begin{aligned} & \concat\red\comp f\tri\sq\tri\comp \concat\red\map\map \comp (f\sq\tri\comp f \tri \map) \tri\map \\ = & \concat\red\comp\concat\red\map\map \comp (\key{f\sq\tri\comp f \tri \map})\sq\tri \comp (\key{f\sq\tri\comp f \tri \map}) \tri\map \\ = & \concat\red\comp\concat\red\map\map \comp f\sq\tri\sq\tri \comp f \tri\map\sq\tri \comp f\sq\tri\tri\map \comp f\tri\map\tri\map \\ \end{aligned}$

在 FFT 当中， $f = n\fk$ 。我们假设 $f\sq\tri\sq\tri$ 应用到的列表是一个 $a\times b\times c\times d$ 的四维列表，且 $bd = n$ 。那么

$\begin{aligned} f\sq\tri\sq\tri & = n\fk \sq\tri\sq\tri \\ & = (n\fk)^d \map\tri\sq\tri \\ & = ((n\fk)^d \map\tri)^b\map\tri \\ & = (n\fk)^{db}\map\tri\map\tri \\ & = 1\fk\map\tri\map\tri \\ & = 1\fk\map\map\map\map \\ \end{aligned}$

相当于啥都不干。

推导 FFT

设 $n = p\times q$ ， $\mathcal{F} \comp\concat\red$ 作用于一个 $p\times q$ 的二维列表。那么

$\begin{aligned} \mathcal{F}_n\comp \concat\red & = +\red\map \comp n\fk\tri\tri \comp \key{n\repeat \circ \concat\red} \\ & = +\red\map \comp n\fk\tri\tri \comp \concat\red \comp \key{q\repeat \comp p\repeat \comp\concat\red} \\ & = +\red\map \comp \key{n\fk\tri\tri \comp \concat\red \comp\concat\red\map\map} \comp q\repeat \comp p\repeat \\ \end{aligned}$

这里我们设 $K = n\fk\sq\tri\sq\tri \comp n\fk \tri\map\sq\tri \comp n\fk\sq\tri\tri\map \comp n\fk\tri\map\tri\map$ ：

$\begin{aligned} & +\red\map \comp \key{n\fk\tri\tri \comp \concat\red \comp\concat\red\map\map} \comp q\repeat \comp p\repeat \\ = & \key{+\red\map \comp \concat\red} \comp\concat\red\map\map \comp K \comp q\repeat \comp p\repeat \\ = & \concat\red \comp \key{+\red\map\map \comp \concat\red\map\map} \comp K \comp q\repeat \comp p\repeat \\ = & \concat\red \comp (+\red\circ \concat\red)\map\map \comp K \comp q\repeat \comp p\repeat \\ = & \concat\red \comp (+\red\circ +\red\map)\map\map \comp K \comp q\repeat \comp p\repeat \\ \end{aligned}$

最后一步用到了加法结合律。然后我们推导一下后面的 $K \comp q\repeat \comp p\repeat$ 。注意到 $q\repeat\comp p\repeat$ 作用于 $p\times q$ 的结果是得到一个 $q\times p \times p\times q$ 的列表，因此

$\begin{aligned} K & = pq\fk\key{\sq}\tri\key{\sq}\tri \comp pq\fk \tri\map\key{\sq}\tri \comp pq\fk\key{\sq}\tri\tri\map \comp n\fk\tri\map\tri\map \\ & = \key{1\fk\map\tri\map\tri} \comp q\fk \tri\map\map\tri \comp p\fk\map\tri\tri\map \comp n\fk\tri\map\tri\map \\ & = q\fk \tri\map\map\tri \comp \key{p\fk\map\tri\tri\map} \comp \key{n\fk\tri\map\tri\map} \\ & = q\fk \tri\map\map\tri \comp n\fk\tri\map\tri\map \comp p\fk\map\tri\tri\map \\ \end{aligned}$

$K$ 的类型为 $q\times p\times p\times q\to q\times p\times p\times q$ 。接下来来推导 $\mathcal{F}_n\comp \concat\red\comp\Tau$ ：

$\begin{aligned} &\mathcal{F}_n \comp\concat\red \comp \Tau\\ = &\concat\red \comp (+\red\circ +\red\map)\map\map\\ & \comp q\fk \tri\map\map\tri \comp n\fk\tri\map\tri\map \comp p\fk\map\tri\tri\map \\ & \comp q\repeat \comp p \repeat \comp \key{\Tau}\\ = &\concat\red \comp (+\red\circ +\red\map)\map\map\\ & \comp \key{q\fk \tri\map}\map\tri \comp \key{n\fk\tri\map}\tri\map \comp \key{p\fk\map\tri}\tri\map \comp \key{\Tau}\map\map \\ &\comp q\repeat \comp p \repeat\\ = &\concat\red \comp (+\red\circ +\red\map)\map\map \comp \key{\Tau\map\map}\\ & \comp q\fk \map\tri\map\tri \comp n\fk\map\tri\tri\map \comp p\fk\tri\map\tri\map \\ & \comp q\repeat \comp p \repeat\\ = &\concat\red \comp (+\red\circ +\red\map\comp \Tau)\map\map\\ & \comp q\fk \map\tri\map\tri \comp n\fk\map\tri\tri\map \comp p\fk\tri\map\tri\map \\ & \comp q\repeat \comp p \repeat\\ = &\concat\red \comp (+\red\circ +\red\map)\map\map\\ & \comp q\fk \map\tri\map\tri \comp n\fk\map\tri\tri\map \comp p\fk\tri\map\tri\map \\ & \comp q\repeat \comp p \repeat\\ \end{aligned}$

最后一步用到了加法交换律。注意到

$k\repeat\comp f = f \map \comp k \repeat$
分配律： $+\red \comp k\fk \map =k\fk\comp+\red$

因此

$\begin{aligned} &\concat\red \comp (\key{+\red\circ +\red\map})\map\map\\ & \comp q\fk \map\tri\map\tri \comp n\fk\map\tri\tri\map \comp p\fk\tri\map\tri\map \\ & \comp q\repeat \comp p \repeat\\ = &\concat\red \comp +\red\map\map \comp \key{+\red}\map\map\map \\ & \comp \key{q\fk\map}\tri\map\tri \comp \key{n\fk\map}\tri\tri\map \comp p\fk\tri\map\tri\map \\ & \comp q\repeat \comp p \repeat\\ = &\concat\red \comp +\red\map\map \comp q\fk\tri\map\tri \\ & \comp n\fk\tri\tri\map \comp +\red\map\map\map \comp p\fk\tri\map\tri\map \\ & \comp \key{q\repeat} \comp p \repeat\\ = &\concat\red \\ &\comp +\red\map\map \comp q\fk\tri\map\tri \comp q \repeat \\ &\comp n\fk\tri\tri \\ &\comp +\red\map\map \comp p\fk\tri\map\tri \comp p \repeat \\ \end{aligned}$

设 $G(k) = +\red\map\map\comp k\fk\tri\map\tri\comp k\repeat$ 。

$\begin{aligned} G(k) & = +\red\map\map\comp k\fk\tri\map\tri\comp \key{k\repeat} \\ & = +\red\map\map\comp \key{k\fk\tri}\map\tri\comp \key\Tau \comp k\repeat\map \\ & = \Tau\comp+\red\map\map\comp k\fk\tri\tri\map \comp k\repeat\map\\ & = \Tau\comp(+\red\map\comp k\fk\tri\tri \comp k\repeat)\map \\ & = \Tau\comp\mathcal{F}_k\map \end{aligned}$

因此

$\begin{aligned} \mathcal{F}_n\comp \concat\red & = (\mathcal{F}_n\comp \concat\red\comp\Tau) \comp \Tau \\ & = (\concat\red \comp \Tau\comp\mathcal{F}_q\map \comp n\fk\tri\tri \comp \Tau\comp\mathcal{F}_p\map) \comp \Tau \end{aligned}$

记 $(\concat\red)^{-1}_{p\times q}$ 表示将一个长度为 $p\times q$ 的序列分成一个 $p\times q$ 的二维列表，我们就得到了 FFT 的最终表示：

$\begin{aligned} \mathcal{F}_n & = \concat\red \comp \Tau\comp\mathcal{F}_q\map \comp n\fk\tri\tri \comp \Tau\comp\mathcal{F}_p\map \comp \Tau\comp (\concat\red)^{-1}_{p\times q} \end{aligned}$

平时我们写的 FFT 是取 $p = 2$ 或者 $q = 2$ （二分治），这时 $\mathcal{F}_2$ 就会退化为 $x_0, x_1\to x_0 + x_1, x_0 - x_1$ （不需要乘法）。

参考文献

Deriving the fast Fourier Algorithm