FFT笔记

Posted on 2019-02-03 Views:

~~上次更新是在十月份？我到底摸了多久~~

卷积

对于两个序列 $\{x\}, \{y\}$ ，定义它们的卷积 $\{x \ast y\}$ 为：

$(x \ast y)_k = \sum_{j = 0}^k x_jy_{k- j}$

单位根

定义复平面单位圆上的次单位根

$\omega_{n} = e^{\frac{2\pi i}{n}} = \cos\left(\frac{2\pi}{n}\right) + i\sin\left(\frac{2\pi}{n}\right)$

易证其有如下的性质：

$\omega_{2n}^{2k} = \omega_n^k$
$\omega_n^k = \omega_n^{k +n}$
$\omega_n^k = - \omega_{n}^{k + n / 2}$

DFT

对于一长度为的序列 $\{x\}$ ，定义其离散傅里叶变换（DFT） $\mathcal{F}\{x\} = \{X\}$ 为：

$\mathcal{F}\{x\}_k = X_k = \sum_{j = 0}^{n - 1} x_i\omega_n^{jk}$

则有其离散傅里叶逆变换（IDFT）为：

$x_k = \frac{1}{n}\sum_{j = 0}^{n - 1} X_i\omega_n^{-jk}$

（标准定义里面 $\omega_n^{jk}$ 和 $\omega_n^{-jk}$ 应该互换，但是不影响）

暴力证明：考虑 $\mathcal{F}\{x_n\}$ 的计算式：

$\begin{pmatrix} 1 &1 &1 &\cdots &1 \\ 1 &\omega_n &\omega_n^2 &\cdots &\omega_n^{n - 1} \\ 1 &\omega_n^2 &\omega_n^4 &\cdots &\omega_n^{2(n-1)} \\ \vdots &\vdots &\vdots &\ddots &\vdots \\ 1 &\omega_n^{n - 1} &\omega_{n}^{2(n-1)} &\cdots &\omega_n^{(n-1)^2} \end{pmatrix} \begin{pmatrix} x_1 \\x_2 \\x_3 \\ \vdots \\x_{n - 1} \end{pmatrix} = \begin{pmatrix} X_1 \\X_2 \\X_3 \\ \vdots \\X_{n - 1} \end{pmatrix}$

设左边矩阵为

，接下来从充分性上证明 $V^{-1}_{jk} = \frac{\omega_n^{-jk}}{n}$ ：

$\left(VV^{-1}\right)_{jk} = \sum_{l = 0}^{n - 1}\omega_n^{jl}\cdot\frac{\omega_n^{-lk}}{n} = \frac{1}{n}\sum_{l = 0}^{n - 1}\omega_n^{l(j - k)}$

显然当

时有 $\left(VV^{-1}\right)_{jk} = 1$ ，而 $j \neq k$ 时，由等比数列求和公式：

$\begin{aligned} \frac{1}{n}\sum_{l = 0}^{n - 1}\omega_n^{l(j - k)} &= \frac{1}{n} \cdot \frac{1 - \omega_n^{n(j - k)}}{1 - \omega_n^{j - k}} = 0 \end{aligned}$

因此 $VV^{-1} = I$ ，得证。

上面的DFT有如下意义：

将 $\{x_n\}$ 作为一组正交基 $\{\omega_n^0, \omega_n^1, \cdots\}, \{\omega_n^0, \omega_n^{2}, \cdots\}$ 的系数，计算和向量。（在此意义下逆变换就是将向量进行分解，因此逆变换的意义显然，无需上述证明）
将 $\{x_n\}$ 作为以个多项式的系数，计算该多项式在 $\omega_n^0, \omega_n^1,\cdots$ 点上的取值。

DFT是线性的，即和的DFT等于DFT的和，且常数倍的DFT是DFT的常数倍。

同时DFT的另一个重要性质被称为卷积性，即

$\mathcal{F}\{x\}_k\mathcal{F}\{y\}_k = \mathcal{F}\{x \ast y\}_k$

证明：把式子写出来：

$\left(\sum_{j = 0}^{n - 1} x_j \omega_n^{jk}\right)\left(\sum_{j = 0}^{n - 1} y_j \omega_n^{jk}\right) \overset{?}{=} \sum_{j = 1}^{n - 1}\left(\sum_{l = 0}^jx_ly_{j - l}\right)\omega_n^{jk}$

显然左边式子展开对应的任意一项 $x_ay_b\omega_n^{(a + b)k}$ 都唯一地出现在右式

内层求和的展开式中，因此两个式子总是等价的。

DFT将卷积的操作转化为的操作，因此可以用于加速卷积。

FFT

因为DFT对应点相乘这一操作本身是 O(n) 的，因此我们只要找到快速的DFT/IDFT算法，就可以避免 O(n^2) 地计算卷积。FFT就是这样一种基于分治的快速DFT算法，时间复杂度为 $O(n\log n)$ 。

首先假设我们需要计算一个序列 $\{x\}$ 的DFT（假设是的整数次幂），我们按照下标的奇偶性将这个序列分为 $\{a\}, \{b\}$ 两部分，其中 $a_j = x_{2j}$ 且 $b_j = x_{2j + 1}$ 。那我们有：

$\begin{aligned} \mathcal{F}\{x\}_k &= \sum_{j = 0}^{n - 1} x_j\omega_n^{jk} \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} \left[(x_{2j}\omega_n^{2jk} + x_{2j + 1}\omega_n^{(2j + 1)k}\right] \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} x_{2j}\omega_n^{2jk} + \omega_n^k \sum_{j = 0}^{\frac{n}{2} - 1} x_{2j + 1}\omega_n^{2jk} \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} a_j\omega_{n / 2}^{jk} + \omega_n^k \sum_{j = 0}^{\frac{n}{2} - 1} b_j\omega_{n / 2}^{jk} \\ &= \mathcal{F}\{a\}_k + \omega_{n}^k\mathcal{F}\{b\}_k \end{aligned}$

~~看起来好像没有问题？~~我们最后一步推导其实是有一点小问题的，因为当 $\mathcal{F}\{a\}$ 依据我们的定义只有 $\frac{n}{2}$ 项，因此上面的推导只能让我们计算 $0 \le k < \frac{n}{2}$ 的值。对于 $k \ge \frac{n}{2}$ 的情况要如何处理呢？

首先，注意到：

$\omega_n^{2jk} = \omega_n^{2jk - jn} = \omega{_n^{2j(k - n / 2)}}$

因此我们稍微变一下：

$\begin{aligned} \mathcal{F}\{x\}_k &= \sum_{j = 0}^{n - 1} x_j\omega_n^{jk} \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} \left[x_{2j}\omega_n^{2j(k - n /2)} + x_{2j + 1}\omega_n^{2j(k - n /2) + k}\right] \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} x_{2j}\omega_n^{2j(k - n /2)} + \omega_n^k \sum_{j = 0}^{\frac{n}{2} - 1} x_{2j + 1}\omega_n^{2j(k - n /2)} \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} a_j\omega_{n / 2}^{j(k - n /2)} + \omega_n^k \sum_{j = 0}^{\frac{n}{2} - 1} b_j\omega_{n / 2}^{j(k - n /2)} \\ &= \mathcal{F}\{a\}_{k - n / 2} - \omega_{n}^{k - n / 2}\mathcal{F}\{b\}_{k - n / 2} \end{aligned}$

（其实本质上就是DFT的序列其实是可以周期拓展的）

这使得我们得以计算 $\frac{n}{2} \le k < n$ 时DFT的值。因此，最为朴素的FFT实现就是基于递归不断将 $\{x_n\}$ 进行分治，并通过上式对于分治结果进行合并。时间复杂度自然就是 $O(n\log n)$ 。

然而事实一定没有那么简单。递归写法的FFT常数很大，很大程度上是因为合并本身不是in-place的，需要额外申请空间。更快的in-place FFT实现，Cooley–Tukey算法，是基于这样一个思路的：

假设每次分出来的 $\{a\}$ 向左走， $\{b\}$ 向右走，那递归版的FFT就会画出一棵很漂亮的二叉树，其中叶子节点对应的下标和序列本身的下标顺序是有很大差别的。我们可以预先快速计算出叶子节点的顺序，然后直接从叶子开始就地一层一层地合并上去。

那叶子节点的顺序遵循怎样的规律呢？其实我们依据奇偶分组本质上是在依据二进制的末位分组，在这样的思想下，我们~~脑补一下~~发现这就是叶子节点的下标把原来序列的下标的二进制翻转一下（低位变高位，高位变低位）（即 $\{00, 01,10, 11\}$ 到叶子节点变成 $\{00, 10, 01, 11\}$ ）（可以结合画图理解）。这样的话我们维护一个从高位到低位的二进制加法器就行了。

代码如下：

#include <cstdio>
#include <complex>
#include <cmath>

using namespace std;
typedef complex<double> cplx;
const double PI = 3.1415926535897932384626;

void fft(int len, cplx *x, cplx *X, int dir = 1) {
    for (int i = 0, k = 0; i < len; i++) {
        X[i] = x[k];
        int y = len >> 1;
        while (y & k) k ^= y, y >>= 1; // 维护倒序加法器
        k |= y; 
    }
    for (int s = 2; s <= len; s <<= 1) {
        cplx w_(cos(2 * PI / s), dir * sin(2 * PI / s));
        for (int l = 0; l < len; l += s) {
            int mid = l + (s >> 1);
            cplx w = 1;
            for (int i = 0; i < s >> 1; i++) { // 逐层向上合并（一次蝴蝶操作）
                cplx t = w * X[mid + i];
                X[mid + i] = X[l + i] - t;
                X[l + i] = X[l + i] + t;
                w *= w_; 
            }
        }
    }
}