一、连续傅里叶变换 和 离散傅里叶变换
假设$f$是复平面上的勒贝格可积的函数，则定义$f$的连续傅里叶变换$F$是
对于任意实数$w$,
$$
F(\omega) =
\int_{-\infty} ^{\infty} f(t) e^ { -i \omega t }dt
\tag{1}
$$
在这里，我们解释$\omega$为角频率，$F(\omega)$是复数，并且是信号在该频率成分处的幅度和相位。
傅里叶逆变换是
$$
f(t) = \frac {1} {2 \pi}
\int_{-\infty} ^{\infty} F(\omega)e^{i \omega t}dw
\tag{2}
$$
怎么理解傅里叶变换这2个公式呢。
傅里叶提出，任何连续周期函数，都可以由不同频率不同相位不同幅度的正弦波组合而成。
上述的傅里叶变换公式表达的是，如何求给定的频率，它在信号$f$中的幅度和相位。
在计算机中，我们无法处理连续的积分行为，因此只能进行离散化的采样了。
对于原始信号$f(t)$的离散采样$f(n)$, 我们定义，它的离散傅里叶变换（DFT）为
$$
F[k] = \sum_{n=0} ^{N-1} {f[n]e^{-\frac {2 \pi } {N} i k n}}
\tag{3}
$$
离散傅里叶逆变换（IDFT）为
$$
f[n] = \frac {1}{N} \sum_{k=0} ^{N-1} F[k] e^{\frac{2 \pi} {N} i k n}
\tag{4}
$$
这里k表达的是离散的频率，n是离散的信号采样点，$0 \le k \lt N$, $0 \le n \lt N$。
二、离散傅里叶变换的快速算法
快速傅里叶变换利用的是分治法的思想。
2.1 奇偶分割
假设$N=2^L$, 我们将$f[n]$按n的奇偶分为2组：
令 $r \in [0, 1... N/2-1]$，
则
$$
\begin{aligned}
F[k] &= \sum_{n=0} ^{N-1} {f[n]e^{-\frac {2 \pi } {N} i k n}}
\ &= \sum_{r=0} ^{\frac{N}{2}-1} {f[2r]e^{-\frac {2 \pi } {N} i k {(2r)}}} +
\sum_{r=0} ^{\frac{N}{2}-1} {f[2r+1]e^{-\frac {2 \pi } {N} i k {(2r+1)}}}
\end{aligned}
$$
我们可以进行一些函数代换，假设$f_1[r]=f[2r], f_2[r]=f[2r+1]$, $F_1[k]$ 和 $F_2[k]$分别是$f_1[n], f_2[n]$的傅里叶变换, 其中$r, k \in [0,1,2... \frac N 2-1]$。
则上式等于
$$
\begin{aligned}
F[k] &= \sum_{r=0} ^{\frac{N}{2}-1} {f[2r]e^{-\frac {2 \pi } {N/2} i k {r}}} +
\sum_{r=0} ^{\frac{N}{2}-1} {f[2r+1]e^{-\frac {2 \pi } {N/2} i k {r}}} * e^{-\frac {2 \pi } {N} i k}
\ &= \sum_{r=0} ^{\frac{N}{2}-1} {{f_1}[r]e^{-\frac {2 \pi } {N/2} i k {r}}} +
e^{-\frac {2 \pi } {N} i k} \sum_{r=0} ^{\frac{N}{2}-1} {f_2[r]e^{-\frac {2 \pi } {N/2} i k {r}}}
\ &= F_1[k] + e^{-\frac {2 \pi } {N} i k} F_2[k]
\end{aligned}
$$
上述式子中，我们成功的将$F[k]$转换成了2个子傅里叶变换之和，我们可以了解到上述式中k的范围在$[0,1,2... \frac N 2-1]$中，因为两个子傅里叶变换的k的区间只有这么大。
我们求出了离散傅里叶变换分解式
$$
F[k] = F_1[k] + e^{-\frac {2 \pi } {N} i k} F_2[k]， k \in [0,1,2...\frac N 2 - 1] \tag{5}
$$
，但是k的区间还不够，还差一半。
2.2 求下一半的公式
下面将描述如何求取$F[k], k > N/2 -1$。
由
$$
e^{-\frac {2 \pi } {N/2} i {(k + N/2)} {r}} = e^{-\frac {2 \pi} {N/2} ikr} * e^{-2 \pi ir}
$$
由欧拉方程$e^{i\theta} = cos\theta + i sin\theta$，另外$r$属于正整数，可知
$$
e^{-2 \pi ir} = cos(-2\pi r) + i sin(-2 \pi r) = 1
$$
因此
$$
e^{-\frac {2 \pi } {N/2} i {(k + N/2)} {r}} = e^{-\frac {2 \pi} {N/2} ikr}
$$
因此
$$
\begin{aligned}
F_1[\frac N 2 + k]
&= \sum_{r=0} ^{\frac{N}{2}-1} {{f_1}[r]e^{-\frac {2 \pi } {N/2} i (k + N/2) {r}}}
\ &= \sum_{r=0} ^{\frac{N}{2}-1} {{f_1}[r]e^{-\frac {2 \pi } {N/2} i k {r}}} \ &= F_1[k]
\end{aligned}
$$
同理，$F_2[\frac N 2 + k] = F_2[k]$,
因此
$$
\begin{aligned}
F[\frac N 2 + k]
&= F_1[\frac N 2 + k] + e^{-\frac {2 \pi } {N} i ({\frac N 2 + k})} F_2[\frac N 2 + k]
\ &= F_1[\frac N 2 + k] + e^{-\pi i } * e^{-\frac {2 \pi } {N} i k} F_2[\frac N 2 + k]
\ &= F_1[\frac N 2 + k] - e^{-\frac {2 \pi } {N} i k} F_2[\frac N 2 + k]
\ &= F_1[k] - e^{-\frac {2 \pi } {N} i k} F_2[k]
\ &= F_1[k] + e^{-\frac {2 \pi } {N} i {(k + \frac N 2)}} F_2[k]
\end{aligned}
$$
离散傅里叶变换的另一半的分解式
$$
\begin{aligned}
F[\frac N 2 + k] &= F_1[k] - e^{-\frac {2 \pi } {N} i k} F_2[k]
\ &= F_1[k] + e^{-\frac {2 \pi } {N} i {(k + \frac N 2)}} F_2[k]
\ &k \in [0,1,2...\frac N 2 - 1] \tag{6}
\end{aligned}
$$
注意到，这里公式的2种方式都是可以滴，有的代码实现就是用第二个子式实现的。
2.3 蝶形运算
为了简化下面的一些符号的表示，我们令
$$
W^k_N = e^{-\frac {2 \pi } {N} i k}
$$
我们定义，蝶形运算公式
$$
\begin{aligned}
F[k] &= F_1[k] + W^k_N F_2[k]
\
F[k + \frac N 2] &= F_1[k] - W^k_N F_2[k] = F_1[k] + W^{k+\frac N 2}_N F_2[k]
\end{aligned}
\tag{7}
$$
它的图示为
上图中，$F_1$和$F_2$分别是$F$的奇数部分和偶数部分的子傅里叶变换的结果，他们的大小为$\frac N 2$。
2.4 两个DFT4利用蝶形变换生成DFT8
如下图所示，我们利用2个黑盒生成的DFT4生成DFT8.
首先，我们先取DFT8的原数据的奇数和偶数的信号，然后用某种黑盒过程生成DFT4，然后再利用蝶形变换生成DFT8.
2.5 利用4个DFT2生成2个DFT4，再生成DFT8
上述说的黑盒过程，可以进一步利用蝶形变换分解。下图的绿色部分标记出了新增的蝶形变换部分。 注意到，上面的序号发生了变化，每个4个DFT的奇数取出来。
2.6 将DFT2直接转换为蝶形变换
通过上述三个递归步骤，可以了解到FFT的详细做法了。
利用分治思想，一直2分，直到只有N=2时，即可直接变为蝶形变换。
2.7 序列输入顺序
我们可以注意到，最原始的DFT8在完全分解成蝶形变换时，$f[t],\space \space t \in [0, 8)$，序列发生了变化，因为除了N=2之外，每一次分解为蝶形变换时，因为取奇偶导致顺序发生了错位。
这种看起来没啥规律，但是很明显，这种规律交换的蝶形运算，最后这个序号肯定有规律的。规律刚好是二进制的位反 (bitreverse)。
如下图所示
三、FFT具体实现
我们注意到，上述过程
FFT的具体蝶形如何做，有很多种方案，下面这个链接的PDF里总结了好多种。
下面我着重讲述2个最常见的类型，第一个是最原版的先换顺序然后进行蝶形变换。第二个是GPU中最常见的FFT实现。
首先，为了简化细节，我直接将第二章描述的蝶形变换内部细节标记为一个黑盒，如下图。
3.1 经典快速傅里叶变换算法(Cooley-Tukey FFT)
经典FFT算法，就是先进行一次序号的位反计算出新的序号，然后按照鲽形图进行递归或者非递归计算。
下图是Cooley-Tukey FFT的DFT8的详细流程图
我们可以注意到，在Cooley-Tukey FFT的实现中，除了一开始需要映射表的序号之外，所有的蝶形变换中，输入的两个元素的序号和输出的2个元素的序号没有发生变化。
例如DFT2->DFT4中，4和6元素的下标分别是1, 3, 他们进行蝶形变换后输出的下标也是1, 3。
3.2 适合GPU的FFT(Stockham FFT)
上述的FFT，输入数据要进行一次序号重映射行为，在GPU端，这种映射可通过查表实现。但是，有一种Stockham FFT，它提出了一种巧妙的映射方式，使得id能直接计算映射，这在GPU中相比查表会更加快捷，对于大型FFT来说，构建表传输表都比较慢。
下图来自其他地方，这两个图简单明了的描述了区别。
下图是Stockham FFT的详细流程图。
下面，我们解析一下Stockham FFT的流程图
我们从第0轮开始分析，
首先，数据被分成8组，每个数据单独成组。然后，(0,4)(1,5)(2,6)(3,7)分别是DFT2的奇偶子项，他们分别做蝶形变换，完成了第0轮。
第一轮中，数据被分成了4组，$(0,4)(1,5)(2,6)(3,7)$, 其中$[(0,4)(2,6)]$与$[(1,5)(3,7)]$分别是DFT4的奇偶子项，然后他们分别做蝶形变换，完成第一轮。
以此类推。。。
下面，我们解析以下如何计算Stockham FFT的映射。
3.2.1 Stockham蝶形变换的输入下标的求取
我们观察可以发现以下特征：
- 所有蝶形变换的2个输入对的下标都符合$x, x + \frac N 2$的形式，$N$是数据的长度。例如上面的图里N=8.
- 所有蝶形变换的第一个输入的下标组合起来刚好是$[0, 1, 2, ... \frac N 2)$的形式。
根据这2点，我们可以在代码里确定蝶形变换的输入的下标了。
假设某个蝶形变换左边输入的下标为$x$，则第二个输入的下标是$x + \frac N 2$, 其中$x \in [0,1...\frac N 2]$。
注意到每一轮FFT合并过程中，都执行了$N/2$次蝶形变换，因此我们可以定义：
假设某个蝶形变换左边输入的下标为$x$，则称这个蝶形变换为第$x$个蝶形变换。
然后我们来求，蝶形变换的输出的下标。
3.2.2 Stockham蝶形变换的输出下标的求取
可以观察到：
- 每个蝶形变换的2个输出的下标间隔会随着每一轮蝶形变换而呈指数扩张。例如 DFT1 -> DFT2中，2个输出的下标间隔是1. DFT2->DFT4中，2个输出下标的间隔是2。 DFT4->DFT8中，2个输出下标间隔是4.
- 每一轮内，每组DFT内的蝶形变换，都是左边输出排在右边输出前面。然后，每一轮内不同组之间延后摆放。
因此，我们可以求出以下结论：
第m轮(m从0开始)，该轮内的子DFT的大小是$p = 2^m$, 输出的DFT的大小是$2^{m+1}$。
第m轮(m从0开始)，第$x$个蝶形变换，它在某个子DFT内的下标是$k = x mod p$。
第m轮(m从0开始)，第$x$个蝶形变换的左边输出的下标为$$
floor(\frac x {2^m}) * 2^{m+1} + x mod(2^m) =(x-k)*2 + k$$
右边输出的下标为
$$
floor(\frac x {2^m}) * 2^{m+1} + x mod(2^m) + 2^m =(x-k) * 2 + k + p$$
3.2.3 Stockham根据输出下标求取输入的下标
在我浏览的一些实现中，GPU是通过输出的下标进行分发任务的，也就是说，每个GPU函数中输出的下标确定了，这时候我们需要反向求出蝶形变换的2个输入下标，以及该输出在蝶形变换中的正负。
我们通过计算和分析可得出
假设第m轮(m从0开始)，第x个输出下标所属的
- 蝶形变换的序号是 $floor(\frac x {2^{m+1}})*2^m + x mod(2^m)$
- 根据蝶形变换下标的定义可知 左边的输入的下标就是蝶形变换的序号
- 右边的输入的下标是 $floor(\frac x {2^{m+1}})*2^m + x mod(2^m) + \frac N 2$
- 输出的正负性是 $x mod 2^{m+1}<2^m$ 为正，否则为负。(一般这种实现，不需要利用正负性，直接利用公式$(7)$的第二种形式即可)