快速傅立叶变换

已知 $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 ki}{n}},k=0,1,\dots ,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,\dots ,n-1$ 。

矩阵 $M_n$ 的定义如下：

$$M_n=\left[\begin{array}{cccc}{w}_n^{0,0}& w_n^{0,1}& \dots & w_n^{0,n-1}\\ w_n^{1,0}& w_n^{1,1}& \dots & w_n^{1,n-1}\\ & & ⋮& \\ w_n^{n-1,0}& w_n^{n-1,1}& \dots & w_n^{n-1,n-1}\end{array}\right]$$

下文假设 $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,u=\left(\begin{array}{c}{w}_n^0\\ w_n^1\\ ⋮\\ w_n^{n-1}\end{array}\right)=\left(\begin{array}{c}{w}_8^0\\ w_8^1\\ ⋮\\ w_8^7\end{array}\right),$$
其中 $\times$ 代表把 $M_8$ 中的每一列，与 $u$ 按元素逐个相乘，即
$$M_8\times u=\left[\begin{array}{cccc}{w}_8^{0,0}\cdot w_n^0& w_8^{0,1}\cdot w_n^0& \dots & w_8^{0,7}\cdot w_n^0\\ w_8^{1,0}\cdot w_n^1& w_8^{1,1}\cdot w_n^1& \dots & w_8^{1,7}\cdot w_n^1\\ & & ⋮& \\ w_8^{7,0}\cdot w_n^7& w_8^{7,1}\cdot w_n^7& \dots & w_8^{7,7}\cdot w_n^7\end{array}\right]$$
3、紫色矩阵与红色矩阵符号相反，即 $-M_8\times u$。 即换句话说，它相当于把红色向量顺时针旋转180度。

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

先把 $M_n$ 的列重排，即，偶数列放在一起，奇数列放在一起。新的矩阵记作 ${\bar{M}}_n$，我们有
$${\bar{M}}_n=\left[\begin{array}{cc}{M}_{n/2}& M_{n/2}\times u\\ M_{n/2}& -M_{n/2}\times u\end{array}\right].$$

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

新的向量记作 ${\bar{v}}_n$，我们有
$$\bar{v}=\left[\begin{array}{c}{v}_{\text{even}}\\ v_{\text{odd}}\end{array}\right],$$

其中
$$v_{\text{even}}=\left(\begin{array}{c}{v}_0\\ v_2\\ ⋮\\ v_{n-2}\end{array}\right),v_{\text{odd}}=\left(\begin{array}{c}{v}_1\\ v_3\\ ⋮\\ v_{n-1}\end{array}\right).$$

我们有
$$\begin{array}{rl}\mathcal{F}(v)& =M_{n}v={\bar{M}}_n\bar{v}\\ & =\left[\begin{array}{cc}{M}_{n/2}& M_{n/2}\times u\\ M_{n/2}& -M_{n/2}\times u\end{array}\right]\left[\begin{array}{c}{v}_{\text{even}}\\ v_{\text{odd}}\end{array}\right]\\ & =\left[\begin{array}{c}{M}_{n/2}{v}_{\text{even}}+M_{n/2}{v}_{\text{odd}}\times u\\ M_{n/2}{v}_{\text{even}}-M_{n/2}{v}_{\text{odd}}\times u\end{array}\right]\end{array}$$
即，
$$\mathcal{F}(v)=\left[\begin{array}{c}\mathcal{F}(v_{\text{even}})+\mathcal{F}(v_{\text{odd}})\times u\\ \mathcal{F}(v_{\text{even}})-\mathcal{F}(v_{\text{odd}})\times u\end{array}\right].$$
这就是快速傅立叶变换算法，可以用递归的方式实现（代码见下文）。

下面分析计算复杂度。

令 $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 = [v[i] for i in range(0, n, 2)]
    v_odd = [v[i] for i in range(1, n, 2)]
    q1 = fft(v_even)
    q2 = fft(v_odd)
    s = [np.e ** (-2 * np.pi * complex(0, 1) * j / n) for j in range(n//2)]
    part1 = (np.array(q1) + np.array(q2) * s).tolist()
    part2 = (np.array(q1) - np.array(q2) * s).tolist()

    return part1 + part2

逆变换

再看看它的逆变换。

计算 $M_n\cdot M_n$，我们得到
$$M_n\cdot M_n=\left[\begin{array}{ccccc}1& 0& \cdots & 0& 0\\ 0& \cdots & 0& 0& 1\\ 0& \cdots & 0& 1& 0\\ ⋮& & ⋮& & \\ 0& 1& 0& \cdots & 0\end{array}\right]\times n$$
我们有，
$$\frac{1}{n}{M}_n{M}_{n}v=\left(\begin{array}{c}{v}_0\\ v_{n-1}\\ v_{n-2}\\ ⋮\\ v_1\end{array}\right)=:{\sigma }_{n-1}(v).$$
其中 ${\sigma }_{n-1}(v)$ 代表对向量 $v$ 的后 $n-1$ 个元素按下标逆序排列。

注意到 ${\sigma }_{n-1}({\sigma }_{n-1}(v))=v$，我们有
$$\begin{array}{r}\frac{1}{n}{\sigma }_{n-1}[M_n(M_{n}v)]=v\end{array}$$
即，
$${\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）的思想，对每个子问题递归求解，从而节省计算时间。