4、FFT原理与频谱分析：离散傅里叶变换（DFT）到FFT的推导、频谱分辨率、幅值谱与功率谱
在TinyML振动异常检测中，FFT（快速傅里叶变换）是将时域振动信号转换为频域特征的核心数学工具。无论是轴承磨损、齿轮断齿还是电机不平衡，其故障特征往往以特定频率成分的幅值变化形式隐藏在原始波形中。理解从DFT到FFT的推导逻辑、掌握频谱分辨率与窗函数的关系、区分幅值谱与功率谱的物理含义，是后续进行窗函数选择和参数调优的理论基石。本章将系统拆解FFT的数学原理、工程参数及其在资源受限MCU上的实现约束。
4.1 从DFT到FFT：数学推导与计算效率革命
核心原理：离散傅里叶变换（DFT）将长度为N的有限长序列x[n]映射为N个频域复数值X[k]，其定义为：
X[k] = Σ_{n=0}^{N-1} x[n] · e^(-j·2π·k·n/N), k = 0, 1, ..., N-1直接计算每个X[k]需要N次复数乘法和N-1次复数加法，计算全部N个频点则总复杂度为O(N²)。当N=1024时，约需100万次复数运算，这对TinyML设备上的实时处理是不可接受的。FFT（快速傅里叶变换）的核心思想是利用旋转因子W_N = e^(-j·2π/N)的周期性和对称性，将N点DFT分解为多个小规模DFT的组合。以基-2时间抽取（DIT）算法为例，将序列按奇偶位置拆分为两个N/2点子序列，递归执行直到2点蝶形运算。每次蝶形运算仅需1次复数乘法和2次复数加法，总计算量降为O(N·log₂N)。当N=1024时，运算量约为10240次，比直接DFT减少约99%。
4.2 频谱分辨率：窗函数与FFT长度的耦合关系
核心原理：频谱分辨率Δf定义为相邻两个频点之间的频率间隔，计算公式为Δf = fs / N。其中fs为采样率，N为FFT点数。分辨率决定了频谱中能区分两个相邻频率分量的最小间隔。例如，当fs=1000Hz、N=256时，Δf≈3.9Hz；若两个故障特征频率相差小于3.9Hz，则它们在频谱中会合并为一个峰，无法区分。
判断方法与操作流程：
- 确定目标频率间隔：根据设备转速和故障机理，计算需要区分的最小频率差（如轴承外圈故障特征频率与转频的间隔）。
- 选择FFT点数：在采样率fs固定的前提下，N必须满足N ≥ fs / Δf_required。例如需要0.5Hz分辨率且fs=2000Hz，则N≥4000，实际取N=4096。
- 考虑窗函数主瓣宽度：矩形窗主瓣宽度为2Δf，Hann窗为4Δf，Blackman窗为6Δf。加窗后实际可分辨的最小频率间隔约为窗函数主瓣宽度对应的频率范围。因此，若使用Hann窗，则实际分辨率约为4·fs/N。
- 权衡时间长度与实时性：N越大，需要采集的时域数据长度T=N/fs越长，频谱更新周期越长。在TinyML场景中，需在分辨率与响应速度之间折中。
| FFT点数N | 采样率fs (Hz) | Δf (Hz) | Hann窗实际分辨能力 (Hz) | 数据时长T (s) |
|---|---|---|---|---|
| 256 | 1000 | 3.91 | 15.6 | 0.256 |
| 512 | 1000 | 1.95 | 7.8 | 0.512 |
| 1024 | 1000 | 0.98 | 3.9 | 1.024 |
| 2048 | 2000 | 0.98 | 3.9 | 1.024 |
4.3 幅值谱与功率谱：物理含义与工程选择
核心原理：幅值谱（Amplitude Spectrum）直接反映各频率分量的幅度大小，其计算为|X[k]| = √(Re(X[k])² + Im(X[k])²)。对于实信号，单边幅值谱的幅值等于时域信号中该频率分量的实际振幅（需乘以2/N进行归一化）。功率谱（Power Spectrum）定义为P[k] = |X[k]|² / N，表示信号功率在各频率上的分布密度，其物理单位是信号幅值的平方。功率谱对幅值较小的频率分量有压缩作用，更突出能量集中的频段。
判断方法与操作流程：
- 故障特征提取：当需要检测特定频率的幅值变化（如齿轮啮合频率的幅值增大），使用幅值谱更直观，可直接与阈值比较。
- 能量分布分析：当需要评估振动信号在不同频带的能量占比（如低频段能量占比），使用功率谱更合适，因为能量具有可加性。
- 归一化处理：在TinyML模型中，幅值谱通常需要除以最大值进行归一化，而功率谱常取对数（dB单位）以压缩动态范围。
// CMSIS-DSP 幅值谱计算示例（伪代码）
arm_cfft_f32(&fft_instance, input, 0, 1); // 执行FFT
for (k = 0; k < N/2; k++) {
float re = input[2*k];
float im = input[2*k+1];
amplitude[k] = sqrtf(re*re + im*im) * 2.0f / N; // 单边幅值归一化
power[k] = (re*re + im*im) / N; // 功率谱
}4.4 案例拆解：电机轴承故障的频谱分析实战
场景与约束：某工厂一台额定转速为1470rpm（转频fr=24.5Hz）的异步电机，配置加速度传感器，采样率fs=2000Hz。目标是在TinyML设备（Cortex-M4，128KB RAM）上实现轴承外圈故障的实时检测。已知轴承外圈故障特征频率BPFO≈3.5×fr≈85.75Hz。
初始现象与目标：电机正常运行时振动幅值较小，但频谱中在85Hz附近出现微弱峰值（幅值约0.02g）。需要设计FFT参数，使该特征频率能被稳定提取，并区分于转频的2倍频（49Hz）和3倍频（73.5Hz）。
关键参数与判断过程：
- 频率间隔分析：BPFO=85.75Hz与3倍转频73.5Hz相差12.25Hz，与4倍转频98Hz相差12.25Hz。因此分辨率需优于12Hz，取Δf≤5Hz较为安全。
- FFT点数选择：N ≥ fs/Δf = 2000/5 = 400，取N=512。此时Δf=3.9Hz，Hann窗实际分辨能力约15.6Hz，可区分12.25Hz的间隔。
- 窗函数选择：由于需要抑制频谱泄漏（转频幅值较大），选用Hann窗，主瓣宽度适中，旁瓣衰减31dB。
- 内存估算：512点FFT需要512×4字节×2（实部+虚部）=4KB，加上窗函数和输入缓冲，总内存约6KB，满足128KB RAM约束。
实际操作步骤：
// 1. 采集512点数据（时长0.256s）
// 2. 加Hann窗
for (n = 0; n < 512; n++) {
windowed[n] = raw[n] * 0.5f * (1.0f - cosf(2*PI*n/511));
}
// 3. 执行FFT（CMSIS-DSP）
arm_cfft_f32(&fft_instance, windowed, 0, 1);
// 4. 计算幅值谱并提取85.75Hz附近峰值
float target_amp = amplitude[round(85.75/3.9)]; // 约第22个频点结果对比与复盘：使用N=512+Hann窗后，频谱中85.75Hz处幅值稳定在0.025g±0.003g，而73.5Hz处幅值为0.018g，两者可清晰区分。若改用N=256（Δf=7.8Hz），则85.75Hz与73.5Hz的间隔仅1.6个频点，加窗后频谱泄漏导致两峰合并，无法识别。若使用矩形窗，虽然分辨率更高，但转频的旁瓣泄漏会淹没85Hz附近的微弱信号。最终方案选择N=512+Hann窗，在分辨率和泄漏抑制之间取得平衡。
4.5 常见误区与总结复盘
误区一：FFT点数越大越好。虽然增大N可提高分辨率，但会延长数据采集时间，降低频谱更新率，且增加内存和计算开销。在TinyML场景中，应根据故障特征频率的最小间隔和系统实时性要求综合选择N。
误区二：窗函数只影响幅值精度。窗函数不仅影响幅值精度（相干增益），更关键的是影响频谱泄漏程度和主瓣宽度。矩形窗泄漏最大但分辨率最高，Hann窗泄漏较小但主瓣展宽，Blackman窗旁瓣衰减最大但主瓣最宽。选择窗函数需在泄漏抑制与分辨率之间权衡。
误区三：幅值谱和功率谱可以互换使用。两者物理含义不同，幅值谱保留幅度信息，功率谱强调能量分布。在特征提取时，应根据模型输入需求选择。若使用功率谱，建议取对数或开方以压缩动态范围，避免小幅值特征被淹没。
误区四：忽略FFT的归一化。不同库函数对FFT输出的缩放不同。CMSIS-DSP的arm_cfft_f32输出未归一化，需手动乘以2/N（单边谱）才能得到真实幅值。若直接使用未归一化的幅值进行阈值判断，会导致误判。