波形发生识别系统开发日志 —— FIR 数字滤波

这篇写采集链路上那一级 FIR 低通. 怎么定参数, 代码怎么写, 效果怎么验

工程叫 Wave2, 滤波相关代码在 Core/Src/fir.c 和 Core/Inc/fir.h

为什么要加这一级

F401 的识别链路原本是裸的

1
ADC 采 500 点 -> 提取 5 个特征 -> 顺序判定 -> 多数表决 -> 上屏

原始码值直接喂给 WaveId_Extract(). PC 上跑模型时这样没问题, 因为模拟信号是干净数学生成的.

接上实物之后就露出问题了. 链路上有两个噪声源

  • G431 的 DAC 输出端没有模拟滤波, 方波和锯齿的陡沿后面拖着明显振铃
  • 两板之间连的是一根普通杜邦线, 没有屏蔽

振铃对特征的影响很直接. 方波的一个陡沿过去, 后面跟着一串高频抖动, 这些点的相邻差超过了 WV_FLAT_FRAC（3% Vpp）的判据, 被判成”在动”. flat 因此被拉低, 波形看起来就没那么”方”了.

锯齿波的情况更糟, 因为它只有一个陡沿, 振铃全压在那一个沿上.

所以加了一级低通, 把边沿的高频毛刺磨掉, 让特征提取拿到的是一份干净的数据.

设计约束

截止频率不能随便定, 有两条硬约束夹着.

基波必须无损

采集是 100 kHz, 一帧 500 点, 信号基频 1 kHz.

1
2
每周期点数 = 100000 / 1000 = 100 点
一帧 500 点 = 正好 5 个完整周期

如果低通把 1 kHz 削掉一部分, 后面所有特征算的都是假数据. 所以 1 kHz 必须落在通带里, 而且衰减要小到可以忽略.

调幅波的载波要保住

这条更苛刻. 调幅波在 G431 那边的生成方式是

1
2
/* Wave1/Core/Src/wavegen.c —— 调幅波 */
s = sinf(wc) * (1.0f + m * cosf(we)) / (1.0f + m);

载波 10 kHz, 包络 1 kHz, 调制深度 m = 4/5.

10 kHz 在 100 kHz 采样率下一周期只有 10 个点. 这个密度下, 任何靠近 10 kHz 的衰减都会直接削掉峰值, 让 env_ratio 失真, 调幅波就认不出来了.

实测不同截止频率对调幅波 Vpp 的影响

1
2
3
4
taps  fc       AM Vpp 偏差
9     15 kHz   -19.4%      载波被削烂
9     25 kHz    -4.6%
5     40 kHz    -2.3%      可以接受

所以截止频率压不下去, 40 kHz 附近是比较合适的位置.

为什么选 FIR 不选 IIR

IIR 的相位是非线性的, 不同频率的延迟不一样. 波形过一遍 IIR, 形状会被”扭”——方波的上升沿和下降沿会呈现不对称的过冲.

这个作品识别看的就是形状, 相位失真会直接污染 flat 和 peak_share.

FIR 可以做成严格线性相位, 条件很简单: 系数对称.

1
2
3
h[0] == h[N-1]
h[1] == h[N-2]
...

系数对称时, 所有频率的延迟都等于 (N-1)/2 个采样点, 波形只是整体平移, 形状一点不变.

代价是同样抑制强度下阶数要高不少, 但在这个尺度上算力不是瓶颈, 后面具体算.

Gibbs 现象的权衡

阶数不能无脑往上堆, 这里有个反直觉的地方.

方波理论上是无限带宽的. 任何低通截断它的高频分量, 都会在边沿两侧留下过冲, 这就是 Gibbs 现象. 阶数越高, 过渡带越陡, 过冲越明显.

实测 15 阶、截止 15 kHz（30 kHz 处抑制 -57 dB, 理论上很强）的结果

1
2
3
4
5
6
SINE        -0.11%
SQUARE      +8.83%   <- 方波 Vpp 被抬高
TRIANGLE    -2.17%
SAW         +1.04%
STAIR       +5.56%
AM         -12.69%

方波 Vpp 被抬高了近 9%. 而 Vpp 在这个作品里是承重结构

  • WV_MIN_VPP 判断信号是否存在
  • flat 的分子分母都以 Vpp 为基准
  • peak_share 的直方图分桶按 Vpp 划分

Vpp 一变, 整套特征就跟着歪.

扫了一遍阶数和截止频率的组合

1
2
3
4
5
6
7
8
taps  fc(kHz)   最大 Vpp 偏差
3     20         2.1%
5     40         2.5%     <- 选这个
9     15         0.8%     但调幅波被削 19.4%
9     30        11.5%
15    15         7.8%
15    35        18.9%
21    25        16.3%

规律很清楚: 阶数越高过冲越大. 这不是参数没调好, 是 Gibbs 现象的必然结果.

最终参数

1
2
3
4
/* Core/Src/fir.c */
static const float s_fir_taps[FIR_TAPS] = {
  -0.012383558f,  0.103321705f,  0.818123706f,  0.103321705f, -0.012383558f,
};

5 阶, 截止 40 kHz, Hamming 窗. 系数严格对称, 直流增益精确为 1（系数和 = 1.000000000000）.

频率响应

1
2
3
4
5
6
 1 kHz   -0.00 dB    基波无损
 5 kHz   -0.05 dB
10 kHz   -0.20 dB    调幅波载波, 几乎无损失
20 kHz   -0.90 dB
30 kHz   -2.22 dB
40 kHz   -3.83 dB

六种波形的 Vpp 保真

1
2
3
4
5
6
SINE        -0.0%
SQUARE      +2.5%
TRIANGLE    -0.5%
SAW         +0.5%
STAIR       +2.5%
AM          -2.3%

全部在 ±2.5% 以内.

这个指标值得单独强调一下. 滤波器的意义是”让波形变干净”, 但前提是不能把波形本身改掉. Vpp 偏差全部压在 ±2.5% 以内, 说明这一级滤波只动了高频毛刺, 没碰波形的骨架. 后面的特征提取才能放心用.

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
/* Core/Src/fir.c */
#define FIR_TAPS        5U
#define FIR_GROUP_DELAY (FIR_TAPS / 2U)   /* = 2 */

void FIR_LowPass(uint16_t *buf, uint32_t len)
{
  uint32_t i;

  /* ！！必须从后往前算！！ */
  for (i = len; i-- > 0U; )
  {
    float acc = 0.0f;
    uint32_t k;

    for (k = 0U; k < FIR_TAPS; k++)
    {
      int32_t idx = (int32_t)i + (int32_t)k - (int32_t)FIR_GROUP_DELAY;

      /* 边界用重复首尾点, 不用补零 */
      if (idx < 0)              { idx = 0; }
      if (idx >= (int32_t)len)  { idx = (int32_t)len - 1; }

      acc += (float)buf[idx] * s_fir_taps[k];
    }

    if (acc < 0.0f)         { acc = 0.0f; }
    else if (acc > 4095.0f) { acc = 4095.0f; }

    buf[i] = (uint16_t)(acc + 0.5f);
  }
}

三个实现细节都值得说.

就地逆序迭代

这是就地滤波的正确写法.

如果从前往后算, 第 i 个点的输出会立刻覆盖掉原始数据, 后面算 i+1 时窗口里读到的就是被污染过的输入. 那样实际实现的就不是 FIR, 而是一个奇怪的 IIR, 结果完全不可控.

倒着算就没这个问题: 算第 i 点时, 窗口覆盖 [i-2, i+2], 而这些位置的值要么还是原始数据, 要么是已经算好但不再被读取的位置.

边界重复首尾点

帧头帧尾的窗口会越界. 两种处理方式:

补零的话, 帧头第一个点会被当成”从 0 跳变到实际值”, 制造出一个假的半 Vpp 台阶. flat 会把这个台阶数成”运动点”, 直接污染特征.

重复首尾点没有真实的延拓信息, 但至少不会凭空造出一个突变. 对一帧 500 点来说, 只有头尾各 2 个点受影响, 占比 0.8%, 可以忽略.

浮点累加后夹紧

acc 用 float 累加, 最后夹到 [0, 4095] 再转回 uint16_t.

不夹的话, 负的累加结果转无符号会回绕成一个巨大的值（比如 -5 变成 65531）, 画到屏幕上就是一根戳出天际的线.

加 0.5f 是为了四舍五入, 直接截断会引入固定的 -0.5 码值偏移, 累积起来会啃掉一点 Vpp.

群延迟

5 阶 FIR 的群延迟是 (5-1)/2 = 2 个采样点.

1
2 / 100 kHz = 20 微秒

相对 1 ms 的信号周期, 这个延迟是 2%, 而且它是同一帧内的整体平移, 不影响波形形状, 也不影响特征提取（特征都是相对量, 和绝对时间无关）.

算力

1
500 点 x 5 抽头 = 2500 次乘加

F401 是 Cortex-M4 @84 MHz, 带单周期浮点单元.

1
2
2500 x 3 周期/抽头 = 7500 周期
7500 / 84 MHz = 89 微秒

采集周期是 5 ms, 所以占 1.8%.

对比一下旁边那个开销: I2C 整屏刷新一次要 10 ms. 滤波连刷屏的 1% 都不到, 完全不是瓶颈.

居中对齐、逆序迭代、边界判断这些开销都算进去了, 即使按 4 周期/抽头估也只有 0.6 ms, 占 12%, 还有大量余量.

接进识别链路

1
2
3
4
5
6
7
/* Core/Src/app.c —— app_process_frame() */
  ADC_GetFrame(s_frame);          /* 拷出 500 点 */

  /* 数字滤波：在识别之前把边沿毛刺压掉 */
  FIR_LowPass(s_frame, WAVEID_FRAME_LEN);

  WaveId_Classify(s_frame, &raw); /* ★ 识别入口 */

位置很关键, 必须卡在 ADC_GetFrame 之后、WaveId_Classify 之前.

这样做的效果是”看”和”认”用的是同一份数据: 屏幕画的是滤波后的波形, 特征提取算的也是滤波后的波形. 屏幕上看到什么, 算法就用什么.

自检路径也走了同一套滤波, 保证自检结果和实时识别对得上.

效果: 波形变干净了

用串口接 G431 看原始波形

想直接看到滤波前后的差别, 最方便的办法是把 G431 的串口接到电脑, 用 VOFA+ 看它输出的原始波形.

接线

1
2
G431 PC4  --(USB-TTL 的 RXD)-->  电脑
G431 GND  --(USB-TTL 的 GND)-->  电脑

参数

1
2
3
4
5
波特率   921600
数据位   8
校验     无
停止位   1
协议     justfloat

G431 的发送格式是二进制 justfloat, 每个点 8 字节: 4 字节小端 float 加 4 字节帧尾 00 00 80 7F. 一帧发 200 点, 也就是 1600 字节.

注意 USB-TTL 模块的 TXD 要悬空, 只接 RXD 和 GND. G431 的 PC4 是它的发送脚, 接错成对接会两边对着推

打开 VOFA+ 之后, 切波形看得很清楚:

  • 方波和锯齿的陡沿上会看到毛刺和过冲
  • 三角波和正弦相对干净
  • 调幅波的载波上有明显的高频抖动

这就是 F401 收到的东西.

滤波之后

F401 侧的波形画在 OLED 上, 直接就能和上面串口看到的原始波形对比:

  • 方波边沿的高频毛刺被磨平, 沿变得更”利落”, 但过渡带宽度只有几个采样点, 肉眼看不出圆角
  • 锯齿的陡降沿干净了, 不再有沿后抖动
  • 正弦和三角几乎看不出变化（它们本来就没有陡沿, 高频成分少）
  • 调幅波的载波保留完整, 包络仍然清楚

关键在于形状没有被改: 方波还是方的, 锯齿还是斜的, 阶梯的台阶还是台阶. 滤波只带走了边沿上那些不该有的高频抖动.

实测结果

六种波形在板子上全部能认, 偶尔有个别帧抖动, 靠 5 帧多数表决抹平.

板子上的观感和 PC 模型不完全一样, 因为模型里只有噪声, 而实物链路上还有 DAC 的非线性和线路上的容性负载, 这些在模型里模拟不出来.

早期在 PC 上做过一版纯数学的仿真, 里面用理想方波加人为注入的振铃, 结论偏保守. 实物的真实波形没有仿真里那么极端, 滤波后的实际效果比那个仿真好. 以板子上跑出来的为准.

关键参数

项目 值 出处
阶数 5 FIR_TAPS
截止频率 40 kHz 设计值
窗函数 Hamming —
群延迟 2 采样点（20 µs） FIR_GROUP_DELAY
系数和 1.000000000000 归一化
单帧耗时 ~89 µs 占采集周期 1.8%
开关 APP_USE_FIR app.c

加进 Keil 工程

fir.c 是后加的文件, 需要手动拖进工程

  1. 左侧 Project 窗口, 右键 app.c 所在的那个组
  2. Add Existing Files to Group …
  3. 选 Core\Src\fir.c, Add, Close
  4. 重新编译

头文件不用管, Core\Inc 已经在 Include Paths 里了.

漏了这步会报 Undefined symbol FIR_LowPass.

todo