标签：even end 变换 快速 Mn bmatrix text 傅立叶 odd

已知 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)=Mn​v.

根据矩阵乘法公式，把 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,0​wn1,0​wnn−1,0​​wn0,1​wn1,1​wnn−1,1​​……⋮…​wn0,n−1​wn1,n−1​wnn−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=⎝⎜⎜⎜⎛​wn0​wn1​⋮wnn−1​​⎠⎟⎟⎟⎞​=⎝⎜⎜⎜⎛​w80​w81​⋮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​⋅wn0​w81,0​⋅wn1​w87,0​⋅wn7​​w80,1​⋅wn0​w81,1​⋅wn1​w87,1​⋅wn7​​……⋮…​w80,7​⋅wn0​w81,7​⋅wn1​w87,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/2​Mn/2​​Mn/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ˉ=[veven​vodd​​],

其中
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​=⎝⎜⎜⎜⎛​v0​v2​⋮vn−2​​⎠⎟⎟⎟⎞​,vodd​=⎝⎜⎜⎜⎛​v1​v3​⋮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)​=Mn​v=Mˉn​vˉ=[Mn/2​Mn/2​​Mn/2​×u−Mn/2​×u​][veven​vodd​​]=[Mn/2​veven​+Mn/2​vodd​×uMn/2​veven​−Mn/2​vodd​×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⋮0​0⋯⋯1​⋯00⋮0​001⋯​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). n1​Mn​Mn​v=⎝⎜⎜⎜⎜⎜⎛​v0​vn−1​vn−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​(Mn​v)]=v​
即，
F − 1 = 1 n σ n − 1 F . \mathcal{F}^{-1} = \frac{1}{n}\sigma_{n-1}\mathcal{F}. F−1=n1​σn−1​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 ) 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）的思想，对每个子问题递归求解，从而节省计算时间。

标签：even,end,变换,快速,Mn,bmatrix,text,傅立叶,odd
来源： https://blog.csdn.net/qx3501332/article/details/106397339

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。