其他算法-FFT

FFT算法

快速傅里叶变换 (Fast Fourier Transform)， 即利用计算机计算离散傅里叶变换DFT的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。
学习FFT前需要先掌握复数，单位根和多项式

复数

复数z=a+bi，a和b都是实数，$i=\sqrt{-1}$，z可以写成$z=r{e}^{i\theta }$，r是z的模，$\theta$为z的幅角，其中$e^{i\theta }=cos\theta +isin\theta$，任意的a+bi对应复平面上的(a,b)，a被称为实部，b被称为虚部

单位根

对于$z=r{e}^{i\theta }$，有：

单位根的定义是$z^n=1$在复数域内的n个根，即$r^n{e}^{i\theta n}=1$，容易看出，$r=1,\theta =\frac{2k\pi }{n}$，所以$z^n=1$的n个根为：
$$w_n^k=e^{i\frac{2k\pi }{n}}=cos\frac{2k\pi }{n}+isin\frac{2k\pi }{n}$$
其中，$k=0,1,2,...,n-1$，$w_n^1=e^{i\frac{2\pi }{n}}$称为主单位根

单位根的3个引理

• 消去引理：$w_{dn}^{dk}=w_n^k$，比如$w_n^{\frac{n}{2}}=w_2^1=-1$
  证明：
  $$w_{dn}^{dk}=(e^{i\frac{2\pi }{dn}}{)}^{dk}=e^{i\frac{2\pi k}{n}}=w_n^k$$
• 折半引理----FFT的核心
  $$(w_n^{k+\frac{n}{2}}{)}^2=(w_n^k{)}^2=w_{\frac{n}{2}}^k$$
  证明：
  $$w_n^{k+\frac{n}{2}}=w_n^k{w}_n^{\frac{n}{2}}=-w_n^k$$
  $$(-w_n^k{)}^2=w_n^{2k}=w_{\frac{n}{2}}^k$$
• 求和引理
  $$\sum _{i=0}^{n-1}(w_n^k{)}^i=0$$
  证明：
  根据等比数列求和公式有：
  $$\sum _{i=0}^{n-1}(w_n^k{)}^i=\frac{(w_n^k{)}^n-1}{w_n^k-1}$$
  $$=\frac{(w_n^n{)}^k-1}{w_n^k-1}=\frac{(1{)}^k-1}{w_n^k-1}=0$$

多项式

在实际应用中，FFT也常用于处理多项式乘法；假设有一个多项式$F(x)=a_0+a_{1}x+...+a_n{x}^n={\sum }_{i=0}^n{a}_i{x}^i$，则次数为$degree(F)=n$，次数界为大于次数的整数，因此次数界$>n$；

多项式加法

现在有两个多项式$A(x)={\sum }_{i=0}^n{a}_i{x}^i$和$B(x)={\sum }_{i=0}^n{b}_i{x}^i$；
则$A(x)+B(x)$为：
$$C(x)=\sum _{i=0}^n{c}_i{x}^i,c_i=a_i+b_i$$

多项式乘法

先假设有多项式M：$6x^3+7x^2-10x+9$；
多项式N：$-2x^3+4x-5$；
多项式M乘多项式N，需要用多项式N的每一项乘多项式M，得到3个多项式，结果则是这3个多项式相加，经过观察发现结果的次数为$degree(M)+degree(N)$；
如果推广到同次多项式$A(x)$和$B(x)$，则有：
$$C(x)=A(x)B(x)=\sum _{i=0}^{2n}{c}_i{x}^i$$
而$c_i={\sum }_{j=0}^i{a}_j{b}_{i-j}$，如果记$c$为$c_i$构成的序列，$a$为$a_i$构成的序列，$b$为$b_i$构成的序列，则$c$为$a$与$b$进行卷积：
$$c=a\otimes b$$
对于degree不同的两个多项式，可通过添加0系数上升至同一次数n

多项式的表示

系数表示
用多项式的系数表示多项式：
$$a=[a_0,a_1,...,a_n{]}^T$$
可见，加法的时间复杂度为$O(n)$，乘法的时间复杂度为$O(n^2)$；
点值表示
用至少n个多项式上的点表示，即：
{ ( x 0 , A ( x 0 ) , ( x 1 , A ( x 1 ) ) , . . . , ( x n , A ( x n ) } \left \{ (x_{0},A(x_{0}),(x_{1},A(x_{1})),...,(x_{n},A(x_{n}) \right \} {(x0​,A(x0​),(x1​,A(x1​)),...,(xn​,A(xn​)}
$A(x_0)$只是多项式$A(x)$在$x_0$处的取值；
点值表示的多项式运算应确保$x_i$对应：
加法：$C(x_i)=A(x_i)+B(x_i)$，时间复杂度为$O(n)$；
乘法：$C(x_i)=A(x_i)B(x_i)$，时间复杂度为$O(n)$；

---

插值
已知点值，给出一个$x$求值$A(x)$叫插值，插值一般由拉格朗日插值公式计算：
$$A(x)=\sum _{i=0}^{n}A(x_i)\frac{{\prod }_{j\ne i}^{}(x-x_j)}{{\prod }_{j\ne i}^{}(x_i-x_j)}$$
时间复杂度为$O(n^2)$；
比如多项式有两个点值 { ( x 1 , y 1 ) , ( x 2 , y 2 ) } \left \{ (x_{1},y_{1}),(x_{2},y_{2}) \right \} {(x1​,y1​),(x2​,y2​)}，插值得到$(x,y)$为：
$$y=y_1\frac{x-x_2}{x_1-x_2}+y_2\frac{x-x_1}{x_2-x_1}$$

---

DFT离散傅里叶变换

DFT即Discrete Fourier Transform，DFT是对一个序列进行变换，假设有多项式$A(x)={\sum }_{i=0}^{n-1}{a}_i{x}^i$，现在把单位根$w_n^0,w_n^1,...,w_n^{n-1}$代入多项式得到：
$$y=[y_0,y_1,...,y_{n-1}]$$
其中：
$$y_i=A(w_n^i)=\sum _{j=0}^{n-1}{a}_j(w_n^i{)}^j$$
记$y=DF{T}_n(a)$；
可见，由于单位根是已知的，只要给出序列a，就能得到y，而计算一次$y_i$的时间复杂度为$O(n)$，$y_i$共n个，所以时间复杂度为$O(n^2)$

IDFT离散傅里叶逆变换

DFT中，$y_i={\sum }_{j=0}^{n-1}{a}_j(w_n^i{)}^j={\sum }_{j=0}^{n-1}{a}_j{w}_n^{ij}$，如果用矩阵表示可以写出：
$$y=V_{n}a$$
$$(V_n{)}_{i,j}=w_n^{ij}$$
因此有：
$$a=V_n^{-1}y$$
即：
$$a_i=\frac{1}{n}\sum _{j=0}^{n-1}{w}_n^{-ij}{y}_j$$

FFT快速傅里叶变换

FFT即Fast Fourier Transform，通过单位根的引理将时间复杂度从DFT的$O(n^2)$降为$O(nlogn)$

证明

在DFT中，有序列a：
$$a=[a_0,a_1,...,a_{n-1}]$$
将a拆开为：
$$a^{[0]}=[a_0,a_2,...,a_{n-2}]$$
$$a^{[1]}=[a_1,a_3,...,a_{n-1}]$$
偶数项$a^{[0]}$和奇数项$a^{[1]}$分别对应多项式$A^{[0]}(x),A^{[1]}(x)$；
现在写出：
$$A(x)=a_0+a_{1}x+...+a_{n-1}{x}^{n-1}$$
$$A^{[0]}(x)=a_0+a_{2}x+a_4{x}^2+...+a_{n-2}{x}^{\frac{n}{2}-1}$$
$$A^{[1]}(x)=a_1+a_{3}x+a_5{x}^2+...+a_{n-1}{x}^{\frac{n}{2}-1}$$
如何关联三者？
先写出：
$$A^{[0]}(x^2)=a_0+a_2{x}^2+a_4{x}^4+...+a_{n-2}{x}^{n-2}$$
$$x{A}^{[1]}(x^2)=(a_1+a_3{x}^2+a_5{x}^4+...+a_{n-1}{x}^{n-2})x$$
$$=a_1+a_3{x}^3+a_5{x}^5+...+a_{n-1}{x}^{n-1}$$
所以有：
$$A(x)=A^{[0]}(x^2)+x{A}^{[1]}(x^2)$$
此时，$y_k$为：
$$y_k=A(w_n^k)=A^{[0]}((w_n^k{)}^2)+w_n^k{A}^{[1]}((w_n^k{)}^2)$$
同理，对于：
$$A(w_n^{k+\frac{n}{2}})=A^{[0]}((w_n^{k+\frac{n}{2}}{)}^2)+w_n^{k+\frac{n}{2}}{A}^{[1]}((w_n^{k+\frac{n}{2}}{)}^2)$$
根据折半引理$(w_n^{k+\frac{n}{2}}{)}^2=w_{\frac{n}{2}}^k$有：
$$A(w_n^{k+\frac{n}{2}})=A^{[0]}(w_{\frac{n}{2}}^k)+w_n^{k+\frac{n}{2}}{A}^{[1]}(w_{\frac{n}{2}}^k)$$
在证明折半引理时，有一个中间结论：$w_n^{k+\frac{n}{2}}=-w_n^k$，因此可以化简为：
$$A(w_n^{k+\frac{n}{2}})=A(-w_n^k)=A^{[0]}(w_{\frac{n}{2}}^k)-w_n^k{A}^{[1]}(w_{\frac{n}{2}}^k)$$
反过来得到：
$$y_k=A(w_n^k)=A^{[0]}(w_{\frac{n}{2}}^k)+w_n^k{A}^{[1]}(w_{\frac{n}{2}}^k)$$
此时，分别在$a^{[0]}$和$a^{[1]}$上计算$\frac{n}{2}$次多项式即可；
FFT的高效之处
在前面已经知道：
$$y_k=A(w_n^k)=A^{[0]}(w_{\frac{n}{2}}^k)+w_n^k{A}^{[1]}(w_{\frac{n}{2}}^k)$$
另外：
$$A(w_n^{k+\frac{n}{2}})=A^{[0]}(w_{\frac{n}{2}}^k)-w_n^k{A}^{[1]}(w_{\frac{n}{2}}^k)$$
注意到：$y_{k+\frac{n}{2}}=A(w_n^{k+\frac{n}{2}})$，也就是说，$y_{k+\frac{n}{2}}$和$y_k$在计算上仅相差$w_n^k$的正负，所以可以保存$A^{[0]}(w_{\frac{n}{2}}^k)$和$A^{[1]}(w_{\frac{n}{2}}^k)$，即空间换时间，因此$w_n^k$的计算规模缩小为从0到$\frac{n}{2}$，即$DF{T}_n(a)$可以转移为计算$DF{T}_{\frac{n}{2}}(a^{[0]})$和$DF{T}_{\frac{n}{2}}(a^{[1]})$
时间复杂度
基于FFT，原来的$DF{T}_n(a)$可以转移为$DF{T}_{\frac{n}{2}}(a^{[0]})$和$DF{T}_{\frac{n}{2}}(a^{[1]})$：

而$DF{T}_{\frac{n}{2}}(a^{[0]})$或者$DF{T}_{\frac{n}{2}}(a^{[1]})$又可以分别进一步缩小规模(递归)：
比如$DF{T}_{\frac{n}{2}}(a^{[0]})$转移为$DF{T}_{\frac{n}{4}}(a^{[0][0]})$和$DF{T}_{\frac{n}{4}}(a^{[0][1]})$；
因此时间复杂度$T(n)$应表示为，总的时间等于二倍子问题时间加上$n$个结果的合并时间：
$$T(n)=2T(\frac{n}{2})+O(n)$$

---

详细解释：$n$个结果的合并时间
这里递归的目的是返回{$A^{[0]}(w_{\frac{n}{2}}^k),A^{[1]}(w_{\frac{n}{2}}^k)$}，$k=0,1,..,\frac{n}{2}-1$；
获取完整的DFT还需要在$A^{[1]}(w_{\frac{n}{2}}^k)$前乘$w_n^k$，在FFT的高效之处提过$w_n^k$的系数：
前$\frac{n}{2}$的$w_n^k$系数为+，后$\frac{n}{2}$的$w_n^k$系数为-，所以完整的结果$y$还需要做n次合并处理：
$$A^{[0]}(w_{\frac{n}{2}}^k)+w_n^k{A}^{[1]}(w_{\frac{n}{2}}^k)$$
和：
$$A^{[0]}(w_{\frac{n}{2}}^k)-w_n^k{A}^{[1]}(w_{\frac{n}{2}}^k)$$

---

时间复杂度为：
$$T(n)=2T(\frac{n}{2})+O(n)$$
重复代入T(n)化简：
$$T(n)=2(2T(\frac{n}{4})+O(\frac{n}{2}))+O(n)=4T(\frac{n}{4})+2O(n)$$
归纳为：
$$T(n)=2^{k}T(\frac{n}{2^k})+kO(n)$$
令$n=2^k$得到：
$$T(n)=nT(\frac{n}{n})+lo{g}_{2}nO(n)=n+nlogn$$
所以时间复杂度为$O(nlogn)$

FFT的应用

多项式部分提到，多项式相乘实际上是两个系数向量作卷积，信号系统里有时域的卷积等于频域相乘，根据DFT和IDFT的内容，可以得到：
$$a\otimes b=IDF{T}_{2n}(DF{T}_{2n}(a)\circ DF{T}_{2n}(b))$$
符号$\circ$代表两序列的元素对应相乘，可以理解为python内置高阶函数中的map；使用FFT去作DFT，将可以节省时间的开销