dramforever

coding, thoughts and stuff irl

6 步搞懂 FFT

2019-09-15

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{cases} n & k \bmod n = 0 \\ 0 & k \bmod n \neq0 \end{cases} \]

第二步:DFT

DFT 是一个 \(n\) 元序列到 \(n\) 元序列的映射,如下:

\[ \mathsf{DFT}(a)_k = \sum_p a_p \omega_n^{pk} \]

其逆如下:

\[ \mathsf{IDFT}(a)_k = {1 \over n} \sum_p a_p \omega_n^{-pk} \]

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

\[ A=\mathsf{DFT}(a) \]

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

考虑从 \(a\)\(b\) 序列 \(c\)(用 \(\circ\) 表示序列按位乘,如 \(\langle 1,2,3 \rangle \circ \langle 4,5,6 \rangle = \langle 4, 10, 18 \rangle\)):

\[ c = \mathsf{IDFT}(A \circ B) \]

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

\[ \begin{align} c_k & = {1 \over n} (\sum_p (\sum_q a_q \omega_n^{qp})(\sum_r b_r \omega_n^{rp})\ \omega_n^{-pk} ) \\ & = {1 \over n} (\sum_q \sum_r a_q b_r (\sum_p \omega_n^{(q+r-k) p})) \end{align} \]

根据单位根基本性质,将最里面的求和式化简:(方括号表示真取 \(1\) 假取 \(0\)

\[ \sum_p \omega_n^{(q+r-k) p} = n \cdot [(q+r-k) \bmod n=0] \]

故有:

\[ c_k = \sum_q \sum_r a_q b_r [(q+r-k) \bmod n=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\),也就是说相当于 \(\bmod n\) 不需要考虑,那么循环卷积特化为多项式乘法:

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

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

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

第四步:DFT 的两个性质

线性性(此处加法是按位加),比较容易得到:

\[ \mathsf{DFT}(A+B) = \mathsf{DFT}(A) + \mathsf{DFT}(B) \]

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

\[ x = \langle a_0, 0, a_1, 0, \dots,a_{n-1},0 \rangle \\ y = \langle 0, b_0, 0, b_1, \dots,0, b_{n-1} \rangle \\ \]

则对 \(x\),利用其中只有一半元素非零,以及单位根的性质可得: \[ \begin{align} 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_{k \bmod n} \end{align} \]

类似的,对 \(y\)

\[ \begin{align} 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_{k \bmod n} \end{align} \]

(取模的原因是 \(k\) 可能超过 \(n-1\),但是 DFT 原式代入 \(k\)\(k + n\) 是一样的)

第五步:FFT

取上面的 \(x + y\) 得:

\[ z = x+y=\langle a_0, b_0, a_1, b_1, \dots , a_{n-1}, b_{n-1} \rangle \]

利用线性性有:

\[ Z_k = X_k + Y_k = A_{k \bmod n} + \omega_{2n}^k B_{k \bmod n} \]

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

这个方法叫做 FFT。

第六步:二进制重排

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

如对于 \(n = 8\) 按照分治的顺序排好是:

\[ \langle 0, 4, 2, 6, 1, 5, 3, 7 \rangle \]

然后迭代合并:

\[ \begin{align} & \langle 0, 4, 2, 6, 1, 5, 3, 7 \rangle \\ \to\ & \langle 04, 26, 15, 37 \rangle \\ \to\ & \langle 0246, 1357 \rangle \\ \to\ & \langle 01234567 \rangle \end{align} \]

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

完工撒花!

本文使用 PandocMarkdown 原始文件生成,参数为:

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