3. 快速傅里叶变换（FFT）：从数学公式到工程利器

说实话，我刚入行做音频处理那会儿，第一次看到 DFT 的公式，心里是有点发怵的。N 点 DFT 要算 N² 次复数乘法和加法，采样率 44.1kHz 的 1 秒音频，N=44100，算一次就要将近 20 亿次运算。这在嵌入式设备上根本跑不动。后来我接触到 FFT，才真正体会到什么叫「数学之美」——它把计算量从 N² 降到了 N log₂N，这个差距有多大？我们后面用数据说话。

3.1 蝶形运算：FFT 的核心思想

FFT 说白了，就是利用 DFT 中旋转因子的对称性和周期性，把一个大问题拆成两个小问题，再递归地拆下去。我习惯把它理解成「分而治之」的典范。

以最常用的基-2 时间抽取 FFT（DIT-FFT）为例，它的核心操作就是「蝶形运算」。为什么叫蝶形？你画出来看看就知道了——两个输入节点、两个输出节点，交叉连接，形状像蝴蝶的翅膀。

X[k] = x[n] + W_N^k · x[n+N/2]
X[k+N/2] = x[n] - W_N^k · x[n+N/2]

嗯，这里要注意：每次蝶形运算只需要 1 次复数乘法和 2 次复数加法。而直接 DFT 的一个点就需要 N 次复数乘法。差距就是这么拉开的。

我曾经在优化一个实时语音降噪算法时，发现 FFT 的蝶形运算顺序搞反了，结果频域数据全乱了。排查了整整两天，最后发现是旋转因子的索引写错了。从那以后，我每次写 FFT 都会先在纸上画一遍 8 点的蝶形图，确认数据流方向。

3.2 计算复杂度对比：N² vs N log₂N

你想想看，如果 N=1024，直接 DFT 需要 1024² ≈ 104 万次复数乘法。而 FFT 只需要 1024 × log₂1024 = 1024 × 10 = 10240 次。差了整整 100 倍。N 越大，这个差距越恐怖。

3.3 Python 实现与可视化

我习惯用 NumPy 的 FFT 接口，但为了让你理解底层原理，我们先手写一个基-2 DIT-FFT。嗯，这里要注意：输入长度必须是 2 的幂次，否则需要补零。

import numpy as np
import matplotlib.pyplot as plt

def dit_fft(x):
    """基-2 时间抽取 FFT（递归实现）"""
    N = len(x)
    if N <= 1:
        return x
    
    # 奇偶分离
    even = dit_fft(x[0::2])
    odd = dit_fft(x[1::2])
    
    # 蝶形运算
    W = np.exp(-2j * np.pi * np.arange(N//2) / N)
    return np.concatenate([
        even + W * odd,
        even - W * odd
    ])

# 生成测试信号：50Hz + 120Hz 正弦波叠加
fs = 1000  # 采样率 1kHz
t = np.arange(0, 1, 1/fs)
x = 0.7 * np.sin(2 * np.pi * 50 * t) + 1.2 * np.sin(2 * np.pi * 120 * t)

# 计算 FFT
N = 1024
X = dit_fft(x[:N])
freqs = np.fft.fftfreq(N, 1/fs)

# 可视化
plt.figure(figsize=(10, 6))
plt.plot(freqs[:N//2], np.abs(X[:N//2]))
plt.title('FFT 频谱（手写实现）')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.grid(True, alpha=0.3)
plt.show()

实际工程中，我建议直接用 NumPy 的 np.fft.fft()，它底层调用的是 FFTPACK 或 MKL，经过高度优化。但理解蝶形运算的原理，能帮你更好地理解频谱泄漏、加窗、补零这些后续概念。

3.4 工程实践中的注意事项

• 补零策略：如果信号长度不是 2 的幂次，补零到下一个 2 的幂次。补零不会增加频率分辨率，但能让 FFT 点数更规整，计算更快。
• 实数 FFT（RFFT）：对于实数信号，可以利用共轭对称性只算一半，计算量再减半。NumPy 的 np.fft.rfft() 就是干这个的。
• 位反转排序：DIT-FFT 的输入需要按位反转顺序排列。我早期写代码时忘了这一步，结果频谱全乱了。后来我写了个小工具函数专门做位反转。

def bit_reverse(x):
    """位反转排序（用于 DIT-FFT）"""
    N = len(x)
    bits = int(np.log2(N))
    indices = np.arange(N)
    rev_indices = np.zeros(N, dtype=int)
    for i in range(N):
        rev = 0
        for b in range(bits):
            if i & (1 << b):
                rev |= 1 << (bits - 1 - b)
        rev_indices[i] = rev
    return x[rev_indices]