SciPy 滤波函数详解:butter、sosfilt、sosfiltfilt 与 gaussian_filter1d
把 RespiraScope 中实际使用的第三方滤波函数逐个讲清楚:它们接收什么参数、输出什么结果、为什么实时和离线要用不同函数。
背景问题
很多工程项目里“使用 SciPy 滤波”只写成一句话,但真正调参时必须知道第三方函数在做什么。RespiraScope 当前的滤波核心在 filter_strategies.py:
sos = signal.butter(order, [low_wn, high_wn], btype="bandpass", output="sos")
if realtime:
return signal.sosfilt(sos, data)
return signal.sosfiltfilt(sos, data)
随后项目会调用:
process_data = gaussian_smooth(process_data, sigma=config.gaussian_sigma)
这一篇专门解释这些第三方函数的设计意义。
设计方案对比
| 方案 | 实现方式 | 优点 | 问题 |
|---|---|---|---|
| 手写滤波器 | 自己实现差分方程 | 可控性最高 | 容易写错,数值稳定性难保证 |
使用 lfilter | 传入 b/a 系数滤波 | 传统接口,资料多 | 高阶滤波时数值稳定性不如 SOS |
| 使用 SOS 滤波 | butter(..., output="sos") + sosfilt | 数值更稳,SciPy 推荐 | 需要理解 SOS 结构 |
| 使用双向 SOS | sosfiltfilt | 零相位,更适合离线复盘 | 不能用于实时流式数据 |
项目选择 SOS 是一个工程稳健性选择。滤波器阶数、截止频率和采样率组合不当时,普通 b/a 系数可能更容易出现数值问题;SOS 把高阶滤波器拆成多个二阶节,稳定性更好。
signal.butter()
signal.butter() 的作用是设计 Butterworth 滤波器。它不直接滤波,而是返回滤波器系数。项目调用时传了四个关键参数:
| 参数 | 项目值 | 说明 |
|---|---|---|
N | order | 滤波器阶数,阶数越高过渡越陡 |
Wn | [low_wn, high_wn] | 归一化截止频率,范围相对 Nyquist 频率 |
btype | "bandpass" | 带通滤波,只保留一个频段 |
output | "sos" | 输出二阶节系数 |
low_wn 和 high_wn 不是直接的 Hz,而是相对 Nyquist 频率的归一化值:
nyquist = sampling_rate / 2
low_wn = cutoff_low / nyquist
high_wn = cutoff_high / nyquist
当采样率是 50 Hz 时,Nyquist 频率是 25 Hz。呼吸频段 0.10 到 0.67 Hz,在归一化后就是很低的一段。这也说明为什么采样率必须参与滤波器设计:同样的 BPM 范围,在不同采样率下对应不同的归一化频率。
signal.sosfilt()
sosfilt(sos, data) 使用 SOS 系数对数据进行因果滤波。因果的意思是:输出当前点时只能使用当前和过去的数据。这正适合实时流式系统。
RespiraScope 在实时模式下使用 sosfilt,原因是 Monitor 页面要尽快看到新点。代价是相位延迟,尤其在峰谷位置上会表现为“标记稍晚”。项目不把这种延迟当 bug,而是通过 Record End 后的离线复算来获得更适合复盘的结果。
适合 sosfilt 的场景:
- 实时显示。
- 低延迟优先。
- 数据不断到来,不能等待完整片段。
- 允许峰谷位置有轻微延迟。
signal.sosfiltfilt()
sosfiltfilt(sos, data) 会先正向滤波,再反向滤波。这样可以抵消相位延迟,得到零相位滤波效果。它更适合离线数据,但不能用于真正实时的当前点输出,因为反向滤波需要知道未来数据。
项目在 OfflineFilterStrategy 中使用它。Record End 后,系统已经拿到了完整记录片段和前后辅助点,此时再追求低延迟没有意义,应该优先波形质量和峰谷位置。
| 函数 | 是否可实时 | 相位延迟 | 项目用途 |
|---|---|---|---|
sosfilt | 可以 | 有 | Monitor 实时波形 |
sosfiltfilt | 不可以 | 基本抵消 | Record 离线复盘 |
gaussian_filter1d()
带通滤波后,波形可能仍有细小毛刺。gaussian_filter1d(data, sigma) 会用一维高斯核做平滑。sigma 越大,平滑越强,细节越少。
RespiraScope 把 Gaussian 平滑放在 Butterworth 之后,是因为它只想对已经限制频段后的波形做细节修整,而不是用高斯平滑替代频段选择。
NumPy 在这里做什么
项目没有把 NumPy 当复杂矩阵库使用,而是用它做基础数值处理:
| 函数 | 用法 | 为什么需要 |
|---|---|---|
np.median(data) | 计算原始数据中位数 | restore_baseline 时恢复显示基线 |
np.round(data, 1) | 保留 1 位小数 | 降低输出噪声和 JSON 体积 |
np.asarray(data) | 转成数组 | 方便 find_peaks 和向量计算 |
np.percentile(data, 95) | 估计高分位 | 峰谷自适应阈值使用信号跨度 |
np.diff(data) | 计算相邻差分 | 估计噪声和 BPM 间隔 |
代价和边界
SciPy 函数成熟,但不会替你判断参数是否有医学意义。比如 high_bpm 设置过高,会把较快噪声也纳入有效频段;sigma 设置过高,会让浅呼吸峰谷变得不明显;order 设置更高,过渡更陡,但也可能带来振铃或更强的边界效应。
离线 sosfiltfilt 还要求数据足够长。项目用 FILTER_STARTUP_DELAY = 300 防止短片段进入滤波,这不是随意限制,而是为了避免数据太短导致滤波结果没有解释价值。
关键函数表
| 函数 | 参数重点 | 输出 | 项目位置 |
|---|---|---|---|
signal.butter() | order、Wn、btype、output="sos" | SOS 滤波器系数 | butterworth_filter() |
signal.sosfilt() | sos、data | 实时滤波数组 | RealtimeFilterStrategy |
signal.sosfiltfilt() | sos、data | 零相位离线滤波数组 | OfflineFilterStrategy |
gaussian_filter1d() | data、sigma | 平滑后的数组 | gaussian_smooth() |
np.median() | 原始数据窗口 | 基线值 | _filter_common() |
np.round() | 滤波结果 | 保留 1 位小数的输出 | _filter_common() |
读源码建议
重点读 filter_strategies.py。先看 butterworth_filter(),确认它如何把业务参数变成 SciPy 参数;再看 FilterStrategy._filter_common(),确认滤波链路的固定顺序。最后回到 signal_processor.py,看实时模式如何只取最后 step_size 个新滤波点输出。