在软件上可以轻易地计算上百万点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. 物理直觉
假设原始信号（复信号）$x[n]$长度为$N$，采样率为$f_s$。将信号序列理解为按行填充的矩阵，形状为$(N_1 \times N_2)$，那么矩阵的第$i$行表示$N_2$段序列中的第$i$段，矩阵的第$j$列表示$N_2$段序列中第$j$个数据构成的集合。
我们观察一下矩阵的列，它们实际上是在原序列中每隔$N_2$个点取值得到的，相当于对原序列降采样，将采样率由$f_s$降低为$f_s / N_2$。降采样会造成频谱混叠，原始频谱会折叠在一起。如下图所示，红色为非混叠频谱，蓝色、绿色、紫色则表示降采样导致的原始频谱的延拓，黄色区域表示降采样后的数据的FFT所能表示的频率范围，在该频率范围内发生了频谱混叠（后文称延拓后的频谱为混叠频谱）。
对矩阵的列做FFT，得到的就是降采样后的频谱。我们把降采样后的频谱画出来，如下图所示，在每一个频率点上的复向量，由$N_2$个复向量（我这里的例子$N_2 = 4$）复合而成，它们分别代表延拓后的每一个频谱（混叠频谱）在当前频率点上的复向量。
对矩阵的每一列做FFT之后，得到一列列频谱，这些频谱之间的差别在于，每一列时域数据在做FFT之前其初始相位（初相）都不同。矩阵的某一行（对应某个频率点）表示的是降采样后的多份频谱在同一个频率点处的复向量，但是由于每一列降采样时域序列的初相不同，不同的降采样频谱的复向量，在该频率点处存在差异。
如下图所示，3幅图分别表示降采样频谱中的3个不同频率点$k_1～k_3$的复向量，其中每一幅图都是一个复平面，用于展示该频率点上的复向量。在复平面中编号为1-4的复向量分别表示4个不同初相（均匀时间间隔），对应矩阵中该频率点所在的一行。
我们首先只考虑非混叠分量，记作$\vec V_1$，矩阵在按列做完FFT之后，频率点和初相索引分别为$[k, n]$的复向量表示为
$$ \vec V_1 (k, n) = \vec V_1 (k, 0) \cdot e^{j 2\pi (\frac k{N_1} \cdot \frac {f_s}{N_2})\cdot n \Delta t} $$
其中，$\frac {f_s}{N_2}$为降采样后的等效采样率（原始信号矩阵中某一列数据的采样率），$\frac k{N_1}$为第个频率点的频率因子（用频率因子乘以等效采样率，就是该点的频率），$\Delta t = \frac 1{f_s}$为时域数据序列的采样时间间隔，化简后可得
$$ \vec V_1 (k, n) = \vec V_1 (k, 0) \cdot e^{j \frac {2\pi}{N_1 N_2} nk} $$
我们把第$k$行所有的复向量都消除掉初相造成的相位旋转，也就说把$\vec V (k, n)$的初相与$\vec V(k, 0)$的初相对齐，那么就要给每一个元素都乘以相位旋转因子$W_{N_1 N_2}^{nk}$
$$ W_{N_1 N_2}^{nk} = (e^{-j \frac {2\pi}{N_1 N_2}})^{nk} $$
这样一来，对于第$k$行的任意一个复向量，它的频率为$k \frac {f_s}{N_2}$的分量都会被旋转到与当前行第一个元素的该频率分量完全相同的初相，成为不随时间变化的直流分量。
现在我们考虑其他混叠分量。矩阵的第$[k, n]$个复向量实际上是混叠和非混叠分量的和
$$ \vec V (k, n) = \sum_i \vec V_i (k_i, n) = \sum_i \vec V_i (k_i, 0) e^{j \frac {2\pi}{N_1 N_2} n k_i} $$
其中，$i = 1$表示非混叠分量，正如上文提到的；$i \ne 1$表示混叠分量，根据混叠的数学原理，混叠分量的频率为非混叠分量的倍频，即$k_i = i \cdot k$ ，在我们的例子中，$i = 1, 2, 3, 4$，那么该复向量经过相位旋转后，结果为
$$ \vec V (k, n) \cdot W_{N_1 N_2}^{nk} = \vec V_1 (k, n) e^0 + \vec V_2 (2k, n) e^{j \frac {2\pi}{N_1 N_2} nk} + \vec V_3 (3k, n) e^{j \frac {2\pi}{N_1 N_2} n2k} + \vec V_4 (4k, n) e^{j \frac {2\pi}{N_1 N_2} n3k} $$
在乘以相位旋转因子后，混叠分量虽然也会被旋转，但因为它们的频率比非混叠分量更高，同样的旋转不足以让他们的初相对齐，所以它们旋转后仍然具有随时间线性变化的相位分布。如下图所示，左图表示第$k$行第一个元素的所有频率分量（红色表示非混叠频率，蓝色、绿色和紫色表示混叠频率），右图表示改行的其他某个元素的所有频率分量，虚线标注了左图每个分量的相位，作为参考相位。对右图所有分量乘以相位旋转因子（基于非混叠频率计算出的旋转因子）后，可以发现只有非混叠分量的相位与参考相位（左图）完成了对齐，但其他混叠分量经过相同的相位旋转后并未与各自的参考相位对齐，所以混叠分量仍然存在随时间变化的波动项（时域信息）。
为了把矩阵里同一行中的混叠在一起的不同频率分量区分开，我们对每一行再做一次FFT，结果就是所有混叠分量被完美地分离出来，其中非混叠分量作为直流分量保存在了每一行的第一列，第$n$个混叠分量位于每一行的第$n$列。根据$\vec V (k, n)$的表达式和FFT的数学定义，可以严格证明这个结果。如下图所示，我们把完成第二次FFT之后的矩阵，进行转置，以行优先的顺序读取矩阵的每个元素，重新排列为一个行向量，最终就得到了我们在最开始要的原始信号的FFT频谱。
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，可以发现它们是一致的。
5. 推广
5.1 基2-FFT
如果我们把长度为$2N$的时域序列分解为$(2 \times N)$的矩阵，那经过以下过程就可以完成FFT。
- 对每列计算2点FFT，即蝶形运算。
- 对矩阵的每个元素乘以相位旋转因子。
- 对每行进行N点FFT。
- 转置矩阵，重排为行向量。
如果$N$本身也是偶数，就可以继续分解。
我们可以很自然地想到一个推广：一个长度为$N$（$N$是2的幂）的时域序列，可以进行$log_2^N$次递归分解，直到分解为最基本的蝶形运算。这个就是基2-FFT的数学原理，是Cooley-Tukey分解的特殊形式。
5.2 基N-FFT
更进一步推广，一个长度为$N = n^x$的时域序列，可以进行$log_n^N = x$次递归运算，完成FFT。我们可以把它称为基N-FFT。