在软件上可以轻易地计算上百万点FFT，因为PC或服务器的内存足够大，但是在硬件电路（ASIC或FPGA）中，RAM非常稀有且昂贵，直接做如此大规模的FFT是不现实的。我们可以把大规模的FFT分解为两个维度上的小规模FFT，从而解决RAM不足的问题，这个算法叫做Cooley-Tukey算法。
1. 数学原理
对长度为$N$的原始信号（复信号）$x[n]$直接做离散傅里叶变换，结果如下
$$ X[f] = \sum_{n = 0}^{N-1} x[n] e^{-j\frac {2\pi}N fn} $$
将原始信号的数据序列长度分解为$N = N_1 \cdot N_2$，按列优先填充的顺序，将数据序列重排为$N_1 \times N_2$的矩阵。我们要计算的离散傅里叶变换序列，则是按行优先的顺序，排列为同样纬度$N_1 \times N_2$的矩阵。（为什么要这样定义？见后文物理直觉章节）
$$ \begin{matrix} \begin{bmatrix} x[0] & \dots & x[N - 1] \end{bmatrix} & \Rightarrow & \begin{bmatrix} x[0] & \dots & x[(N_2 - 1)N_1] \\ \vdots & & \vdots \\ x[N_1 - 1] & \dots & x[N_2 N_1 - 1] \\ \end{bmatrix} & \Rightarrow & \begin{bmatrix} X[0] & \dots & X[N_2 - 1] \\ \vdots & & \vdots \\ X[(N_1 - 1)N_2] & \dots & X[N_2 N_1 - 1] \\ \end{bmatrix} \\ time\ domain(N) & & time\ domain(N_1 \times N_2) & & frequency\ domain(N_1 \times N_2) \end{matrix} $$
信号数据序列的索引$n$可以表示为$n = p + q N_1$ ，傅里叶变换频谱的索引$f$可以表示为$f = rN_2 + s$，其中整数$p, r \in [0, N_1)$和$q, s \in [0, N_2)$分别表示矩阵的行和列。将展开后的索引代入到离散傅里叶变换中，
$$ X[rN_2 + s] = \sum_{p = 0}^{N_1 - 1} \sum_{q = 0}^{N_2 - 1} x[p + qN_1] e^{-j \frac {2\pi}{N_1 N_2} (pr N_2 + qs N_1 + ps + qr N_1 N_2)} $$
将指数中的常数和索引分开，适配求和次序，并化简，可以得到
$$ X[rN_2 + s] = \sum_{p = 0}^{N_1 - 1} \left\{ \left\{ \sum_{q = 0}^{N_2 - 1} x[p + qN_1] e^{-j \frac {2\pi}{N_2} qs} \right\} \cdot e^{-j \frac {2\pi}{N_1 N_2} ps} \right\} e^{-j \frac {2\pi}{N_1} pr} $$
其中，最内层的求和是对矩阵的行做DFT，定义原始矩阵按行做DFT之后的矩阵为$F$，$F[p, s]$表示其第$p$行第$s$列的值
$$ F[p, s] = \sum_{q = 0}^{N_2 - 1} x[p + qN_1] e^{-j \frac {2\pi}{N_2} qs} $$
定义相位旋转因子底数$W_N$为
$$ W_N = e^{-j \frac {2\pi}{N_1 N_2}} $$
那么原始DFT可化简为
$$ X[rN_2 + s] = \sum_{p = 0}^{N_1 - 1} ( F[p, s] \cdot W_N^{ps} ) e^{-j \frac {2\pi}{N_1} pr} $$
定义矩阵$F$的每个元素乘以相位旋转因子之后的矩阵为$G$，$G[s, p]$表示其第$p$行第$s$列的值
$$ G[p, s] = F[p, s] \cdot W_N^{ps} $$
代入到原始DFT中，可得
$$ X[rN_2 + s] = \sum_{p = 0}^{N_1 - 1} G[p, s] \cdot e^{-j \frac {2\pi}{N_1} pr} $$
而这其实是对矩阵$G$的每一列做DFT。将最终的矩阵$X$按行优先填充的方式还原为长度为$N = N_1 \times N_2$的序列，就得到了原始信号数据序列的DFT频谱。
$$ \begin{matrix} \begin{bmatrix} X[0] & \dots & X[N_2 - 1] \\ \vdots & & \vdots \\ X[(N_1 - 1)N_2] & \dots & X[N_2 N_1 - 1] \\ \end{bmatrix} & \Rightarrow & \begin{bmatrix} X[0] & \dots & X[N - 1] \end{bmatrix} \\ frequency\ domain(N_1 \times N_2) & & frequency\ domain(N) \end{matrix} $$
2. 物理直觉
TODO
3. 算法总结
- 将长度为N（N = R x C）的数据序列分解为R段长度为C的字序列，按照行优先的顺序排列成R x C的矩阵。
- 按照列优先的顺序对矩阵的每一列做R点FFT。
- 矩阵的每个元素（索引为r, c）乘以相位旋转因子${W_N}^{rc}$。
- 按照行优先的顺序对矩阵的每一行做C点FFT。
- 矩阵转置。
- 将矩阵的每一行拼接起来，还原出长度为N点频谱。
- 若上述过程中出现的FFT点数仍然过大，则递归进行二维分解，直至完成当前的FFT计算。
4. 代码示例
以完成2048点FFT为例，将其排列为为32 x 64的矩阵。
import numpy as np
def decompose_fft_demonstration():
fs = 2048.0
N = 2048
R, C = 32, 64
# signal
t = np.arange(N) / fs
signal = 1.0 * np.exp(1j * 2 * np.pi * 50 * t) + \
0.5 * np.exp(1j * 2 * np.pi * 120 * t) + \
0.1 * (np.random.randn(N) + 1j * np.random.randn(N))
# ==========================================
# just FFT
# ==========================================
fft_ref = np.fft.fft(signal)
# ==========================================
# Cooley-Tukey
# ==========================================
matrix = signal.reshape(R, C)
matrix_step1 = np.fft.fft(matrix, axis=0)
k_r = np.arange(R).reshape(-1, 1)
c = np.arange(C).reshape(1, -1)
twiddle = np.exp(-1j * 2 * np.pi * k_r * c / N)
matrix_step2 = matrix_step1 * twiddle
matrix_step3 = np.fft.fft(matrix_step2, axis=1)
fft_decomposed = matrix_step3.T.reshape(-1)
# ==========================================
# compare
# ==========================================
error = np.linalg.norm(fft_ref - fft_decomposed)
print(f"Total error: {error:.10e}")
is_close = np.allclose(fft_ref, fft_decomposed)
print(f"close? {'yes' if is_close else 'no'}")
return fft_ref, fft_decomposed
fft_ref, fft_dec = decompose_fft_demonstration()可以对比2个频谱fft_ref和fft_decomposed，可以发现它们是一致的。