返回系列

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 结构
使用双向 SOSsosfiltfilt零相位,更适合离线复盘不能用于实时流式数据

项目选择 SOS 是一个工程稳健性选择。滤波器阶数、截止频率和采样率组合不当时,普通 b/a 系数可能更容易出现数值问题;SOS 把高阶滤波器拆成多个二阶节,稳定性更好。

signal.butter()

signal.butter() 的作用是设计 Butterworth 滤波器。它不直接滤波,而是返回滤波器系数。项目调用时传了四个关键参数:

参数项目值说明
Norder滤波器阶数,阶数越高过渡越陡
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 范围,在不同采样率下对应不同的归一化频率。

BPM 边界 6 到 40 BPM 换算 Hz 0.10 到 0.67 Hz 归一化 除以 Nyquist SOS 系数 signal.butter butter 负责“设计滤波器”,真正滤波由 sosfilt 或 sosfiltfilt 完成

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 越大,平滑越强,细节越少。

滤波后仍有细小抖动 Gaussian 平滑后的趋势 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 个新滤波点输出。