快速傅立叶变换FFT学习笔记

什么是FFT？

FFT（Fast Fourier Transformation） 是离散傅氏变换（DFT）的快速算法，即快速傅氏变换。FFT使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。FFT可以将多项式乘法的复杂度从$O(n^2)$降到$O(nlogn)$。

下图是FFT的整体计算流程，FFT变换的复杂度为$O(nlogn)$，FFT域上的pointwise乘法的复杂度为$O(n)$，逆FFT变换的复杂度为$O(nlogn)$，总体复杂度为$O(nlogn)$。

多项式的表示方法

方法1：系数表示

从多项式函数的定义，将所有系数视为系数向量，而由全部系数组成的向量$a$叫做该多项式的系数表达:

$$f(x)=\sum _{i=0}^{n?1}{a}_i{x}^i$$

$$a=(a_0,a_2,\dots ,a_{n?1})$$

举个简单的例子：$f(x)=5x^0+6x^1+7x^2$的系数表示为 { 5 , 6 , 7 } \{5, 6, 7\} {5,6,7}。反之， { 5 , 6 , 7 } \{5, 6, 7\} {5,6,7}的系数编码结果为$f(x)=5x^0+6x^1+7x^2$。

系数表示特点是对多项式加法友好，时间复杂度是$O(n)$。但是对多项式乘法不友好，采用多项式逐项相乘，时间复杂度仍为$O(n^2)$。
注：系数表示的乘法表示卷积操作。

方法2：点值表示

任意选取$n$个不同的自变量$x$带入多项式函数$f(x)$进行求值运算，将得到$n$个不同的结果。于是，多项式的点值表达就是由这$n$个数值点组成的集合：

$$\{(x_0,f(x_0)),(x_1,f(x_1)),\dots ,(x_{n?1},f(x_{n?1}))\}$$

举个简单的例子：$f(x)=5x^0+6x^1+7x^2$的点值表示为$\{(0,f(0)),(1,f(1)),(2,f(2))\}$。反之， { 5 , 6 , 7 } \{5, 6, 7\} {5,6,7}的点值编码应该满足$f(0)=5,f(1)=6,f(2)=7$

点值表示的特点是对多项式乘法友好，时间复杂度是$O(n)$【因为可以做element-wise乘法（EWMM）】，但是多项式加法不友好。

复数

复数的定义：设$a,b$为实数，则$z=a+bi$的数称为复数，其中$a$称为实部，$b$称为虚部。
复数的模为：$|z|=\sqrt{a^2+b^2}$。一个复数的共轭复数为：$\overline{z}=a?bi$，即改变虚部的符号。

单位复数根
对于任意一个复数$\omega$，其$n$次幂的结果为1，就称复数$\omega$是$n$次单位复数根，即
$${\omega }^n=1$$

可以看到，$n$次单位复数根有$n$个，其几何意义为：$n$个单位复数根均匀地分布在以复平面原点为圆心的单位圆上。

在几何意义的单位圆中，我们将圆周角$2\pi$均分成$n$份，则$\frac{2\pi }{n}$叫做单位根的幅角。
由欧拉公式：
$$e^{i\frac{2\pi }{n}}=cos\frac{2\pi }{n}+isin\frac{2\pi }{n}$$

定义${\omega }_n$表示一个$n$次单位根：
$${\omega }_n=e^{i\frac{2\pi }{n}}$$

${\omega }_n$也是主$n$次单位根，其余$w_n^1、w_n^2$等叫做$n$次单位根的幂次，记为：
$${\omega }_n^k=e^{i\frac{2k\pi }{n}},k=0,1,\dots ,n?1$$

于是，很容易知道：
$${\omega }_n^0={\omega }_n^n=1,{\omega }_n^{\frac{n}{2}}=?1$$

单位复数根的性质1：消去引理
$${\omega }_{dn}^{dk}=e^{i\frac{2dk\pi }{dn}}=e^{i\frac{2k\pi }{n}}=w_n^k$$

单位复数根的性质1：折半引理
$${\omega }_n^{k+\frac{n}{2}}={\omega }_n^k{\omega }_n^{\frac{n}{2}}=?{\omega }_n^k$$
于是也可以得到：
$$({\omega }_n^{k+\frac{n}{2}}{)}^2=(?{\omega }_n^k{)}^2=({\omega }_n^k{)}^2={\omega }_n^{2k}={\omega }_{\frac{n}{2}}^k$$
好处是将$n$降到了原来的一半。

DFT：离散傅立叶变换

假设多项式：
$$A(x)=\sum _{i=0}^{n?1}{a}_i{x}^i$$

把$n$次单位根的幂次$x={\omega }_n^k$分别代入多项式：
$$y_k=A({\omega }_n^k)=\sum _{i=0}^{n?1}{a}_i{\omega }_n^{ki},k=0,1,\dots ,n?1$$

记$y=(y_0,y_1,\dots ,y_{n?1})$是系数向量$a=(a_0,a_1,\dots ,a_{n?1})$的离散傅立叶变换，即DFT。

IDFT：离散傅立叶逆变换
$$a_j=\frac{1}{n}\sum _{k=0}^{n?1}{y}_k{\omega }_n^{?ki}$$

DFT对应多项式求值
IDFT对应插值，求多项式系数

值得注意的是，DFT的复杂度仍然是$O(n^2)$。

FFT和蝶形计算

FFT的原理是将多项式分解成奇偶两部分，并用分治的思想依次计算下去。
$$A(x)=a_0+a_1{x}^1+?+a_{n?1}{x}^{n?1}$$

$$A_0(x)=a_0+a_2{x}^1+?+a_{n?2}{x}^{\frac{n?2}{2}}$$

$$A_1(x)=a_1+a_3{x}^1+?+a_{n?1}{x}^{\frac{n?2}{2}}$$

$$A(x)=A_0(x^2)+x{A}_1(x^2)$$

证明：
$$\begin{gathered} A(x)=A_0(x^2)+x{A}_1(x^2)=a_0+a_2{x}^2+a_4{x}^4+?+a_{n?2}{x}^{n?2}+ \\ {a}_1{x}^1+a_3{x}^3+a_5{x}^5+?+a_{n?1}{x}^{n?1} \end{gathered}$$
得证！

将$x={\omega }_n^k$代入$A(x)=A_0(x^2)+x{A}_1(x^2)$中，得到：
$$A({\omega }_n^k)=A_0({\omega }_n^{2k})+{\omega }_n^k{A}_1({\omega }_n^{2k})=A_0({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{A}_1({\omega }_{\frac{n}{2}}^k)$$

将$x={\omega }_n^{k+\frac{n}{2}}$代入$A(x)=A_0(x^2)+x{A}_1(x^2)$中，得到：
$$A({\omega }_n^{k+\frac{n}{2}})=A_0({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}{A}_1({\omega }_n^{2k+n})=A_0({\omega }_n^{2k})?{\omega }_n^k{A}_1(w_n^{2k})=A_0({\omega }_{\frac{n}{2}}^k)?{\omega }_n^k{A}_1(w_{\frac{n}{2}}^k)$$

我们发现，$A({\omega }_n^k)$和$A({\omega }_n^{k+\frac{n}{2}})$的第一项完全相同，仅第二项为相反数。因此，如果知道$A_0({\omega }_{\frac{n}{2}}^k)$和$A_1({\omega }_{\frac{n}{2}}^k)$的值，我们就可以同时知道$A({\omega }_n^k)$和$A({\omega }_n^{k+\frac{n}{2}})$，所以可以用分治思想计算FFT，原问题的规模缩减了一半。

总结一下，FFT的计算如下：
$$A({\omega }_n^k)=A_0({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{A}_1({\omega }_{\frac{n}{2}}^k)$$

$$A({\omega }_n^{k+\frac{n}{2}})=A_0({\omega }_{\frac{n}{2}}^k)?{\omega }_n^k{A}_1(w_{\frac{n}{2}}^k)$$

可以通过这样的方式将一个多项式一直分解下去，如下图是对16点输入的分解：

在计算FFT时，需要成对的点做蝶形运算，这里成对的点就是0和8、4和12等，这个分组的过程可以用bit reverse实现。

8点FFT计算图示：

每一对数的蝶形运算：

Bit Reverse确定蝶形运算对

从下面这个8点FFT可以很清楚地看到，FFT蝶形运算时打乱了输入的顺序（倒位序），倒位序是由bit reverse操作得到的。

FFT的输入为倒位序，输出为自然顺序。

Bit reverse的原理其实并不复杂，从上文中16点输入的奇偶分解那个图就很容易看出来。

RFFT和FFT

RFFT中的R是实数的意思，RFFT是FFT的特殊版本，为实数输入设计。RFFT利用了实数的傅立叶变换为共轭对称这个事实，因此RFFT只需要计算一半的傅立叶变换系数。所以RFFT效率明显高于FFT，并且也只有一半的存储开销。

因此，当我们的输入为实数时（比如图像卷积任务），我们就可以利用实数的傅立叶变换为共轭对称这个特性，用RFFT替换FFT来提高计算效率。

实数的FFT为什么是共轭对称？

我们直接看DFT的计算公式（把上文中的索引i改成了t，方便和复数i区分开）：
$$A({\omega }_n^k)=\sum _{t=0}^{n?1}{a}_t{\omega }_n^{kt}$$

其中,
$${\omega }_n^k=e^{i\frac{2k\pi }{n}}$$

于是，代入得到：
$$A({\omega }_n^k)=\sum _{t=0}^{n?1}{a}_t{e}^{i\frac{2kt\pi }{n}}(1)$$

同时，我们可以计算出与上面点对称的点：
$${\omega }_n^{n?k}=e^{i\frac{2(n?k)\pi }{n}}$$

同样，代入得到：
$$A({\omega }_n^{n?k})=\sum _{t=0}^{n?1}{a}_t{\omega }_n^{(n?k)t}$$

其中，
$${\omega }_n^{(n?k)t}={\omega }_n^{nt?kt}={\omega }_n^{nt}/{\omega }_n^{kt}={\omega }_n^{?kt}$$

于是，
$$A({\omega }_n^{n?k})=\sum _{t=0}^{n?1}{a}_t{\omega }_n^{?kt}=\sum _{t=0}^{n?1}{a}_t{e}^{?i\frac{2kt\pi }{n}}(2)$$

容易发现，式（1）和（2）只是在$e$的指数上为相反数关系！
根据欧拉公式，对于（1）：
$$e^{i\frac{2kt\pi }{n}}=cos\frac{2kt\pi }{n}+isin\frac{2kt\pi }{n}$$

对于（2）：
$$e^{?i\frac{2kt\pi }{n}}=cos\frac{2kt\pi }{n}?isin\frac{2kt\pi }{n}$$

证毕！

一大堆参考文献

• 零知识证明 - 理解FFT的蝶形运算
• 十分简明易懂的FFT（快速傅里叶变换）
• 大数乘法—多项式与快速傅里叶变换
• 快速傅立叶变换（Fast Fourier Transform）
• fft海面模拟（二）
• 彻底搞懂快速傅里叶变换FFT–旋转因子
• 快速傅里叶变换（FFT）之一：Radix-2 DIT FFT
• Rfft 和 FFT 有什么区别？