3. 音频信号处理核心链路

说到音频DSP,FFT绝对是绕不开的核心。不管是降噪、回声消除还是音频编码,背后都站着这个算法。我个人习惯把FFT比作「音频世界的显微镜」——它能把时域信号拆成频域,让你看清每个频率分量长什么样。

但问题是,嵌入式系统资源有限。一个浮点FFT跑在Cortex-M4上,可能吃掉你大半的MIPS。所以这一章,咱们重点聊怎么把FFT做快、做省、做准。

本章核心目标:理解FFT的蝶形运算本质,掌握位反转技巧,学会定点FFT的实现方法,以及那些我在项目中踩过的加速坑。

3.1 FFT原理快速回顾

DFT的公式大家都见过:

X[k] = Σ x[n] · W_N^(nk),  n=0..N-1

直接算的话,复杂度是O(N²)。N=1024时,要算一百多万次复数乘法——这在实时音频里根本没法用。FFT的核心思想就是利用旋转因子的对称性和周期性,把大DFT拆成小DFT。

说白了就是分治。把序列按奇偶分成两组,递归计算,最后合并。复杂度降到O(N log₂N)。N=1024时,只需要约一万次操作,差了两个数量级。

旋转因子(Twiddle Factor):

W_N^k = cos(2πk/N) - j·sin(2πk/N)

嗯,这里要注意。旋转因子是复数,而且有对称性:W_N^(k+N/2) = -W_N^k。这个性质在蝶形运算里会反复用到。

3.2 蝶形运算——FFT的基本单元

蝶形运算长什么样?看这个图就明白了:

基2蝶形运算单元 A B A + B·W A - B·W × W 蝶形运算:一次复数乘法和两次复数加法

一个蝶形运算包含:

  • 1次复数乘法(B × W)
  • 2次复数加法(A + 结果,A - 结果)

整个FFT就是由log₂N级蝶形运算组成的。每级有N/2个蝶形,总共(N/2)·log₂N个蝶形。

我的经验:在ARM Cortex-M上,一次复数乘法大约需要4个周期。如果N=1024,光乘法就要约2万个周期。加上内存访问开销,实际跑下来大概要5-8万周期。所以优化空间很大。

3.3 位反转——数据重排的关键

FFT的输入顺序和输出顺序不一样。时域抽取(DIT)的FFT,输入是自然顺序,输出是位反转顺序。频域抽取(DIF)则反过来。

什么叫位反转?举个例子:N=8时,索引3的二进制是011,反转后是110,也就是6。所以输入的第3个点,在输出里排在第6位。

位反转的查表法实现:

// 预计算位反转表,N=1024
uint16_t bit_rev[N];
for (int i = 0; i < N; i++) {
    uint16_t rev = 0;
    for (int j = 0; j < 10; j++) {
        rev = (rev << 1) | ((i >> j) & 1);
    }
    bit_rev[i] = rev;
}

// 重排数据
for (int i = 0; i < N; i++) {
    if (i < bit_rev[i]) {
        swap(&x[i], &x[bit_rev[i]]);
    }
}

注意那个if判断。如果不加,每个点会被交换两次,等于白干。我曾经在这个坑里浪费了半天——调试时发现输出全乱,最后发现是交换逻辑写错了。

避坑指南:位反转表可以提前算好存成常量数组,运行时直接查表。别在实时处理循环里现场算位反转,那会吃掉大量CPU时间。我见过有人每帧都算一遍位反转,结果FFT还没跑完,音频已经卡顿了。

3.4 定点FFT实现

浮点FFT在PC上跑得很欢,但到了嵌入式平台就尴尬了。很多DSP和MCU没有硬件浮点单元(FPU),或者浮点运算慢得离谱。这时候就得用定点数。

Q格式定点数:

  • Q15:1位符号 + 15位小数,范围[-1, 1),精度2⁻¹⁵
  • Q31:1位符号 + 31位小数,范围[-1, 1),精度2⁻³¹

定点FFT的关键是防止溢出。蝶形运算里,A + B·W的结果可能超过1。所以每级运算后需要做缩放。

我常用的策略是「每级右移1位」:

// 定点蝶形运算,Q15格式
int16_t a_re, a_im, b_re, b_im;
int16_t w_re, w_im;

// 计算 B * W
int32_t t_re = (int32_t)b_re * w_re - (int32_t)b_im * w_im;
int32_t t_im = (int32_t)b_re * w_im + (int32_t)b_im * w_re;

// 右移15位回到Q15
t_re >>= 15;
t_im >>= 15;

// 蝶形加减,然后右移1位防溢出
int16_t out0_re = (a_re + (int16_t)t_re) >> 1;
int16_t out0_im = (a_im + (int16_t)t_im) >> 1;
int16_t out1_re = (a_re - (int16_t)t_re) >> 1;
int16_t out1_im = (a_im - (int16_t)t_im) >> 1;

这样每级衰减6dB,N=1024时总共衰减60dB。输出结果需要乘以N才能恢复幅度。不过对于降噪算法来说,我们通常只关心相对幅度,所以这个缩放可以接受。

实现方式 精度 速度(相对) 适用平台
浮点FFT 高(~10⁻⁷) 1x(基准) 有FPU的MCU、PC
Q15定点FFT 中(~10⁻⁴) 3-5x 无FPU的MCU、DSP
Q31定点FFT 较高(~10⁻⁶) 2-3x 32位MCU、DSP
块浮点FFT 高(自适应) 1.5-2x 高端DSP

我的建议:如果平台支持SIMD指令(比如ARM的SMID),可以用Q15格式一次处理两个16位数据。我在Cortex-M4上做过测试,SIMD优化的定点FFT比纯C版本快了将近4倍。

3.5 FFT加速技巧

光有定点实现还不够,实际项目中还得上各种优化手段。下面是我总结的几个实用技巧:

技巧1:旋转因子预计算

别在FFT循环里现场算cos和sin。提前算好存成查找表。N=1024时,只需要N/2=512个复数因子,占用约2KB内存,换来几十倍的加速。

技巧2:循环展开

蝶形运算的内循环可以手动展开。比如一次处理4个蝶形,减少循环开销和分支预测失败。

// 展开4个蝶形
for (int k = 0; k < N/2; k += 4) {
    // 蝶形1
    // 蝶形2
    // 蝶形3
    // 蝶形4
}

技巧3:内存对齐

确保输入输出缓冲区按16字节对齐。这样可以用SIMD加载指令一次读4个16位数据。不对齐的话,有些MCU会触发异常,或者性能骤降。

技巧4:减少内存拷贝

很多FFT实现会先做位反转重排,然后原地计算。但如果你能接受位反转顺序的输出,就可以省掉重排这一步。降噪算法里,频域处理通常不关心顺序,只要知道每个bin对应的频率就行。

我曾经踩过的坑:在某个项目里,我为了省内存,让FFT的输入和输出共用同一个缓冲区。结果发现位反转重排时把还没处理的数据覆盖了。后来改成双缓冲区,虽然多占了一倍内存,但逻辑清晰多了,调试也方便。

3.6 本章知识体系

下面这张图总结了FFT优化的核心脉络:

FFT优化知识体系 FFT原理 蝶形运算 位反转 定点实现 复数乘法优化 查表法 Q格式缩放 防溢出策略 加速技巧:预计算 | 循环展开 | 内存对齐 | SIMD 目标:在嵌入式平台上实现实时、低功耗、高精度的FFT

FFT优化不是一蹴而就的事。你得根据具体平台、实时性要求和精度需求来权衡。我个人习惯是先跑一个浮点版本验证算法,然后逐步定点化、优化,每一步都做性能对比测试。

嗯,这一章的内容就到这。记住:FFT优化的本质是「用空间换时间,用精度换速度」。找到那个平衡点,你的音频DSP系统就能跑得又快又稳。