任意长度FFT

针对长度 $n$ 为合数和质数的情况，分别使用 Cooley-Tukey 和 Rader 算法进行实现。

Cooley-Tukey

---

Cooley-Tukey 是一种求解长度为合数的 FFT 算法。基 2 的 FFT 是一种特殊的 Cooley Tukey，简单回顾一下：

$\begin{aligned} &(E_0,E_1,\dots,E_{n/2-1})=\text{FFT}(x_0,x_2,\dots,x_{n-2})\\ &(O_0,O_1,\dots,O_{n/2-1})=\text{FFT}(x_1,x_3,\dots,x_{n-1})\\ &\text{for}\ k=0,1,\dots,n/2-1:\ X_k=E_k+w_{n}^kO_k\\ &\text{for}\ k=n/2,n/2+1,\dots,n-1:\ X_k=E_{k-n/2}+w_{n}^kO_{k-n/2} \end{aligned}$

简单来说，就是将长度为 $n$ 的 DFT 拆分成两个长度为 $n/2$ 的 DFT，从而使得计算的时间复杂度降低。

Cooley-Tukey 也是使用了这种分而治之的思想，根据长度的因数对序列进行拆分。假设 $x$ 的长度为 $N$，其中 $N=N_1N_2$，那么 DFT 的公式为

$X_k=\sum_{n=0}^{N-1}\exp(-i\frac{2\pi}{N}nk)x_n$

将以下两个式子代入上式

$k=N_2k_1+k_2,\ k_1\in\{0,1,\dots,N_1-1\},\ k_2\in\{0,1,\dots,N_2-1\}\\ n=N_1n_2+n_1,\ n_1\in\{0,1,\dots,N_1-1\},\ n_2\in\{0,1,\dots,N_2-1\}$

可以得到

$\begin{aligned} X_{N_2k_1+k_2}&=\sum_{n_1=0}^{N_1-1}\sum_{n_2=0}^{N_2-1}\exp\left(-i\frac{2\pi}{N_1N_2}(N_1n_2+n_1)(N_2k_1+k_2)\right)x_{N_1n_2+n_1}\\ &=\sum_{n_1=0}^{N_1-1}\exp\left(-i\frac{2\pi}{N_1}n_1k_1\right)\exp(-i\frac{2\pi}{N}n_1k_2)\sum_{n_2=0}^{N_2-1}x_{N_1n_2+n_1}\exp\left(-i\frac{2\pi}{N_2}(n_2k_2)\right) \end{aligned}$

分析一下以上式子，$\sum\limits_{n_2=0}^{N_2-1}x_{N_1n_2+n_1}\exp\left(-i\dfrac{2\pi}{N_2}(n_2k_2)\right)$ 正是长度为 $N_2$ 的 DFT，对应的序列为 $x_{N_1n_2+n_1},\ n_2\in\{0,1,\dots,N_2-1\}$

记 $y_{n_1,k_2}=\sum\limits_{n_2=0}^{N_2-1}x_{N_1n_2+n_1}\exp\left(-i\dfrac{2\pi}{N_2}(n_2k_2)\right),z_{n_1,k_2}=\exp(-i\dfrac{2\pi}{N}n_1k_2)y_{n_1,k_2}$，则有

$X_{N_2k_1+k_2}=\sum_{n_1=0}^{N_1-1}\exp\left(-i\frac{2\pi}{N_1}n_1k_1\right)z_{n_1,k_2}$

这正是长度为 $N_1$ 的 DFT。

根据以上的分析，对于 $N$ 为合数的情况，我们可以先做 $N_1$ 个长度为 $N_2$ 的 DFT，然后每个元素乘上一个系数，再做 $N_2$ 个长度为 $N_1$ 的 DFT。那么时间复杂度可由下式表示：

$T(N)=N_1T(N_2)+N_2T(N_1)+\Theta(N)$

具体的算法步骤如下：

• 将长为 $N$ 序列（向量）按列排列变为 $N_1\times N_2$ 的矩阵
• 对矩阵的每一行分别做 DFT
• 矩阵的每一个元素乘上系数 $\exp(-i\dfrac{2\pi}{N}n_1n_2)$，$n_1,n_2$ 分别表示该元素所处的行和列（从 0 开始）
• 对矩阵的每一列分别做 DFT
• 将 $N_1\times N_2$ 矩阵按行排列变为长为 $N$ 的向量

推导思路如下：

记 $x[n_1][n_2]=x[N_1n_2+n_1]$，对每一行进行 DFT 变换，即固定 $n_1$，得到 $X_1[n_1][k_2]$，乘上系数后得到 $X_2[n_1][k_2]$，再对每一列进行 DFT 变换，固定$k_2$，得到 $X[k_1][k_2]=X[N_2k_1+k_2]$，所以最后按行排列即可得到原序列对应的 DFT。

Rader

---

对于长度为质数的情况，就不能将其简单的拆分了。注意到假如 $N$ 为质数的话，那么 $N-1$ 一定为合数，如果能把一个元素单独抽出来，那么剩下就可以用 Cooley Tukey 进行求解了。

Rader 算法需要用到数论的知识，首先介绍相关的知识。由于相关内容超出这门课的范围，就不详细展开证明，只进行叙述。

假设 $N$ 为质数，定义 $a,b$ 之间的运算 $\cdot$ 为 $(ab)\mod N$，则 $\{0,1,\dots,N-1\}$ 和 $\cdot$ 可以构成一个群，这个群叫做整数模 $N$ 乘法群。根据数论知识，这样的群存在一个整数 $g$，也叫做原根，使得对于任意非零的 $n$ 都有唯一对应的 $q\in\{0,1,\dots,N-2\}$，对应关系为 $n=g^{p}\mod N$。同样对于任意非零的 $k$ 都有唯一对应的 $p\in\{0,1,\dots,N-2\}$，使得 $k=g^{-p}\mod N$。此处，$k=g^{-p}\mod N$ 等价于 $kg^{p}\mod N=1$。

有了以上的结论，就可以将 $\{1,\dots,N-1\}$ 映射到 $\{0,\dots,N-1\}$。 将上述 $n,\ k$ 代入 DFT 的公式中，单独写出第 0 项，可以得到

$\begin{aligned} X_0&=\sum_{n=0}^{N-1}x_n\\ X_{g^{-p}(\text{mod}N)}&=x_0+\sum_{q=0}^{N-2}x_{g^{q}(\text{mod}N)}\exp\left(-i\frac{2\pi}{N}g^{q-p}(\text{mod}N)\right) \end{aligned}$

令 $a_q=x_{g^{q}(\text{mod}N)}$，$b_q=\exp\left(-i\dfrac{2\pi}{N}g^{-q}(\text{mod}N)\right)$，$q\in\{0,\dots,N-2\}$，那么 $\sum\limits_{q=0}^{N-2}x_{g^{q}(\text{mod}N)}\exp\left(-i\dfrac{2\pi}{N}g^{q-p}(\text{mod}N)\right)$ 实际上是 $\{a_q\}$ 和 $\{b_q\}$ 的循环卷积。根据循环卷积定理，有

$a*b=\check{(\hat{a}\cdot\hat{b})}$

这样循环卷积就可以转化为两个 DFT 和 一个 IDFT，这就转化成了我们熟悉的问题。

对于这种情况，为了方便进行 DFT，可以将其填充到 2 的幂次方，使用 基-2 的 FFT 算法进行求解。如果不将其补到 2 的幂次，使用 Rader 算法进行求解的过程中可能还会产生大的质数，这种情况下时间也会很大。

对于长度为 $N$ 的两个序列，按以下方法进行补充：

• 首先找到一个最小的 $m\in\Z$，使得 $M=2^{m}>2N$
• 令 $\tilde a=[a_0,\underbrace{0,\dots,0}_{M-N},a_1,\dots,a_{N-1}],$ $\tilde b=[b_0,\dots,b_{N-1},b_0,\dots]$ (一直重复 $b$，直到长度为 $M$，对超出的部分直接截断)
• 求 $\tilde a*\tilde b$，并取前 $N$ 项，则得到 $a*b$

证明：

$(\tilde a*\tilde b)_k=\sum\limits_{j=0}^{M-1}\tilde{a}_j\tilde b_{k-j}=a_0b_k+\sum\limits_{j=M-N}^{M-1}\tilde a_j\tilde b_{k-j}=a_0b_k+\sum\limits_{j=1}^{N-1}a_j\tilde b_{k+N-j}=\sum\limits_{j=0}^{N-1}a_jb_{k-j}=(a*b)_k$