Lesson 67: 快速傅里叶变换 N=8
练习任务
难度：中
用 C 语言实现 Cooley-Tukey 快速傅里叶变换 (FFT) 算法，对 N=8 点的实数序列 {1, 2, 3, 4, 4, 3, 2, 1} 进行频谱分析。你需要完成 5 个核心函数：
bit_reverse()— 位反转置换：将输入数组按二进制位反转重排butterfly()— 单级蝶形运算：执行一层 Cooley-Tukey 蝶形操作fft()— 主 FFT 流程：调用位反转 + 3 级蝶形，打印中间结果print_complex()— 复数打印：按real+imagi格式输出复数main()— 主流程：构建复数数组，驱动 FFT，打印最终频谱
Complex 结构体、固定输入序列、旋转因子表 W[4] 和常量 N=8 已提供，学员只需实现算法逻辑。验证方式：make test，编译后输出与 expected_output.txt 逐行对比。
提示：蝶形运算的关键陷阱——必须先保存
even的原始值再计算。如果先修改了a[even_idx]，lower = a[even_idx] - T就会用错误的被减数。为什么位反转要从i < j时才交换？去掉这个条件会发生什么？
核心知识点
- DFT 定义与 FFT 动机 — 离散傅里叶变换 O(N²) 不可行，分治策略降至 O(N log N)，加速比可达数万倍
- Cooley-Tukey 分治策略 — 按索引奇偶将 N 点 DFT 分解为两个 N/2 点 DFT，递归至 N=1，再通过蝶形组合
- 蝶形运算 —
upper = even + W·odd, lower = even - W·odd，必须先保存 even 原始值 - 旋转因子 (Twiddle Factor) —
W_N^k = e^{-2πi·k/N} = cos(2πk/N) - i·sin(2πk/N)，复平面单位圆上的点 - 位反转置换 — 高效位运算算法在线生成位反转序，使蝶形对相邻排列，支持原地迭代计算
- N=8 三级蝶形网络 — 逐级追踪 group_size=2/4/8，理解索引规律和 twiddle_idx 计算
- 复数运算基础 — 乘法公式
(a+bi)(c+di) = (ac-bd) + (ad+bc)i，蝶形中 T 的计算 - Hermitian 对称性 — 实输入信信号的共轭对称
X[N-k] = conj(X[k])，节省一半计算量 - DIT vs DIF — 时域抽取(DIT)位反转输入→正常序输出；频域抽取(DIF)正常序输入→位反转输出
- 原地计算 — 蝶形直接在原数组更新，O(N) 空间，无需额外数组
代码框架
#include <stdio.h>
#include <math.h>
#define N 8
typedef struct {
double real;
double imag;
} Complex;
/* 旋转因子 W_N^k = e^(-2πi·k/N), N=8, 只需前 4 个 */
const Complex W[4] = {
{ 1.000000, 0.000000}, // W8^0 = 1
{ 0.707107, -0.707107}, // W8^1 = √2/2 - i√2/2
{ 0.000000, -1.000000}, // W8^2 = -i
{-0.707107, -0.707107}, // W8^3 = -√2/2 - i√2/2
};
/* 位反转置换：使用高效位运算算法 */
void bit_reverse(Complex a[], int n) {
// ① j = 0
// ② for (i = 1; i < n; i++):
// bit = n / 2
// while (j & bit): j ^= bit; bit >>= 1
// j ^= bit
// if (i < j): 交换 a[i] 和 a[j]
}
/* 单级蝶形运算: stage 控制 group_size = 1 << stage */
void butterfly(Complex a[], int n, int stage) {
// ③ 计算 group_size = 1 << stage, half = group_size >> 1, step = n / group_size
// ④ 外层: for (g = 0; g < n; g += group_size) 遍历每组
// ⑤ 内层: for (k = 0; k < half; k++) 遍历组内蝶形对:
// even_idx = g + k, odd_idx = g + k + half
// twiddle_idx = k * step
// 保存 even_saved = a[even_idx]
// T = W[twiddle_idx] × a[odd_idx]
// 更新 a[even_idx] = even_saved + T
// 更新 a[odd_idx] = even_saved - T
}
/* 复数打印：一行实现，%+ 标志处理正负号 */
void print_complex(Complex c) {
// ⑥ 使用 printf 的 %+ 标志实现格式输出
}
/* 主 FFT 流程 */
void fft(Complex a[], int n, const Complex W[]) {
// ⑦ bit_reverse(a, n)
// ⑧ 打印 "Bit-reversed order:" 及每个元素
// ⑨ for (stage = 1; stage <= 3; stage++):
// butterfly(a, n, stage)
// 打印 "Stage %d (group_size=%d):" 及每个元素
}
int main(void) {
// ⑩ 声明 Complex a[N]
// ⑪ 填充: a[i].real = input[i]; a[i].imag = 0.0
// input[] = {1, 2, 3, 4, 4, 3, 2, 1}
// ⑫ 打印输入标题和序列
// ⑬ 调用 fft(a, N, W)
// ⑭ 打印 "Final frequency-domain result:" 及 X[0]~X[7]
return 0;
}阅读骨架后，尝试自己填充 // ① 到 // ⑭ 标记的部分。核心挑战在于：位反转的位运算逐位翻转逻辑如何工作？蝶形中 twiddle_idx = k * step 是如何推导出来的？复数乘法 (ac-bd) + (ad+bc)i 为什么实部是减号？print_complex 的 %+ 标志有什么作用？
TIP
先不要往下翻看参考解答。用纸笔完整追踪 N=8 三级蝶形网络：画 8 个节点的 3 层信号流图，逐级标注旋转因子和中间结果。重点理解 Stage 1(group_size=2)时每组只有一个蝶形(4 组)→Stage 2(group_size=4)时每组 2 个蝶形(2 组)→Stage 3(group_size=8)时整组 4 个蝶形(1 组)的递进关系。
深度讲解
1. 从 DFT 到 FFT——为什么要快？
1.1 DFT 定义与计算量
离散傅里叶变换将时域 N 个采样点变换为频域 N 个频谱分量：
N-1
X[k] = Σ x[n] · e^(-2πi·k·n/N), k = 0, 1, ..., N-1
n=0直接计算：每个 k 需要 N 次复数乘法 + N-1 次复数加法，N 个 k 共需 O(N²) 次复数乘法。当 N=1024 时约 100 万次，N=10⁶ 时达 10¹² 次——这在实际中不可接受。
1.2 Cooley-Tukey 的分治洞察
1965 年，Cooley 和 Tukey 发表了标志性论文，揭示了 DFT 的递归结构：
将输入序列 x[n] 按索引奇偶分成两组:
偶数索引: x[0], x[2], x[4], x[6] → 两个 N/2 点 DFT → E[k]
奇数索引: x[1], x[3], x[5], x[7] → 两个 N/2 点 DFT → O[k]
组合公式 (这就是"蝶形"):
X[k] = E[k] + W_N^k · O[k]
X[k+N/2] = E[k] - W_N^k · O[k]关键洞察：如果 N/2 点 DFT 的复杂度是 T(N/2)，组合需要 N/2 次蝶形运算(每对一次复数乘法+两次加减法)，则 T(N) = 2·T(N/2) + O(N) → T(N) = O(N log N)。
| N | DFT O(N²) | FFT O(N log₂ N) | 加速比 |
|---|---|---|---|
| 8 | 64 | 24 | 2.7× |
| 256 | 65,536 | 2,048 | 32× |
| 1024 | 1,048,576 | 10,240 | 102× |
| 4096 | 16,777,216 | 53,248 | 315× |
| 10⁶ | 10¹² | 2×10⁷ | 50,000× |
N=8 的分治树结构：
N=8 DFT
/ \
N=4 DFT N=4 DFT
/ \ / \
N=2 DFT N=2 DFT N=2 DFT N=2 DFT
/ \
N=1 DFT N=1 DFT ... (共 8 个叶节点)IMPORTANT
FFT 被称为"20 世纪最重要的数值算法"绝非夸张。从 JPEG 图像压缩(DCT 是 FFT 近亲)、MP3 音频编码(MDCT)、4G/5G WiFi 的 OFDM 调制，到医学 MRI 成像、地震数据分析——FFT 的身影无处不在。本课 N=8 是理解它的最小完整实例。
2. 蝶形运算——FFT 的基本计算单元
2.1 蝶形的数学定义
每个蝶形处理两个复数 a（偶数索引）和 b（奇数索引），产生两个输出：
蝶形信号流图:
a (even) ──────┬────── a + W·b (上支路: 和)
│
W (旋转因子)
│
b (odd) ──────┼────── a - W·b (下支路: 差)
│
✕ (交叉)数学表达：upper = a + W·b, lower = a - W·b
2.2 蝶形的 C 实现与关键陷阱
void butterfly(Complex a[], int n, int stage) {
int group_size = 1 << stage; // 2, 4, 8
int half = group_size >> 1; // 1, 2, 4
int step = n / group_size; // 4, 2, 1
for (int g = 0; g < n; g += group_size) { // 外层: 遍历每组
for (int k = 0; k < half; k++) { // 内层: 组内蝶形对
int even_idx = g + k;
int odd_idx = g + k + half;
int twiddle_idx = k * step;
/* 关键: 必须先保存 even 的原始值! */
Complex even_saved = a[even_idx]; // 保存
/* 复数乘法: T = W[twiddle_idx] × a[odd_idx] */
Complex T;
T.real = W[twiddle_idx].real * a[odd_idx].real
- W[twiddle_idx].imag * a[odd_idx].imag;
T.imag = W[twiddle_idx].real * a[odd_idx].imag
+ W[twiddle_idx].imag * a[odd_idx].real;
/* 蝶形更新 */
a[even_idx].real = even_saved.real + T.real;
a[even_idx].imag = even_saved.imag + T.imag;
a[odd_idx].real = even_saved.real - T.real;
a[odd_idx].imag = even_saved.imag - T.imag;
}
}
}CAUTION
最常见的 FFT Bug：不保存 even_saved。如果写成 a[even_idx].real = a[even_idx].real + T.real 然后 a[odd_idx].real = a[even_idx].real - T.real，第二行使用的 a[even_idx] 已经是更新后的值，导致 lower 计算错误。这个 Bug 在 N=8 的中间阶段可能不明显，但在大 N 下会导致完全错误的频谱。
2.3 twiddle_idx = k * step 的推导
为什么旋转因子索引不是直接用 k，而是 k * step？关键在于不同 stage 的组间间隔：
Stage 1: group_size=2, step=4
4 组, 每组 1 个蝶形(k=0)
每组蝶形对应原始序列中 step=4 个位置的旋转
→ W 为 0, 0, 0, 0 (所有蝶形用 W8^0)
Stage 2: group_size=4, step=2
2 组, 每组 2 个蝶形(k=0,1)
→ W 为 0,2, 0,2 (W8^0, W8^4 交替)
Stage 3: group_size=8, step=1
1 组, 每组 4 个蝶形(k=0,1,2,3)
→ W 为 0,1,2,3 (依次使用 W8^0,W8^1,W8^2,W8^3)twiddle_idx = k * step 精确地表达了这一规律——step 控制旋转因子在单位圆上的"行进"速度。
3. 位反转置换——为什么需要它？
3.1 Cooley-Tukey 分治的自然结果
如果直接按分治树处理输入序列，需要复杂的索引计算。位反转提供了一种优雅的替代方案：先按位反转重排输入，然后迭代蝶形运算——无需递归，且输出天然是正常序。
原始输入: [x0, x1, x2, x3, x4, x5, x6, x7]
索引二进制: 000 001 010 011 100 101 110 111
↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓
位反转: 000 100 010 110 001 101 011 111
↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓
位反转序: [x0, x4, x2, x6, x1, x5, x3, x7]3.2 高效位运算算法
不使用查表，通过位运算在线计算每个索引的位反转：
void bit_reverse(Complex a[], int n) {
int j = 0;
for (int i = 1; i < n; i++) {
int bit = n >> 1; // 从最高位开始
while (j & bit) { // 如果 j 的当前位是 1
j ^= bit; // 清除它
bit >>= 1; // 移到下一位
}
j ^= bit; // 翻转当前位
if (i < j) { // 只交换一次！
Complex tmp = a[i];
a[i] = a[j];
a[j] = tmp;
}
}
}以 N=8 为例追踪：
i=1: j=0→bit=4, j&4=0→j^=4→j=4, 1<4 → swap(a[1],a[4]) (1↔4)
i=2: j=4→bit=4, j&4=4→j^=4→j=0, bit=2→j&2=0→j^=2→j=2, 2<2 → 不交换
i=3: j=2→bit=4, j&4=0→j^=4→j=6, 3<6 → swap(a[3],a[6]) (3↔6)
i=4: j=6→bit=4, j&4=4→j^=4→j=2, bit=2→j&2=2→j^=2→j=0,
bit=1→j&1=0→j^=1→j=1, 4<1 → 不交换
i=5: j=1→bit=4, j&4=0→j^=4→j=5, 5<5 → 不交换
i=6: j=5→bit=4, j&4=4→j^=4→j=1, bit=2→j&2=0→j^=2→j=3, 6<3 → 不交换
i=7: j=3→bit=4, j&4=0→j^=4→j=7, 7<7 → 不交换结果：[0,4,2,6,1,5,3,7]——恰好是位反转序。
WARNING
if (i < j) 不是可选的优化，而是正确性保证。去掉这个条件，每对元素会被交换两次，恢复原状——位反转白做了。
4. 旋转因子——复平面单位圆上的舞蹈
4.1 定义与几何意义
旋转因子 W_N^k = e^{-2πi·k/N} = cos(2πk/N) - i·sin(2πk/N) 是复平面单位圆上的点：
Im
↑
W8^2=-i │
W8^3 │ W8^1
(-√2/2-√2/2i) │ (√2/2-√2/2i)
\ │ /
\ 135° │ 45° /
\ │ /
\ │ /
\ │ /
\ │ /
\ │ /
\ │ /
\│/
──────────────●──────────────→ Re
/│\ W8^0=1
/ │ \
/ │ \
/ │ \
/ │ \
/ │ \
/ 225°│ 315° \
/ │ \
W8^5=-W8^1 │ W8^7=-W8^3
W8^6=-W8^2=iN=8 旋转因子表：
| W8^k | 直角坐标 | 角度 |
|---|---|---|
| W8^0 | 1.000 + 0.000i | 0° |
| W8^1 | 0.707 - 0.707i | -45° |
| W8^2 | 0.000 - 1.000i | -90° |
| W8^3 | -0.707 - 0.707i | -135° |
| W8^4 | -1.000 + 0.000i | -180° (= -W8^0) |
| W8^5 | -0.707 + 0.707i | -225° (= -W8^1) |
| W8^6 | 0.000 + 1.000i | -270° (= -W8^2) |
| W8^7 | 0.707 + 0.707i | -315° (= -W8^3) |
NOTE
本题只用到 W8^0 ~ W8^3 (W[0] ~ W[3])。由于对称性 W_N^{k+N/2} = -W_N^k，其余可由取负得到。例如 W8^4 = -W8^0, W8^5 = -W8^1。
4.2 旋转因子的对称性
W_N^{k+N/2} = -W_N^k (半周对称)
W_N^{N-k} = conj(W_N^k) (共轭对称)
W_N^0 = 1 (恒等)
W_N^{N/4} = -i (90°旋转)这些对称性使得只需存储 N/4 个旋转因子，其余可通过取负/取共轭得到——这是 FFT 优化的重要基础。
5. 复数运算基础
5.1 复数乘法公式
/* 复数乘法: (a+bi)(c+di) = (ac-bd) + (ad+bc)i
*
* 推导:
* (a+bi)(c+di) = ac + adi + bci + bdi²
* = ac + (ad+bc)i + bd(-1)
* = (ac - bd) + (ad + bc)i
*/
void complex_mul(Complex *result, Complex x, Complex y) {
result->real = x.real * y.real - x.imag * y.imag; // ac - bd
result->imag = x.real * y.imag + x.imag * y.real; // ad + bc
}
/* 复数加法 */
void complex_add(Complex *result, Complex x, Complex y) {
result->real = x.real + y.real;
result->imag = x.imag + y.imag;
}
/* 复数减法 */
void complex_sub(Complex *result, Complex x, Complex y) {
result->real = x.real - y.real;
result->imag = x.imag - y.imag;
}CAUTION
复数乘法公式中实部是 ac - bd（减号），虚部是 ad + bc（加号）。许多初学者把实部也写成加号——ac + bd——这会导致蝶形运算全盘错误。
6. N=8 三级蝶形完整追踪
以输入 {1, 2, 3, 4, 4, 3, 2, 1} 的位反转后序列 {1, 4, 3, 2, 2, 3, 4, 1} 为起点：
Stage 1 (group_size=2, half=1, step=4)：4 组，每组 1 个蝶形，全部使用 W8^0 = 1。
组0: (a[0],a[1]) = (1,4): T=1×4=4, upper=5, lower=-3
组1: (a[2],a[3]) = (3,2): T=1×2=2, upper=5, lower=1
组2: (a[4],a[5]) = (2,3): T=1×3=3, upper=5, lower=-1
组3: (a[6],a[7]) = (4,1): T=1×1=1, upper=5, lower=3
→ [5, -3, 5, 1, 5, -1, 5, 3]Stage 2 (group_size=4, half=2, step=2)：2 组，每组 2 个蝶形，k=0 用 W8^0=1，k=1 用 W8^2=-i。
组0(索引0-3):
k=0: (a[0],a[2]) = (5,5): T=1×5=5, upper=10, lower=0
k=1: (a[1],a[3]) = (-3,1): T=(-i)×1=-i, upper=-3-i, lower=-3+i
组1(索引4-7):
k=0: (a[4],a[6]) = (5,5): T=1×5=5, upper=10, lower=0
k=1: (a[5],a[7]) = (-1,3): T=(-i)×3=-3i, upper=-1-3i, lower=-1+3i
→ [10, -3-i, 0, -3+i, 10, -1-3i, 0, -1+3i]Stage 3 (group_size=8, half=4, step=1)：1 组，4 个蝶形，k=0~3 对应 W8^0~W8^3。
k=0: (a[0],a[4]) = (10,10): T=1×10=10, upper=20, lower=0
k=1: (a[1],a[5]) = (-3-i, -1-3i):
T = W8^1 × (-1-3i) = (0.707-0.707i)(-1-3i)
= (-0.707-2.121) + i(0.707-2.121) ≈ -2.828 - 1.414i
upper = (-3-i) + T ≈ -5.828 - 2.414i → X[1]
lower = (-3-i) - T ≈ -0.172 - 0.414i → X[3]
k=2: (a[2],a[6]) = (0,0): T=0, upper=0, lower=0
k=3: (a[3],a[7]) = (-3+i, -1+3i):
T = W8^3 × (-1+3i) = (-0.707-0.707i)(-1+3i)
= (0.707+2.121) + i(-2.121+0.707) ≈ 2.828 - 1.414i
upper = (-3+i) + T ≈ -0.172 - 0.414i → X[2]
lower = (-3+i) - T ≈ -5.828 + 2.414i → X[7]最终频域结果：
X[0] = 20.000000 + 0.000000i (DC 分量, Σx[n] = 20)
X[1] = -5.828427 - 2.414214i (基频分量)
X[2] = 0.000000 + 0.000000i (2 次谐波 → 零!)
X[3] = -0.171573 - 0.414214i (3 次谐波)
X[4] = 0.000000 + 0.000000i (奈奎斯特频率 → 零!)
X[5] = -0.171573 + 0.414214i (X[3] 的共轭)
X[6] = 0.000000 + 0.000000i (X[2] 的共轭)
X[7] = -5.828427 + 2.414214i (X[1] 的共轭)7. 实输入的共轭对称性
对于实数输入序列，DFT 满足 Hermitian 对称性：
X[N-k] = conj(X[k]) 即 X[N-k].real = X[k].real
X[N-k].imag = -X[k].imag验证本题：X[7] = -5.828 + 2.414i，conj(X[1]) = conj(-5.828 - 2.414i) = -5.828 + 2.414i = X[7] ✓
直观理解：实数信号没有"虚部来源"，频域中的虚部必须成对出现、相互抵消——所以频谱具有对称性。这一性质的实际价值：对于 N 点实输入，只需计算 X[0] 到 X[N/2]，后半部分由对称性自动得到，节省约一半计算量。
参考解答
练习1: bit_reverse — 位反转置换
void bit_reverse(Complex a[], int n) {
int j = 0;
for (int i = 1; i < n; i++) {
int bit = n >> 1;
while (j & bit) {
j ^= bit; /* 清除当前位 */
bit >>= 1;
}
j ^= bit; /* 翻转当前位 */
if (i < j) { /* 只交换一次 */
Complex tmp = a[i];
a[i] = a[j];
a[j] = tmp;
}
}
}要点：bit 从最高位向低位移动，内层 while 处理进位链，外层 if (i < j) 确保每对元素只交换一次。
练习2-3: butterfly — 单级蝶形运算
void butterfly(Complex a[], int n, int stage, const Complex W[]) {
int group_size = 1 << stage; /* 2, 4, 8 */
int half = group_size >> 1; /* 1, 2, 4 */
int step = n / group_size; /* 4, 2, 1 */
for (int g = 0; g < n; g += group_size) {
for (int k = 0; k < half; k++) {
int even_idx = g + k;
int odd_idx = g + k + half;
int twiddle_idx = k * step;
Complex even_saved = a[even_idx]; /* 关键: 保存原始值 */
/* T = W[twiddle_idx] × a[odd_idx] */
Complex T;
T.real = W[twiddle_idx].real * a[odd_idx].real
- W[twiddle_idx].imag * a[odd_idx].imag;
T.imag = W[twiddle_idx].real * a[odd_idx].imag
+ W[twiddle_idx].imag * a[odd_idx].real;
/* 蝶形更新 */
a[even_idx].real = even_saved.real + T.real;
a[even_idx].imag = even_saved.imag + T.imag;
a[odd_idx].real = even_saved.real - T.real;
a[odd_idx].imag = even_saved.imag - T.imag;
}
}
}核心逻辑：先保存 even_saved，再计算 T = W × odd，最后用 even_saved 计算上支路(和)和下支路(差)。
练习4: print_complex — 复数打印
void print_complex(Complex c) {
printf("%10.6f%+10.6fi", c.real, c.imag);
}%+10.6f 中的 + 标志使正数前显示 + 号，负数自动显示 - 号。10 指定字段宽度以实现对齐。
练习5: fft + main — 完整程序
#include <stdio.h>
#include <math.h>
#define N 8
typedef struct {
double real;
double imag;
} Complex;
const Complex W[4] = {
{ 1.000000, 0.000000},
{ 0.707107, -0.707107},
{ 0.000000, -1.000000},
{-0.707107, -0.707107},
};
void bit_reverse(Complex a[], int n);
void butterfly(Complex a[], int n, int stage, const Complex W[]);
void print_complex(Complex c);
void fft(Complex a[], int n, const Complex W[]) {
bit_reverse(a, n);
printf("Bit-reversed order:\n");
for (int i = 0; i < n; i++) {
printf(" [%d]: ", i);
print_complex(a[i]);
printf("\n");
}
for (int stage = 1; stage <= 3; stage++) {
butterfly(a, n, stage, W);
printf("Stage %d (group_size=%d):\n", stage, 1 << stage);
for (int i = 0; i < n; i++) {
printf(" [%d]: ", i);
print_complex(a[i]);
printf("\n");
}
}
}
int main(void) {
Complex a[N];
double input[N] = {1, 2, 3, 4, 4, 3, 2, 1};
for (int i = 0; i < N; i++) {
a[i].real = input[i];
a[i].imag = 0.0;
}
printf("=== Cooley-Tukey FFT (N=%d) ===\n\n", N);
printf("Input (time domain):\n");
printf(" x = [1, 2, 3, 4, 4, 3, 2, 1]\n\n");
fft(a, N, W);
printf("\nFinal frequency-domain result:\n");
for (int i = 0; i < N; i++) {
printf(" X[%d] = ", i);
print_complex(a[i]);
printf("\n");
}
return 0;
}输出验证：X[0] = 20.00（DC 分量 = Σx[n]），X[2] 和 X[4] 应为零，X[1] 和 X[7] 应互为共轭。
对照检查：bit_reverse 中用了
if (i < j)吗？butterfly 中先保存even_saved了吗？复数乘法实部是ac - bd不是ac + bd吗？twiddle_idx = k * step不是k吗？print_complex用了%+标志吗？
课堂讨论
- 为什么 X[2] 和 X[4] 为零？这与输入信号的什么特性有关？
- 位反转置换中，为什么需要
if (i < j)条件？如果去掉会发生什么？ - 如果 N 不是 2 的幂（例如 N=6），Cooley-Tukey 算法还能直接使用吗？有哪些处理方法？
- 在蝶形运算中，旋转因子索引
twiddle_idx = k * step是如何推导出来的？为什么 Stage 1 中所有蝶形都用 W8^0？ - FFT 结果中的共轭对称性有什么实际用途？可以节省多少计算量？
- DIT(时域抽取)和 DIF(频域抽取)有什么区别？各自的输入输出顺序是怎样的？
讨论问题
Q1: 为什么 X[2] 和 X[4] 为零？
输入 {1,2,3,4,4,3,2,1} 是一个循环对称的序列：x[n] = x[N-n]（循环意义上）。这种对称性在频域中表现为：只有与对称性匹配的频率分量才有能量。X[2] 对应 2 周/8 采样 = 4 采样/周期的分量，输入序列以 8 为循环，偶数谐波恰好与信号的对称中心对齐抵消。X[4] 对应奈奎斯特频率（最高可分辨频率），在这个对称信号中也恰好为零。
更具体地说：如果把输入视为周期延拓的信号，偶数频率分量的正负半周贡献相互抵消，所以 X[2]=X[4]=0。
Q2: 位反转中 if (i < j) 的作用
去掉 if (i < j) 会导致每对元素被交换两次——第一次 i→j 交换，第二次 j→i 时再交换一次，恢复原状。例如 i=1 时 a[1] 与 a[4] 交换，i=4 时 a[4] 又与 a[1] 交换回来。最终数组没有发生任何重排。
if (i < j) 确保每对元素只在"i 小于 j"时交换一次。这是一个经典的"避免重复交换"模式，也出现在求逆序对等算法中。
Q3: N 不是 2 的幂怎么办？
Cooley-Tukey 算法要求 N 是 2 的幂。对于 N 不是 2 的幂的情况，常用处理方法：
补零 (Zero-padding)：将输入补零到最近的 2 的幂。例如 N=6 补零到 N=8。这会改变频谱分辨率，但不改变频率成分——补零在频域中是一种插值。
混合基 FFT (Mixed-radix FFT)：如果 N = p·q（p、q 互质），可以用 Cooley-Tukey 的推广形式分解。例如 N=6 = 2×3，可以用基-2 和基-3 蝶形组合。
素因子 FFT (Prime Factor Algorithm)：N 可分解为互质因子时使用。
Bluestein 算法 (Chirp-Z)：可以用 O(N log N) 时间计算任意 N 的 DFT，但常数因子较大。
Q4: twiddle_idx = k * step 的推导
这个公式源自 Cooley-Tukey 的数学推导。对于 N 点 FFT，当分解为 group_size 大小的子问题后，组内第 k 个蝶形使用的旋转因子应该是 W_N^{k·N/group_size}。而 step = N / group_size，所以 twiddle_idx = k * step。
Stage 1 中 step = N/2 = 4，所以 twiddle_idx = k*4。由于 k 只有 0（group_size=2, half=1），所有蝶形都使用 W_N^0 = 1。这是合理的——Stage 1 的蝶形是首次合并长度为 1 的 DFT 对，它们之间还没有相位差。
Q5: 共轭对称性的实际用途
对于实数输入，只需计算 X[0] 到 X[N/2]（共 N/2+1 个值），后半部分由 X[N-k] = conj(X[k]) 自动得到。节省了约 N/2 - 1 个频域值的计算——接近 50% 的计算量。
实际 FFT 库（如 FFTW）专门提供了 r2c（real-to-complex）接口，利用这一性质：
- 输入：N 个实数
- 输出：N/2+1 个复数（非冗余部分）
- 计算量：约为复数 FFT 的一半
对于 N=1024，复数 FFT 需要约 5120 次蝶形运算；实数 FFT 只需约 2560 次。
Q6: DIT vs DIF 两种实现
| 方式 | 输入顺序 | 输出顺序 | 蝶形操作 |
|---|---|---|---|
| DIT (时域抽取) | 位反转 | 正常序 | 先乘后加减 |
| DIF (频域抽取) | 正常序 | 位反转 | 先加减后乘 |
DIT（本题采用）：将时域序列按奇偶分解，因此需要位反转输入。蝶形公式为 upper = a + W×b, lower = a - W×b。
DIF：将频域序列按前后分解，正常序输入、位反转输出。蝶形公式为 upper = a + b, lower = (a - b)×W——注意乘法的位置不同。
两者计算量完全相同，选择取决于应用场景：如果需要正常序输出（如实时频谱显示），选 DIT；如果需要正常序输入（如流式数据），选 DIF。
课后练习
逆 FFT (IFFT)。实现逆变换：将旋转因子取共轭（
W_N^k → conj(W_N^k)），每一步蝶形后除以 2（或最终除以 N）。用 IFFT 从频域结果恢复原始时域信号。知识点提示：IFFT 公式为
x[n] = (1/N)·Σ X[k]·e^{+2πi·k·n/N}。只需将 W8^k 替换为 conj(W8^k) 并在最后除以 N 即可。验证：对 FFT 输出做 IFFT 应还原原序列。参考解答
c/* 逆 FFT: 将旋转因子取共轭，最后除以 N */ void ifft(Complex a[], int n, const Complex W[]) { /* 构造共轭旋转因子 */ Complex WC[4]; for (int i = 0; i < 4; i++) { WC[i].real = W[i].real; WC[i].imag = -W[i].imag; /* 取共轭 */ } bit_reverse(a, n); for (int stage = 1; stage <= 3; stage++) butterfly(a, n, stage, WC); /* 除以 N */ for (int i = 0; i < n; i++) { a[i].real /= n; a[i].imag /= n; } }扩展 N=16。将 N 扩展到 16，增加旋转因子表到 8 个条目，stage 循环从 3 变为 4。测试序列
{1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0}（交替方波）的频谱。知识点提示：N=16 时 stage=1..4，group_size 为 2,4,8,16。旋转因子 W16^k = cos(2πk/16)-i·sin(2πk/16)。方波的频谱仅在奇次谐波处非零——这正是"方波 = 奇次谐波叠加"的数学证明。
参考解答
c#define N 16 const Complex W[N/2] = { { 1.000000, 0.000000}, // W16^0 { 0.923880,-0.382683}, // W16^1 { 0.707107,-0.707107}, // W16^2 { 0.382683,-0.923880}, // W16^3 { 0.000000,-1.000000}, // W16^4 {-0.382683,-0.923880}, // W16^5 {-0.707107,-0.707107}, // W16^6 {-0.923880,-0.382683}, // W16^7 }; /* fft() 中 stage 循环改为 for (stage=1; stage<=4; stage++) */ /* 其余代码完全相同——这就是 FFT 的优雅之处 */方波
{1,0,1,0...}的频谱：X[4]=8（DC），X[2]=0, X[6]=0, X[10]=0, X[14]=0（偶次谐波为零），奇次谐波非零。幅度谱计算。编写
magnitude_spectrum(Complex X[], int n)函数，计算每个频率分量的幅度|X[k]| = sqrt(real² + imag²)并以柱状图形式打印（用#字符）。知识点提示：幅度谱揭示了信号的能量分布。对于本题
{1,2,3,4,4,3,2,1}，|X[0]| = 20.00 最大（DC 分量），|X[1]| = |X[7]| ≈ 6.31（基频），其他谐波很小或为零。柱状图用#表示幅度，每个#代表 1 单位幅度。参考解答
c#include <math.h> void magnitude_spectrum(Complex X[], int n) { double max_mag = 0.0; double mags[N]; for (int i = 0; i < n; i++) { mags[i] = sqrt(X[i].real * X[i].real + X[i].imag * X[i].imag); if (mags[i] > max_mag) max_mag = mags[i]; } printf("\nMagnitude Spectrum:\n"); for (int i = 0; i < n; i++) { printf(" |X[%d]| = %8.3f ", i, mags[i]); int bars = (int)(mags[i] / max_mag * 40); for (int j = 0; j < bars; j++) printf("#"); printf("\n"); } }
参考资料
- Cooley, J.W. & Tukey, J.W. (1965). "An Algorithm for the Machine Calculation of Complex Fourier Series". Mathematics of Computation, Vol. 19, pp. 297-301. — FFT 诞生的标志性论文
- Brigham, E.O. (1988). The Fast Fourier Transform and Its Applications. Prentice Hall. — FFT 经典教材
- Smith, S.W. (1997). The Scientist and Engineer's Guide to Digital Signal Processing. Chapter 12. — 面向工程师的 DSP 指南
- Numerical Recipes in C, 2nd Edition. Chapter 12: Fast Fourier Transform. — 实用 FFT 实现参考
- FFTW 官网：https://www.fftw.org/ — 工业级 FFT 库，支持自适应优化
"The FFT is perhaps the most important numerical algorithm of our lifetime." — Gilbert Strang