快速傅立叶变换

2023-10-17 11:17:04

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

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

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

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

离散傅立叶变换

下面介绍矩阵 M n M_n Mn 的定义。

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

矩阵 M n M_n Mn 的定义如下：

M n = [ w n 0 , 0 w n 0 , 1 … w n 0 , n − 1 w n 1 , 0 w n 1 , 1 … w n 1 , n − 1 ⋮ w n n − 1 , 0 w n n − 1 , 1 … w n n − 1 , n − 1 ] M_n = \begin{bmatrix} w_n^{0,0} & w_n^{0,1} & \ldots & w_n^{0,n-1}\\ w_n^{1,0} & w_n^{1,1} & \ldots & w_n^{1,n-1}\\ & & \vdots &\\ w_n^{n-1,0} & w_n^{n-1,1} & \ldots & w_n^{n-1,n-1} \end{bmatrix} Mn=⎣⎢⎢⎢⎡wn0,0wn1,0wnn−1,0wn0,1wn1,1wnn−1,1……⋮…wn0,n−1wn1,n−1wnn−1,n−1⎦⎥⎥⎥⎤

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

快速傅立叶变换

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

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

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

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

仔细观察，我们发现：

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

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

设黑色元素为 w n i k w_n^{ik} wnik，那么它右边的红色元素为 w n i k ⋅ w n i w_n^{ik}\cdot w_n^i wnik⋅wni。这样一来，红色矩阵可以表示成
M 8 × u , u = ( w n 0 w n 1 ⋮ w n n − 1 ) = ( w 8 0 w 8 1 ⋮ w 8 7 ) , M_8 \times u, \quad u= \begin{pmatrix} w_n^0\\ w_n^1\\ \vdots\\ w_n^{n-1} \end{pmatrix} = \begin{pmatrix} w^0_8\\ w^1_8\\ \vdots\\ w^7_8 \end{pmatrix}, M8×u,u=⎝⎜⎜⎜⎛wn0wn1⋮wnn−1⎠⎟⎟⎟⎞=⎝⎜⎜⎜⎛w80w81⋮w87⎠⎟⎟⎟⎞,
其中 × \times × 代表把 M 8 M_8 M8 中的每一列，与 u u u 按元素逐个相乘，即
M 8 × u = [ w 8 0 , 0 ⋅ w n 0 w 8 0 , 1 ⋅ w n 0 … w 8 0 , 7 ⋅ w n 0 w 8 1 , 0 ⋅ w n 1 w 8 1 , 1 ⋅ w n 1 … w 8 1 , 7 ⋅ w n 1 ⋮ w 8 7 , 0 ⋅ w n 7 w 8 7 , 1 ⋅ w n 7 … w 8 7 , 7 ⋅ w n 7 ] 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\\ 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\\ & & \vdots &\\ 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} M8×u=⎣⎢⎢⎢⎡w80,0⋅wn0w81,0⋅wn1w87,0⋅wn7w80,1⋅wn0w81,1⋅wn1w87,1⋅wn7……⋮…w80,7⋅wn0w81,7⋅wn1w87,7⋅wn7⎦⎥⎥⎥⎤
3、紫色矩阵与红色矩阵符号相反，即 − M 8 × u -M_8\times u −M8×u。即换句话说，它相当于把红色向量顺时针旋转180度。

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

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

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

新的向量记作 v ˉ n \bar{v}_n vˉn，我们有
v ˉ = [ v even v odd ] , \bar{v} = \begin{bmatrix} v_{\text{even}}\\ v_{\text{odd}} \end{bmatrix}, vˉ=[vevenvodd],

其中
v even = ( v 0 v 2 ⋮ v n − 2 ) , v odd = ( v 1 v 3 ⋮ v n − 1 ) . v_{\text{even}} = \begin{pmatrix} v_0\\ v_2\\ \vdots\\ v_{n-2} \end{pmatrix}, \quad v_{\text{odd}} = \begin{pmatrix} v_1\\ v_3\\ \vdots\\ v_{n-1} \end{pmatrix}. veven=⎝⎜⎜⎜⎛v0v2⋮vn−2⎠⎟⎟⎟⎞,vodd=⎝⎜⎜⎜⎛v1v3⋮vn−1⎠⎟⎟⎟⎞.

我们有
F ( v ) = M n v = M ˉ n v ˉ = [ M n / 2 M n / 2 × u M n / 2 − M n / 2 × u ] [ v even v odd ] = [ M n / 2 v even + M n / 2 v odd × u M n / 2 v even − M n / 2 v odd × u ] \begin{aligned} \mathcal{F}(v) & = M_n v = \bar{M}_n \bar{v}\\ & = \begin{bmatrix} M_{n/2} & M_{n/2}\times u\\ M_{n/2} & -M_{n/2}\times u \end{bmatrix}\begin{bmatrix} v_{\text{even}}\\ v_{\text{odd}} \end{bmatrix}\\ & = \begin{bmatrix}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{bmatrix} \end{aligned} F(v)=Mnv=Mˉnvˉ=[Mn/2Mn/2Mn/2×u−Mn/2×u][vevenvodd]=[Mn/2veven+Mn/2vodd×uMn/2veven−Mn/2vodd×u]
即，
F ( v ) = [ F ( v even ) + F ( v odd ) × u F ( v even ) − F ( v odd ) × u ] . \mathcal{F}(v) = \begin{bmatrix} \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{bmatrix}. F(v)=[F(veven)+F(vodd)×uF(veven)−F(vodd)×u].
这就是快速傅立叶变换算法，可以用递归的方式实现（代码见下文）。

下面分析计算复杂度。

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

算法实现

下面用 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 ⋅ M n M_n \cdot M_n Mn⋅Mn，我们得到
M n ⋅ M n = [ 1 0 ⋯ 0 0 0 ⋯ 0 0 1 0 ⋯ 0 1 0 ⋮ ⋮ 0 1 0 ⋯ 0 ] × 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 Mn⋅Mn=⎣⎢⎢⎢⎢⎢⎡100⋮00⋯⋯1⋯00⋮0001⋯0100⎦⎥⎥⎥⎥⎥⎤×n
我们有，
1 n M n M n v = ( v 0 v n − 1 v n − 2 ⋮ v 1 ) = : σ n − 1 ( v ) . \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). n1MnMnv=⎝⎜⎜⎜⎜⎜⎛v0vn−1vn−2⋮v1⎠⎟⎟⎟⎟⎟⎞=:σn−1(v).
其中 σ n − 1 ( v ) \sigma_{n-1}(v) σn−1(v) 代表对向量 v v v 的后 n − 1 n-1 n−1 个元素按下标逆序排列。

注意到 σ n − 1 ( σ n − 1 ( v ) ) = v \sigma_{n-1}(\sigma_{n-1}(v)) = v σn−1(σn−1(v))=v，我们有
1 n σ n − 1 [ M n ( M n v ) ] = v \begin{aligned} \frac{1}{n}\sigma_{n-1}[M_n(M_n v)] = v \end{aligned} n1σn−1[Mn(Mnv)]=v
即，
F − 1 = 1 n σ n − 1 F . \mathcal{F}^{-1} = \frac{1}{n}\sigma_{n-1}\mathcal{F}. F−1=n1σn−1F.
下面是 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 ) O(n\log n) O(nlogn)，比直接用矩阵乘更高效，原因在于矩阵 M n M_n Mn 具有特殊结构。

通过观察发现， M n M_n Mn 可以被拆解成四个与 M n / 2 M_{n/2} Mn/2 有关联的子矩阵。基于分治法（Divide and Conquer）的思想，对每个子问题递归求解，从而节省计算时间。

码农公寓

离散傅立叶变换

快速傅立叶变换

算法实现

逆变换

总结

相关文章