信号预处理：肌电噪声去除的实战滤波方案
各位同学，今天我们来聊聊信号预处理中一个非常实际的问题——如何去除肌电噪声。说实话，我在医疗贴片项目里被这个噪声折磨过好几次。肌电信号本身频率范围宽，幅度变化大，跟我们要采集的生理信号经常混在一起。你想想看，如果不去掉它，后面算法再牛也是白搭。
这一节我们重点讲三种滤波器：巴特沃斯低通、移动平均滤波、中值滤波。它们各有各的脾气，用对了地方就是神器。
一、巴特沃斯低通滤波器：阶数与截止频率的博弈
巴特沃斯滤波器最大的特点就是通带内最平坦。什么意思？就是你要保留的信号频率段，它几乎不衰减。我刚开始做嵌入式滤波时，总喜欢用高阶数，觉得越陡峭越好。结果呢？相位失真严重，信号都变形了。
1. 阶数选择：不是越高越好
阶数决定了滤波器滚降的陡峭程度。阶数越高，过渡带越窄，但代价也大：
- 计算量翻倍：每增加一阶，多一个乘加运算。在Cortex-M4上跑4阶还行，8阶就开始吃力了
- 相位失真加剧：高阶滤波器会引入明显的群延迟，实时处理时你会看到信号"滞后"
- 数值稳定性变差：我记得有一次用10阶巴特沃斯，结果定点化后直接震荡了
我个人习惯：医疗贴片场景下，2阶到4阶就够了。肌电噪声主要能量在20-500Hz，我们关心的生理信号（比如心电）在0.5-100Hz。用4阶巴特沃斯，截止频率设在40Hz，效果已经很好了。
经验法则：阶数 = 所需衰减量(dB) / (6 × 过渡带宽度(octave))。但别死磕公式，实际调参时从2阶开始试。
2. 截止频率计算：采样率是关键
截止频率不是随便定的。你得先知道采样率。比如贴片ADC采样率是1000Hz，那奈奎斯特频率就是500Hz。我们要保留的信号最高频率是100Hz，那截止频率设在100Hz？别急，这里有个坑。
巴特沃斯滤波器的截止频率是-3dB点。如果你把截止频率设在100Hz，那100Hz的信号已经被衰减了3dB。我建议：截止频率 = 信号最高频率 × 1.2 ~ 1.5。比如信号最高100Hz，截止频率设在120-150Hz比较合理。
避坑指南：我曾经在项目里把截止频率设得太靠近信号频率，结果心电信号的R波峰值被削掉了。后来改成1.3倍，问题解决。
下面给个实际代码，用Python设计一个4阶巴特沃斯低通滤波器：
import numpy as np
from scipy.signal import butter, filtfilt
def design_butterworth_lowpass(fs, cutoff, order=4):
"""
fs: 采样率 (Hz)
cutoff: 截止频率 (Hz)
order: 滤波器阶数
"""
nyquist = 0.5 * fs
normal_cutoff = cutoff / nyquist
b, a = butter(order, normal_cutoff, btype='low', analog=False)
return b, a
# 示例：采样率1000Hz，截止频率120Hz，4阶
fs = 1000
cutoff = 120
b, a = design_butterworth_lowpass(fs, cutoff, order=4)
# 对信号滤波（注意用filtfilt实现零相位）
filtered_signal = filtfilt(b, a, raw_signal)
注意这里用了filtfilt而不是lfilter。为什么？因为filtfilt是双向滤波，相位延迟抵消了。但代价是不能实时处理——它需要整段数据。如果你做的是离线分析，用filtfilt；如果是实时流式处理，只能用lfilter，然后接受一点相位延迟。
二、移动平均滤波：简单但有效
移动平均滤波，说白了就是取最近N个点的平均值。它是最简单的低通滤波器，没有之一。但简单不代表没用。
移动平均的传递函数是：
y[n] = (x[n] + x[n-1] + ... + x[n-N+1]) / N它的-3dB截止频率大约是：fc ≈ 0.443 × fs / N。比如采样率1000Hz，窗口长度N=10，那截止频率约44.3Hz。
优点：计算量极小，适合嵌入式实时处理。没有乘除法，只有累加和移位（如果N是2的幂）。
缺点：阻带衰减慢，对高频噪声抑制效果一般。而且会引入群延迟 = (N-1)/2 个采样点。
我一般在贴片项目里这样用：先用移动平均做初步平滑，再用巴特沃斯做精细滤波。两者配合，效果比单独用任何一种都好。
嵌入式C代码实现：
// 移动平均滤波，窗口长度N=8（2的幂，可用移位加速）
#define N 8
int16_t buffer[N] = {0};
uint8_t index = 0;
int32_t sum = 0;
int16_t moving_average(int16_t sample) {
sum -= buffer[index];
buffer[index] = sample;
sum += buffer[index];
index = (index + 1) & (N - 1); // 循环索引
return sum >> 3; // 除以8
}
这段代码我用了很多年。注意index = (index + 1) & (N - 1)这个技巧——当N是2的幂时，用位与代替取模，效率高很多。
三、中值滤波：对付脉冲噪声的利器
中值滤波跟前面两种完全不同。它不是线性滤波器，而是排序滤波器。原理很简单：取窗口内所有值的中位数作为输出。
为什么需要它？因为肌电噪声里经常有尖峰脉冲——比如肌肉突然收缩一下。这种脉冲用低通滤波去不掉，反而会被"抹开"成一段小波动。中值滤波就擅长处理这种：它直接把异常值替换成周围正常值的中位数。
注意：中值滤波对高斯噪声效果不好。高斯噪声是连续分布的，中值滤波会保留很多噪声成分。所以它只适合去除脉冲噪声，不适合做平滑。
窗口长度怎么选？我建议：
- 3点窗口：去除孤立脉冲，对信号影响最小
- 5点窗口：去除连续2-3个脉冲，信号边缘会稍微变钝
- 7点以上：慎用。信号细节会被严重破坏
嵌入式实现中值滤波有个经典技巧——用排序网络代替完整排序。比如3点中值：
// 3点中值滤波，用比较交换实现
int16_t median_3(int16_t a, int16_t b, int16_t c) {
// 排序网络：a <= b <= c
if (a > b) { int16_t t = a; a = b; b = t; }
if (b > c) { int16_t t = b; b = c; c = t; }
if (a > b) { int16_t t = a; a = b; b = t; }
return b; // 中位数
}
这段代码只有3次比较交换，没有循环，没有内存分配。在STM32上跑，一个采样点只需要几十个时钟周期。
四、三种滤波器的实战组合策略
在实际的医疗贴片项目中，我很少只用一种滤波器。通常的组合方式是：
| 处理阶段 | 滤波器 | 作用 |
|---|---|---|
| 第一级 | 中值滤波（3点） | 去除运动伪迹和脉冲噪声 |
| 第二级 | 移动平均（N=4或8） | 初步平滑，降低高频抖动 |
| 第三级 | 巴特沃斯低通（2-4阶） | 精细滤波，保留目标信号频段 |
这个顺序是有讲究的。中值滤波放在最前面，因为脉冲噪声如果不先去掉，后面移动平均和巴特沃斯都会把它"扩散"开。移动平均放在中间，因为它计算快，能减轻后面巴特沃斯的负担。巴特沃斯放在最后，做精细的频率选择。
我的经验：在Cortex-M4上，这个三级流水线处理一个采样点大约需要5-10微秒。如果采样率是1000Hz，那CPU占用率只有0.5%-1%，完全够用。
最后提醒一句：滤波器的参数一定要根据实际信号调整。别指望一套参数打天下。我每次拿到新的贴片数据，都会先画个频谱图看看噪声分布，再决定截止频率和窗口长度。嗯，这个习惯帮我避免了很多次返工。
下一节我们会讲特征提取，到时候这些滤波后的信号就是我们的"原材料"了。