3、FFT频谱分析入门:傅里叶变换原理、FFT算法简介、频谱图解读(幅值谱、功率谱)、频率分辨率与窗函数的作用

在TinyML振动检测系统中,误报的根源往往不在于模型本身,而在于我们喂给模型的“特征”是否真实反映了物理世界的振动状态。FFT(快速傅里叶变换)是将时域振动信号转换为频域特征的核心工具。本章将深入剖析傅里叶变换的数学直觉、FFT算法的工程实现、频谱图的正确解读方法,以及频率分辨率与窗函数这两个最容易被忽视却直接影响误报率的关键参数。理解这些内容,是后续诊断“频谱泄漏导致误报”“分辨率不足掩盖故障特征”等问题的理论基础。

FFT频谱分析核心链路与参数影响 时域振动信号 采样率 fs, 点数 N 加速度/速度/位移 预处理 去直流、去趋势项 加窗函数 w[n] FFT计算 X[k] = Σ x[n]·W_N^(nk) 复数输出: 实部+虚部 频谱后处理 幅值谱/功率谱 频率轴映射 关键参数影响 频率分辨率 Δf = fs / N Δf 越小 → 频率分辨能力越强 N 固定时,降低 fs 可提高分辨率 窗函数选择 矩形窗: 分辨率高, 泄漏大 汉宁窗: 泄漏小, 主瓣变宽 频谱泄漏与栅栏效应 非整周期截断 → 能量扩散 峰值偏移 → 幅值误差 TinyML特征提取输出 频段能量 (低频/中频/高频) 峰值频率与幅值 | 频谱平坦度 | 边带能量比 → 输入分类器/异常检测模型 参数调整反馈

3.1 傅里叶变换原理:从时域到频域的桥梁

傅里叶变换的核心思想是:任何周期性信号都可以分解为一系列不同频率、不同幅值和不同相位的正弦波之和。对于振动信号而言,时域波形描述的是“振动幅度随时间的变化”,而频域描述的是“振动能量在不同频率上的分布”。连续傅里叶变换的数学表达式为:

X(f) = ∫ x(t) · e^(-j2πft) dt   (从 -∞ 到 +∞)

在实际嵌入式系统中,我们处理的是离散采样信号,因此使用离散傅里叶变换(DFT):

X[k] = Σ x[n] · e^(-j2πkn/N)   (n = 0, 1, ..., N-1)

其中,x[n] 是时域采样序列,N 是采样点数,X[k] 是第 k 个频率分量的复数结果。DFT 的输出包含实部和虚部,幅值 |X[k]| = √(Re² + Im²),相位 φ = arctan(Im/Re)。对于TinyML振动检测,我们通常只关心幅值信息,因为幅值谱直接反映了各频率成分的振动强度。

工程直觉:想象一个电机以 30Hz 旋转,同时存在 120Hz 的轴承故障特征频率。时域波形看起来像一条“扭曲的波浪”,很难直接区分这两个频率成分。但经过傅里叶变换后,频谱图上会清晰出现两个独立的峰值——一个在 30Hz,一个在 120Hz。这就是频域分析的核心价值:将叠加的复杂信号分离为独立的频率分量。

3.2 FFT算法简介:从DFT到快速计算的工程跨越

直接计算 N 点 DFT 需要 O(N²) 次复数乘法,当 N=1024 时,需要约 100 万次乘法运算,这对MCU来说是不可接受的。FFT(快速傅里叶变换)利用旋转因子的周期性和对称性,将计算复杂度降低到 O(N log₂N)。以 N=1024 为例,FFT 只需约 10240 次运算,比 DFT 快约 100 倍。

最常用的 FFT 实现是基-2 时间抽取算法,其核心步骤为:

  1. 位反转排序:将输入序列按二进制位反转重新排列。
  2. 蝶形运算:迭代执行“蝶形”计算单元,每级合并两个子序列。
  3. 多级分解:将 N 点序列逐级分解为 2 点 DFT,共 log₂N 级。

在TinyML平台上(如STM32、ESP32、Cortex-M系列),通常使用CMSIS-DSP库或KissFFT等轻量级实现。关键参数包括:

  • FFT点数 N:必须为 2 的幂次方(如 256、512、1024),直接影响频率分辨率和计算时间。
  • 采样率 fs:决定可分析的最高频率(奈奎斯特频率 = fs/2)。
  • 数据类型:16位定点(q15)或32位浮点(f32),定点运算更快但动态范围受限。
FFT点数 N频率分辨率 Δf (fs=1000Hz)计算时间 (Cortex-M4 @100MHz)适用场景
2563.906 Hz~0.5 ms快速诊断,粗粒度特征
5121.953 Hz~1.2 ms平衡分辨率与速度
10240.977 Hz~2.8 ms精细频谱分析,低速轴承故障
20480.488 Hz~6.5 ms高精度诊断,但内存占用大

3.3 频谱图解读:幅值谱与功率谱的工程选择

FFT 输出的复数数据需要转换为可用的频谱形式。两种最常见的表示方法是幅值谱和功率谱:

  • 幅值谱(Amplitude Spectrum):A[k] = 2·|X[k]| / N(单边谱),单位与原始信号一致(如 m/s²)。物理意义直观,表示该频率分量的振动幅值。
  • 功率谱(Power Spectrum):P[k] = A[k]²,单位是幅值的平方。功率谱对大幅值分量更敏感,适合突出主导频率,但会压缩小幅值成分的可见性。

在TinyML振动检测中,推荐使用幅值谱作为特征输入,原因有三:其一,幅值谱的线性特性保留了不同频率成分间的真实比例关系;其二,幅值谱对噪声的敏感度低于功率谱(功率谱会放大噪声的平方效应);其三,幅值谱的数值范围更稳定,便于后续归一化处理。

频谱解读要点:读取频谱图时,应关注三个维度——峰值频率(对应故障特征频率)、峰值幅值(反映振动能量大小)、边带结构(峰值两侧的调制边带,常指示轴承磨损或齿轮啮合问题)。不要只看最高峰,低频区的微小峰值往往隐藏着早期故障信息。

3.4 频率分辨率与窗函数:误报的两大隐形杀手

频率分辨率定义为 Δf = fs / N,表示频谱中相邻两条谱线之间的频率间隔。Δf 越小,频率分辨能力越强。例如,当 fs=1000Hz、N=512 时,Δf≈1.95Hz,这意味着两个频率差小于 1.95Hz 的信号在频谱上无法区分。对于低速旋转机械(如 10Hz 转频),其故障特征频率可能仅相差 0.5Hz,此时必须提高频率分辨率。

提高分辨率的方法有两种:增加 N(延长采样时间)或降低 fs(但需保证 fs > 2×最高关注频率)。在TinyML场景中,N 受限于内存和计算时间,fs 受限于奈奎斯特采样定理,因此需要权衡。

窗函数的作用是抑制频谱泄漏。当采样时间不是信号周期的整数倍时,FFT 会假设信号在截断边界处周期性延拓,导致能量从真实频率“泄漏”到相邻频率上。矩形窗(即不加窗)的泄漏最严重,主瓣窄但旁瓣高;汉宁窗(Hann)旁瓣衰减快,但主瓣宽度增加约一倍;布莱克曼窗(Blackman)旁瓣更低,但主瓣更宽。

窗函数主瓣宽度旁瓣衰减幅值精度适用场景
矩形窗2Δf-13 dB±1 dB瞬态信号、整周期采样
汉宁窗4Δf-31 dB±0.1 dB一般振动信号(推荐)
布莱克曼窗6Δf-58 dB±0.01 dB强干扰下的微弱信号
平顶窗~5Δf-70 dB±0.001 dB幅值校准场景
常见误区:很多工程师为了追求高分辨率,盲目增加 FFT 点数 N,却忽略了窗函数的影响。实际上,当使用汉宁窗时,主瓣宽度为 4Δf,这意味着即使 Δf 很小,两个频率差小于 4Δf 的信号仍然无法分辨。正确的做法是先确定所需频率分辨率,再反推 N 和窗函数类型。例如,需要分辨 0.5Hz 的频率差,使用汉宁窗时,Δf 应设为 0.125Hz,即 N = fs / 0.125 = 8000,取 N=8192。

3.5 案例拆解:风机轴承故障的频谱误报诊断

背景与约束:某工厂的离心风机配备TinyML振动监测节点,采样率 fs=2000Hz,FFT点数 N=512,使用矩形窗。近期系统频繁报警“轴承故障”,但现场检查轴承并无异常。

初始现象:频谱图在 50Hz 附近出现宽峰,幅值超过报警阈值。工程师怀疑是电源工频干扰,但检查接地后问题依旧。

关键数据:风机转频 fr=24.5Hz,轴承外圈故障特征频率 BPFO≈3.5×fr≈85.75Hz。频谱图中 50Hz 处峰值幅值 0.42 m/s²,85.75Hz 处幅值仅 0.08 m/s²。

分析判断过程:

  1. 首先检查频率分辨率:Δf = 2000/512 = 3.906Hz。这意味着 85.75Hz 与 89.65Hz(相邻谱线)无法区分,BPFO 特征可能被“淹没”在相邻谱线中。
  2. 50Hz 宽峰的实际中心频率为 48.8Hz(第 12 条谱线),并非精确的 50Hz。由于矩形窗的旁瓣泄漏,24.5Hz 转频的谐波(2×24.5=49Hz)能量泄漏到了相邻谱线,形成“假峰值”。
  3. 改用汉宁窗后,49Hz 处的旁瓣被抑制,50Hz 附近峰值幅值降至 0.15 m/s²,不再超阈值。
  4. 同时将 N 增加到 1024,Δf 降至 1.953Hz,85.75Hz 处出现明显峰值(0.21 m/s²),确认轴承早期磨损。

实际操作步骤:

// 伪代码示例:窗函数与FFT参数调整
#define FS 2000
#define N 1024
float hann[N];
for (int i = 0; i < N; i++) {
    hann[i] = 0.5 * (1 - cos(2 * PI * i / (N - 1)));
}
// 对采样数据加窗
for (int i = 0; i < N; i++) {
    windowed_data[i] = raw_data[i] * hann[i];
}
// 执行FFT(CMSIS-DSP)
arm_rfft_fast_f32(&fft_instance, windowed_data, fft_output, 0);
// 计算幅值谱(单边)
for (int k = 0; k < N/2; k++) {
    magnitude[k] = 2 * sqrt(fft_output[2*k]^2 + fft_output[2*k+1]^2) / N;
}

结果对比:调整后,系统在 85.75Hz 处稳定检测到峰值,且幅值随运行时间缓慢增长(从 0.18 到 0.25 m/s²),成功预警了轴承磨损。误报率从每周 3 次降至 0。

风险与代价:增加 N 导致采样时间从 256ms 延长至 512ms,检测响应时间增加约 300ms,但对于风机这类稳态设备可接受。汉宁窗导致主瓣变宽,两个相距 5Hz 以内的频率成分可能无法区分,需结合具体设备特征频率评估。

复盘经验:① 频谱误报的首要排查项是窗函数与频率分辨率的匹配性;② 矩形窗只适用于整周期采样或瞬态分析,稳态振动监测应默认使用汉宁窗;③ 报警阈值应基于加窗后的幅值谱重新标定,不能直接沿用矩形窗的阈值。

3.6 常见误区与总结复盘

  • 误区一:FFT点数越大越好。实际上,N 增大意味着采样时间变长,对于非平稳信号(如变转速设备),过长的采样窗口会模糊时变特征。应优先保证采样时间内信号近似平稳。
  • 误区二:忽略直流分量。FFT 输出的第 0 条谱线是直流分量(信号平均值),如果不做去直流处理,其泄漏会污染低频区。建议在 FFT 前减去信号均值。
  • 误区三:混淆单边谱与双边谱。FFT 输出包含正负频率,实际使用时应取单边谱(0~fs/2),并将幅值乘以 2(直流分量除外)。
  • 误区四:窗函数只影响幅值精度。实际上,窗函数对频率分辨率和旁瓣抑制的影响更为关键,直接决定了能否分辨相邻故障特征频率。
本章核心要点:① 傅里叶变换将时域振动信号分解为频率分量,是TinyML振动特征提取的数学基础;② FFT 将 DFT 计算复杂度从 O(N²) 降至 O(N log₂N),使嵌入式实时分析成为可能;③ 幅值谱比功率谱更适合作为模型输入,物理意义明确且抗噪性更好;④ 频率分辨率 Δf=fs/N 和窗函数选择是影响误报率的两大关键参数,必须联合设计;⑤ 实际工程中,应优先使用汉宁窗,并根据目标频率间隔反推 FFT 点数。