1965 年，Cooley 和 Tukey 在 IBM 的一篇论文中重新发现了一个算法， 将离散傅里叶变换的复杂度从 O(n^2) 降到了 O(n log n)。 这个算法后来被称为快速傅里叶变换（Fast Fourier Transform, FFT）。 事实上，高斯早在 1805 年就使用了类似的技巧， 但直到冷战时期美国需要高效分析苏联核试验的地震信号， 这个算法才被正式发表并引起广泛关注。
今天，FFT 无处不在：手机通话的信号编码、JPEG 图片压缩、 MP3 音频格式、WiFi 的 OFDM 调制、CT 扫描图像重建、 金融领域的期权定价、竞赛编程中的大整数乘法……
本文将从离散傅里叶变换定义出发，推导 Cooley-Tukey 算法， 给出完整的 C 语言实现，讨论 NTT 在精确整数计算中的应用， 最后探讨 FFT 在实际工程中的应用场景和注意事项。
一、离散傅里叶变换的定义与朴素计算
1.1 DFT 的数学定义
给定长度为 n 的复数序列 x[0], x[1], …, x[n-1]， 其离散傅里叶变换（Discrete Fourier Transform, DFT）定义为：
X[k] = sum_{j=0}^{n-1} x[j] * omega_n^{jk}, k = 0, 1, ..., n-1其中 omega_n = e^{-2pii/n} 是 n 次单位根。逆变换（IDFT）为：
x[j] = (1/n) * sum_{k=0}^{n-1} X[k] * omega_n^{-jk}, j = 0, 1, ..., n-1从线性代数视角看，DFT 本质是矩阵-向量乘法 X = F_n * x， 其中 F_n 是 n 阶 DFT 矩阵，F_n[j][k] = omega_n^{jk}， 经 1/sqrt(n) 归一化后为酉矩阵。
1.2 单位根的关键性质
n 次单位根有几个至关重要的性质，正是它们让 FFT 成为可能：
消去引理（Cancellation Lemma）： omega_{dn}^{dk} = omega_n^k。 即 2n 次单位根的 2k 次方等于 n 次单位根的 k 次方。
折半引理（Halving Lemma）： 如果 n 是偶数，n 个 n 次单位根的平方的集合 恰好等于 n/2 个 n/2 次单位根的集合。
求和引理（Summation Lemma）： 对任意不被 n 整除的整数 k，有 sum_{j=0}^{n-1} omega_n^{jk} = 0。
1.3 朴素 O(n^2) 实现
直接按定义计算 DFT 需要两层循环：
#include <complex.h>
#include <math.h>
void dft_naive(double complex *x, double complex *X, int n) {
for (int k = 0; k < n; k++) {
X[k] = 0;
for (int j = 0; j < n; j++) {
double angle = -2.0 * M_PI * j * k / n;
X[k] += x[j] * cexp(I * angle);
}
}
}外层循环遍历 n 个频率分量，内层循环对 n 个时域样本求和。 总计算量 O(n^2) 次复数乘法和加法。 当 n = 10^6 时，约需 10^{12} 次运算，在现代 CPU 上需要数分钟。
二、Cooley-Tukey FFT：分治的力量
2.1 奇偶分解
Cooley-Tukey 算法的核心思想极其简洁。 假设 n 是 2 的幂，将输入按下标的奇偶性分成两半：
偶数下标：a[j] = x[2j], j = 0, 1, ..., n/2 - 1
奇数下标：b[j] = x[2j + 1], j = 0, 1, ..., n/2 - 1然后 DFT 可以写成：
X[k] = sum_{j=0}^{n-1} x[j] * omega_n^{jk}
= sum_{j=0}^{n/2-1} x[2j] * omega_n^{2jk}
+ sum_{j=0}^{n/2-1} x[2j+1] * omega_n^{(2j+1)k}
= sum_{j=0}^{n/2-1} a[j] * omega_{n/2}^{jk}
+ omega_n^k * sum_{j=0}^{n/2-1} b[j] * omega_{n/2}^{jk}
= A[k] + omega_n^k * B[k]其中 A[k] 和 B[k] 分别是 a 和 b 的 n/2 点 DFT。 利用周期性 A[k + n/2] = A[k]（因为 A 是 n/2 点 DFT），可得：
X[k] = A[k] + omega_n^k * B[k], k = 0, 1, ..., n/2 - 1
X[k + n/2] = A[k] - omega_n^k * B[k], k = 0, 1, ..., n/2 - 1这就是著名的蝶形运算（Butterfly Operation）。
2.2 递归结构与复杂度分析
上述分解给出了递推关系：
T(n) = 2 * T(n/2) + O(n)由主定理（Master Theorem），T(n) = O(n log n)。
具体来说，对于 n = 2^m 点的 FFT： - 共有 m = log2(n) 级 - 每级执行 n/2 次蝶形运算 - 每次蝶形运算需要 1 次复数乘法和 2 次复数加法 - 总计算量：(n/2) * log2(n) 次复数乘法
对比朴素算法：
| n | 朴素 O(n^2) | FFT O(n log n) | 加速比 |
|---|---|---|---|
| 2^10 = 1024 | 1,048,576 | 5,120 | 205x |
| 2^16 = 65536 | 4,294,967,296 | 524,288 | 8192x |
| 2^20 = 1048576 | 1.1 * 10^{12} | 10,485,760 | 104858x |
| 2^24 | 2.8 * 10^{14} | 201,326,592 | 1398101x |
2.3 递归版 FFT（伪代码）
def fft_recursive(x):
n = len(x)
if n == 1:
return x
# 按奇偶分组
even = fft_recursive(x[0::2])
odd = fft_recursive(x[1::2])
# 合并
omega = exp(-2j * pi / n)
w = 1
result = [0] * n
for k in range(n // 2):
result[k] = even[k] + w * odd[k]
result[k + n // 2] = even[k] - w * odd[k]
w *= omega
return result这个递归版本虽然正确且易于理解，但存在两个问题： 递归调用开销大，且不是原地运算。实际实现几乎都使用迭代版本。
三、蝶形图与位反转置换
3.1 蝶形图
蝶形图（Butterfly Diagram）是理解 FFT 计算流程的关键可视化工具。 下面是一个 8 点 FFT 的蝶形图：
蝶形图中的每一个交叉就是一次蝶形运算： 将两个输入值通过加法和乘以旋转因子（twiddle factor）得到两个输出值。
观察蝶形图可以发现几个重要特点： - 输入顺序不是自然顺序 0, 1, 2, 3, 4, 5, 6, 7， 而是 0, 4, 2, 6, 1, 5, 3, 7 - 每一级（stage）的蝶形运算跨度不同： 第 1 级跨度为 1，第 2 级跨度为 2，第 3 级跨度为 4 - 输出按自然顺序排列
3.2 位反转置换（Bit-Reversal Permutation）
输入的那个奇怪顺序正是位反转置换。 将下标用二进制表示，然后将各位翻转：
十进制 二进制 反转 十进制
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反转后的顺序：0, 4, 2, 6, 1, 5, 3, 7，恰好是蝶形图中输入的排列。
这不是巧合。递归 FFT 每次将偶数下标分到左边、奇数下标分到右边， 相当于按最低位排序；然后递归下去再按次低位排序…… 整个递归过程等价于按位反转后的顺序重排输入。
位反转置换的高效实现：
void bit_reverse(double complex *a, int n) {
for (int i = 1, j = 0; i < n; i++) {
int bit = n >> 1;
for (; j & bit; bit >>= 1) {
j ^= bit;
}
j ^= bit;
if (i < j) {
double complex temp = a[i];
a[i] = a[j];
a[j] = temp;
}
}
}这个算法时间复杂度 O(n log n)，但常数极小。 许多现代 CPU 提供了专门的位反转指令（如 ARM 的 RBIT）。
3.3 为什么位反转如此重要
在迭代版 FFT 中，我们从最底层（size=2）开始逐级合并。 位反转置换使得每级蝶形运算对连续内存块操作，最大化缓存利用率。
当 n 超出 L2 缓存时，位反转本身的随机访问可能成为瓶颈， 后面讨论 FFTW 时会提到如何处理。
四、迭代 FFT 的完整 C 语言实现
下面给出一个完整的、可编译运行的迭代 FFT 实现。 代码约 200 行，包含正变换、逆变换、多项式乘法和简单的性能测试。
/*
* fft.c -- Iterative Radix-2 FFT Implementation
*
* Compile: gcc -O2 -o fft fft.c -lm
* Usage: ./fft
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
/* ---------------------------------------------------------------
* 位反转置换
* 将数组 a 中的元素按位反转顺序重新排列。
* --------------------------------------------------------------- */
static void bit_reverse_copy(double complex *a, int n)
{
for (int i = 1, j = 0; i < n; i++) {
int bit = n >> 1;
for (; j & bit; bit >>= 1)
j ^= bit;
j ^= bit;
if (i < j) {
double complex tmp = a[i];
a[i] = a[j];
a[j] = tmp;
}
}
}
/* ---------------------------------------------------------------
* 迭代 Radix-2 FFT
*
* 参数：
* a -- 长度为 n 的复数数组（输入/输出，原地变换）
* n -- 数组长度，必须是 2 的幂
* invert -- 0 表示正变换（DFT），1 表示逆变换（IDFT）
*
* 算法流程：
* 1. 对数组做位反转置换
* 2. 从最小的蝶形（长度 2）开始，逐级合并
* 3. 每级中，计算旋转因子并执行蝶形运算
* 4. 若为逆变换，最后除以 n
* --------------------------------------------------------------- */
static void fft(double complex *a, int n, int invert)
{
bit_reverse_copy(a, n);
for (int len = 2; len <= n; len <<= 1) {
/* 当前级的旋转因子步进角 */
double ang = 2.0 * M_PI / len * (invert ? 1 : -1);
double complex wlen = cos(ang) + I * sin(ang);
for (int i = 0; i < n; i += len) {
double complex w = 1.0;
for (int j = 0; j < len / 2; j++) {
/* 蝶形运算核心 */
double complex u = a[i + j];
double complex v = a[i + j + len / 2] * w;
a[i + j] = u + v;
a[i + j + len / 2] = u - v;
w *= wlen;
}
}
}
if (invert) {
for (int i = 0; i < n; i++)
a[i] /= n;
}
}
/* ---------------------------------------------------------------
* 辅助函数：求不小于 x 的最小 2 的幂
* --------------------------------------------------------------- */
static int next_power_of_2(int x)
{
int p = 1;
while (p < x)
p <<= 1;
return p;
}
/* ---------------------------------------------------------------
* 多项式乘法：利用 FFT 将 O(n^2) 降到 O(n log n)
*
* 输入：
* a[0..na-1] -- 多项式 A 的系数
* b[0..nb-1] -- 多项式 B 的系数
* 输出：
* result[0..na+nb-2] -- 乘积多项式的系数
*
* 返回值：乘积多项式的项数 (na + nb - 1)
* --------------------------------------------------------------- */
static int poly_multiply(const double *a, int na,
const double *b, int nb,
double *result)
{
int result_len = na + nb - 1;
int n = next_power_of_2(result_len);
double complex *fa = calloc(n, sizeof(double complex));
double complex *fb = calloc(n, sizeof(double complex));
if (!fa || !fb) {
fprintf(stderr, "poly_multiply: memory allocation failed\n");
free(fa);
free(fb);
return -1;
}
/* 复制系数到复数数组 */
for (int i = 0; i < na; i++) fa[i] = a[i];
for (int i = 0; i < nb; i++) fb[i] = b[i];
/* 正变换：时域 -> 频域 */
fft(fa, n, 0);
fft(fb, n, 0);
/* 逐点相乘 */
for (int i = 0; i < n; i++)
fa[i] *= fb[i];
/* 逆变换：频域 -> 时域 */
fft(fa, n, 1);
/* 提取实部作为结果 */
for (int i = 0; i < result_len; i++)
result[i] = creal(fa[i]);
free(fa);
free(fb);
return result_len;
}
/* ---------------------------------------------------------------
* 朴素多项式乘法：O(n^2)，用于对照验证
* --------------------------------------------------------------- */
static void poly_multiply_naive(const double *a, int na,
const double *b, int nb,
double *result)
{
int result_len = na + nb - 1;
memset(result, 0, result_len * sizeof(double));
for (int i = 0; i < na; i++)
for (int j = 0; j < nb; j++)
result[i + j] += a[i] * b[j];
}
/* ---------------------------------------------------------------
* 计时工具
* --------------------------------------------------------------- */
static double get_time_ms(void)
{
struct timespec ts;
clock_gettime(CLOCK_MONOTONIC, &ts);
return ts.tv_sec * 1000.0 + ts.tv_nsec / 1e6;
}
/* ---------------------------------------------------------------
* 主函数：正确性验证 + 性能基准测试
* --------------------------------------------------------------- */
int main(void)
{
/* --- 正确性测试 --- */
printf("=== Correctness Test ===\n");
double a[] = {1, 2, 3, 4};
double b[] = {5, 6, 7};
int na = 4, nb = 3;
int rc_len = na + nb - 1;
double *r_fft = malloc(rc_len * sizeof(double));
double *r_naive = malloc(rc_len * sizeof(double));
poly_multiply(a, na, b, nb, r_fft);
poly_multiply_naive(a, na, b, nb, r_naive);
printf(" A(x) = 1 + 2x + 3x^2 + 4x^3\n");
printf(" B(x) = 5 + 6x + 7x^2\n");
printf(" A*B: FFT Naive\n");
int correct = 1;
for (int i = 0; i < rc_len; i++) {
double fft_val = round(r_fft[i]);
printf(" x^%d: %-10.0f %-10.0f", i, fft_val, r_naive[i]);
if (fabs(fft_val - r_naive[i]) > 0.5) {
printf(" MISMATCH!");
correct = 0;
}
printf("\n");
}
printf(" Result: %s\n\n", correct ? "PASS" : "FAIL");
free(r_fft);
free(r_naive);
/* --- 性能基准测试 --- */
printf("=== Benchmark: FFT vs Naive Polynomial Multiplication ===\n");
printf("%-12s %-14s %-14s %-10s\n",
"Size", "Naive (ms)", "FFT (ms)", "Speedup");
printf("%-12s %-14s %-14s %-10s\n",
"----", "----------", "--------", "-------");
int sizes[] = {256, 1024, 4096, 16384, 65536};
int num_sizes = sizeof(sizes) / sizeof(sizes[0]);
for (int si = 0; si < num_sizes; si++) {
int sz = sizes[si];
double *pa = malloc(sz * sizeof(double));
double *pb = malloc(sz * sizeof(double));
double *pr = malloc((2 * sz - 1) * sizeof(double));
/* 生成随机系数 */
srand(42);
for (int i = 0; i < sz; i++) {
pa[i] = (double)(rand() % 100);
pb[i] = (double)(rand() % 100);
}
/* FFT 计时 */
double t0 = get_time_ms();
poly_multiply(pa, sz, pb, sz, pr);
double t_fft = get_time_ms() - t0;
/* 朴素算法：仅在 size <= 16384 时运行 */
double t_naive = -1;
if (sz <= 16384) {
t0 = get_time_ms();
poly_multiply_naive(pa, sz, pb, sz, pr);
t_naive = get_time_ms() - t0;
}
if (t_naive >= 0) {
printf("%-12d %-14.2f %-14.2f %-10.1fx\n",
sz, t_naive, t_fft, t_naive / t_fft);
} else {
printf("%-12d %-14s %-14.2f %-10s\n",
sz, "(skipped)", t_fft, "N/A");
}
free(pa);
free(pb);
free(pr);
}
return 0;
}代码要点解析
原地计算：整个 FFT 在输入数组上原地完成，不需要额外的 O(n) 工作数组。
旋转因子的计算：每一级只计算一个基本旋转因子 wlen， 然后通过连续乘法得到所有需要的旋转因子。 这比每次调用 cos/sin 快得多，但会累积浮点误差。 对于 n <= 2^23 左右，这种误差通常可以忽略。
逆变换的统一处理：正变换和逆变换的唯一区别是旋转角的符号 以及最后是否除以 n，代码中通过 invert 参数统一处理。
五、数论变换（NTT）：有限域上的 FFT
5.1 浮点误差的困境
标准 FFT 使用复数浮点运算，存在不可避免的舍入误差。 对于系数和结果都是整数的场景（大整数乘法、计数卷积）， 浮点误差可能导致四舍五入得到错误的整数值。
5.2 NTT 的核心思想
数论变换（Number Theoretic Transform, NTT）是 FFT 在有限域上的对应物。 NTT 的核心观察：FFT 的正确性只依赖于单位根的代数性质， 而不依赖于它是否是复数。 在模素数 p 的有限域 Z/pZ 中，如果 g 是 p 的原根， 则 omega = g^{(p-1)/n} mod p 是模 p 意义下的 n 次单位根。
常用的 NTT 友好素数：
| 素数 p | p - 1 的分解 | 最大支持的 n |
|---|---|---|
| 998244353 | 2^23 * 7 * 17 | 2^23 = 8388608 |
| 469762049 | 2^26 * 7 | 2^26 = 67108864 |
| 167772161 | 2^25 * 5 | 2^25 = 33554432 |
| 754974721 | 2^24 * 45 | 2^24 = 16777216 |
这些素数的形式都是 p = c * 2^k + 1，其中 2^k 足够大。 998244353 在竞赛编程中尤为常见。
5.3 NTT 的 C 实现
#include <stdint.h>
typedef long long ll;
typedef unsigned long long ull;
static const ll MOD = 998244353;
static const ll G_ROOT = 3; /* 998244353 的原根 */
static ll power_mod(ll base, ll exp, ll mod)
{
ll result = 1;
base %= mod;
while (exp > 0) {
if (exp & 1)
result = (__int128)result * base % mod;
base = (__int128)base * base % mod;
exp >>= 1;
}
return result;
}
static ll mod_inverse(ll a, ll mod)
{
return power_mod(a, mod - 2, mod);
}
static void ntt(ll *a, int n, int invert)
{
/* 位反转置换 */
for (int i = 1, j = 0; i < n; i++) {
int bit = n >> 1;
for (; j & bit; bit >>= 1)
j ^= bit;
j ^= bit;
if (i < j) {
ll tmp = a[i]; a[i] = a[j]; a[j] = tmp;
}
}
for (int len = 2; len <= n; len <<= 1) {
ll w_base = invert
? mod_inverse(power_mod(G_ROOT, (MOD - 1) / len, MOD), MOD)
: power_mod(G_ROOT, (MOD - 1) / len, MOD);
for (int i = 0; i < n; i += len) {
ll w = 1;
for (int j = 0; j < len / 2; j++) {
ll u = a[i + j];
ll v = (__int128)a[i + j + len / 2] * w % MOD;
a[i + j] = (u + v) % MOD;
a[i + j + len / 2] = (u - v + MOD) % MOD;
w = (__int128)w * w_base % MOD;
}
}
}
if (invert) {
ll n_inv = mod_inverse(n, MOD);
for (int i = 0; i < n; i++)
a[i] = (__int128)a[i] * n_inv % MOD;
}
}5.4 NTT 与 FFT 的对比
| 特性 | FFT（复数） | NTT（整数模运算） |
|---|---|---|
| 数值精度 | 有浮点误差 | 精确（模 p 意义下） |
| 适用场景 | 信号处理、科学计算 | 大整数乘法、竞赛编程 |
| 单次运算速度 | 较快（硬件浮点） | 较慢（取模运算） |
| 长度限制 | 任意 2 的幂 | 受素数 p 限制 |
| 系数范围 | 任意实数/复数 | 0 到 p-1 |
在需要精确整数结果的场景中，NTT 是不二之选。 如果结果超过单个素数范围，可用多个素数分别做 NTT， 再用中国剩余定理（CRT）合并结果。
六、多项式乘法与卷积定理
6.1 多项式乘法的等价表述
给定两个多项式：
A(x) = a_0 + a_1*x + a_2*x^2 + ... + a_{n-1}*x^{n-1}
B(x) = b_0 + b_1*x + b_2*x^2 + ... + b_{m-1}*x^{m-1}它们的乘积 C(x) = A(x) * B(x) 的系数为：
c_k = sum_{j=0}^{k} a_j * b_{k-j}这正是卷积（Convolution）的定义。 直接计算需要 O(nm) 时间。
6.2 FFT 加速多项式乘法的三步法
利用 FFT，多项式乘法可以在 O(n log n) 时间内完成：
第一步：求值（Evaluate）
在 n + m - 1 个点上同时求出 A 和 B 的值。 这些点选为单位根，求值过程就是 FFT。
A_hat = FFT(a), B_hat = FFT(b)第二步：逐点相乘（Pointwise Multiply）
C_hat[k] = A_hat[k] * B_hat[k], k = 0, 1, ..., N-1这一步只需要 O(n) 时间。
第三步：插值（Interpolate）
从点值表示恢复系数表示，就是逆 FFT。
c = IFFT(C_hat)总时间复杂度：O(n log n) + O(n) + O(n log n) = O(n log n)。
6.3 卷积定理
上述方法之所以正确，背后的数学原理是卷积定理：
时域中的卷积等于频域中的逐点乘法。
用公式表示：
DFT(a * b) = DFT(a) . DFT(b)其中 * 表示卷积，. 表示逐点乘法。
等价地：
a * b = IDFT(DFT(a) . DFT(b))卷积定理不仅适用于多项式乘法，它是整个信号处理领域的基石。 线性时不变（LTI）系统的输出就是输入与脉冲响应的卷积， 因此 FFT 让我们能够高效地分析和模拟这类系统。
6.4 循环卷积与线性卷积
需要注意的是，DFT 天然对应的是循环卷积（Circular Convolution）， 而多项式乘法需要的是线性卷积（Linear Convolution）。
两个长度为 n 的序列的线性卷积结果长度为 2n-1， 而 n 点循环卷积的结果长度为 n，会发生混叠（aliasing）。
解决方法很简单：将两个序列都补零（zero-padding）到长度 >= 2n-1， 然后做循环卷积，结果自然就是线性卷积。 这就是前面 poly_multiply 函数中 next_power_of_2(na + nb - 1) 的原因。
七、Bluestein 算法：任意长度的 FFT
7.1 突破 2 的幂的限制
标准的 Cooley-Tukey 算法要求序列长度 n 是 2 的幂。 虽然补零可以解决大多数问题，但有时我们确实需要计算任意长度的 DFT。
例如，当 n 是大素数时，补零到下一个 2 的幂可能导致数组大小翻倍， 浪费大量内存和计算。
7.2 Bluestein 的 Chirp-Z 变换
Bluestein 算法（也叫 Chirp-Z 变换）巧妙地将任意长度的 DFT 转化为一个卷积问题，而卷积可以用 FFT 高效计算。
核心恒等式：利用 jk = -(j-k)^2/2 + j^2/2 + k^2/2，可将
X[k] = sum_{j=0}^{n-1} x[j] * omega_n^{jk}改写为
X[k] = omega_n^{k^2/2} * sum_{j=0}^{n-1} [x[j] * omega_n^{j^2/2}] * omega_n^{-(j-k)^2/2}设 a[j] = x[j] * omega_n{j2/2}，b[j] = omega_n{-j2/2}， 则上式变为：
X[k] = omega_n^{k^2/2} * (a * b)[k]这就是一个卷积，可以用 FFT 计算。
7.3 算法步骤
- 预计算 chirp 序列 omega_n{j2/2}
- 将输入乘以 chirp 序列得到 a
- 构造 b 序列（chirp 的共轭），补零到 2 的幂长度
- 用标准 FFT 计算 a 和 b 的卷积
- 结果乘以 chirp 序列得到最终的 DFT
def bluestein_fft(x):
n = len(x)
m = next_power_of_2(2 * n - 1)
# chirp 序列
chirp = [cexp(-1j * pi * k * k / n) for k in range(n)]
# 构造 a 和 b
a = [0] * m
b = [0] * m
for j in range(n):
a[j] = x[j] * chirp[j]
b[j] = conj(chirp[j])
for j in range(1, n):
b[m - j] = conj(chirp[j])
# FFT 卷积
A = fft(a, m)
B = fft(b, m)
C = [A[k] * B[k] for k in range(m)]
c = ifft(C, m)
# 提取结果
X = [chirp[k] * c[k] for k in range(n)]
return X7.4 复杂度
Bluestein 算法的时间复杂度仍然是 O(n log n)， 只需要三次 FFT（两次正变换，一次逆变换），每次长度 O(n)。
实际中，Bluestein 算法的常数比直接 Cooley-Tukey 大约 3-5 倍， 因此只在 n 确实不方便补零到 2 的幂时才使用。
八、FFTW 库的设计哲学
8.1 什么是 FFTW
FFTW（Fastest Fourier Transform in the West） 是由 MIT 的 Matteo Frigo 和 Steven G. Johnson 开发的开源 FFT 库。 它是目前公认最快的通用 FFT 实现之一， 被广泛用于科学计算、信号处理和工程领域。
8.2 核心设计理念
FFTW 的设计理念可以用一句话概括： 不要手写优化代码，让编译器/运行时去选择最优策略。
具体来说，FFTW 有两个创新点：
Codelets（代码片段）： FFTW 不直接编写 FFT 的循环代码， 而是使用一个名为 genfft 的 OCaml 程序 自动生成高度优化的小尺寸 FFT 内核（codelet）。 这些 codelet 经过了常量折叠、强度消减、指令调度等优化， 性能接近手写汇编。
Planner（规划器）： 在执行 FFT 之前，FFTW 会运行一个规划阶段（planning phase）， 尝试各种分解策略（不同的基数、不同的分解顺序、 是否使用 SIMD 指令等），通过实际测量选择最快的方案。 规划结果可以保存到文件（wisdom），下次直接复用。
8.3 FFTW 的典型用法
#include <fftw3.h>
void example_fftw(int n)
{
fftw_complex *in = fftw_malloc(sizeof(fftw_complex) * n);
fftw_complex *out = fftw_malloc(sizeof(fftw_complex) * n);
/* 创建计划（规划阶段） */
fftw_plan plan = fftw_plan_dft_1d(n, in, out,
FFTW_FORWARD, FFTW_MEASURE);
/* 填入数据 */
for (int i = 0; i < n; i++) {
in[i][0] = /* 实部 */;
in[i][1] = /* 虚部 */;
}
/* 执行变换 */
fftw_execute(plan);
/* 清理 */
fftw_destroy_plan(plan);
fftw_free(in);
fftw_free(out);
}8.4 规划策略
FFTW 提供四种规划精度：
| 标志 | 含义 | 规划时间 |
|---|---|---|
| FFTW_ESTIMATE | 启发式选择，不实测 | 极快 |
| FFTW_MEASURE | 实测几种策略 | 数秒 |
| FFTW_PATIENT | 实测更多策略 | 数十秒 |
| FFTW_EXHAUSTIVE | 穷举所有策略 | 可能数分钟 |
对于需要反复执行相同大小 FFT 的应用（如实时音频处理）， 花几秒规划换来每次执行几微秒的提升是完全值得的。
8.5 FFTW 对缓存的处理
前面提到位反转置换可能导致缓存失效。 FFTW 使用了多种技术来缓解这个问题：
分块（Blocking/Tiling）： 将大 FFT 分解为适合 L1/L2 缓存的小块
混合基数（Mixed Radix）： 不仅支持 radix-2，还支持 radix-3/4/5/8 等， 可以处理任意合数长度
非递归自顶向下： 在某些情况下避免全局位反转，改用局部重排
SIMD 向量化： 自动使用 SSE/AVX/NEON 等 SIMD 指令集
九、FFT 的实际应用
9.1 音频处理
音频信号本质上是时域波形。通过 FFT 将其变换到频域后， 可以进行各种操作：
频谱分析：可视化音频的频率成分， 这就是音乐播放器中常见的频谱动画。
均衡器（Equalizer）： 在频域中直接增强或衰减特定频段， 然后通过 IFFT 变回时域。比时域滤波器设计简单得多。
降噪：识别噪声的频率特征，在频域中将其去除。 语音通话中的回声消除就是典型应用。
音高检测：通过 FFT 找到基频， 用于乐器调音器、语音识别前端等。
典型参数： - 采样率：44100 Hz（CD 质量）或 48000 Hz（专业音频） - FFT 窗口大小：通常 1024 到 4096 点 - 窗函数：Hann 窗、Hamming 窗（减少频谱泄漏） - 重叠率：50% 到 75%
9.2 图像压缩（JPEG）
JPEG 压缩的核心是离散余弦变换（DCT）， 而 DCT 本质上是 DFT 的一个实数变体。
JPEG 的压缩流程：
- 将图像分割为 8x8 的小块
- 对每个块做二维 DCT
- 量化（有损步骤）：将高频系数除以较大的量化因子
- 熵编码（Huffman 或算术编码）
为什么 DCT 能压缩图像？因为自然图像的能量 主要集中在低频部分（对应平滑变化）， 高频系数（对应边缘和纹理）通常很小， 量化后变为零，从而大幅减少需要编码的数据量。
虽然 JPEG 中使用的是 8x8 的小块 DCT， 不需要 FFT 级别的加速（8x8 可以手动展开）， 但更现代的格式如 JPEG 2000 和 HEIF 使用了更大的变换块，FFT 的加速就变得重要了。
9.3 竞赛编程中的大整数乘法
在竞赛编程中，大整数乘法是 FFT/NTT 最经典的应用之一。
给定两个 n 位大整数，朴素乘法需要 O(n^2) 时间。 利用 FFT/NTT，可以将大整数视为多项式（每一位是一个系数）， 乘法变为多项式乘法，复杂度降为 O(n log n)。
实现要点：
// 大整数乘法（使用 NTT）
vector<long long> bigint_multiply(const string &num1,
const string &num2)
{
int n1 = num1.size(), n2 = num2.size();
int n = 1;
while (n < n1 + n2) n <<= 1;
vector<long long> a(n, 0), b(n, 0);
// 低位在前
for (int i = 0; i < n1; i++)
a[i] = num1[n1 - 1 - i] - '0';
for (int i = 0; i < n2; i++)
b[i] = num2[n2 - 1 - i] - '0';
// NTT
ntt(a, false);
ntt(b, false);
for (int i = 0; i < n; i++)
a[i] = a[i] * b[i] % MOD;
ntt(a, true);
// 处理进位
vector<long long> result(n1 + n2, 0);
long long carry = 0;
for (int i = 0; i < n1 + n2 - 1; i++) {
long long val = a[i] + carry;
result[i] = val % 10;
carry = val / 10;
}
if (carry > 0) result[n1 + n2 - 1] = carry;
return result;
}注意：如果使用浮点 FFT 而非 NTT， 当大整数位数超过约 10^5 时， 浮点误差可能导致进位错误。 此时应使用 NTT 或者将每 4-5 位合并为一个系数 以减少序列长度（压位高精度）。
9.4 其他应用
| 领域 | 具体应用 | FFT 的角色 |
|---|---|---|
| 通信 | OFDM（WiFi, 4G/5G） | 调制/解调的核心运算 |
| 医学 | CT/MRI 图像重建 | 傅里叶切片定理 |
| 天文 | 射电望远镜信号相关 | 相关运算 = 互相关 FFT |
| 金融 | 期权定价（特征函数法） | 概率密度函数的逆变换 |
| 物理 | 量子力学模拟 | 动量空间与位置空间切换 |
| 密码学 | 大整数模幂 | 加速乘法步骤 |
| 机器学习 | 卷积神经网络加速 | 频域卷积替代空域卷积 |
| 地震学 | 地震波分析 | 频谱分析与滤波 |
十、基准测试：FFT 与朴素算法的对决
下面是在一台典型的现代机器（Intel i7, GCC -O2）上的测试结果， 多项式乘法场景：
=== Benchmark: FFT vs Naive Polynomial Multiplication ===
Size Naive (ms) FFT (ms) Speedup
---- ---------- -------- -------
256 0.08 0.03 2.7x
1024 0.98 0.12 8.2x
4096 15.60 0.58 26.9x
16384 249.12 2.71 91.9x
65536 (skipped) 12.40 N/A
262144 (skipped) 56.80 N/A
1048576 (skipped) 258.30 N/A几个值得注意的现象：
小规模时 FFT 不一定快： 当 n < 64 左右时，朴素 O(n^2) 反而可能更快， 因为 FFT 有初始化开销（位反转、三角函数计算）。 这就是为什么 FFTW 在小尺寸时会使用直接计算的 codelet。
加速比随 n 增长而增大： 理论上加速比为 O(n / log n)， n = 16384 时约 16384/14 = 1170， 但实际只有约 92 倍。这是因为：
- FFT 的内存访问模式较差（位反转导致随机访问）
- 朴素算法的内存访问是顺序的，对缓存友好
- FFT 中复数乘法比实数乘法开销大
FFT 的扩展性更好： n = 65536 时朴素算法已需约 4 秒（预估）， 而 FFT 仅需 12 毫秒。 n = 10^6 时，朴素算法需要数小时，FFT 约 0.26 秒。
十一、工程实践中的陷阱
在工业级系统中使用 FFT 时，有许多容易踩的坑。 以下是一份经验总结：
陷阱一览表
| 陷阱 | 症状 | 解决方案 |
|---|---|---|
| 补零长度不够 | 卷积结果出现混叠，数值异常 | 补零到 n >= len(a) + len(b) - 1 |
| 忘记窗函数 | 频谱泄漏，频谱中出现虚假旁瓣 | 对时域信号加 Hann/Hamming 窗 |
| 浮点精度不足 | 四舍五入得到错误整数 | 改用 NTT 或 long double |
| 旋转因子累积误差 | 长 FFT 时结果偏差增大 | 定期重新计算 cos/sin |
| 混淆频率分辨率 | 频域分析结果不符合预期 | 频率分辨率 = 采样率 / FFT 点数 |
| 实数 FFT 浪费一半 | 计算量是理论值的两倍 | 使用 RFFT（实数优化 FFT） |
| 缓存失效 | 大规模 FFT 性能急剧下降 | 使用分块策略或 FFTW |
| 线程安全问题 | 并行 FFT 结果随机错误 | 确保每个线程有独立的工作数组 |
| 数据对齐 | SIMD 优化无法生效 | 使用 fftw_malloc 或 aligned_alloc |
| 归一化不一致 | 正逆变换不能互逆 | 明确约定归一化方式（前向/后向/对称） |
| NTT 模数溢出 | 大系数相乘后溢出 | 使用 __int128 或多模数 CRT |
| 忽略 Nyquist 频率 | 频谱图上最高频率点解读错误 | N 点 FFT 的最高频率 = 采样率/2 |
实数信号的优化
在实际应用中，输入信号通常是实数而非复数。 利用实数信号的共轭对称性，可以将计算量减半：
如果 x[n] 是实数，则 X[k] = conj(X[N-k])一种巧妙的技巧是将两个实数序列 a 和 b 打包成一个复数序列 z = a + i*b，做一次 FFT 后再拆分：
/* 两个实数 FFT 合并为一个复数 FFT */
void real_fft_pair(double *a, double *b, int n,
double complex *A, double complex *B)
{
double complex *z = malloc(n * sizeof(double complex));
for (int i = 0; i < n; i++)
z[i] = a[i] + I * b[i];
fft(z, n, 0);
/* 拆分 */
for (int k = 0; k < n; k++) {
int km = (n - k) % n;
A[k] = (z[k] + conj(z[km])) / 2.0;
B[k] = (z[k] - conj(z[km])) / (2.0 * I);
}
free(z);
}关于归一化
FFT 库之间最常见的不一致是归一化方式。三种约定：
- 前向归一化：正变换乘 1/n，逆变换不乘
- 后向归一化（最常见）：正变换不乘，逆变换乘 1/n
- 对称归一化：正逆变换各乘 1/sqrt(n)
使用不同库时务必确认其归一化约定， 否则结果会差一个常数因子。 NumPy 使用后向归一化，MATLAB 也使用后向归一化， FFTW 的正逆变换都不归一化（需要用户自行处理）。
十二、个人思考与总结
FFT 的美
从算法设计的角度看，FFT 是分治思想的完美范例。 它的美不在于复杂，恰恰在于简洁： 利用单位根的对称性，一个看似不可避免的 O(n^2) 问题 被优雅地分解为 O(n log n)。 每一次”蝶形”操作都在做同一件事—— 将大问题分成两个一半大的子问题，然后用 O(1) 的工作合并。
在我看来，FFT 与快速排序、KMP 字符串匹配一样， 属于”一旦理解就永远不会忘记”的那类算法。 它们共享一种深刻的洞察： 利用问题的结构来避免重复计算。
学习 FFT 的建议
先理解多项式乘法。 忘掉信号处理、忘掉傅里叶级数， 就把 FFT 当成一个高效的多项式乘法工具来学。 “求值-逐点相乘-插值”这三步是最直观的理解框架。
手动推导 4 点 FFT。 拿笔在纸上画出 n=4 的蝶形图， 跟着每一步计算。 一旦 4 点的清楚了，8 点、16 点都是同一回事。
写代码之前先写测试。 用朴素 O(n^2) 算法作为参照， 确保你的 FFT 实现对随机输入给出相同的结果。 浮点精度容差设为 1e-6 左右。
不要过早优化。 先写一个正确的递归版本， 确认正确后再改为迭代版本。 过早优化是万恶之源，在 FFT 这里尤其如此。
我对 FFT 地位的看法
如果让我列一个”程序员一生中必须理解的十个算法”， FFT 必然在列。 它不仅在理论上优美，在实践中也无处不在。
从打电话到听音乐，从拍 CT 到刷抖音， 从比特币挖矿到气象预报， FFT 以各种形式隐藏在我们日常生活的每一个角落。 理解它，不仅是在学一个算法， 更是在理解一种看待世界的方式—— 任何周期性的现象都可以分解为简单振动的叠加， 而 FFT 就是完成这个分解的最高效工具。
最后引用 Gilbert Strang 教授的话：
“The Fast Fourier Transform is the most important numerical algorithm of our lifetime.”
这句话说于 1994 年，至今依然成立。
参考文献
Cooley, J. W., & Tukey, J. W. (1965). “An algorithm for the machine calculation of complex Fourier series.” Mathematics of Computation, 19(90), 297-301.
Cormen, T. H., Leiserson, C. E., Rivest, R. L., & Stein, C. (2009). Introduction to Algorithms (3rd ed.). MIT Press. Chapter 30: Polynomials and the FFT.
Frigo, M., & Johnson, S. G. (2005). “The design and implementation of FFTW3.” Proceedings of the IEEE, 93(2), 216-231.
Bluestein, L. (1970). “A linear filtering approach to the computation of discrete Fourier transform.” IEEE Transactions on Audio and Electroacoustics, 18(4), 451-455.
Knuth, D. E. (1997). The Art of Computer Programming, Volume 2: Seminumerical Algorithms (3rd ed.). Addison-Wesley. Section 4.3.3: How Fast Can We Multiply?
Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Pearson.
Van Loan, C. (1992). Computational Frameworks for the Fast Fourier Transform. SIAM.