FFT，启动！

第一步：单位根

如下形式的复数称为 $n$ 次单位根（其中 $i$ 是虚数单位）：

$${\omega }_n^k=e^{2\pi ik/n}$$

（这里上标是幂次）

基本性质有：

$${\omega }_n^k={\omega }_n^{k+xn}={\omega }_{xn}^{xk}$$

$${\omega }_n^a{\omega }_n^b={\omega }_n^{a+b}$$

$$\sum _{p=0}^{n-1}{\omega }_n^{pk}=\{\begin{array}{cc}n& kmodn=0\\ 0& kmodn\ne 0\end{array}$$

第二步：DFT

DFT 是一个 $n$ 元序列到 $n$ 元序列的映射，如下：

$$DFT(a{)}_k=\sum _p{a}_p{\omega }_n^{pk}$$

其逆如下：

$$IDFT(a{)}_k=\frac{1}{n}\sum _p{a}_p{\omega }_n^{-pk}$$

以下默认：大写字母表示对应小写字母的 DFT，如：

$$A=DFT\left(a\right)$$

第三歩：循环卷积和多项式乘法

考虑从 $a$ 和 $b$ 构造的序列 $c$（用 $\circ$ 表示序列按位乘，如 $⟨1,2,3⟩\circ ⟨4,5,6⟩=⟨4,10,18⟩$）：

$$c=IDFT(A\circ B)$$

展开定义将各项重新组合可得：

$$\begin{array}{cc}{c}_k& =\frac{1}{n}\left(\sum _p\right(\sum _q{a}_q{\omega }_n^{qp}\left)\right(\sum _r{b}_r{\omega }_n^{rp}\left){\omega }_n^{-pk}\right)\\ & =\frac{1}{n}\left(\sum _q\sum _r{a}_q{b}_r\right(\sum _p{\omega }_n^{(q+r-k)p}\left)\right)\end{array}$$

根据单位根基本性质，将最里面的求和式化简：（方括号表示真取 $1$ 假取 $0$）

$$\sum _p{\omega }_n^{(q+r-k)p}=n\cdot \left[\right(q+r-k)modn=0]$$

故有：

$$c_k=\sum _q\sum _r{a}_q{b}_r\left[\right(q+r-k)modn=0]$$

（用类似的思路可以证明 DFT 和 IDFT 互逆，此处略去）

也就是说，等式左侧实现了一个叫做循环卷积的过程，伪代码如下：

c = new array of n zeros
for q:
    for r:
        c[(q + r) mod n] += a[q] * b[r]

特别的，考虑需要计算两个多项式的乘法：我们有系数序列 $a$ 和 $b$，$a$ 有 $u$ 项，$b$ 有 $v$ 项。假设系数存的是从常数项在下标 $0$ 开始，我们给 $a$ 后补 $v$ 个 $0$，给 $b$ 后补 $u$ 个 $0$，计算循环卷积时所有 $q+r$ “溢出” $n$ 的项都为 $0$，也就是说相当于 $modn$ 不需要考虑，那么循环卷积特化为多项式乘法：

c = new array of n zeros
for q:
    for r:
        c[q + r] += a[q] * b[r]

注意到大整数乘法可以化为多项式乘法加处理进位。

直接计算是 $O\left(n^2\right)$ 的，我们通过在 $O\left(n\log n\right)$ 时间计算 DFT 和 IDFT，配合 $O\left(n\right)$ 时间按位乘，在 $O\left(n\log n\right)$ 时间完成这个过程。

第四步：DFT 的两个性质

线性性（此处加法是按位加），比较容易得到：

$$DFT(A+B)=DFT\left(A\right)+DFT\left(B\right)$$

还有一个没那么明显：从 $a$ 和 $b$ 构造两个交替为原序列元素和 $0$ 的序列：

$$x=⟨a_0,0,a_1,0,\dots ,a_{n-1},0⟩y=⟨0,b_0,0,b_1,\dots ,0,b_{n-1}⟩$$

则对 $x$，利用其中只有一半元素非零，以及单位根的性质可得： $$\begin{array}{cc}{X}_k& =\sum _{p=0}^{2n-1}{x}_p{\omega }_{2n}^{pk}\\ & =\sum _{p=0}^{n-1}{a}_p{\omega }_{2n}^{2pk}=\sum _{p=0}^{n-1}{a}_p{\omega }_n^{pk}\\ & =A_{kmodn}\end{array}$$

类似的，对 $y$：

$$\begin{array}{cc}{Y}_k& =\sum _{p=0}^{2n-1}{y}_p{\omega }_{2n}^{pk}\\ & =\sum _{p=0}^{n-1}{b}_p{\omega }_{2n}^{(2p+1)k}=\sum _{p=0}^{n-1}{b}_p{\omega }_n^{pk}{\omega }_{2n}^k\\ & ={\omega }_{2n}^k{B}_{kmodn}\end{array}$$

（取模的原因是 $k$ 可能超过 $n-1$，但是 DFT 原式代入 $k$ 和 $k+n$ 是一样的）

第五步：FFT

取上面的 $x+y$ 得：

$$z=x+y=⟨a_0,b_0,a_1,b_1,\dots ,a_{n-1},b_{n-1}⟩$$

利用线性性有：

$$Z_k=X_k+Y_k=A_{kmodn}+{\omega }_{2n}^k{B}_{kmodn}$$

也就是说，我们可以将 $2n$ 长度的 $z$ 的 DFT 拆成两个长度为 $n$ 的 $a$ 和 $b$ 的 DFT，和一个 $O\left(n\right)$ 的合并过程（是奇偶拆，不是中间分开）。不断递归进行这个过程（单元素序列的 DFT 是其自身），可以在 $O\left(n\log n\right)$ 时间完成一个长度为 $n$ 的 DFT，这里 $n$ 需要是 $2$ 的幂。

这个方法叫做 FFT。

第六步：二进制重排

考虑分治的“分”过程，我们首先将偶数位置和奇数位置分开，然后各自再拆偶数位和奇数位，依此类推。在二进制中，就是先按最低位 0 1 排序，然后两半分别按照次低位排序，依此类推，也就是按照二进制反转排序（如 000 100 010 ... 001 101 ...）

如对于 $n=8$ 按照分治的顺序排好是：

$$⟨0,4,2,6,1,5,3,7⟩$$

然后迭代合并：

$$\begin{array}{cc}& ⟨0,4,2,6,1,5,3,7⟩\\ \to & ⟨04,26,15,37⟩\\ \to & ⟨0246,1357⟩\\ \to & ⟨01234567⟩\end{array}$$

这样我们可以使用两个数组来保存本次迭代和上次迭代的结果，每次一对一对合并，使用 $O\left(n\right)$ 的内存空间完成 FFT（直接递归需要 $O\left(n\log n\right)$），而且一般更快。

完工撒花！

注

• 在不同的实现中 DFT 的参数可能不同，如 numpy.fft.fft 的实现中，单位根上的的幂次是我们这里所用的取负，$pk$ 变为 $-pk$（或者说单位根取共轭），ifft 相反。

• 类似的方法可以构造拆成大于 $2$ 份的 FFT，在此略去。

本文使用 Pandoc 从 Markdown 原始文件生成，参数为：

$ pandoc --mathjax -f markdown-auto_identifiers -t html < fft.md > index.html