12/20/2018 重写。

概述

快速傅里叶变换(Fast Fourier Transform,FFT)是一种可在 $\Theta(n \log n)$ 时间内完成的离散傅里叶变换(Discrete Fourier transform,DFT)算法。

这篇文章将简单地讲解它的工作原理,以及使用 C++ 的代码实现。

基本概念

多项式

在数学中,由若干个单项式相加组成的代数式叫做多项式。多项式中的每个单项式叫做多项式的项,这些单项式中的最高项次数,就是这个多项式的次数。

系数表示

设 $A(x)$ 表示一个 $n-1$ 次多项式,所有项系数组成的 $n$ 维向量 $(a_0,a_1,\cdots,a_{n-1})$ 唯一确定了这个多项式。

$$\begin{aligned} A(x)=& \sum_{i=1}^{n-1}a_ix^i \\ =& a_0+a_1x+a_2x^2\cdots+a_{n-1}x^{n-1}\end{aligned}$$

点值表示

将一组互不相同的 $n$ 个 $x$ 带入多项式,得到的 $n$ 个值。设它们组成的 n 维向量  $(x_0,x_1,\cdots,x_{n-1})$,$(y_0,y_1,\cdots,y_{n-1})$ 唯一确定了这个多项式。

$$\begin{aligned} y_i=& A(x_i) \\ =&\sum_{j=1}^{n-1}a_j x_i^j\end{aligned}$$

多项式乘法

令 $A(x)=\sum_{i=0}^{n-1}a_i x^i$,$B(x)=\sum_{i=0}^{n-1}b_i x^i$,有

$$C(x)=A(x)\times B(x) = \sum_{k=0}^{2n-2}(\sum_{k=i+j}a_i b_j)x^k$$

两个 $n-1$ 次多项式相乘,得到的是一个 $2n-2$ 次多项式,朴素算法的时间复杂度为 $\Theta(n^2)$。

如果用两个多项式在 $2n-1$ 个点取得的点值表示,有

$$C_{y_i}=(\sum_{j=0}^{2n-1}a_j x_i^j)\times(\sum_{j=0}^{2n-1}b_j x_i^j)=A_{y_i}\times B_{y_i}$$

已知 $A$ 和 $B$ 点值表示,求出 $C$ 的点值表示,时间复杂度为 $\Theta(n)$。

复数

复数,为实数的延伸,它使任一多项式方程都有根。复数当中有个“虚数单位” $i$,它是 $-1$ 的一个平方根,即 $i^{2}=-1$。任一复数都可表达为 $x+y_i$,其中 $x$ 及 $y$ 皆为实数,分别称为复数之“实部”和“虚部”。

复平面

数学中,复平面(complex plane)是用水平的实轴与垂直的虚轴建立起来的复数的几何表示。它可视为一个具有特定代数结构笛卡儿平面(实平面),一个复数的实部用沿着 x-轴的位移表示,虚部用沿着 y-轴的位移表示。

复平面的想法提供了一个复数的几何解释。在加法下,它们像向量一样相加;两个复数的乘法在极坐标下的表示最简单——乘积的长度或模长是两个绝对值或模长的乘积,乘积的角度或辐角是两个角度或辐角的和。特别地,用一个模长为 1 的复数相乘即为一个旋转。

单位根

在复平面上,以原点为圆心,1 为半径作得的圆称为单位圆。

数学上,$n$ 次单位根是 $n$ 次幂为 $1$ 的复数。它们位于复平面的单位圆上,构成正 $n$ 边形的顶点,其中一个顶点是 $1$。设所得的幅角为正且最小的单位根为 $\omega_n$,其余的记为 $\omega_n^k$,有

$$\omega_n^k = \cos k \dfrac{2\pi}{n}+i \sin k \dfrac{2\pi}{n}$$

其中 $\omega_n^0 = \omega_n^n = 1$。

性质

$$\begin{aligned} \omega_{2n}^{2k}=&\omega_n^k \\ \omega_n^{k+\frac{n}{2}}=& -\omega_n^k \end{aligned}$$

欧拉公式

$$e^{ix}=\cos x+i\sin x$$

单位根的代码实现

由欧拉公式,求 $n$ 次单位根 $\omega_n$,有

其中,$\text{exp}$ 用于计算 $e$ 的幂,$\text{complex}$ 是复数类。

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

考虑多项式 $A(x)$ 的表示。将 $n$ 次单位根的 $0$ 到 $n – 1$ 次幂带入多项式的系数表示,所得点值向量 $(A(\omega_n^0),A(\omega_n^1),\cdots,A(\omega_n^{n-1}))$ 称为其系数向量 $(a_0,a_1,\cdots,a_{n-1})$ 的离散傅里叶变换。

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

考虑将多项式按照系数下标的奇偶分成两部分

$$A(x)=(a_0 + a_2 x^2 + a_4 x^4 + \cdots) + (a_1 x + a_3 x^3 + a_5 x^5 + \cdots) $$

$$\begin{aligned}A_1(x)=&a_0+a_2x+a_4x^2+\cdots \\ A_2(x)=& a_1 + a_3 x + a_5 x^2 + \cdots \end{aligned}$$

则有

$$A(x)=A_1(x^2)+A_2(x^2)$$

设 $k< \dfrac{n}{2}$,有

$$\begin{aligned} A(\omega_n ^ k) &= A_1(\omega_n ^ {2k}) + \omega_n ^ {k} A_2(\omega_n ^ {2k}) \\ &= A_1(\omega_{\frac{n}{2}} ^ {k}) + \omega_n ^ {k} A_2(\omega_{\frac{n}{2}} ^ {k}) \ \end{aligned} $$

$$\begin{aligned} A(\omega_n ^ {k + \frac{n}{2}}) &= A_1(\omega_n ^ {2k + n}) + \omega_n ^ {k + \frac{n}{2}} A_2(\omega_n ^ {2k + n}) \\ &= A_1(\omega_n ^ {2k} \times \omega_n ^ n)-\omega_n ^ {k} A_2(\omega_n ^ {2k} \times \omega_n ^ n) \\ &= A_1(\omega_n ^ {2k})-\omega_n ^ {k} A_2(\omega_n ^ {2k}) \\ &= A_1(\omega_{\frac{n}{2}} ^ {k})-\omega_n ^ {k} A_2(\omega_{\frac{n}{2}} ^ {k}) \ \end{aligned}$$

也就是说,已知 $A_1(x)$ 和 $A_2(x)$ 在 $\omega_{ \frac{n}{2} } ^ 0,\ \omega_{ \frac{n}{2} } ^ 1,\ \dots,\ \omega_{ \frac{n}{2} } ^ { \frac{n}{2} – 1 }$ 的点值表示,就可以在 $\Theta(n)$ 的时间内求得 $A(x)$ 在 $ \omega_n ^ 0,\ \omega_n ^ 1,\ \dots,\ \omega_n ^ {n – 1} $ 处的点值表示。而关于 $A_1(x)$ 和 $A_2(x)$ 的问题都是相对于 $A(x)$ 规模缩小至 $\dfrac{2}{1}$ 的子问题,分治的边界为常数 $a_i$。这个分治算法在 $\Theta(n\log n)$ 的时间内完成了 DDF。

傅里叶逆变换 (Inverse Fourier Transform)

将点值表示的多项式转化为系数表示,同样可以使用快速傅里叶变换,这个过程叫做傅里叶逆变换

设 $ (y_0,\ y_1,\ y_2,\ \dots,\ y_{n – 1}) $ 为 $ (a_0,\ a_1,\ a_2,\ \dots,\ a_{n – 1}) $ 的傅里叶变换。考虑另一个向量 $ (c_0,\ c_1,\ c_2,\ \dots,\ c_{n – 1}) $ 满足

$$ c_k = \sum\limits_{i = 0} ^ {n – 1} y_i (\omega_n ^ {-k}) ^ i $$

即多项式 $ B(x) = y_0 + y_1 x + y_2 x ^ 2 + \dots + y_{n – 1} x ^ {n – 1} $ 在 $\omega_n ^ 0,\ \omega_n ^ {-1},\ \omega_n ^ {-2},\ \dots,\ \omega_n ^ {-(n – 1)} $ 处的点值表示。

将上式展开,得

$$\begin{aligned} c_k &= \sum\limits_{i = 0} ^ {n – 1} y_i (\omega_n ^ {-k}) ^ i \\ &= \sum\limits_{i = 0} ^ {n – 1} (\sum\limits_{j = 0} ^ {n – 1} a_j (\omega_n ^ i) ^ j) (\omega_n ^ {-k}) ^ i \\ &= \sum\limits_{i = 0} ^ {n – 1} (\sum\limits_{j = 0} ^ {n – 1} a_j (\omega_n ^ j) ^ i) (\omega_n ^ {-k}) ^ i \\ &= \sum\limits_{i = 0} ^ {n – 1} (\sum\limits_{j = 0} ^ {n – 1} a_j (\omega_n ^ j) ^ i (\omega_n ^ {-k}) ^ i) \\ &= \sum\limits_{i = 0} ^ {n – 1} \sum\limits_{j = 0} ^ {n – 1} a_j (\omega_n ^ {j – k}) ^ i \\ &= \sum\limits_{j = 0} ^ {n – 1} a_j (\sum\limits_{i = 0} ^ {n – 1} (\omega_n ^ {j – k}) ^ i) \end{aligned} $$

对于 $ \sum\limits_{i = 0} ^ {n – 1} (\omega_n ^ {j – k}) ^ i $,当 $k \ne j$ 时,由单位根性质可知(或者用等比数列求和),它的值为 $0$;当 $k = j$ 时,显然为 $n$。

所以,当且仅当 $j=k$ 时,$a_j$ 对 $c_k$ 有 $na_j$ 的贡献,即

$$ \begin{aligned} c_i &= n a_i \\ a_i &= \frac{1}{n} c_i \end{aligned} $$

所以,使用单位根的倒数代替单位根,做一次类似快速傅里叶变换的过程,再将结果每个数除以 $n$,即为傅里叶逆变换的结果。

代码实现

预处理下标

直接分治计算的效率较低,我们使用迭代完成。

具体地,我们以 $n=7$ 为例,寻求分治过程中下标的变化

抵达分治边界时的下标,恰好是原下标的二进制翻转。我们可以预处理出下标的对应关系。

蝴蝶操作

考虑合并两个子问题的过程,假设 $A_1(\omega_{ \frac{n}{2} } ^k)$ 和 $A_2(\omega_{ \frac{n}{2} } ^k)$ 分别存在 $a(k)$ 和 $a(\frac{n}{2} + k)$ 中,$A(\omega_n ^ {k})$ 和 $A(\omega_n ^ {k + \frac{n}{2} })$ 将要被存放在 $b(k)$ 和 $b(\frac{n}{2} + k)$ 中,合并的操作可表示为

$$\begin{aligned} b(k) & \leftarrow a(k) + \omega_n ^ k \times a(\frac{n}{2} + k) \\ b(\frac{n}{2} + k) & \leftarrow a(k)-\omega_n ^ k \times a(\frac{n}{2} + k) \ \end{aligned} $$

蝴蝶操作的意义是,引入临时变量,使得合并在原地进行,无需新增一个数组

$$ \begin{aligned} t & \leftarrow \omega_n ^ k \times a(\frac{n}{2} + k) \\ a(\frac{n}{2}+k) & \leftarrow a(k)-t \\ a(k) & \leftarrow a(k)+t \end{aligned} $$

代码

这份代码里,我们利用 $FFT$ 求出了 $A(x)$ 和 $B(x)$ 在 $s$ 次单位根处的点值表示,$\Theta(n)$ 求出卷积的点值表示,然后再次 $FFT$ 求得卷积的系数表示。

感谢 MenciWikipedia 的珍贵资料。

  1. oyyj603470138
    Sep 11, 2018

    您应该从群上调和分析(DFT)讲起。。

    Reply
    • xgzepto
      Sep 11, 2018

      所以我坑了

      Reply