FFT离散的快速傅里叶变换及其实现

FFT快速傅里叶变换（离散型）

先说个比较慢的版本便于理解，下次给出迭代的FFT，基于一个叫蝴蝶变换的东西。

读者所需预备知识

1. 会写出多项式及其相关概念
2. 知道什么是多项式的系数表示与点值表示
3. 知道啥叫复数

多项式的表示

• 多项式：$P(x)={\sum }_{i=0}^{n-1}{a}_i\ast x^i$

• 多项式的次数界：大于最高次数的数字，例如说上面的多项式的次数界为$n$

• 霍纳法则：$P(x)={\sum }_{i=0}^{n-1}{a}_i{x}_i=a_0+x(a_1+x(a_2+...+x(a_{n-1})))$，有了霍纳法则我们可以在$\theta (n)$的时间内求出多项式的值。

• 点值表示: { ( x 0 , y 0 ) , ( x 1 , y 1 ) , … , ( x n − 1 , y n − 1 ) } \{(x_0,y_0),(x_1,y_1),\dots,(x_{n-1},y_{n-1})\} {(x0​,y0​),(x1​,y1​),…,(xn−1​,yn−1​)}或 { ( x 0 , P ( x 0 ) ) , ( x 1 , P ( x 1 ) ) , … , ( x n − 1 , P ( x n − 1 ) ) } \{(x_0,P(x_0)),(x_1,P(x_1)),\dots,(x_{n-1},P(x_{n-1}))\} {(x0​,P(x0​)),(x1​,P(x1​)),…,(xn−1​,P(xn−1​))}

• 我们能求出一个多项式的$n$​个点值表达，那么我们能够唯一对应一个多项式系数表达式证明，反之同理（插值）。（即点值和系数存在一一对应关系，可以互相转换）

DFT与FFT

在聊DFT之前先看个很重要的东西，必须要理解。

单位复数根

我们断言，如果给多项式$C(x)$​​和$Q(x)$​​​代入$n$​​个特殊的点，这些点叫做n次单位复数根或简称n次单位根，那么可以在$\theta (nlgn)$​的时间内求出$P(x)$​的点值表式。

啥是n次单位根？

• $n$​​​次单位根是满足${\omega }^n=1$​​的复数$\omega$​​​​。注意$\omega$是一个复数，也就是$a+ib$的形式的东西。

• 复数要表示在复平面上，如图这就是一个3次单位根（每个箭头都是一个），其中${\omega }_3^0={\omega }_3^3=1$​​，上边这个圆叫单位圆
• 回忆一下复数的乘法（满足平行四边形法则），是不是辐角相加，模长相乘。在单位圆上由于模长都是1，所以两个单位根相乘，就相当于在把那个箭头做旋转。
  • 例如${\omega }_3^2\ast {\omega }_3^1={\omega }_3^3$,就相当于把第二象限的那个箭头旋转到$\mathbb{R}$轴正方向上
• 引入欧拉公式如何理解：$e^{ix}=\cos (x)+i\ast \sin (x)$,将$x=2\pi$（也就是一个圆周）代入得到$e^{2\pi i}=1$
• 将圆周分成$n$​​份，发现这个玩意恰好和单位根定义一样，所以有$\omega =e^{\frac{2\pi i}{n}}$​,为了方便区分我们把圆分成了几份我们这样记$w_n=w$​

证明几个性质

1. (消去引理) ${\omega }_{dn}^{dk}={\omega }_n^k$，可以理解为把圆分成了$dn$份取$dk$份和圆分成$n$份取$k$份取得的面积不变，证明：${\omega }_{dn}^{dk}=(e^{\frac{2\pi i}{dn}}{)}^{dk}=(e^{\frac{2\pi i}{n}}{)}^k={\omega }_n^k$
2. (负方向) ${\omega }_n^{n/2}=-1$​，理解为$\frac{n}{2}$​恰好旋转到$x$​轴的负方向。
3. (周期性) ${\omega }_n^{2k+n}={\omega }_n^{2k}$，理解为我们是在单位圆上旋转$+n$相当于多旋转一圈。
4. (折半引理) ${\omega }_n^{k+n/2}=-{\omega }_n^k$，证明：${\omega }_n^{k+n/2}={\omega }_n^{n/2}\ast {\omega }_n^k=-1\ast {\omega }_n^k=-{\omega }_n^k$
   • 推论：$({\omega }_n^{k+n/2}{)}^2={\omega }_n^{2k}={\omega }_{n/2}^k$，证明：$({\omega }_n^{k+n/2}{)}^2={\omega }_n^{2k+n}={\omega }_n^{2k}={\omega }_{n/2}^k$
5. (求和引理) ${\sum }_{j=0}^{n-1}{\omega }_n^{kj}={\sum }_{j=0}^{n-1}({\omega }_n^k{)}^j=0$，证明如下：${\sum }_{j=0}^{n-1}({\omega }_n^k{)}^j=\frac{({\omega }_n^k{)}^n-1}{{\omega }_n^k-1}=\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1}=\frac{1^k-1}{{\omega }_n^k-1}=0$

DFT（系数转点值）

如果我们把${\omega }_n^0,{\omega }_n^1,{\omega }_n^2,\cdots ,{\omega }_n^{n-1}$作为$x_0,x1,x_2,\cdots ,x_{n-1}$代入$P(x)$,定义$y_k=P({\omega }_n^k)={\sum }_{j=0}^{n-1}{a}_j\ast {\omega }_n^{kj}$，

那么向量$y=(y_0,y_1,y_2,\cdots ,y_{n-1})$,就是系数向量$a=(a_0,a_1,a_2,\cdots ,a_{n-1})$的离散傅里叶变换（DFT）。

记$y=DF{T}_n(a)$，表示一个$n$次数界多项式的离散傅里叶变换。

FFT

如果我们把多项式$P(x)$拆分成两个多项式，$P_{odd}(x)$和$P_{even}(x)$，分别表示奇数系数项和偶数系数项

• $P_{odd}(x)=a_1+a_3{x}^2+a_5{x}^4+\cdots +a_{n-1}{x}^{n-2}\stackrel{u=x^2}{=}{a}_1+a_{3}u+a_5{u}^2+\cdots +a_{n-1}{u}^{n/2-1}$

• $P_{even}(x)=a_0+a_2{x}^2+a_4{x}^4+\cdots +a_{n-2}{x}^{n-2}\stackrel{u=x^2}{=}{a}_0+a_{2}u+a_4{u}^2+\cdots +a_{n-2}{u}^{n/2-1}$

那么可以发现$P(x)=P_{even}(u)+x{P}_{odd}(u)=P_{even}(x^2)+x{P}_{odd}(x^2)$

总结：$P(x)=P_{even}(x^2)+x{P}_{odd}(x^2)$，这表示我们只要求两个次数界为$\frac{n}{2}-1$的多项式在$({\omega }_n^0{)}^2,({\omega }_n^1{)}^2,({\omega }_n^2{)}^2,\cdots ,({\omega }_n^{n-1}{)}^2$处的值，即可得到原多项式的值。由折半引理和它的推论可以得到：
$$\begin{gathered} P({\omega }_n^k)=P_{even}({\omega }_{n/2}^k)+{\omega }_n^k{P}_{odd}({\omega }_{n/2}^k) \\ P({\omega }_n^{k+n/2})=P_{even}({\omega }_{n/2}^k)-{\omega }_n^k{P}_{odd}({\omega }_{n/2}^k) \end{gathered}$$
给出FFT的伪代码

def FFT(A):
    # 假设n是2的幂，如果不是向后补位到2的幂
    n = A.length
    if n == 1: return
    # divide
    Peven = (A0,A2,..,An-2)
    Podd  = (A1,A3,..,An-1)
    Yeven = FFT(Peven)
    Yodd  = FFT(Podd)
    # conquer
    Wn = e^(2pi)/n = cos(2pi/n) + sin(2pi/n)
    w = 1
    Y = list()
    # 折半引理和推论
    for k in range(n/2-1):
    	Y[k]     = Yeven[k] + w * Yodd[k]
        Y[k+n/2] = Yeven[k] - w * Yood[k]
        w = w * Wn # 旋转一个单位，相当于辐角+(1/n圆)
        
    return Y

• 我们通过$FFT(C)，FFT(Q)$将两个多项式$C(x),Q(x)$的系数表达式转变为了点值表达式，也就是求得两多项式的$DFT$

IDFT插值(点值转系数)

先抛出结论：证明
$$IDFT(y)=\frac{1}{n}\sum _{k=0}^{n-1}{y}_k\ast {\omega }_n^{-kj}$$

• 此时我们发现上式其实也是在做$DFT$，只不过里边的${\omega }_n^{kj}变为{\omega }_n^{-kj}$，而前边也多了个$\frac{1}{n}$。

  • ${\omega }_n^{-kj}=e^{-\frac{2k\pi i}{n}}=\cos (-\frac{2k\pi }{n})+i\sin (-\frac{2k\pi }{n})\stackrel{奇偶性}{=}\cos (\frac{2k\pi }{n})-i\sin (\frac{2k\pi }{n})$
  • 这说明${\omega }_n^{-kj}$相当于将${\omega }_n^{kj}$的虚部取反，实部不变。

• 所以我们考虑改造一下我们的$FFT$算法

# inv=-1表示IDFT，inv=1表示DFT
 def FFT(A, inv=1):
     n = A.length
     if n == 1: return
     
     Peven = (A0,A2,..,An-2)
     Podd  = (A1,A3,..,An-1)
     Yeven = FFT(Peven, inv)
     Yodd  = FFT(Podd,  inv)
    
 	 # 注意这里变化了，如果是逆变换虚部要取反
     Wn = e^(2pi)/n = cos(2pi/n) + inv * sin(2pi/n)
     w = 1
     Y = list()
   
     for k in range(n/2-1):
     	 Y[k]     = Yeven[k] + w * Yodd[k]
         Y[k+n/2] = Yeven[k] - w * Yood[k]
         w = w * Wn
         
     return Y

如何调用：

def solve():
    C = (a0, a1, ... , a2)
    Q = (b0, b1, ... , b2)
    # 如果他两的长度不是2的幂，补齐到2的幂
    n = power2(max(C.length,Q.length)) 
    # DFT转点值
    Yc = FFT(C, 1)
    Yq = FFT(Q, 1)
    Y = list()
    # 点值乘法 P(i) = C(i) * Q(i)
    for i in range(n):
        Y[i] = Yc[i] * Yq[i]
    # a = IDFT(Y)
    a = FFT(Y, -1)
    return a

思路总结：

要求$P(x)=C(x)\ast Q(x)$

1. 如果我们能通过选取一系列特殊的点使得求出点值表达式的时间加快,具体来说时间复杂度为$\theta (nlgn)$​
2. 对点值进乘法运算，使用$\theta (n)$的时间
   • 可以得到$P(x_i)=C(x_i)\ast Q(x_i)$
   • 可以的到P的点值表示形式： P ( x ) = { ( x 0 , C ( x 0 ) ∗ Q ( x 0 ) ) , ⋯   , ( x n − 1 , C ( x n − 1 ) ∗ Q ( x n − 1 ) ) } P(x)=\{(x_0, C(x_0)*Q(x_0)),\cdots,(x_{n-1}, C(x_{n-1})*Q(x_{n-1}))\} P(x)={(x0​,C(x0​)∗Q(x0​)),⋯,(xn−1​,C(xn−1​)∗Q(xn−1​))}
3. 将得到的P的点值表达式通过$IDFT$转换为系数表达式，时间复杂度$\theta (nlgn)$
4. 总之就是如下图

一些需要的证明

• 证明(多项式插值的唯一性)

  • 把第一个矩阵每一行看成一个多项式每一项代入的未知量，第二个矩阵看成多项式的系数，
    那么由矩阵乘法可以推出所有的多项式的值(用第一个矩阵的每一行乘第二个矩阵的每一列)。

  $$\left[\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^{n-1}\\ 1& x_1& x_1^2& \cdots & x_1^{n-1}\\ 1& x_2& x_2^2& \cdots & x_2^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{n-1}& x_{n-1}^2& \cdots & x_{n-1}^{n-1}\end{array}\right]\ast \left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ \cdots \\ a_{n-2}\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ \cdots \\ y_{n-2}\\ y_{n-1}\end{array}\right]$$

  • 第一个矩阵叫做范德蒙(Vandermonde)矩阵，由线性代数的知识我们知道他是可逆的。
  • 那么给等号右边的值矩阵乘以范德蒙矩阵的逆矩阵就可以推出系数矩阵，这代表两个矩阵可以互相转化。

• 证明(IDFT的推导)

  • 写出将单位根代入多项式后的矩阵相乘的形式

  $$\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }_n& {\omega }_n^2& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& \cdots & {\omega }_n^{2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& \cdots & {\omega }_n^{(n-1)(n-1)}\end{array}\right]\ast \left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ \cdots \\ a_{n-2}\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ \cdots \\ y_{n-2}\\ y_{n-1}\end{array}\right]$$

  • 记左1矩阵为：$V_n$（表示n阶范德蒙矩阵)，左2矩阵为：$a$ （表示多项式的系数矩阵），右1矩阵为：$y$ （表示多项式的值矩阵）

  • 那么我们就是要求$a$，也就是给等号两边同时乘以$V_n$的逆矩阵$V_n^{-1}$，那么这个逆矩阵该如何表示呢？

  • 定理：$V_n^{-1}$对于$j,k=0,1,\cdots ,n-1$的行标为$j$，列标为$k$处的元素为${\omega }_n^{-kj}/n$。

    要证上边的定理成立，就是要证明$V_n^{-1}{V}_n=I_n$，其中$I_n$是$n$阶单位阵。

    考虑$(i,j)$处的元素：
    $$[V_n^{-1}{V}_n{]}_{ij}=\sum _{k=0}^{n-1}(\frac{{\omega }_n^{-ki}}{n}){\omega }_n^{kj})=\sum _{k=0}^{n-1}{\omega }_n^{k(j-i)}$$
    如果$i=j$则上式的和为1(主对角线上元素为1)，否则根据求和引理该和为0，那么乘积的结果为：
    $$\left[\begin{array}{ccccc}1& 0& 0& \cdots & 0\\ 0& 1& 0& \cdots & 0\\ 0& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 1\end{array}\right]$$
    这恰好是单位矩阵$I_n$,得证。

  • 那么就有下式成立：
    $$a=V_n^{-1}y=\sum _{k=0}^{n-1}\frac{{\omega }_n^{-kj}}{n}\ast y_k=\frac{1}{n}\sum _{k=0}^{n-1}{y}_k\ast {\omega }_n^{-kj}$$

• (简单理解欧拉公式)
  • $e^{ix}=\cos (x)+i\sin (x)$,其实就是把一个复数在复平面上对着$\mathbb{R}$轴和$I$轴做投影，$\mathbb{R}$轴投影叫做$\cos (x)$,$I$轴投影叫做$i\sin (x)$，然后数学家恰好发现这些东西等于$e^{ix}$（可以用泰勒展开验证，详见下边链接）。
  • 深入理解请参考：欧拉公式