快速傅立叶变换

已知 $n$ 维向量 $v\in\mathbb{C}^n$，它的 离散傅立叶变换 $\mathcal{F}(v)$ 是一个 $n\times n$ 矩阵 $M_n$ 与它的乘积： $$ \mathcal{F}(v) = M_n v. $$

根据矩阵乘法公式，把 $M_n$ 的每一行与 $v$ 分别做向量积，即可得到结果。每一行做向量积耗时 $O(n)$，一共 $n$ 行，所以用矩阵乘法计算的时间复杂度为 $O(n^2)$。

快速傅立叶变换算法可以达到 $O(n\log n)$ 的时间复杂度（用 Python 实现，仅需十几行代码）。

快速傅立叶变换之所以比矩阵乘法高效，是因为 $M_n$ 有特殊的结构。这种结构，可以采用分而治之的思想（分治法），来提升计算效率。

离散傅立叶变换

下面介绍矩阵 $M_n$ 的定义。

考虑复数域上的 $n$ 次单位根 $w$，即 $w^n = 1$，我们有 $$ w = e^{\frac{2\pi k i}{n}}, \quad k = 0, 1, \ldots, n-1. $$ 令 $w_n$ 代表第 $n$ 个单位根，即 $$ w_n = e^{\frac{2\pi(n-1)i}{n}} = e^{-\frac{2\pi i}{n}}. $$ 令 $w_n^{j,k}$ 代表 $w_n$ 的 $j\times k$ 次方，其中 $j, k$ 分别代表行和列的下标。 需要注意，下标从0开始计数，即 $j, k = 0,1,\ldots, n-1$ 。

矩阵 $M_n$ 的定义如下：

$$ M_n = \begin{bmatrix} w_n^{0,0} & w_n^{0,1} & \ldots & w_n^{0,n-1}\newline w_n^{1,0} & w_n^{1,1} & \ldots & w_n^{1,n-1}\newline & & \vdots &\newline w_n^{n-1,0} & w_n^{n-1,1} & \ldots & w_n^{n-1,n-1} \end{bmatrix} $$

下文假设 $n$ 是2的幂次方，即存在整数 $k > 0$ 使得 $n = 2^k$。

快速傅立叶变换

我们知道 $$ w_n = e^{-\frac{2\pi i}{n}} = \cos(\frac{2\pi}{n}) - i\sin(\frac{2\pi}{n}), $$ 它是复平面上的一个单位向量（如下所示）：

令 $\theta = \frac{2\pi}{n}$，那么 $w_n^k$ 相当于把 $w_n$ 沿着顺时针方向旋转 $(k-1)\cdot\theta$ （如下所示）：

有了这样的认知，我们可以对 $M_n$ 进行可视化。以 $M_{16}$ 为例，如下图所示：

一咋看有些混乱，下面把它分成四个子矩阵，按颜色区分：

仔细观察，我们发现：

1、黑色矩阵和蓝色矩阵相同，它们实际上是 $M_8$；

2、每个红色向量，相当于把它“左边”的黑色向量，顺时针方向旋转，其中旋转角度的大小取决于向量所在的行。

设黑色元素为 $w_n^{ik}$，那么它右边的红色元素为 $w_n^{ik}\cdot w_n^i$。这样一来，红色矩阵可以表示成 $$ M_8 \times u, \quad u= \begin{pmatrix} w_n^0\newline w_n^1\newline \vdots\newline w_n^{n-1} \end{pmatrix} = \begin{pmatrix} w^0_8\newline w^1_8\newline \vdots\newline w^7_8 \end{pmatrix}, $$ 其中 $\times$ 代表把 $M_8$ 中的每一列，与 $u$ 按元素逐个相乘，即 $$ M_8\times u = \begin{bmatrix} w_8^{0,0}\cdot w_n^0 & w_8^{0,1}\cdot w_n^0 & \ldots & w_8^{0,7}\cdot w_n^0\newline w_8^{1,0}\cdot w_n^1 & w_8^{1,1} \cdot w_n^1 & \ldots & w_8^{1,7}\cdot w_n^1\newline & & \vdots &\newline w_8^{7,0} \cdot w_n^7& w_8^{7,1} \cdot w_n^7& \ldots & w_8^{7,7}\cdot w_n^7 \end{bmatrix} $$ 3、紫色矩阵与红色矩阵符号相反，即 $-M_8\times u$。 即换句话说，它相当于把红色向量顺时针旋转180度。

有了上面的观察，我们即将得到快速傅立叶变换的计算方法。

第一步，先把 $M_n$ 的列重排，即，偶数列放在一起，奇数列放在一起。新的矩阵记作 $\bar{M}_n$，我们有

$$ \bar{M_n} = \begin{bmatrix} M_{n/2} & M_{n/2}\times u \newline M_{n/2} & -M_{n/2}\times u \end{bmatrix}. $$

第二步，为了不改变计算结果，还要把 $v$ 中的元素按下标重新排列，即，偶数位和奇数位分别单独放一起。

新的向量记作 $\bar{v}_n$，我们有

$$ \bar{v} = \begin{bmatrix} v_{\text{even}} \newline v_{\text{odd}} \end{bmatrix}, $$

其中 $$ v_{\text{even}} = \begin{pmatrix} v_0\newline v_2\newline \vdots\newline v_{n-2} \end{pmatrix}, \quad v_{\text{odd}} = \begin{pmatrix} v_1\newline v_3\newline \vdots\newline v_{n-1} \end{pmatrix}. $$

第三步，计算 $\bar{M}_n$ 和 $\bar{v}$ 的乘积。

$$ \begin{aligned} \mathcal{F}(v) & = M_n v = \bar{M_n} \bar{v} \newline & = \begin{bmatrix} M_{n/2} & M_{n/2}\times u\newline M_{n/2} & -M_{n/2}\times u \end{bmatrix}\begin{bmatrix} v_{\text{even}}\newline v_{\text{odd}} \end{bmatrix}\newline & = \begin{bmatrix}M_{n/2}v_{\text{even}} + M_{n/2}v_{\text{odd}}\times u \newline M_{n/2}v_{\text{even}} -M_{n/2}v_{\text{odd}}\times u \end{bmatrix} \end{aligned} $$

于是得到下面的递归表达式。 $$ \mathcal{F}(v) = \begin{bmatrix} \mathcal{F}(v_{\text{even}}) + \mathcal{F}({v_{\text{odd}}})\times u \newline \mathcal{F}(v_{\text{even}}) - \mathcal{F}({v_{\text{odd}}})\times u \newline \end{bmatrix}. $$ 这就是快速傅立叶变换算法，可以用递归的方式实现（代码见下文）。

下面分析计算复杂度。

令 $T(n)$ 代表向量长度为 $n$ 的离散傅立叶变换的计算时间。根据上述公式，计算 $\mathcal{F}(v_{\text{even}})$ 和 $\mathcal{F}(v_{\text{odd}})$ 分别耗时 $T(n/2)$，拼接耗时 $O(n)$。我们有 $$ T(n) = 2T(n/2) + O(n). $$ 根据主定理（Master theorem），我们得到 $T(n) = O(n\log n)$。

算法实现

下面用 Python 实现快速傅立叶变换。

def FFT(v):
    """ 快速傅立叶变换
    注意：输入向量 v 的长度是2的幂次方。
    :param v: array like
    :return: array like（复向量）
    """
    n = len(v)
    if n == 1:
        return v
    v_even = FFT(v[0:n:2])
    v_odd = FFT(v[1:n:2])
    u = [w_n**i for i in range(n//2)]
    part1 = [v_even[i] + u[i]*v_odd[i] for i in range(n//2)]
    part2 = [v_even[i] - u[i]*v_odd[i] for i in range(n//2)]
    return part1 + part2

逆变换

再看看傅立叶变换的逆变换。

计算 $M_n \cdot M_n$，我们得到

$$ M_n \cdot M_n = \begin{bmatrix} 1 & 0 & \cdots & 0 & 0\ 0 & \cdots & 0 & 0 & 1 \ 0 & \cdots & 0 & 1 & 0 \ \vdots & & \vdots & & \ 0 & 1 & 0 & \cdots & 0 \end{bmatrix} \times n $$ 我们有，

$$ \frac{1}{n}M_nM_n v = \begin{pmatrix} v_0\ v_{n-1}\ v_{n-2}\ \vdots\ v_1 \end{pmatrix} =: \sigma_{n-1}(v). $$

其中 $\sigma_{n-1}(v)$ 代表对向量 $v$ 的后 $n-1$ 个元素按下标逆序排列。

注意到 $\sigma_{n-1}(\sigma_{n-1}(v)) = v$，我们有 $$ \begin{aligned} \frac{1}{n}\sigma_{n-1}[M_n(M_n v)] = v \end{aligned} $$ 即，

$$ \mathcal{F}^{-1} = \frac{1}{n}\sigma_{n-1}\mathcal{F}. $$ 下面是 python 实现。

def iFFT(v):
    """
    傅立叶变换的逆变换
    注意：输入向量 v 的长度是2的幂次方。
    :param v: array like
    :return: array like（复向量）
    """
    n = len(v)
    u = FFT(v)
    x = [i/n for i in u]  # u 除以 n
    # 对后n-1个元素逆序排列
    y = x[1:]
    y.reverse()  
    return x[0:1] + y

快速傅立叶变换的时间复杂度是 $O(n\log n)$，比直接用矩阵乘更高效，原因在于矩阵 $M_n$ 具有特殊结构。通过观察发现，$M_n$ 可以被拆解成四个与 $M_{n/2}$ 有关联的子矩阵。基于分治法（Divide and Conquer）的思想，对每个子问题递归求解，从而节省计算时间。