快速傅里叶变换(FFT)算法

By Long Luo

之前的文章 傅里叶变换(Fourier Transform) 详细介绍了傅里叶变换 $(\textit{Fourier Transform})$ 是什么和做什么，相信你已经感受到了傅里叶变换的强大之处。

但理论要联系实际，今天我们就来学习 快速傅里叶变换 $(\textit{Fast Fourier Transform, FFT})$ ^[1]的实际运用：多项式乘法。

快速傅里叶变换(Fast Fourier Transform, FFT)

快速傅里叶变换 $(\textit{Fast Fourier Transform, FFT})$ ^[2] 是一种可在 $O(nlogn)$ 时间内完成的 离散傅里叶变换^[3] $(\textit{Discrete Fourier transform, DFT})$ 算法。

FFT可以做什么？

傅里叶变换 $(\textit{Fourier Transform})$ 本质上是信号与三角函数进行卷积运算，而快速傅里叶变换 $(\textit{FFT})$ 就是提高卷积的计算效率，时间复杂度从 $O(n^2)$ 降低到 $O(nlogn)$。

$\textit{FFT}$ 在算法中的运用主要是用来加速多项式乘法或大数乘法。

多项式乘法

正如之前的文章 卷积(Convolution) 所说，多项式乘法也是一种卷积运算。

在计算中，泰勒级数^[4] 可以使用多项式函数来逼近一个函数，所以计算中多项式乘法非常重要。

大数乘法

超大数字乘法，可以参考 超大数字的四则运算是如何实现的呢？ 。朴素的算法就是列竖式做乘法，算法时间复杂度为 $O(n^2)$ ，如果数字太大的话，效率也不够高，如果应用 $(\textit{FFT})$ 则可以使算法时间复杂度降至 $O(nlogn)$。

不妨设十进制数字 $num=123456789$ ，很容易知道：

$$123456789 = 1 \times 10^8 + 2 \times 10^7 + 3 \times 10^6 + 4 \times 10^5 + 5 \times 10^4 + 6 \times 10^3 + 7 \times 10^2 + 8 \times 10^1 + 9 \times 10^0$$

令 $x = 10$，则可以转化为：

$$1 \times x^8 + 2 \times x^7 + 3 \times x^6 + 4 \times x^5 + 5 \times x^4 + 6 \times x^3 + 7 \times x^2 + 8 \times x^1 + 9 \times x^0$$

所以大数乘法就是 $x = 10$ 情况下的多项式乘法！

那下面我们就以多项式乘法的为例来学习快速傅里叶变换 $(\textit{FFT})$ 具体是如何做的。

多项式

在学习多项式乘法之前，我们需要先学习下有关多项式的知识。

多项式有两种表示方法: 系数表示法与点值表示法。

系数表示法

设多项式 $A(x)$ 为一个 $n$ 项 $n - 1$ 次的多项式，显然，所有项的系数组成的系数向量 $(a_0, a_1, a_2, \dots, a_{n-1})$ 唯一确定了这个多项式。

$$A(x) = \sum_{i=0}^{n-1}a_i \cdot x^i = a_0 + a_1x^1 + a_2x^2 + \cdots + a_{n-1}x^{n-1} \iff A(x) = {a_0, a_1, \dots, a_{n-1}}$$

点值表示法

点值表示法是把这个多项式看成一个函数，从其中选取 $n$ 个不同的点，从而利用这 $n$ 个点来唯一地表示这个函数。

设

$$\begin{array}{c} A(x_0) = y_0 = a_0 + a_1x_0+a_2x_0^2+a_3x_0^3+ \cdots + a_{n-1}x_0^{n-1} \\ A(x_1) = y_1 = a_0 + a_1x_1+a_2x_1^2+a_3x_1^3+ \cdots + a_{n-1}x_1^{n-1} \\ A(x_2) = y_2 = a_0 + a_1x_2+a_2x_2^2+a_3x_2^3+ \cdots + a_{n-1}x_2^{n-1} \\ \vdots \\ A(x_{n-1}) = y_{n-1} = a_0 + a_1x_{n-1}+a_2x_{n-1}^2+a_3x_{n-1}^3+ \cdots + a_{n-1}x_{n-1}^{n-1} \end{array}$$

那么用点值表示法表示 $A(x)$ 如下：

$$A(x) = a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1} \iff A(x) = {(x_0,y_0), (x_1,y_1), \cdots,(x_{n-1},y_{n-1})}$$

为什么用 $n$ 个不同点就能唯一地表示一个 $n-1$ 次函数？

证明如下：

• Proof $1$ :

两点确定一条直线。再来一个点，能确定这个直线中的另一个参数，那么也就是说 $n$ 个点能确定 $n-1$ 个参数（不考虑倍数点之类的没用点）。

• Proof $2$^[5] :

假设原命题不成立，则存在两个不同的 $n-1$ 次多项式函数 $A(x)$ ，$B(x)$ ，那么 $A(x)$ ，$B(x)$ 有 $n-1$ 个交点，即任何 $i \in [0, n-1]$，有 $A(x_i) = B(x_i)$ 。

令 $C(x) = A(x) - B(x)$ ，则 $C(x)$ 也是一个 $n-1$ 次多项式。对于任何 $i \in [0, n-1]$，都有 $C(x_i) = 0$ 。

即 $C(x)$ 有 $n$ 个根，这与代数基本定理（一个 $n-1$ 次多项式在复数域上有且仅有 $n-1$ 个根）相矛盾，故 $C(x)$ 并不是一个 $n-1$ 次多项式，推导矛盾。

故原命题成立。

多项式乘法

考虑两个多项式 $A(x)$ ，$B(x)$ ，其乘积 $C(x) = A(x) \cdot B(x)$ 。

假设 $A(x)$ 的项数为 $n$ ，其系数构成的 $n$ 维向量为 $(a_0, a_1, a_2, \cdots, a_{n-1})$ ；
$B(x)$ 的项数为 $m$ ，其系数构成的 $m$ 维向量为 $(b_0, b_1, b_2, \cdots, b_{m-1})$^[6] 。

我们要求 $C(x)$ 的系数构成的 $n+m-1$ 维的向量，先考虑暴力做法：

public int[] mutply(int[] A, int[] B) {
    int n = A.length;
    int m = B.length;

    int[] C = new int[n + m];

    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            C[i + j] += A[i] * B[j];
        }
    }

    return C;
}

可见时间复杂度是 $O(n^2)$ 。

如何加速多项式乘法？

实际运用中多项式的项数非常多，比如 $10^5$ 这种级别，那么有没有什么方法可以加速运算呢？

已知在一组插值节点 $(x_0, x_1, x_2, \cdots, x_n)$ 中 $A(x)$ ， $B(x)$ (假设多项式的项数相同，没有则视为 $0$ 。)的点值向量分别为 $(y_{a0}, y_{a1}, y_{a_2}, \cdots, y_{an})$ ，$(y_{b0}, y_{b1}, y_{b_2}, \cdots, y_{bn})$ ，则：

$$\begin{aligned} A(x)&={(x_0, y_{a0}), (x_1, y_{a1}), (x_2, y_{a2}), \cdots, (x_n, y_{an})} \\ B(x)&={(x_0, y_{b0}), (x_1, y_{b1}), (x_2, y_{b_2}), \cdots, (x_n, y_{bn})} \end{aligned}$$

那么 $C(x)=A(x) \cdot B(x)$ ，那么其点值表示法可以在 $O(n)$ 的时间内求出：

$$C(x) = {(x_0, y_{a0} \cdot y_{b0}), (x_1, y_{a1} \cdot y_{b1}), (x_2, y_{a_2} \cdot y_{b2}), \cdots, (x_n, y_{an} \cdot y_{bn})}$$

多项式乘法的系数表示法的时间复杂度是 $O(n^2)$ ，而点值表示法的时间复杂度是 $O(n)$ 。

其实可以观察到，多项式乘法的系数表示法就是做卷积运算，而点值表示法是做乘法运算。卷积定理^[7] 告诉我们在一个域中的卷积相当于另一个域中的乘积。

还记得傅立叶变换的本质就是从另外一个维度看待世界吗？

在这里，系数表示法就是时域上表达，很复杂，而点值表示法就是频域上表示，就非常简单。那么我们将其转换为频域上表示，就可以大大降低其复杂度！

但在得到了 $C(x)$ 的点值表达式之后，还是不行，因为我们要的是系数表达式。

接下来问题变成了从点值回到系数。如果我们带入到高斯消元法的方程组中去，会把复杂度变得非常高。光是计算 $x^i(0 \leq i \leq n)$ 就是 $n$ 项，这就已经 $O(n^2)$ 了，更别说还要把 $n + 1$ 个方程进行消元…

不过我们暂时忽略如何将点值表达式变成系数表达式，关注点放在乘法运算上。因为系数表达式和点值表达式运算复杂度差距如此之大，如果能在运算之前把多项式变成点值表示法，做完乘法之后，再将点值表示法变成系数表示法不就可以大大提高效率了吗？

点值表示法 和 系数表示法 如何转换？

想象有个魔法黑匣子，可以实现多项式的点值表示法和系数表示法的互换。

每次运算前，我们先向黑匣子里输入 $A(x)$ ，$B(x)$ 系数表达式，黑匣子内部先将 $A(x)$ ， $B(x)$ 都变成点值表达式，黑匣子进行乘法运算之后，再在内部转换为系数表达式，返回给我们，不就行了吗？

所以这个黑匣子里面是什么？能不能实现这个黑匣子呢？

其实这个魔法黑匣子的本质就是我们今天要研究的快速傅里叶变换 $(\textit{FFT})$ ，其思路就是由系数表达式到点值表达式，生成结果的点值表达式，再将点值表达式转换为结果的系数表达式。

离散傅里叶变换(Discrete Fourier Transform)

上一章讲到了一个魔法黑匣子可以实现多项式的点值表达式和系数表达式的转换，这一章我们就来研究如何将多项式系数表达式转换为点值表达式。

DFT

考虑一个 $n$ 项( $n = 2^x$ )的多项式 $A(x)$ ，其系数向量为 $(a_0, a_1, a_2, \cdots, a_{n-1})$ ：

$$A(x) = a_0 x^0 + a_1 x^1 + \cdots + a_{n-1} x^{n-1}$$

将 $n$ 次单位根的 $[0, n-1]$ 次幂分别带入 $A(x)$ 得到其点值向量 $(A(w_n^{0}), A(w_n^{1}), A(w_n^{2}), \cdots, A(w_n^{n-1}))$ 。

这个过程称为离散傅里叶变换 $(\textit{Discrete Fourier Transform})$ 。

如果朴素代入计算，时间复杂度也是 $O(n^2)$ ，所以我们必须要利用到单位根 $\omega$ 的特殊性质以减少运算。

奇偶次幂分组

对于 $A(x) = a_0 x^0 + a_1 x^1 + a_2 x^2 + a_3 x^3 + \cdots + a_{n-1} x^{n-1}$ ，

将其按照奇偶次幂分组：

$$\begin{aligned} A(x) &= (a_0 x^0 + a_2 x^2 + a_{4} x^{4} + \cdots + a_{n-2} x^{n-2})+(a_1 x^1 + a_3 x^3 + a_{5} x^{5} + \cdots + a_{n-1} x^{n-1}) \\ &=(a_0 x^0 + a_2 x^2 + a_{4} x^{4} + \cdots + a_{n-2} x^{n-2}) + x \cdot(a_1 + a_3 x^2 + a_{5} x^{4} + \cdots + a_{n-1} x^{n-2}) \end{aligned}$$

令 $A_e(x)=(a_0 x^0 + a_2 x + a_{4} x^{2} + \cdots + a_{n-2} x^{\frac{n-2}{2}})$ ，$A_o(x)=(a_1 + a_3 x + a_{5} x^{2} + \cdots + a_{n-1} x^{\frac{n-2}{2}})$ ，

那么易得： $A(x) = A_e(x^2) + xA_o(x^2)$ 。

分类讨论

1. 当 $0 \leq k \leq \frac{n}{2}-1, k \in Z$ ，代入单位根 $\omega_{n}^{k}$ ，则：

$$A(\omega_{n}^{k}) = A_e(\omega_{n}^{2k}) + \omega_{n}^{k} A_o(\omega_{n}^{2k})$$

由上文提到的折半引理

$$A(\omega_{n}^{k}) = A_e(\omega_{\frac{n}{2}}^{k}) + \omega_{n}^{k} \cdot A_o(\omega_{\frac{n}{2}}^{k})$$

2. 当$\frac{n}{2} \leq k + \frac{n}{2} \leq n-1$时，

$$A(\omega_{n}^{k+\frac{n}{2}}) = A_e(\omega_{n}^{2k+n}) + \omega_{n}^{k + \frac{n}{2}} \cdot A_o(\omega_{n}^{2k+n})$$

其中 $\omega_{n}^{2k+n} = \omega_{n}^{2k} \cdot \omega_{n}^{n} = \omega_{n}^{2k} = \omega_{\frac{n}{2}}^{k}$ 。

由消去引理 $\omega_{n}^{k+\frac{n}{2}} = -\omega_{n}^{k}$ ，

那么：

$$A(\omega_{n}^{k+\frac{n}{2}}) = A_e(\omega_{\frac{n}{2}}^{k}) - \omega_{n}^{k} \cdot A_o(\omega_{\frac{n}{2}}^{k})$$

注意： $k$ 与 $k + \frac{n}{2}$ 取遍了 $[0, n-1]$ 中的 $n$ 个整数，保证了可以由这 $n$ 个点值反推解出系数。

DFT小结

综合这两个式子，如果知道 $A_e(x)$ ，$A_o(x)$ 分别在 $\omega_{\frac{n}{2}}^{0}, \omega_{\frac{n}{2}}^{1}, \cdots, \omega_{\frac{n}{2}}^{\frac{n}{2}-1}$ 处的点值，就可以 $O(n)$ 的时间内求出 $A(x)$ 在 $\omega_{n}^{0}, \omega_{n}^{1}, \cdots, \omega_{n}^{n-1}$ 处的点值。

$A_e(x)$ ，$A_o(x)$ 都是 $A(x)$ 一半的规模，可以转化为子问题递归求解。

所以时间复杂度：

$$T(n) = 2T(\frac{n}{2}) + O(n) = O(nlogn)$$

所以明白在一开始的时候 $n = 2^x$ 吧？

分治 $\textit{DFT}$ 能处理的多项式长度只能是 $2^m(m \in N^ \ast )$，否则在分治的时候左右不一样长。

离散傅里叶逆变换(Inverse Discrete Fourier Transform)

上一节我们讲了离散傅里叶变换 $(\textit{Discrete Fourier Transform})$，即实现了魔法黑匣子的一半，实现了多项式系数表达式转换成点值表达式，那么这一章我们就来实现魔法黑匣子的另一半，将多项式点值表达式转化为系数表达式。

将点值表达式的多项式转化为系数表达式，这个过程叫做离散傅里叶逆变换 $(\textit{Inverse Discrete Fourier Transform})$。

IDFT求解过程

问题：

对于多项式 $A(x)=a_0 x^0 + a_1 x^1 + a_2 x^2 + a_3 x^3 + \cdots + a_{n-1} x^{n-1}$ ，已知 $n$ 个点，其 $n$ 维点值向量为 $(A(x_0), A(x_1), \cdots, A(x_{n-1}))$，请求解其 $n$ 维系数向量 $(a_0, a_1, \cdots, a_{n-1})$？^[6:1]

1. 设 $(d_0, d_1, \cdots, d_{n-1})$ 为 $(a_0, a_1, \cdots, a_{n-1})$ 得到的离散傅里叶变换的结果。

构造一个多项式：

$$F(x) = d_0 + d_1 x + d_2 x^2 + \cdots + d_{n-1} x^{n-1}$$

2. 设向量 $(c_0, c_1, \dots, c_{n-1})$，其中 $c_k$ 为 $F(x)$ 在 $x = \omega_{n}^{-k}$ 的点值表示，即 $c_k = \sum_{i=0}^{n-1} d_i \cdot (\omega_{n}^{-k})^i$。

如何得到$d_i$？

因为

$$c_k = \sum_{i=0}^{n-1}[\sum_{j=0}^{n-1}a_j \cdot (\omega_{n}^{i})^j] \cdot (\omega_{n}^{-k})^i$$

由和式的性质

$$c_k = \sum_{j=0}^{n-1}a_j \sum_{i=0}^{n-1}(\omega_{n}^{i})^j \cdot (\omega_{n}^{-k})^i = \sum_{j=0}^{n-1}a_j \sum_{i=0}^{n-1}(\omega_{n}^{i})^{j-k}$$

令 $S(j, k) = \sum_{i=0}^{n-1}(\omega_{n}^{i})^{j-k}$，对其进行化简：

设 $j-k = \delta$，则：

$$S(j, k) = \omega_{n}^{0} + \omega_{n}^{\delta} + \omega_{n}^{2\delta} + \cdots + \omega_{n}^{(n-1)\delta}$$

可见 $\omega_{n}^{k}$ 构成等比数列，其公比为 $\omega_{n}^{\delta}$ 。

• 当 $\omega_{n}^{\delta} = 1$ 即 $\delta=0$ 时，$S(j,k) = n$ 此时 $\delta = 0 \Rightarrow j-k=0 \Rightarrow j = k$ ；

• 当 $\omega_{n}^{\delta} \ne 1$ 即 $\delta \ne 0$ 时，由等比数列求和公式：

$$S(j,k) = \frac{\omega_{n}^{0}[(\omega_{n}^{\delta})^{n}-1]}{\omega_{n}^{\delta}-1} = \frac{1 [(\omega_{n}^{n})^{\delta}-1]}{\omega_{n}^{\delta}-1} = \frac{[(1)^{\delta}-1]}{\omega_{n}^{\delta}-1} = \frac{0}{\omega_{n}^{\delta}-1}=0$$

，此时 $j \ne k$。

综合可得：$S(j, k) = [j=k] \cdot n$ 。

求出$a_k$

将 $S(j, k)$ 带入原式：

$$c_k = \sum_{j=0}^{n-1}a_j \cdot S(j,k) = \sum_{j=0}^{n-1}a_j \cdot [j=k] \cdot n = a_k \cdot n$$

所以：$a_k = \frac{c_k}{n}$ ，最终我们得到了原多项式 $A(x)$ 的系数向量 $(a_0, a_1, \cdots, a_n)$ 中的 $a_k$。

小结

对于多项式 $A(x)$ 由插值节点 $(\omega_{n}^{0}, \omega_{n}^{1}, \omega_{n}^{2}, \cdots, \omega_{n}^{n-1})$ 做离散傅里叶变换得到的点值向量 $(d_0, d_1, \cdots, d_{n-1})$。

将 $(\omega_{n}^{0}, \omega_{n}^{-1}, \omega_{n}^{-2}, \cdots, \omega_{n}^{-(n-1)})$ 作为插值节点，$(d_0, d_1, \cdots, d_{n-1})$ 作为系数向量，做一次离散傅里叶变换得到的向量每一项都除以 $n$ 之后得到的 $(\frac{c_0}{n}, \frac{c_1}{n}, \cdots, \frac{c_{n-1}}{n})$ 就是多项式的系数向量 $(a_0, a_1, \cdots, a_{n-1})$。

注意：$\omega_{n}^{-k}$ 是 $\omega_{n}^{k}$ 的共轭复数。

通过这个过程我们就实现了将点值转换为系数表示。

FFT的本质

通过上面两节的 $\textit{DFT}$ 和 $\textit{IDFT}$，我们实现了魔法黑匣子的功能，也就是实现了多项式的点值表达式和系数表达式的互相转换。

下面进行总结^[8]：

从时域到频域：DFT

$\textit{DFT}$：

$$\textit{DFT}(a_0, a_1, \dots, a_{n-1}) = (y_0, y_1, \dots, y_{n-1}) = (A(w_n^0), A(w_n^1), \dots, A(w_n^{n-1})) = (A(w_n^0), A(w_n^1), \dots, A(w_n^{n-1}))$$

从时域到频域：IDFT

$\textit{IDFT}$ ：

$$\textit{InverseDFT}(y_0, y_1, \dots, y_{n-1}) = (a_0, a_1, \dots, a_{n-1})$$

矩阵运算

我们可以将 $A(x)$ 写成矩阵乘法的形式:

$$\begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{pmatrix} = \begin{pmatrix} w_n^0 & w_n^0 & w_n^0 & w_n^0 & \cdots & w_n^0 \\ w_n^0 & w_n^1 & w_n^2 & w_n^3 & \cdots & w_n^{n-1} \\ w_n^0 & w_n^2 & w_n^4 & w_n^6 & \cdots & w_n^{2(n-1)} \\ w_n^0 & w_n^3 & w_n^6 & w_n^9 & \cdots & w_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w_n^0 & w_n^{n-1} & w_n^{2(n-1)} & w_n^{3(n-1)} & \cdots & w_n^{(n-1)(n-1)} \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{pmatrix}$$

其中的 $W$ 矩阵叫做 Vandermonde Matrix 。

已知 $A(x)$ 点值表达式，也就是已知 $Y=(y_0, y_1, \cdots, y_{n-1})$ 和 $W$ ，求 $A(x)$ 的系数向量 $(a_0, a_1, \cdots, a_{n-1})$ ，$A = \frac{Y}{W} = Y \times W^{-1}$ 。

$$\begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} w_n^0 & w_n^0 & w_n^0 & w_n^0 & \cdots & w_n^0 \\ w_n^0 & w_n^1 & w_n^2 & w_n^3 & \cdots & w_n^{n-1} \\ w_n^0 & w_n^2 & w_n^4 & w_n^6 & \cdots & w_n^{2(n-1)} \\ w_n^0 & w_n^3 & w_n^6 & w_n^9 & \cdots & w_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w_n^0 & w_n^{n-1} & w_n^{2(n-1)} & w_n^{3(n-1)} & \cdots & w_n^{(n-1)(n-1)} \end{pmatrix}^{-1} \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{pmatrix}$$

$W$ 的 逆矩阵(Inverse Matrix) $W^{-1}$ 为：

$$\frac{1}{n} \begin{pmatrix} w_n^0 & w_n^0 & w_n^0 & w_n^0 & \cdots & w_n^0 \\ w_n^0 & w_n^{-1} & w_n^{-2} & w_n^{-3} & \cdots & w_n^{-(n-1)} \\ w_n^0 & w_n^{-2} & w_n^{-4} & w_n^{-6} & \cdots & w_n^{-2(n-1)} \\ w_n^0 & w_n^{-3} & w_n^{-6} & w_n^{-9} & \cdots & w_n^{-3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w_n^0 & w_n^{-(n-1)} & w_n^{-2(n-1)} & w_n^{-3(n-1)} & \cdots & w_n^{-(n-1)(n-1)} \end{pmatrix}$$

因此我们可以得到以下公式：

$$a_k = \frac{1}{n} \sum_{j=0}^{n-1} y_j w_n^{-k j}$$

$y_k$的获取：

$$y_k = \sum_{j=0}^{n-1} a_j w_n^{k j}$$

注意到我们只要对矩阵 $W$ 中每一个元素取共轭复数并除以 $n$ ，就得到了其 逆矩阵 。

FFT 代码实现

上面讲了这么多，下面开始编码实现吧:-)

Complex Number

复数可以使用C++ STL自带的 $\texttt{std::complex<T>}$ ，依照精度要求 $\texttt{T}$ 一般为 $\texttt{double}$ ，也可以自己封装实现^[9]：

struct Complex {
  double x, y;

  Complex(double _x = 0.0, double _y = 0.0) {
    x = _x;
    y = _y;
  }

  Complex operator-(const Complex &b) const {
    return Complex(x - b.x, y - b.y);
  }

  Complex operator+(const Complex &b) const {
    return Complex(x + b.x, y + b.y);
  }

  Complex operator*(const Complex &b) const {
    return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
  }
};

FFT 递归写法

我们可以按照我们刚才的推导，得到递归形式写法，代码如下所示：

const double PI = acos(-1.0);    // PI = arccos(-1)

/**
 * FFT 实现
 *
 * @param a
 * @param invert true means IFFT, else FFT
 * @return y
 */
vector<complex<double>> FFT(vector<complex<double>> &a, bool invert) {
    //第一个参数为一个多项式的系数, 以次数从小到大的顺序, 向量中每一项的实部为该项系数
    int n = a.size();

    // 如果当前多项式仅有常数项时直接返回多项式的值
    if (n == 1) {
        return a;
    }

    vector<complex<double>> Pe(n / 2), Po(n / 2); // 文中的Pe与Po的系数表示法

    for (int i = 0; 2 * i < n; i++) {
        Pe[i] = a[2 * i];
        Po[i] = a[2 * i + 1];
    }

    // Divide 分治
    // 递归求 ye = Pe(xi), yo = Po(xi)
    vector<complex<double>> ye = FFT(Pe, invert);
    vector<complex<double>> yo = FFT(Po, invert);

    // Combine
    vector<complex<double>> y(n);

    // Root of Units
    double ang = 2 * PI / n * (invert ? -1 : 1);
    complex<double> omega(cos(ang), sin(ang)); // omega为第一个n次复根,
    complex<double> curRoot(1, 0);   // curr为第零0个n次复根, 即为 1

    for (int i = 0; i < n / 2; i++) {
        y[i] = ye[i] + curRoot * yo[i];  // 求出P(xi)
        y[i + n / 2] = ye[i] - curRoot * yo[i]; // 由单位复根的性质可知第k个根与第k + n/2个根互为相反数
        curRoot *= omega;   // cur * omega得到下一个复根
    }

    return y;  // 返回最终的系数
}

由于 $\textit{FFT}$ 和 $\textit{IFFT}$ 操作流都是系数不一致，所以我们可以将其写成一个函数。值得注意的是，$\textit{IFFT}$ 之后需要都系数除于 $n$ 才是最终的结果。

复杂度分析

• 时间复杂度：$O(nlogn)$，其中 $n$ 为数组长度，每次需要 $O(n)$ ，总共 $O(logn)$ 递归。

• 空间复杂度：$O(n)$，我们需要额外的 $n$ 长度用于存储数据，函数递归栈需要空间是 $O(logn)$，所以总空间复杂度为：$O(n)$ + $O(logn)$ = $O(n)$。

总结

跨度一年，终于把这个 $\textit{FFT}$ 的坑填完了，虽然大学时学《信号与系统》时知道了 $\textit{FFT}$ ，但真正理解 $\textit{FFT}$ 还是要等到写完了这几篇文章。

下一篇讲聚焦于 $\textit{FFT}$ 的优化！

文章如有错误之处，敬请批评指正。

参考资料

---

1. Wiki: Fast Fourier Transform ↩︎

2. Wiki: 快速傅里叶变换 ↩︎

3. Wiki: 离散傅里叶变换 ↩︎

4. Wiki: 泰勒级数 ↩︎

5. 快速傅里叶变换算法 ↩︎

6. 一小时学会快速傅里叶变换(Fast Fourier Transform) ↩︎ ↩︎

7. Wiki: Convolution theorem ↩︎

8. Algorithm: Fast Fourier Transform ↩︎

9. 学习笔记 - 快速傅里叶变换 ↩︎