3. DSP处理器架构

傅里叶变换，说白了就是把信号从时域搬到频域。你想想看，一段音频波形在时间轴上起起伏伏，我们很难直接看出它包含哪些频率成分。但一旦做了FFT，每个频率的能量分布就一目了然了。

我在音频项目里几乎天天跟FFT打交道。从均衡器到降噪算法，从频谱显示到音高检测，FFT都是绕不开的核心。今天我们就把它彻底讲透。

3.1 从DFT到FFT：为什么需要快速算法？

离散傅里叶变换（DFT）的公式其实很简单：

X[k] = Σ x[n] · W_N^(kn),  n=0..N-1,  k=0..N-1
其中 W_N = e^(-j·2π/N)  —— 我们叫它旋转因子

直接算这个公式，复杂度是 O(N²)。N=1024 时，需要大约 100 万次复数乘加。这在实时音频处理里根本跑不动。

3.2 基2-FFT算法：分而治之

基2-FFT 要求序列长度 N 是 2 的整数次幂。为什么？因为这样才能不断对半拆分。

核心思路就一句话：把大DFT拆成两个小DFT，递归下去。

具体来说，把 x[n] 按奇偶索引分成两半：

x_even = [x[0], x[2], x[4], ...]  长度 N/2
x_odd  = [x[1], x[3], x[5], ...]  长度 N/2

然后：

X[k] = X_even[k] + W_N^k · X_odd[k]
X[k + N/2] = X_even[k] - W_N^k · X_odd[k]

这就是蝶形运算的核心公式。你看，一个 N 点 DFT 变成了两个 N/2 点 DFT，再加上一次蝶形组合。递归下去，直到 N=2 时直接计算。

3.3 蝶形运算：FFT的基本单元

蝶形运算长什么样？我画个图你就明白了。

你看这个图，输入 A 和 B，经过一次蝶形，输出 A+B·W 和 A-B·W。形状像不像蝴蝶展开的翅膀？这就是「蝶形」这个名字的由来。

整个 FFT 就是由 log₂N 级这样的蝶形运算组成的。每一级有 N/2 个蝶形单元。

FFT 点数 N	级数 log₂N	每级蝶形数 N/2	总蝶形数
64	6	32	192
256	8	128	1024
1024	10	512	5120
4096	12	2048	24576

3.4 位反转排序：输入数据的重排

这里有个坑。基2-FFT 要求输入数据按「位反转」顺序排列，而不是自然顺序。

什么意思？举个例子。N=8 时，自然索引是 0,1,2,3,4,5,6,7。写成 3 位二进制：

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。

位反转排序的实现很简单：

// 位反转排序，N 必须是 2 的幂
void bit_reverse(float* real, float* imag, int N) {
    int j = 0;
    for (int i = 0; i < N - 1; i++) {
        if (i < j) {
            // 交换 real[i] 和 real[j]
            float tmp = real[i];
            real[i] = real[j];
            real[j] = tmp;
            // 交换 imag[i] 和 imag[j]
            tmp = imag[i];
            imag[i] = imag[j];
            imag[j] = tmp;
        }
        int k = N >> 1;
        while (k <= j) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
}

3.5 FFT实现优化：让代码跑得更快

理论讲完了，来点实战的。我在音频项目里总结了几条 FFT 优化经验。

3.5.1 预计算旋转因子

每次蝶形运算都要计算 W_N^k = cos(2πk/N) - j·sin(2πk/N)。如果每次都现场算三角函数，那性能就崩了。

我的做法：初始化时一次性算好所有旋转因子，存到查找表里。运行时直接查表。

// 预计算旋转因子
float* cos_table = (float*)malloc(N * sizeof(float));
float* sin_table = (float*)malloc(N * sizeof(float));
for (int i = 0; i < N; i++) {
    cos_table[i] = cos(2 * M_PI * i / N);
    sin_table[i] = -sin(2 * M_PI * i / N);
}

3.5.2 原地计算

FFT 可以原地计算——输入输出共用同一块内存。这样省了一半的内存带宽。对于嵌入式系统来说，这很关键。

3.5.3 循环展开与SIMD

现代 CPU 的 SIMD 指令集（ARM 的 NEON、x86 的 SSE/AVX）可以一次处理多个浮点运算。把蝶形运算写成 SIMD 友好的形式，性能能翻几倍。

3.5.4 定点数实现

在低端 DSP 或 MCU 上，浮点运算可能很慢。这时候可以用定点数（Q15、Q31 格式）实现 FFT。精度会损失一些，但速度能快好几倍。

优化手段	加速比（估算）	适用场景
预计算旋转因子	2x ~ 5x	所有平台
原地计算	1.5x ~ 2x	内存受限系统
SIMD 向量化	3x ~ 8x	ARM NEON / x86 SSE
定点数实现	2x ~ 4x	无 FPU 的 MCU
查表 + 循环展开	1.5x ~ 3x	所有平台

3.6 完整FFT实现（基2、原地、查表）

最后，给你一个完整的 C 语言实现。这段代码我在好几个项目里用过，稳定可靠。

void fft(float* real, float* imag, int N) {
    // 1. 位反转排序
    bit_reverse(real, imag, N);

    // 2. 蝶形运算
    for (int len = 2; len <= N; len <<= 1) {
        int half = len >> 1;
        float w_step_cos = cos_table[N / len];
        float w_step_sin = sin_table[N / len];

        for (int i = 0; i < N; i += len) {
            float w_cos = 1.0f;
            float w_sin = 0.0f;

            for (int j = 0; j < half; j++) {
                int idx1 = i + j;
                int idx2 = i + j + half;

                // 蝶形运算
                float t_real = w_cos * real[idx2] - w_sin * imag[idx2];
                float t_imag = w_cos * imag[idx2] + w_sin * real[idx2];

                real[idx2] = real[idx1] - t_real;
                imag[idx2] = imag[idx1] - t_imag;
                real[idx1] = real[idx1] + t_real;
                imag[idx1] = imag[idx1] + t_imag;

                // 更新旋转因子
                float tmp_cos = w_cos * w_step_cos - w_sin * w_step_sin;
                float tmp_sin = w_cos * w_step_sin + w_sin * w_step_cos;
                w_cos = tmp_cos;
                w_sin = tmp_sin;
            }
        }
    }
}