4、快速傅里叶变换（FFT）：FFT算法原理、计算复杂度优势及常见变体（基2、基4等）
在TinyML的振动故障诊断场景中，传感器采集到的时域信号往往混杂着多种频率成分。直接对时域波形进行统计特征提取（如均值、峰值、峭度）虽然计算简单，却难以区分不同故障源对应的频率特征。快速傅里叶变换（Fast Fourier Transform, FFT）作为离散傅里叶变换（DFT）的高效实现，能够在极短的时间内将时域信号转换为频域表示，是频带能量分析和故障特征提取的基石。本章将深入剖析FFT的数学原理、计算复杂度优势，并对比基2、基4等常见算法变体，帮助学习者在资源受限的嵌入式设备上做出合理的算法选型。
4.1 FFT算法原理：从DFT到分治策略
离散傅里叶变换（DFT）将长度为N的时域序列 x[n] 转换为频域序列 X[k]，其数学定义为：
X[k] = Σ_{n=0}^{N-1} x[n] · W_N^{kn}，其中 W_N = e^{-j2π/N}直接计算DFT需要 N² 次复数乘法和 N(N-1) 次复数加法。当N=1024时，约需100万次复数乘法，这对TinyML设备上的MCU（如Cortex-M4，主频80MHz）而言，单次变换耗时可能超过50ms，无法满足实时性要求。FFT的核心突破在于利用旋转因子 W_N^{kn} 的两个重要性质：
- 对称性：W_N^{k+N/2} = -W_N^k，即旋转因子在复平面上关于原点对称。
- 周期性：W_N^{k+N} = W_N^k，旋转因子以N为周期重复。
基于上述性质，Cooley-Tukey算法将N点DFT分解为两个N/2点DFT的递归组合。以基2时间抽取（DIT）为例，将输入序列按奇偶位置分为两组：
x_even[n] = x[2n]，x_odd[n] = x[2n+1]，n = 0,1,...,N/2-1则N点DFT可表示为：
X[k] = X_even[k] + W_N^k · X_odd[k]X[k+N/2] = X_even[k] - W_N^k · X_odd[k]其中 X_even[k] 和 X_odd[k] 分别是偶数序列和奇数序列的 N/2 点DFT。这一分解过程递归进行，直到分解为2点DFT（即蝶形运算的基本单元）。整个计算流程形成 log₂N 级流水，每级包含 N/2 个蝶形运算，每个蝶形仅需1次复数乘法和2次复数加法。
蝶形运算单元：基2蝶形接收两个输入 a 和 b，输出 A = a + W·b 和 B = a - W·b。在嵌入式实现中，复数乘法通常被拆分为4次实数乘法和2次实数加法，因此优化重点在于减少乘法次数。
4.2 计算复杂度优势与工程意义
FFT将DFT的计算复杂度从 O(N²) 降至 O(N log₂N)。下表对比了不同N值下的乘法次数差异：
| 变换点数 N | DFT 复数乘法次数 | 基2 FFT 复数乘法次数 | 加速比 | MCU耗时估算（80MHz） |
|---|---|---|---|---|
| 64 | 4,096 | 192 | 21.3× | 约 0.1ms |
| 256 | 65,536 | 1,024 | 64× | 约 0.5ms |
| 1024 | 1,048,576 | 5,120 | 204.8× | 约 2.5ms |
| 4096 | 16,777,216 | 24,576 | 682.7× | 约 12ms |
在TinyML场景中，这一复杂度降低具有三重意义：
- 实时性：1024点FFT在80MHz MCU上仅需2~3ms，满足10ms级控制周期要求。
- 功耗控制：计算时间缩短意味着CPU可以更快进入休眠状态，降低平均功耗。
- 内存占用：FFT可采用原位（in-place）计算，只需N个复数存储单元，而DFT需要额外的N×N中间矩阵。
4.3 常见FFT变体：基2、基4与混合基
不同FFT变体的核心差异在于蝶形运算的基数（radix）和分解策略。选择哪种变体取决于变换点数N的因数分解特性以及硬件资源约束。
4.3.1 基2 FFT（Radix-2）
基2 FFT要求 N = 2ᵐ，是最经典且应用最广泛的变体。其特点包括：
- 结构规整：每级蝶形结构完全相同，便于循环展开和流水线设计。
- 位反转重排：输入序列需按位反转顺序排列，输出则为自然顺序。
- 旋转因子预计算：W_N^k 可提前计算并存储为查找表，避免运行时重复计算三角函数。
在CMSIS-DSP库中，arm_cfft_f32() 函数即采用基2/基4混合实现，针对Cortex-M内核做了指令级优化。对于N=512或N=1024的振动信号帧，基2 FFT是首选方案。
4.3.2 基4 FFT（Radix-4）
基4 FFT要求 N = 4ᵐ，每次分解将N点DFT拆分为4个N/4点DFT。其蝶形运算接收4个输入，输出4个结果，每个蝶形需要3次复数乘法和8次复数加法。相比基2，基4的乘法次数减少约25%，因为旋转因子 W_N^0、W_N^{N/4}、W_N^{N/2}、W_N^{3N/4} 存在大量±1和±j的简化运算。
基4 FFT的级数为 log₄N，比基2少一半。例如N=1024时，基2需要10级，基4仅需5级。但基4蝶形的控制逻辑更复杂，且要求N必须是4的幂次。在工程中，当N=256、1024、4096时，基4 FFT能显著减少计算延迟。
4.3.3 混合基与分裂基 FFT
实际应用中，N往往不是2或4的整数幂。例如N=600（用于特定采样率下的频率分辨率需求）。此时可采用混合基策略：将N分解为多个因子的乘积（如600 = 2³ × 3 × 5²），分别使用不同基数的蝶形组合。分裂基FFT（Split-Radix）则是在基2框架中嵌入基4蝶形，对偶数索引使用基2、对奇数索引使用基4，在N=2ᵐ时达到理论最低乘法次数。
选型建议：在TinyML故障诊断中，优先选择N=2ᵐ（如256、512、1024）以使用基2或基4 FFT。若N无法满足2的幂次，可采用补零（zero-padding）方式将信号延长至2的幂次，但需注意补零不增加频率分辨率，仅使频谱插值更平滑。
4.4 案例拆解：电机轴承故障的FFT特征提取
某工业电机振动监测系统采用STM32F446（Cortex-M4，180MHz）作为主控，加速度传感器以12.8kHz采样率采集振动信号。系统需在每帧数据上提取频带能量特征，用于区分正常状态、外圈故障和内圈故障。
初始参数与约束：
- 采样率 fs = 12.8kHz，每帧时长 80ms，帧长 N = 1024 点。
- 故障特征频率：外圈故障 BPFO ≈ 96.3Hz，内圈故障 BPFI ≈ 149.7Hz。
- MCU可用RAM 128KB，FFT需支持原位计算以节省内存。
- 实时性要求：每帧处理时间 ≤ 10ms。
实现步骤：
- 加窗处理：对1024点时域信号施加汉宁窗，减少频谱泄漏。
- FFT变换：调用CMSIS-DSP的
arm_cfft_f32()，内部自动选择基2/基4混合算法，实测耗时约1.8ms。 - 幅值谱计算：取FFT输出前512个频点的模值，得到幅值谱 |X[k]|。
- 频带能量聚合：将0~6.4kHz划分为8个等宽频带，每个频带宽度800Hz，计算各频带内幅值平方和作为能量特征。
- 特征归一化：将8维能量特征除以总能量，得到归一化频带能量分布向量。
结果对比：
| 状态 | 频带1 (0-800Hz) | 频带2 (800-1600Hz) | 频带3 (1600-2400Hz) | 频带4 (2400-3200Hz) | 诊断结论 |
|---|---|---|---|---|---|
| 正常 | 0.72 | 0.18 | 0.07 | 0.03 | 能量集中于低频 |
| 外圈故障 | 0.41 | 0.33 | 0.18 | 0.08 | 中频段能量上升 |
| 内圈故障 | 0.28 | 0.25 | 0.29 | 0.18 | 高频段显著增加 |
复盘与代价：基2 FFT在该场景下完全满足实时性要求，但若将帧长提升至4096点以获取更高频率分辨率，基4 FFT可将计算时间从约12ms降至约8ms，代价是代码复杂度增加和旋转因子表占用更多Flash（约4KB）。此外，FFT对采样时钟抖动敏感，需确保ADC触发时序稳定，否则频谱会出现噪声基底抬升。
4.5 常见误区与工程总结
误区1：FFT点数越大越好。增大N虽然提高频率分辨率（Δf = fs/N），但计算量和内存占用同步增长。在TinyML中，应根据故障特征频率的最小间隔选择N，而非盲目追求高分辨率。
误区2：忽略窗函数的影响。直接对截断信号做FFT会产生频谱泄漏，导致相邻频带能量串扰。矩形窗旁瓣衰减仅-13dB，而汉宁窗可达-31dB。在频带能量计算中，建议至少使用汉宁窗或汉明窗。
误区3：混淆复数FFT与实数FFT。振动信号是实数序列，使用复数FFT会浪费一半计算量。应使用实数FFT（RFFT）或利用实序列的共轭对称性将N点实数FFT转换为N/2点复数FFT，计算量再减半。
总结复盘：FFT是TinyML频域特征提取的核心引擎。基2 FFT以其结构规整和实现简单成为默认选择；基4 FFT在N为4的幂次时提供约25%的乘法节省；混合基则提供了灵活性。在实际工程中，应结合MCU算力、内存容量、实时性要求和故障特征频率分布，综合确定FFT点数、算法变体和窗函数类型。CMSIS-DSP、Arm-2D等库提供了高度优化的FFT实现，开发者应优先复用这些经过验证的代码，而非从零编写。