傅里叶变换解析
阅读说明
本文不涉及各个应用领域对傅里叶变换的解释，仅将其作为一种特定运算进行分析。
如未经特殊说明，则：
- 向量均视为列向量；
- 下标始终从$0$开始。
离散傅里叶变换
原始定义
令$\vec{x}=(x_0,x_1,…,x_{N-1})$为长度为$N$的向量，定义向量$X$：
$$
X_k = \sum_{n=0}^{N-1}x_n\cdot\omega^{nk}_{N}
$$
其中$\omega_N=\mathrm{e}^{-\mathrm{i}\frac{2\pi}{N}}$，称$X$为$x$的离散傅里叶变换（DFT, Discrete Fourier Transform）。
矩阵形式
令$N\times N$矩阵$W_N=[a_{i,j}]$，其中$a_{i,j}=\omega_N^{ij}$，则有：
$$
\begin{aligned}
X=W_N\vec{x}
\end{aligned}
$$
复杂度分析
不难看出此时计算DFT的时间复杂度与一个矩阵向量乘法操作相当，即为$O(n^2)$.
IDFT
在纯数学的视角中，一般会在上面提到的DFT定义的等式右侧乘一个系数$\frac{1}{\sqrt{N}}$，此时变换矩阵${W’}_N=\frac{W_N}{\sqrt{N}}$是一个酉矩阵，其逆矩阵等于自身的共轭转置${W’}^H$。
我们仍然沿用先前的定义，不难得到变换矩阵$W_N$的逆矩阵为$\frac{1}{N}W_N^H$.
于是我们考察逆离散傅里叶变换，则有
$$
\begin{aligned}
\vec{x}&=\frac{1}{N}W_N^HX\
&=\frac{1}{N}\overline{W_N}X\
&=\frac{1}{N}\overline{W_N\overline{X}}\
&=\frac{1}{N}\overline{\mathrm{DFT(\overline{X})}}
\end{aligned}
$$
于是，IDFT问题被归约为了一个DFT问题。
Cooley-Tukey算法
分治思想
Cooley-Tukey算法是一种计算离散傅里叶变换的算法，由于其复杂度性质比较优秀，一般也被称作快速傅里叶变换。
不妨设$N$为合数，此时有$N=N_1\times N_2$.
此时有：
$$
\begin{aligned}
X_k &= \sum_{n_1=0}^{N-1}x_{n_1}\cdot\omega^{n_1k}{N}\
&= \sum{n_1=0}^{N_1-1}\sum_{n_2=0}^{N_2-1}x_{n_2N_1 + n_1}\cdot\omega_N^{(n_2N_1 +n_1)k}\
&= \sum_{n_1=0}^{N_1-1}\omega_N^{n_1k}\sum_{n_2=0}^{N_2-1}x_{n_2N_1 + n_1}\cdot\omega_N^{n_2N_1k}\
\end{aligned}
$$
令$k=k_1+N_2k_2$，其中$k_1<N_2, k_2<N_1$则
$$
\begin{aligned}
X_{k_1+N_2k_2} &= \sum_{n_1=0}^{N_1-1}\omega_N^{n_1(k_1+N_2k_2)}\sum_{n_2=0}^{N_2-1}x_{n_2N_1 + n_1}\cdot\omega_N^{n_2N_1(k_1+N_2k_2)}\
&= \sum_{n_1=0}^{N_1-1}\omega_N^{n_1k_1}\omega_{N_1}^{n_1k_2}\sum_{n_2=0}^{N_2-1}x_{n_2N_1 + n_1}\cdot\omega_{N_2}^{n_2k_1}\
\end{aligned}
$$
记$Y(n_1,k_1)=\sum_{n_2=0}^{N_2-1}x_{n_2\cdot N_1 + n_1}\cdot\omega_{N_2}^{n_2k_1}$.
与DFT定义对比，不难发现这是对向量$\vec{x}$的一部分$(x_{n_1},x_{N_1+n_1},…,x_{(N_2-1)\cdot N_1 +n_1})$的$N_2$点DFT。
换而言之，即是将向量$\vec{x}$按$N_1$为步长切分为$N_1$个$N_2$长度的子向量并进行DFT.
现在考虑将$Y(n_1)$带入，则有
$$
\begin{aligned}
X_{k1+N_2k_2}&=\sum_{n_1=0}^{N_1-1}\omega_{N_1}^{n_1k_1}\omega_N^{n_1k_2}Y(n_1,k_1)
\end{aligned}
$$
不难注意到此时求和式是一个对$\omega^{n_1k_1}_N(Y(0, k_1),Y(1,k_1),…,Y(N_1-1,k_1))$的$N_1$点DFT.
于是我们可以将大规模的DFT归约为多次规模稍小的DFT，并采用上述公式对子问题进行合并，即可得到原问题的答案。
复杂度分析
不妨令单次DFT问题的时间复杂度记为$T_{\mathrm{DFT}}(N)$，其中$N$为问题规模。
我们将一个规模为$N$的DFT问题分解为规模分别为$N_1$和$N_2$的DFT问题，其中，规模为$N_1$的DFT问题需要计算$N_2$次，规模为$N_2$的DFT问题需要计算$N_1$次。
可以得到，算法的总复杂度为$T_{\mathrm{DFT}}(N)=N_2\times T_{\mathrm{DFT}}(N_1) + N_1 \times T_{\mathrm{DFT}}(N_2)+O(N_1)$
假设$N=2^k$，取$N_1=2$（此时记作基2的快速傅里叶变换），根据主定理有$T_{\mathrm{DFT}}(N)=O(N\log N)$.
Bailey’s FFT
算法思想与复杂度分析
Bailey’s FFT是Cooley-Tukey算法的变种，是一种对矩阵运算单元友好的傅里叶变换算法。
不妨设$N=a^2$，取$N_1=N_2=a$，则有：
$$
\begin{aligned}
T_{\mathrm{DFT}}(N)&=N_2\times T_{\mathrm{DFT}}(N_1) + N_1 \times T_{\mathrm{DFT}}(N_2) + O(N_1)\
&=2 a\times T_{\mathrm{DFT}}(a) + O(a)\
\end{aligned}
$$
不对规模为$a$的DFT进行进一步分治，采用DFT的定义直接计算，则有$T_{\mathrm{DFT}}(a)=O(N)$，故$T_{\mathrm{DFT}(N)}=O(N^{1.5})$.
考虑到$N$一般不是完全平方数，在渐近意义上最优需要取$N_1\approx N_2\approx \sqrt{N}$.
换而言之，Bailey’s FFT实际上就是仅执行一次分解的Cooley-Tukey算法。
为什么只分解一次？
很容易注意到，仅分解一次的FFT算法带来的是算法复杂度的退化。那么这样的分解有什么优势呢？
如果使用基2的FFT算法，会出现步长跨度极大的内存访问，观察下面FFT的分治过程：
- ${𝑥0,𝑥1,𝑥2,𝑥3,𝑥4,𝑥5,𝑥6,𝑥7}$
- ${𝑥0,𝑥2,𝑥4,𝑥6},{𝑥1,𝑥3,𝑥5,𝑥7}$
- $ {𝑥0,𝑥4},{𝑥2,𝑥6},{𝑥1,𝑥5},{𝑥3,𝑥7} $
- ${𝑥0},{𝑥4},{𝑥2},{𝑥6},{𝑥1},{𝑥5},{𝑥3},{𝑥7}$
在算法流程从第二层回归第一层时，需要计算4次2点的DFT操作，而每次DFT操作的两份数据的跨度为4。而这只是问题之一，多次的分解涉及大量的数据重排，这也需要大量的内存访问。在问题规模较大时，算法的内存局部性会比较差。
实践中，我们必须充分考虑现代计算机的性能瓶颈，即内存墙问题。内存墙问题的核心矛盾在于计算单元的增长速度远远超过内存带宽的增长速度，内存访问效率过低会导致计算单元长期处于空置状态。这种步长跨度极大的内存访问，容易造成缓存失效等问题，进一步导致数据供给速率的不足。
因此，我们需要一种内存局部性比较好的算法，以期望提高计算单元的利用率。
Bailey’s FFT将原始问题分解为$N_2$次$N_1$点的DFT和$N_1$次$N_2$点的DFT，分解次数显著减少，数据搬运和重排的操作会被显著地缩减。代价是计算单元的任务量从$O(N\log N)$显著地增加到了$O(N^{1.5})$.
如果计算单元是CPU，从这个角度来说，这样的代价是不可接受的。尽管CPU利用率高了，但是复杂度的退化使得计算效率反而下降。
所以Bailey’s FFT实际上并不是设计给常见的CPU使用的算法，它适用于具有大量并行计算单元的处理器。并行计算单元具有更高的计算能力，天生适合处理计算密度高的算法。而Bailey’s FFT正是用时间复杂度退化换取计算密度的一种算法，尽管这显著增加了计算单元的任务量，但是并行化带来的效率提升抵消了这部分代价的增长，加上显著地减少了数据搬运和重排，这使得任务整体的运算速率被提高了。
具体地，我们可以注意到公式：
$$
X_{k_1+N_2k_2} = \sum_{n_1=0}^{N_1-1}\omega_N^{n_1k_1}\omega_{N_1}^{n_1k_2}\sum_{n_2=0}^{N_2-1}x_{n_2N_1 + n_1}\cdot\omega_{N_2}^{n_2k_1}\
$$
内层$N_1$次$N_2$点的DFT实际上是相互独立的。
如果我们将$\vec{x}$重整为$N_1\times N_2$的矩阵$A=[a_{i,j}]$，其中$a_{i,j}=x_{i+jN_1}$，便可以将内层的DFT操作归约为一个矩阵乘法：
$$
Y=AW_{N_2}
$$
这可以视作对$A$的行向量作DFT。
设$N_1 \times N_2$矩阵$T=[t_{i,j}]$，其中$t_{i,j}=\omega_N^{ij}$，令$Y’=Y\odot T$.
不妨令最终结果为$N_1\times N_2$矩阵$R=[r_{i,j}]$，其中$r_{i,j}=X_{j+N_2i}$，于是有
$$
R=W_{N_1}Y’
$$
我们推导出了矩阵形式的Cooley-Tukey算法，实际上也就是Bailey’s FFT的核心流程，其关键点在于将核心的DFT操作归约为了两个矩阵乘法。在实践中，我们可以使用矩阵运算单元来完成这一算法中的两个矩阵乘法，使用向量处理单元来实现逐元素相乘，从而完成了对FFT的并行化加速。
再回到原始问题：为什么只分解一次？
先前提到Bailey’s FFT减少了分解次数，确保数据重排的开销更小。而实际上，它被发明出来主要是用于解决大内存跨步的问题。
讨论之前，我们需要知道一个前提，W矩阵是一个对称矩阵，转置是0开销的。
假设矩阵在存储系统中以行优先的顺序存储，显然按行读取会非常的高效。但是我们注意到在公式
$$
R=W_{N_1}Y’
$$
为了计算$R$矩阵的一个元素，我们需要读取$Y$的一列，此时就出现了大内存跨步的问题。
为了避免这一问题，我们可以将公式变形为
$$
R^T=Y’^TW_{N_1}^T
$$
通过转置操作，可以避免计算过程中出现大的内存跨步。
故核心在于我们将其归约为矩阵乘法操作，在此基础上利用转置操作来避免大内存跨步的访问。
但是我们注意到，转置操作也需要大内存跨步访问，应该怎么办呢？由于转置操作具有极好的可预测性和规律性，可以针对这一操作做单独的优化，与具体的体系结构相关。与本文主题无关，不深入讨论。
于是，我们确保了计算过程中不存在等待内存访问，使得计算单元被尽可能地填满。
实践中$N_1$和$N_2$的选取需要具体考察机器的并行数，选择合适的数字以确保尽可能提高运算单元的利用率。
Bluestein算法
算法思想
以上讨论的算法建立在$N$是一个可分解的合数的假设上，如果$N$是一个质数，那么显然上述算法会迅速退化成$O(N^2)$的朴素DFT。
于是我们需要引入一种新的算法——Bluestein算法。
考虑DFT原始定义：
$$
X_k = \sum_{n=0}^{N-1}x_n\cdot\omega^{nk}{N}
$$
注意到：
$$
nk=\frac{n^2+k^2-(k-n)^2}{2}
$$
于是有：
$$
\begin{aligned}
\omega_N^{nk}&=\omega_N^{\frac{n^2+k^2-(k-n)^2}{2}}\
&=\mathrm{e}^{-\mathrm{i}\frac{\pi}{N}(n^2+k^2-(k-n)^2)}\
&=\mathrm{e}^{-\mathrm{i}\frac{n^2\pi}{N}}\mathrm{e}^{-\mathrm{i}\frac{k^2\pi}{N}}\mathrm{e}^{\mathrm{i}\frac{(k-n)^2\pi}{N}}\
\end{aligned}
$$
代入DFT定义有：
$$
X_k =\mathrm{e}^{-\mathrm{i}\frac{k^2\pi}{N}}\sum{n=0}^{N-1}(x_n\cdot\mathrm{e}^{-\mathrm{i}\frac{n^2\pi}{N}})\mathrm{e}^{\mathrm{i}\frac{(k-n)^2\pi}{N}}
$$
不妨令$A_n=x_n\cdot\mathrm{e}^{-\mathrm{i}\frac{n^2\pi}{N}}, B_n =e^{-\mathrm{i}\frac{n^2\pi}{N}},C_n=\mathrm{e}^{\mathrm{i}\frac{n^2\pi}{N}}$，则有
$$
X_k =B_k\sum_{n=0}^{N-1}A_nC_{k-n}
$$
定义1：对无限长度的序列$A$和序列$B$，有$C_n=\sum_{k=-\infty}^{+\infty}A_kB_{n-k}$，则称$C$为$A$与$B$的离散线性卷积，记作$C=A\ast B$.
我们用 $((n))_N$ 表示取模运算，即满足 $n = kN + r$ 且 $0 \le r < N$ 的唯一整数 $r$。
定义2：对于长度均为$N$的有限长序列$A$和序列$B$，有$C_n=\sum_{k=0}^{N-1}A_kB_{((n-k))_N}$，则称$C$为$A$与$B$的循环卷积，记作$C = A\circledast B$。
卷积定理：令$C = A\circledast B$，则有$\mathrm{DFT}(C)=\mathrm{DFT}(A)\cdot\mathrm{DFT}(B)$.
证明：
令$A’=\mathrm{DFT}(A), B’=\mathrm{DFT}(B), C’=\mathrm{DFT}(C)$，于是有：
$$
\begin{aligned}
{C’}k&=\sum{n=0}^{N-1}C_n\omega_N^{nk}\
&=\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}A_mB_{((n-m))N}\omega_N^{nk}\
&=\sum{m=0}^{N-1}A_m\sum_{n=0}^{N-1}B_{((n-m))N}\omega_N^{nk}\
&=\sum{m=0}^{N-1}A_m\sum_{n=0}^{N-1}B_{n}\omega_N^{(n+m)k}\
&=\sum_{m=0}^{N-1}A_m\omega_N^{mk}\sum_{n=0}^{N-1}B_{n}\omega_N^{nk}\
&={A’}_k\cdot {B’}_k
\end{aligned}
$$
故有$C’=A’\cdot B’$。
证毕。
回到公式
$$
X_k =B_k\sum_{n=0}^{N-1}A_nC_{k-n}
$$
注意求和式的形式，如果我们令$A_i=0,\forall i\notin [0,N). C_i=0, \forall i\notin[-N+1, N-1]$，则有
$$
X_k =B_k\sum_{n=-\infty}^{+\infty}A_nC_{k-n}
$$
此时公式中的求和式本质上是一个线性卷积。
卷积定理向我们指出了使用DFT计算循环卷积的可行性，我们需要寻求一种方法将线性卷积转化为循环卷积。
定理：
设 $A$ 和 $B$ 为两个有限长序列，其长度分别为 $L$ 和 $M$。序列元素 $A_n$ 和 $B_n$ 分别在 $0 \le n \le L-1$ 和 $0 \le n \le M-1$ 范围内有定义，在定义域外均为 0。
令序列 $C$ 为 $A$ 和 $B$ 的线性卷积，即：
$$
C_n = (A * B)n = \sum{m=-\infty}^{+\infty} A_m B_{n-m}
$$
其非零元素的长度为 $P = L + M - 1$。
选定整数 $N$，使得 $N \ge P$。构造长度为 $N$ 的补零序列 $A’$ 和 $B’$：
$$
A’_n = \begin{cases} A_n & 0 \le n < L \ 0 & L \le n < N \end{cases} , \quad B’_n = \begin{cases} B_n & 0 \le n < M \ 0 & M \le n < N \end{cases}
$$
令序列 $D$ 为 $A’$ 和 $B’$ 的 $N$ 点循环卷积，即：
$$
D_n = (A’ \circledast B’)n = \sum{m=0}^{N-1} A’m B’{((n-m))_N}
$$
对于所有 $0 \le n \le N-1$，有：
$$
D_n = C_n
$$
即循环卷积的结果在主值区间内与线性卷积完全相等。
证明：
为了处理循环卷积中的模运算项 $B’_{((n-m))N}$，我们定义 $B$ 的周期延拓序列 $\tilde{B}$，其周期为 $N$：
$$
\tilde{B}k = \sum{r=-\infty}^{+\infty} B{k - rN}
$$
在 $0 \le k \le N-1$ 区间内，$\tilde{B}_k = B’_k = B_k$。
将此代入循环卷积 $D_n$ 的定义式：
$$
D_n = \sum_{m=0}^{N-1} A’m \tilde{B}{n-m}
$$
由于 $A’m$ 仅在 $0 \le m \le L-1$ 内等于 $A_m$，其余位置为 0，且 $L < N$，我们可以将 $A’m$ 替换为 $A_m$，并将求和范围扩展至 $(-\infty, +\infty)$ 而不改变结果值：
$$
D_n = \sum{m=-\infty}^{+\infty} A_m \tilde{B}{n-m}
$$
将 $\tilde{B}$ 的定义代入上式：
$$
\begin{aligned} D_n &= \sum_{m=-\infty}^{+\infty} A_m \left( \sum_{r=-\infty}^{+\infty} B_{n - m - rN} \right) \ &= \sum_{r=-\infty}^{+\infty} \left( \sum_{m=-\infty}^{+\infty} A_m B_{(n - rN) - m} \right) \end{aligned}
$$
观察括号内的求和式 $\sum_{m} A_m B_{(n - rN) - m}$，根据线性卷积的定义，这正是序列 $C$ 在索引 $n - rN$ 处的值 $C_{n - rN}$。
因此，我们得到了公式：
$$
D_n = \sum_{r=-\infty}^{+\infty} C_{n - rN}
$$
我们需要证明在 $n \in [0, N-1]$ 时，求和式中仅有 $r=0$ 一项非零。
由线性卷积性质可知，仅当 $0 \le k \le L + M - 2$ 时，$C_k \ne 0$。
对不等式$0\leq n-rN\leq L+M-2$，化简得$rN\leq n\leq L+M-2+rN$
显然$r=0$时，$C_{n-rN}\neq0$.
考虑不等式左半边，则有$rN\leq n$，在$r>0$时，该不等式不可能成立，故对于$r>0$的情况，$C_{n-rN}$恒等于$0$.
考虑不等式右半边，则有$n\leq L+M-2+rN$，如要使之不成立，则有$n>L+M-2+rN$始终成立。已知$n\geq0$，故有$0>L+M-2+rN$，此时仅考虑$r<0$的情况，有$N>\frac{2-L-M}{r}$始终成立，代入$r=-1$，有$N>L+M-2$始终成立。
由于 $N \ge L + M - 1$，在主值区间 $0 \le n \le N-1$ 内，所有混叠项 ($r \ne 0$) 均为 0。因此$D_n = C_n$.
证毕。
回到公式
$$
X_k =B_k\sum_{n=-\infty}^{+\infty}A_nC_{k-n}
$$
注意到求和式中 $C_{k-n}$ 的下标可能为负（范围为 $-(N-1)$ 到 $N-1$）。
回顾上述定理的证明过程，其本质是将有限长序列进行周期延拓，并在一个充分大的周期 $M$（$M \ge 2N-1$）内进行计算以避免混叠。
根据周期延拓的性质 $C_n = C_{n + M}$，线性卷积所需的负下标数据 $C_{-k}$，在循环卷积的周期性视窗中等价于正下标 $M-k$ 处的数据，即令$C_{M-k}=C_{-k}$。
因此，我们将 $C_n$ 中负下标的部分平移映射到数组的尾部。此时，周期延拓后的序列在运算区间内与原线性卷积核完全一致，从而保证了定理的成立。
值得注意的是，根据上面的定理，延拓的周期长度应满足$M\geq 3N-2$，然而在本算法的语境下，我们仅需要$M\geq2N-1$.
假设我们按照上述要求将$A$和$B$进行延拓，分别记为序列$A’$和$B’$.
我们注意到
$$
\begin{aligned}
Y_k &=\sum_{n=-\infty}^{+\infty}A_nC_{k-n}\
&=\sum_{n=0}^{M-1}A_nC_{((k-n))N}\
&=\sum{n=0}^{N-1}A_nC_{((k-n))_N}\quad (\forall i\geq N, A_i=0)\
\end{aligned}
$$
不难发现，序列$C$真正对结果产生贡献的索引区间为$[0,N-1]$和$[M-N+1,M-1]$，只要这两个区间交集为空即可，不难解得$M\geq2N-1$
我们对序列$A$和序列$C$进行零延拓，确保满足定理条件后，根据卷积定理有：
$$
X=B\cdot\mathrm{IDFT}(\mathrm{DFT}(A)\cdot\mathrm{DFT}(C))
$$
于是我们可以沿用原本的FFT算法计算任意长度的DFT了。
复杂度分析
显然，由于问题被归约为了三次规模为$M$的DFT操作，其算法复杂度取决于DFT操作的复杂度。
如果采用Radix-2的FFT算法，那么复杂度为$O(n\log n)$.
如果采用Bailey’s FFT，那么复杂度为$O(n^{1.5})$.