FFT 是科学计算的基石
FFT 将 DFT 从 O(m²) 降至 O(m log m)，是 20 世纪最重要的算法之一。掌握 FFT 的原理和优化是并行计算研究者的基本功。
CS267 · 第 20 课
快速傅里叶变换（FFT）是科学计算中最重要的算法之一。它将 DFT 的复杂度从 O(m²) 降至 O(m log m)，广泛应用于信号处理、偏微分方程求解、大整数乘法等领域。本课深入探讨 FFT 的数学原理、串行与并行实现、通信复杂度理论以及现代自动调优框架。
离散傅里叶变换是连续傅里叶变换的离散版本，是分析离散信号频域特性的核心工具。一维 DFT 将一个 m 元素向量 v 变换为另一个 m 元素向量 F·v，其中变换矩阵 F 定义为：
F[j,k] = ω^(jk), 0 ≤ j, k ≤ m-1这里 ω = e^(2πi/m) = cos(2π/m) + i·sin(2π/m) 是 m 次单位根（m-th root of unity），满足 ω^m = 1。矩阵 F 是 Vandermonde 矩阵的特殊形式，具有对称性和正交性。
二维 DFT 作用于 m×m 矩阵 V，定义为 F·V·F——先对每列做一维 DFT，再对每行做一维 DFT。高维 DFT 以此类推，每个维度独立变换。这种可分离性使得高维 FFT 可以通过多次一维 FFT 实现。
当 m=4 时，ω = i，ω² = -1，ω³ = -i，ω⁴ = 1。此时 F[j,k] = i^(jk)，构成 4×4 的对称变换矩阵。DFT 可以理解为将信号从时域（或空域）转换到频域。
逆 DFT（IDFT）与 DFT 几乎相同，仅差一个共轭和归一化：
F^(-1) = complex_conjugate(F) / m这意味着实现一个 FFT 核就能同时支持正逆变换。此外，正弦变换（sin transform）是 F 的虚部，余弦变换（cos transform）是 F 的实部，它们用于求解泊松方程等应用。
快速傅里叶变换之所以成为"最重要的算法"，源于它在众多科学和工程领域的广泛应用：
在信号处理中，FFT 将时域信号转换到频域。通过将频域中幅度较小的分量置零（低通/高通滤波），再逆变换回时域，即可实现噪声去除、信号增强。例如，含噪信号 sin(7t) + 0.5·sin(5t) 在 128 个采样点上，通过将 FFT 系数小于 0.25 的分量置零，可以有效滤除高频噪声。
二维 FFT 将图像从像素空间转换到频域。JPEG 压缩的核心思想类似：保留 FFT 系数中幅度最大的 2.5%（低频分量），舍弃高频分量，从而实现 40:1 的压缩比而视觉质量损失很小。对于 200×320 的图像，仅保留最大的 2.5% 分量即可重建出视觉可接受的图像。
Schönhage-Strassen 算法利用 FFT 实现大整数乘法，复杂度为 O(b·log b·log log b)，其中 b 是位数。这是目前已知最快的渐近大数乘法算法之一，在密码学（如 RSA）中有重要应用。
泊松方程 ∇²u = f 是描述静电势、热传导、流体运动等现象的基本偏微分方程。FFT 提供了近乎最优的求解方法。
一维泊松方程离散化后得到三对角矩阵 L₁，其特征分解为：
L₁ = F · D · F^T其中 F 是正弦变换（FFT 的虚部），F(j,k) = √(2/(n+1))·sin(jkπ/(n+1))，D 是对角特征值矩阵，D(j,j) = 2(1 - cos(jπ/(n+1)))。
二维泊松方程 L₂X = B 可以写成 L₁·X + X·L₁ = B 的形式。代入特征分解后，算法分为三步：
总复杂度为 2 次二维 FFT 加上 n² 次加法和除法，即 O(n² log n)，其中 n×n 是网格大小。
FFT 方法的 O(N log N) 复杂度（N = n² 未知量）与 Multigrid 相当，远优于稠密 LU（O(N³)）和带状 LU（O(N²)）。在 PRAM 模型下，FFT 仅需 O(log N) 时间。
| 算法 | 串行复杂度 | PRAM 复杂度 | 内存 | 处理器数 |
|---|---|---|---|---|
| Dense LU | N³ | N | N² | N² |
| Band LU | N² | N | N^(3/2) | N |
| Jacobi | N² | N | N | N |
| Conjugate Gradient | N^(3/2) | N^(1/2)·log N | N | N |
| FFT | N·log N | log N | N | N |
| Multigrid | N | (log N)² | N | N |
| Lower bound | N | log N | — | — |
直接计算 DFT 需要 O(m²) 次运算。FFT 的核心思想是将求值问题转化为多项式的分治求值。
DFT 的第 j 个输出可以写成：
(F·v)[j] = Σ(k=0..m-1) ω^(jk)·v(k) = V(ω^j)其中 V(x) = Σ(k=0..m-1) x^k · v(k) 是一个 m-1 次多项式。DFT 等价于在 m 个单位根处求多项式的值。
将多项式 V(x) 按奇偶次项拆分：
V(x) = V_even(x²) + x · V_odd(x²)其中 V_even 包含偶次项，V_odd 包含奇次项。两者都是 m/2 次多项式。关键观察是：我们只需要在 m 个点 (ω^j)² 处求值，但由于 (ω^(j+m/2))² = (ω^j)²，这 m 个点只有 m/2 个不同的值！
FFT(v, ω, m): // m 是 2 的幂
if m = 1 return v[0]
else
v_even = FFT(v[0:2:m-2], ω², m/2)
v_odd = FFT(v[1:2:m-1], ω², m/2)
ω_vec = [ω^0, ω^1, ..., ω^(m/2-1)] // 预计算
return [v_even + (ω_vec .* v_odd),
v_even - (ω_vec .* v_odd)]递推关系 T(m) = 2T(m/2) + O(m) 的解为 T(m) = O(m log m)。相比直接 DFT 的 O(m²)，FFT 将复杂度降低了一个数量级以上。当 m = 10⁶ 时，FFT 约需 2000 万次运算，而直接 DFT 需 10¹² 次。
递归 FFT 的调用树是一棵 log m 层的完全二叉树。迭代实现自底向上遍历这棵树，避免了函数调用开销。
迭代算法从最小的子问题（大小为 2 的 DFT）开始，逐层合并。每一层将两个较小的 DFT 合并为一个更大的 DFT。这种实现方式天然适合循环结构，是现代 FFT 库的主流选择。
迭代 FFT 的输出按"位反转"（bit-reversed）顺序排列。例如，当 m=8 时，输出索引的二进制表示为：
| 自然顺序 | 二进制 | 位反转 | 对应输出 |
|---|---|---|---|
| 0 | 000 | 000 | 0 |
| 1 | 001 | 100 | 4 |
| 2 | 010 | 010 | 2 |
| 3 | 011 | 110 | 6 |
| 4 | 100 | 001 | 1 |
| 5 | 101 | 101 | 5 |
| 6 | 110 | 011 | 3 |
| 7 | 111 | 111 | 7 |
许多 FFT 应用（如泊松方程求解）可以按位反转顺序使用结果，无需额外排序。但需要自然顺序输出的应用必须执行一次位反转排列，这类似于小规模的转置操作。
实用 FFT 算法通常结合递归和迭代：顶层使用递归以利用内存层次结构（cache 友好），底层使用迭代以减少函数调用开销。FFTW 等库正是采用这种混合策略。
FFT 的并行化面临数据依赖的挑战——蝶形运算（butterfly）模式决定了哪些计算可以并行执行。
FFT 的数据依赖图是一个 log m 层的蝶形网络。在 PRAM（并行随机存取机）模型下，每层的蝶形运算完全并行，因此仅需 O(log m) 时间。但这忽略了通信成本。
在分布式内存系统中，数据分布策略直接影响通信量：
| 布局 | 局部性 | 通信阶段 |
|---|---|---|
| Block Layout | 前 log(m/p) 步局部 | 后 log p 步需通信 |
| Cyclic Layout | 前 log(m/p) 步局部 | 后 log p 步需通信 |
无论哪种布局，通信都集中在 log p 个阶段中。每个阶段需要传输 m/p 个字。
在 LogP 模型下（参数：延迟 α，带宽 β，每处理器计算速率 1）：
Time(FFT) = 2m·log(m)/p ... 完全并行计算
+ log(p) · α ... log p 个通信阶段
+ m·log(p)/p · β ... 每阶段传输 m/p 字Transpose 算法通过重新组织通信，将 FFT 的通信集中到一次全局转置中，从而优化整体性能。
所有通信集中在转置步骤中。将向量视为 m/p × p 的二维数组，转置就是交换行列索引。
Time(transposeFFT) = 2m·log(m)/p ... 计算
+ (p-1) · α ... 通信延迟
+ m·(p-1)/p² · β ... 通信带宽相比普通 FFT 的 log(p)·α 延迟，Transpose 仅需 (p-1)·α。当通信可以流水线化时，延迟项可进一步降至 α。带宽项从 m·log(p)/p 降至 m·(p-1)/p² ≈ m/p，也显著改善。
FFT 的并行通信复杂度已有严格的理论下界，Transpose 算法达到了这些下界。
在 cache 大小为 M 的串行计算机上，对规模为 m（m > m）的 FFT，必须在主存和 cache 之间传输的数据量为：
W = Θ(m·log(m) / log(M))证明基于每个数据字在 cache 中最多被重用 log M 次的观察。Transpose 算法达到此下界，这也启发了 cache-oblivious 递归算法的设计。
在 p 个处理器之间执行规模为 m 的 FFT，必须传输的数据量为：
W = Θ(m·log(m) / (p·log(m/p)))证明假设无重计算（no recomputation），且与局部内存大小无关（要求局部内存 ≥ m/p）。
当 log(m/p) ≥ log(p)（即 m ≥ p²）时，2 ≥ log(m)/log(m/p) ≥ 1，因此：
W_transpose = Θ(m/p) = Θ(m·log(m)/(p·log(m/p)))Transpose 算法达到了通信复杂度的理论下界！
| 下界 | 提出者 | 年份 | 适用场景 |
|---|---|---|---|
| Θ(m·log m / log M) | Hong-Kung | 1981 | 串行 cache 模型 |
| Θ(m·log m / (p·log(m/p))) | Aggarwal-Chandra-Snir | 1990 | 并行分布式内存 |
二维和三维 FFT 是科学计算（如流体力学、分子动力学）的核心，其并行化策略比一维更复杂。
| 策略 | 描述 | 通信需求 |
|---|---|---|
| (1) 2D Block Layout | 矩阵二维分块，行/列 FFT 均需通信 | 高（行列均需通信） |
| (2) Block Row + Transpose | 行 FFT 局部执行，转置后列 FFT 局部执行 | 一次全局转置 |
| (3) Block Row + Parallel Columns | 行 FFT 并行，列 FFT 也并行 | 两次通信 |
策略 (2) 最优，因为通过一次转置将通信集中，且可以与计算重叠。
FFT 需要至少一次全局转置操作，每个处理器需向其他所有处理器发送 1/p 的数据。Bisection Bandwidth（对分带宽）是网络最窄截面的带宽，限制了转置的性能：
LogGP 模型扩展了 LogP，增加了重叠能力。现代网络（带 DMA）允许通信与计算并行执行。GASNet 等单边通信库利用 RDMA 实现高效 one-sided 通信，相比 MPI 的双边通信，在中型消息上可获得 2.2 倍带宽提升。
三维 FFT 是实际应用中最常见的场景，其并行化涉及全局转置和通信-计算重叠。
将 NX × NY × NZ 的网格在 Z 维度上分为 P 份，每处理器持有 NZ/P 个 NX×NY 平面。前两个维度的 FFT 完全局部，第三维度的 FFT 需要通信。
Step 1: 列方向 FFT（全部局部）
Step 2: 行方向 FFT（全部局部）
Step 3: Z 方向 FFT（需全局转置）全局转置有三种不同的消息粒度：
| 算法 | 消息数 | 单消息大小 | 重叠能力 |
|---|---|---|---|
| Packed Slabs (3D) | p-1 | 大（如 512KB） | 无 |
| Slabs (2D) | NZ/P × (p-1) | 中（如 64KB） | 中等 |
| Pencils (1D) | (NX·NZ/P) × (p-1) | 小（如 8KB） | 最大 |
每处理器将所有目标数据打包成 p-1 个大消息，一次性发送。对于 512³ 网格、64 处理器，每处理器发送 63 个 512KB 的消息。这种方式无法与计算重叠，但实现简单。
将数据分成 NZ/P 个平面，每个平面内按行分成 p 个 slab。FFT 完一个 slab 后立即发送。对于 512³ 网格，每处理器发送 512 个 64KB 的消息。可以在不同平面的计算和通信之间重叠。
最细粒度：每完成一行（pencil）就立即发送。对于 512³ 网格，每处理器发送 4096 个 8KB 的消息。通信事件充分分散在计算中，实现最大重叠。
在 NAS FT 基准测试中，UPC（GASNet）实现的 Pencils 算法比原始 Fortran/MPI 实现快约 1.9 倍。MPI 在处理 Pencils 的大量小消息时效率显著下降（甚至崩溃），而 GASNet 的 one-sided RDMA 能充分利用网络带宽。
GASNet 提供 one-sided put/get 通信，无需远端 CPU 参与。在多个平台上对比：
FFTW（Fastest Fourier Transform in the West）是一个自动调优的 FFT 库，由 Matteo Frigo 和 Steven Johnson 开发，获得了 1999 年 Wilkinson 数值软件奖。
FFTW 实现了多种 FFT 算法（Cooley-Tukey、Prime Factor、Rader、Bluestein 等）。Planner 在实际运行前测量不同算法组合的速度，选择最优方案。这类似于 PHIPAC/ATLAS 在稠密线性代数中的自动调优策略。
传统 FFT 采用广度优先遍历（breadth-first），对整个数组做 log n 遍扫描。FFTW 使用自然递归（深度优先），当子问题小到足以放入 cache 时自动停止划分。这种cache-oblivious 策略对任何 cache 大小都有效，且对所有 cache 层级同时最优。
递归的基例是"codelet"——针对小规模 DFT（通常 ≤ 64）的高度优化代码。Codelet 由专用的离线"编译器"生成，充分利用 SIMD 指令（如 SSE2）和寄存器级优化。
fftw_complex *in, *out;
fftw_plan p;
in = fftw_alloc_complex(n);
out = fftw_alloc_complex(n);
// 创建"plan"——FFTW 自动选择最优算法
p = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_MEASURE);
// 多次执行同一大小的变换
fftw_execute(p);
fftw_execute(p);
...
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);Cache-oblivious 算法不需要知道 cache 大小就能达到 cache 最优。Singleton（1967）和 Frigo 等人（1999）证明，递归分治天然具有 cache-oblivious 特性——当递归子问题小到能放入 cache 时，自动获得高命中率，无论 cache 大小如何。
FFTW 在各种平台上表现优异：
Spiral 项目将 FFT 的自动调优思想推广到更广泛的 DSP（数字信号处理）算法领域。
Spiral 不仅为软件实现自动调优，还支持硬件（FPGA/ASIC）的自动生成。它使用形式化的数学变换规则空间和搜索/学习算法来发现最优实现。
更多信息请访问 www.spiral.net，提供在线代码生成器和相关论文。
Q1：当 m=4 时，单位根 ω 的值是多少？
Q2：分治 FFT 的递推关系 T(m) = 2T(m/2) + O(m) 的解是什么？
Q3：在并行 FFT 的 Transpose 算法中，假设 m ≥ p²，通信数据量达到什么渐近复杂度？
Q4：FFTW 的 cache-oblivious 递归策略为什么对任何 cache 大小都有效？
Q5：在三维 FFT 的三种转置算法中，Pencils 算法相比 Packed Slabs 算法的主要优势是什么？
FFT 将 DFT 从 O(m²) 降至 O(m log m)，是 20 世纪最重要的算法之一。掌握 FFT 的原理和优化是并行计算研究者的基本功。
Transpose 算法将通信集中并达到理论下界。理解 Hong-Kung 和 ACS 下界有助于设计最优并行算法。
FFTW 和 Spiral 通过多算法组合、cache-oblivious 递归和自动代码生成，实现了接近峰值的性能，无需手动调优。
GASNet 的 RDMA put/get 在 FFT 的全局转置中显著优于 MPI，尤其在大量小消息场景下。理解通信模型对高性能 FFT 实现至关重要。