跳转至

06 频域滤波与数字滤波器

前面几章多半是在时域里处理数据：限幅挡尖峰，平均压随机抖动，一阶滞后换稳定显示。频域滤波关心的是另一个问题：哪些频率应该留下，哪些频率应该压掉。

这在 MCU 测量里很常见：

  • 压力传感器上叠着 50 Hz 工频干扰。
  • 电机电流里有固定 PWM 纹波。
  • ADC 过采样后要降采样，同时不能把高频噪声带下来。
  • 想判断一个固定频率的振动、按键音、谐波或 DTMF 频率是否存在。
  • 需要比一阶滞后更明确的截止频率、更陡的衰减，或者更可控的相位延迟。

频域滤波不是“更高级的平滑”。它更像是给信号画一张频率通行证：低频、某个频段、某个固定频率，分别用不同规则处理。

先看选型表

工程目标 优先考虑 代价和边界
普通低通，要求线性相位 FIR 阶数可能较高，RAM 和乘加次数增加
普通低通，MCU 资源紧张 一阶 IIR、二阶 Biquad 相位非线性，系数不当会不稳定
压制固定频率干扰 陷波 Biquad、整周期平均 干扰频率漂移时效果下降
过采样后降采样 CIC、FIR 抽取 CIC 有通带下垂，积分器可能溢出
检测某个频率是否存在 Goertzel 适合少数频点，不适合看完整频谱
平滑同时保留峰形 Savitzky-Golay 需要固定窗口，实时输出会有延迟

选型时先问四个问题：

  1. 采样率 \(f_s\) 是多少？
  2. 要保留的最高有效频率是多少？
  3. 要压掉的干扰频率在哪里？
  4. 能接受多少延迟、多少 RAM、多少乘加次数？

如果这四件事说不清，滤波器设计就只能靠试。

频域滤波的基本约束

数字滤波器看到的频率范围只有 \(0\) 到 \(f_s/2\)。\(f_s/2\) 是 Nyquist 频率。超过它的模拟频率会混叠到低频，进了 ADC 以后就很难靠数字滤波补救。

所以频域滤波有一个前提：采样前的模拟信号已经被适当限带。 如果 ADC 前没有模拟抗混叠措施，数字 FIR、IIR、陷波都不能把已经混叠进来的假低频变回原来的高频。

一个滤波器通常要描述这些指标：

  • 通带：希望基本保留的频率范围。
  • 阻带：希望明显压低的频率范围。
  • 截止频率 \(f_c\)：低通滤波器中幅度下降到约 \(0.707\) 的位置。
  • 过渡带：从通带过渡到阻带的频率范围。过渡带越窄，滤波器通常越复杂。
  • 相位或延迟：输出相对输入晚多少。控制环、姿态估计和相位测量尤其在意这一点。

下面的算法都围绕这些指标展开。

FIR 滤波器

FIR 是有限冲激响应滤波器。输出只由当前输入和有限个历史输入决定，不使用历史输出：

\[ y[n] = b_0 x[n] + b_1 x[n-1] + \cdots + b_M x[n-M] \]

其中 \(b_0, b_1, \ldots, b_M\) 是滤波器系数，\(M+1\) 是抽头数（tap 数）。滑动平均就是最简单的 FIR：所有系数都等于 \(1/N\)。

为什么 FIR 常用

FIR 的几个优点很适合测量系统：

  • 天然稳定。 没有历史输出反馈，不会因为极点跑到单位圆外而发散。
  • 容易做线性相位。 系数左右对称时，滤波器对不同频率引入相同的群延迟，波形形状不容易被相位扭曲。
  • 适合严格保留波形形状。 比如称重、压力曲线、光谱、振动包络等场景。

代价也很直接：阶数越高，每次更新的乘加次数越多，环形缓冲区也越大。

延迟和阶数

如果 FIR 系数对称，抽头数为 \(N\)，群延迟为：

\[ \text{delay} = \frac{N-1}{2 f_s} \]

例如 \(f_s = 1000\;\text{Hz}\)，\(N = 31\)，群延迟就是 \(15\;\text{ms}\)。这对显示值可能没问题，但对控制环可能太慢。

FIR 的核心取舍是：

\[ \text{过渡带越窄，阻带要求越高} \quad \Rightarrow \quad \text{需要更多抽头，延迟和计算量越大} \]

系数从哪里来

工程中不要随手填 FIR 系数。常见来源有三类：

  • 滑动平均。 所有系数相等，适合简单低通和零点对准固定干扰。
  • 窗函数法。 从理想低通的冲激响应出发，再乘窗函数，适合教学和中等要求场景。
  • 工具设计。 用 Python、MATLAB 或滤波器设计工具生成系数，适合明确通带、阻带和纹波要求的项目。

窗函数法的低通 FIR 常写成：

\[ h[n] = 2\frac{f_c}{f_s} \cdot \text{sinc}\left(2\frac{f_c}{f_s}(n-\frac{M}{2})\right) \cdot w[n] \]

其中 \(w[n]\) 是窗函数，\(M=N-1\)。实际使用时还要把系数归一化，让所有系数之和为 1，这样直流输入不会被放大或缩小。

这里的 \(\text{sinc}(u)\) 使用归一化定义：

\[ \text{sinc}(u)=\frac{\sin(\pi u)}{\pi u} \]

初学阶段可以先不手推系数，但必须知道：FIR 系数决定幅频响应，窗口长度决定延迟和计算量。

MCU C 实现

下面是一个固定 32 抽头 FIR。系数数组 coeff[0] 对应最新样本，coeff[31] 对应最旧样本。第一次初始化时用首个输入填满缓冲区，可以避免从 0 开始造成启动假瞬态。

#include <stdint.h>

#define FIR32_TAPS 32u

typedef struct {
    float x[FIR32_TAPS];
    uint8_t index;
    uint8_t init;
} Fir32;

void fir32_init(Fir32 *f, float first)
{
    for (uint8_t i = 0; i < FIR32_TAPS; i++) {
        f->x[i] = first;
    }
    f->index = 0;
    f->init = 1;
}

float fir32_update(Fir32 *f, const float coeff[FIR32_TAPS], float x)
{
    if (!f->init) {
        fir32_init(f, x);
    }

    f->x[f->index] = x;

    float y = 0.0f;
    uint8_t j = f->index;

    for (uint8_t i = 0; i < FIR32_TAPS; i++) {
        y += coeff[i] * f->x[j];
        j = (j == 0u) ? (FIR32_TAPS - 1u) : (uint8_t)(j - 1u);
    }

    f->index++;
    if (f->index >= FIR32_TAPS) {
        f->index = 0;
    }

    return y;
}

代码用显式判断完成索引回绕，重点是说明环形缓冲区中"最新样本"和"历史样本"的对应关系。实际工程中只要保证索引始终落在 0 到 N-1，具体回绕方式可以按项目习惯实现。

代价和边界

FIR 对孤立尖峰没有特殊免疫力。一个尖峰会按系数分布进入未来 \(N\) 个输出点。遇到尖峰干扰，应先用限幅、中值或 Hampel 处理，再进入 FIR。

FIR 的另一个边界是延迟。想要很窄的过渡带，就要更多抽头；如果控制系统不能接受几十毫秒延迟，二阶 IIR 可能更合适。

FIR 阶数意味着什么

FIR 的"阶数"和抽头数直接相关：\(N\) 个抽头，对应阶数 \(N-1\)。抽头越多，过渡带通常可以越窄，阻带可以压得更深。代价是乘加次数增加、RAM 增加、线性相位延迟增加。

对称 FIR（系数左右对称）有一个实用优化：可以只算一半乘法，另一半对称累加。奇数抽头时中心抽头单独计算。这对 MCU 上节省计算量很有效。

一句话：FIR 多阶不危险，但会慢、会占 RAM、会带来固定延迟。

IIR 与 Biquad

IIR 使用历史输出参与计算。它的好处是用较少状态和较少乘加得到较强的滤波效果；代价是相位通常不是线性的，而且系数不当会不稳定。

一阶滞后已经在第 02 章讲过，它就是最常见的一阶 IIR：

\[ y[n] = y[n-1] + \alpha(x[n]-y[n-1]) \]

这一章重点看工程中更常用的二阶 IIR：Biquad（二阶节）。后面的二阶低通、陷波和高阶 IIR 级联，都会围绕这个结构展开。

Biquad 差分方程

常见归一化形式为：

\[ y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] - a_1 y[n-1] - a_2 y[n-2] \]

这里默认 \(a_0 = 1\)。Biquad 可以实现低通、高通、带通、陷波和峰值滤波。很多音频、振动、电机控制和传感器项目都会把多个 Biquad 串起来形成更高阶滤波器。

什么时候用 Biquad

Biquad 适合这些场景：

  • 比一阶滞后需要更陡的衰减。
  • FIR 延迟太大或乘加次数太多。
  • 需要一个明确的二阶低通、带通或陷波。
  • 采样率和截止频率固定，可以预先生成系数。

不适合这些场景：

  • 对相位线性要求很高。
  • 系数需要频繁在线改变，但没有做平滑或状态处理。
  • 使用无 FPU MCU，且 Q 很高、系数量化误差无法接受。

低通 Biquad 系数

常用的二阶低通可以用 \(f_c\) 和 \(Q\) 生成系数。先定义：

\[ \omega_0 = 2\pi \frac{f_c}{f_s} \]
\[ A = \frac{\sin(\omega_0)}{2Q} \]

未归一化系数为：

\[ \begin{aligned} b_0 &= \frac{1-\cos(\omega_0)}{2} \\ b_1 &= 1-\cos(\omega_0) \\ b_2 &= \frac{1-\cos(\omega_0)}{2} \\ a_0 &= 1 + A \\ a_1 &= -2\cos(\omega_0) \\ a_2 &= 1 - A \end{aligned} \]

实际代码中要把所有 \(b_0,b_1,b_2,a_1,a_2\) 都除以 \(a_0\)，使 \(a_0=1\)。

其中 \(Q\) 决定二阶响应的形状。常见 Butterworth 低通取：

\[ Q \approx 0.707 \]

这个值通带比较平坦，没有明显峰值。\(Q\) 太大时截止附近会出现峰值，甚至让系统更容易振荡。

MCU C 实现

下面使用转置直接 II 型（Direct Form II Transposed）。它只需要两个状态变量，适合 MCU。系数应在初始化时生成，或者直接由工具离线生成后写入配置表。

#include <math.h>
#include <stdint.h>

typedef struct {
    float b0;
    float b1;
    float b2;
    float a1;
    float a2;
    float z1;
    float z2;
} Biquad;

void biquad_reset(Biquad *q)
{
    q->z1 = 0.0f;
    q->z2 = 0.0f;
}

float biquad_update(Biquad *q, float x)
{
    float y = q->b0 * x + q->z1;
    q->z1 = q->b1 * x - q->a1 * y + q->z2;
    q->z2 = q->b2 * x - q->a2 * y;
    return y;
}

void biquad_set_lowpass(Biquad *q, float fs, float fc, float quality)
{
    const float pi = 3.14159265358979323846f;

    if (fs <= 0.0f) {
        q->b0 = 1.0f;
        q->b1 = 0.0f;
        q->b2 = 0.0f;
        q->a1 = 0.0f;
        q->a2 = 0.0f;
        biquad_reset(q);
        return;
    }

    if (fc <= 0.0f) fc = 0.001f * fs;
    if (fc >= 0.49f * fs) fc = 0.49f * fs;
    if (quality <= 0.0f) quality = 0.707f;

    float w0 = 2.0f * pi * fc / fs;
    float c = cosf(w0);
    float s = sinf(w0);
    float a = s / (2.0f * quality);

    float b0 = (1.0f - c) * 0.5f;
    float b1 = 1.0f - c;
    float b2 = (1.0f - c) * 0.5f;
    float a0 = 1.0f + a;
    float a1 = -2.0f * c;
    float a2 = 1.0f - a;

    q->b0 = b0 / a0;
    q->b1 = b1 / a0;
    q->b2 = b2 / a0;
    q->a1 = a1 / a0;
    q->a2 = a2 / a0;
    biquad_reset(q);
}

sinf、cosf 在无 FPU MCU 上代价很高。实际项目中更推荐离线算好系数，或者只在参数改变时计算一次，不要在采样中断里计算三角函数。

上面的 biquad_reset() 会把内部状态清零。如果传感器输入不是从 0 附近开始，输出会有一段启动瞬态。测量系统中可以忽略前几拍输出，或者按当前输入预置状态。

稳定性和调试

Biquad 的稳定性由反馈系数决定。工程上不要随手改 a1、a2。如果输出突然变成很大的数、出现持续振荡或 NaN，优先检查：

  • \(f_c\) 是否小于 \(f_s/2\)。
  • \(Q\) 是否过大。
  • 系数是否已经除以 \(a_0\)。
  • 多个 Biquad 串联时是否每节都单独保存状态。
  • 定点实现时系数和中间变量是否溢出。

高阶 IIR 为什么要拆成 Biquad

需要四阶、六阶甚至更高阶的 IIR 时，不建议直接写高阶差分方程。工程上的标准做法是拆成多个二阶 Biquad 串联，每个 Biquad 保存自己的状态。

为什么不直接写高阶？因为：

  • 系数量化误差在高阶直接型中会被放大，极点容易偏移到单位圆外。
  • Q 值过高时，中间变量的动态范围很大，定点实现容易溢出。
  • 级联顺序不同，数值稳定性和内部动态范围不同。工程上应优先使用设计工具给出的 SOS 顺序和缩放系数；如果手动调整，高 Q 或峰值增益大的节要重点检查中间变量是否溢出。

一句话：IIR 多阶省计算，但稳定性和系数精度风险更高，所以用 Biquad 级联，不要手搓高阶直接型。

陷波滤波器

陷波滤波器用于压制一个固定频率，例如 50 Hz 工频、60 Hz 工频、电机某个机械共振频率、PWM 纹波等。

它不是通用低通。陷波的目标是：尽量只挖掉目标频率附近的一小段，不影响其他频率。

参数含义

陷波设计至少需要三个参数：

  • \(f_s\)：采样率。
  • \(f_0\)：要压制的中心频率。
  • \(Q\)：品质因数，近似表示中心频率与带宽的比例。

近似关系为：

\[ Q \approx \frac{f_0}{\text{BW}} \]

其中 \(\text{BW}\) 是陷波带宽。\(Q\) 越大，陷波越窄；\(Q\) 越小，陷波越宽。

Q 怎么选

如果干扰频率非常稳定，比如采样系统和工频同步良好，可以用较高 Q，让陷波只影响很窄范围。

如果干扰频率会漂，比如实际工频在 49.8 Hz 到 50.2 Hz 之间变化，Q 太高反而压不住。此时应降低 Q，让陷波稍宽一点。

典型起点：

场景 Q 起点
工频 50/60 Hz，频率较稳定 20 ~ 50
电机固定转速纹波 10 ~ 30
频率有明显漂移 5 ~ 15

Biquad 陷波系数

二阶陷波常用系数如下：

\[ \omega_0 = 2\pi \frac{f_0}{f_s} \]
\[ A = \frac{\sin(\omega_0)}{2Q} \]
\[ \begin{aligned} b_0 &= 1 \\ b_1 &= -2\cos(\omega_0) \\ b_2 &= 1 \\ a_0 &= 1 + A \\ a_1 &= -2\cos(\omega_0) \\ a_2 &= 1 - A \end{aligned} \]

同样要把所有系数除以 \(a_0\)。

void biquad_set_notch(Biquad *q, float fs, float f0, float quality)
{
    const float pi = 3.14159265358979323846f;

    if (fs <= 0.0f) {
        q->b0 = 1.0f;
        q->b1 = 0.0f;
        q->b2 = 0.0f;
        q->a1 = 0.0f;
        q->a2 = 0.0f;
        biquad_reset(q);
        return;
    }

    if (f0 <= 0.0f) f0 = 0.001f * fs;
    if (f0 >= 0.49f * fs) f0 = 0.49f * fs;
    if (quality <= 0.0f) quality = 10.0f;

    float w0 = 2.0f * pi * f0 / fs;
    float c = cosf(w0);
    float s = sinf(w0);
    float a = s / (2.0f * quality);

    float b0 = 1.0f;
    float b1 = -2.0f * c;
    float b2 = 1.0f;
    float a0 = 1.0f + a;
    float a1 = -2.0f * c;
    float a2 = 1.0f - a;

    q->b0 = b0 / a0;
    q->b1 = b1 / a0;
    q->b2 = b2 / a0;
    q->a1 = a1 / a0;
    q->a2 = a2 / a0;
    biquad_reset(q);
}

陷波的边界

陷波只适合频率明确的干扰。若噪声是宽带随机噪声，陷波不会让整体抖动明显变小。若干扰频率随时间变化，固定陷波会漏掉一部分能量，甚至在目标附近造成相位变化。

验证陷波效果时，不要只看时域曲线“好像稳了”。应记录原始数据和滤波后数据，比较 \(f_0\) 附近频率分量是否下降，同时确认目标信号频率没有被误伤。

CIC 滤波器

CIC（Cascaded Integrator-Comb）常用于过采样和抽取。它的特点是只用加法和减法，不需要乘法，适合 Sigma-Delta ADC、数字抽取和高速采样后降采样。

一个 CIC 通常包含三部分：

高采样率输入 -> 多级积分器 -> 降采样 -> 多级梳状器 -> 低采样率输出

它解决什么问题

假设 ADC 以很高采样率采集数据，但后续算法只需要较低输出率。直接每 \(R\) 个点取一个点会混叠。CIC 在降采样前先做一种特殊低通，再抽取输出。

一阶 CIC 可以理解为“先累加，再每 \(R\) 点做差”，效果接近 \(R\) 点滑动平均。多级 CIC 相当于多个这种平均效果叠加，阻带更强，但通带下垂也更明显。

频率特性

\(N\) 级、抽取率为 \(R\) 的 CIC 频率响应可以写成：

\[ H(z) = \left( \frac{1-z^{-R}}{1-z^{-1}} \right)^N \]

直流增益为：

\[ G = R^N \]

所以 CIC 输出通常要除以 \(R^N\)，或者在定点系统中右移相应位数。

MCU C 实现

下面是一个 2 级 CIC 抽取器。输入每来一个样本调用一次；只有返回 1 时，*out 才是新的低速输出。

#include <stdint.h>

typedef struct {
    int64_t integ1;
    int64_t integ2;
    int64_t comb1_delay;
    int64_t comb2_delay;
    uint16_t decim;
    uint16_t count;
} Cic2Decimator;

void cic2_init(Cic2Decimator *c, uint16_t decim)
{
    c->integ1 = 0;
    c->integ2 = 0;
    c->comb1_delay = 0;
    c->comb2_delay = 0;
    c->decim = (decim == 0u) ? 1u : decim;
    c->count = 0;
}

uint8_t cic2_update(Cic2Decimator *c, int32_t x, int32_t *out)
{
    c->integ1 += x;
    c->integ2 += c->integ1;

    c->count++;
    if (c->count < c->decim) {
        return 0u;
    }
    c->count = 0;

    int64_t comb1 = c->integ2 - c->comb1_delay;
    c->comb1_delay = c->integ2;

    int64_t comb2 = comb1 - c->comb2_delay;
    c->comb2_delay = comb1;

    int64_t gain = (int64_t)c->decim * (int64_t)c->decim;
    *out = (int32_t)(comb2 / gain);
    return 1u;
}

CIC 的积分器会不断累加，必须认真检查位宽。上面用 int64_t 是为了降低溢出风险，但对直流或近直流输入，积分器值会随运行时间持续增长，没有自然回绕机制。硬件 CIC 通常依赖定点模运算回绕（溢出后自动截断高位），C 语言中有符号整数溢出是未定义行为，不能依赖。工程上应估算：\(|x| \cdot R^2\) 不能超过 int64_t 范围，且长时间运行时需要定期复位或检查。

代价和边界

CIC 的优点是省乘法，适合高速数据流。缺点是通带会下垂：靠近通带边缘的有效信号会被压低。如果需要高精度幅度测量，CIC 后面常接一个补偿 FIR。

CIC 不适合所有低通场景。普通传感器 100 Hz、200 Hz 采样，只想让显示值稳定，用滑动平均或一阶 IIR 往往更简单。

Goertzel 算法

Goertzel 用来检测某个指定频率的能量。它可以看作“只算一个 DFT 频点”的算法。

适合：

  • DTMF 双音检测。
  • 电力系统某次谐波检测。
  • 判断某个振动频率是否存在。
  • 只关心一两个固定频率，而不想计算完整 FFT。

不适合：

  • 想观察完整频谱。
  • 目标频率很多。
  • 频率快速漂移且没有跟踪机制。

频点关系

Goertzel 对一段长度为 \(N\) 的数据工作。第 \(k\) 个 DFT 频点对应频率：

\[ f_k = k \frac{f_s}{N} \]

若目标频率是 \(f_0\)，应选：

\[ k \approx \frac{f_0 N}{f_s} \]

如果 \(f_0\) 不能刚好落在某个 DFT bin 上，会出现频谱泄漏，检测能量会分散到邻近频点。工程上可以调整 \(N\)，让 \(f_0\) 尽量对准 bin；也可以加窗，但加窗会改变幅度标定。

递推公式

Goertzel 的递推为：

\[ s[n] = x[n] + 2\cos\left(2\pi\frac{k}{N}\right)s[n-1] - s[n-2] \]

设：

\[ \text{coeff} = 2\cos\left(2\pi\frac{k}{N}\right) \]

最后能量为：

\[ P = s_1^2 + s_2^2 - \text{coeff}\cdot s_1 s_2 \]

其中 \(s_1\) 和 \(s_2\) 是递推结束后的最后两个状态。

MCU C 实现

#include <math.h>
#include <stdint.h>

float goertzel_coeff(uint16_t k, uint16_t n)
{
    const float pi = 3.14159265358979323846f;

    if (n == 0u) {
        return 0.0f;
    }

    return 2.0f * cosf(2.0f * pi * (float)k / (float)n);
}

float goertzel_power(const float *x, uint16_t n, float coeff)
{
    if ((x == 0) || (n == 0u)) {
        return 0.0f;
    }

    float s0 = 0.0f;
    float s1 = 0.0f;
    float s2 = 0.0f;

    for (uint16_t i = 0; i < n; i++) {
        s0 = x[i] + coeff * s1 - s2;
        s2 = s1;
        s1 = s0;
    }

    return s1 * s1 + s2 * s2 - coeff * s1 * s2;
}

goertzel_coeff() 只需要在初始化时调用。采样过程中不要每个样本都算 cosf。

阈值怎么定

Goertzel 输出的是能量，不是简单的“有/没有”。阈值通常要靠数据标定：

  1. 采集没有目标频率时的数据，记录背景能量范围。
  2. 采集有目标频率时的数据，记录目标能量范围。
  3. 选择两者之间有裕量的位置作为阈值。
  4. 增加连续多帧确认，避免偶发噪声误判。

如果输入幅度会变化，最好同时估计总能量，用目标频点能量占比作为判断指标。

Savitzky-Golay 平滑

Savitzky-Golay 滤波不是简单平均，而是在窗口内用低阶多项式拟合数据。它的特点是：平滑噪声的同时，比较能保留峰形和斜率。

适合：

  • 光谱、色谱、峰值分析。
  • 缓慢变化但需要保留局部形状的曲线。
  • 离线数据处理或允许固定延迟的实时处理。

不适合：

  • 尖峰异常值很多的场景。
  • 不能接受窗口延迟的控制环。
  • 信号本身突变很快的场景。

5 点二阶平滑

一个常用的 5 点二阶 SG 平滑系数为：

\[ \frac{1}{35}[-3,\ 12,\ 17,\ 12,\ -3] \]

如果窗口内 5 个样本为 \(x_0,x_1,x_2,x_3,x_4\)，输出对应中心点 \(x_2\)：

\[ y_2 = \frac{-3x_0 + 12x_1 + 17x_2 + 12x_3 - 3x_4}{35} \]

注意：这是中心窗口滤波。要输出 \(y_2\)，必须已经拿到 \(x_3\) 和 \(x_4\)，所以实时使用时至少有 2 个采样周期延迟。

float sg5_smooth_center(const float x[5])
{
    return (-3.0f * x[0] + 12.0f * x[1] + 17.0f * x[2]
            + 12.0f * x[3] - 3.0f * x[4]) / 35.0f;
}

SG 本质上也是 FIR，只是系数来自多项式拟合。它不应该拿来替代异常值剔除；如果窗口里有大尖峰，拟合也会被拉偏。

工程验证方法

频域滤波器不要只靠“看曲线顺眼”验证。至少做四类检查：

  1. 正弦测试。 输入不同频率、相同幅度的正弦波，记录输出幅度，检查是否符合预期幅频响应。
  2. 阶跃测试。 输入阶跃，观察延迟、过冲和稳定时间。IIR 尤其要看有没有振荡。
  3. 实际日志回放。 用真实采集数据离线跑滤波器，对比原始数据和滤波后数据。
  4. 边界测试。 检查启动阶段、参数极限、溢出、NaN、采样率变化和异常输入。

对 MCU 项目，建议至少保存：

timestamp, raw, filtered

如果是频率检测，还应保存：

frame_index, target_power, total_power, decision

这样后面调参数时，才能知道是滤波器真的有效，还是只是某一段数据看起来更平滑。

小结：频域滤波怎么选

算法 最适合解决 关键参数 主要风险
FIR 线性相位低通、波形保真 抽头数、系数、\(f_c\) 阶数高时延迟和计算量大
Biquad IIR 低成本低通/高通/带通 \(f_c\)、\(Q\)、系数 系数错误会不稳定，相位非线性；高阶时用 Biquad 级联
陷波 固定频率干扰 \(f_0\)、\(Q\) 频率漂移会漏掉干扰
CIC 过采样抽取 抽取率 \(R\)、级数 \(N\) 通带下垂、积分器溢出
Goertzel 单频检测 \(N\)、\(k\)、阈值 频点不对齐会泄漏
Savitzky-Golay 保留峰形的平滑 窗口长度、多项式阶数 有固定延迟，不抗尖峰

频域滤波的关键不是把算法名字记住，而是把指标说清楚：采样率、目标频率、截止频率、允许延迟、计算资源和验证方法。说清这些，滤波器才是工程设计；说不清这些，就只是换一种方式调参数。