Jerry_Spike 是我给一个尖峰信号检测工具起的项目代号。Jerry 是当时开发机上的用户名节选Spike 就是信号里那些短促有力的尖峰脉冲起名的时候没想太多。结果这个工具后来陪着我处理了神经电生理峰电位、工业振动冲击波形、音频瞬态噪声三类完全不同源的数据从最初一个几百行的小脚本慢慢长成了带滤波、自适应阈值、波形筛选、事件导出的完整流程。这篇文章直接拿 Jerry_Spike 当例子把尖峰检测从原理到落地完整拆一遍。如果你手里也有一堆波形数据要在杂讯里找脉冲这里应该有一份能直接抄作业的方案。1. 项目概述与设计思路1.1 尖峰检测要解决的真实问题尖峰检测不是超过一条线就算数它指的是在一段连续采样的信号里找出持续时间很短、幅度明显突出的瞬态事件。最常见的场景就是神经电生理数据里的峰电位神经元放电表现为几十到几百微秒的快速脉冲幅度不大却携带着关键的生理信息。换到工业侧机械设备撞击、齿轮啮合冲击、轴承早期故障的信号也往往以尖峰的形式呈现。再比如音频里的爆音、瞬态噪声本质上也是尖峰。这类信号最让人头疼的地方在于目标形态不稳定。同样是放电尖峰不同的神经元放电波形差别很大同样是轴承故障不同转速下的冲击脉冲宽度也不一样。再加上背景噪声、基线漂移、传感器本身的低频晃动单纯靠肉眼在波形上框选尖峰短数据还好数据一长就完全不可行了。我在做这个工具之前也试过直接调商用软件里的固定阈值检测器效果非常勉强。阈值调高了小幅值的真实尖峰全部漏掉调低了噪声上下一抖就触发一堆假事件。真正让我下决心自己写一套检测链路的原因是我需要同时控制误检和漏检并且希望整个检测流程可以被批量运行、可以重复复现这样后续处理不同批次的数据时标准才能统一。Jerry_Spike 的定位就是一套可配置的批量尖峰检测链路输入原始波形和采样率输出一个事件时间戳列表以及每个事件对应的幅度、宽度等特征。1.2 固定阈值为什么不可靠很多刚接触尖峰检测的人第一反应是设一条固定阈值线超过阈值就标记成尖峰。这个方法在实验室理想信噪比条件下也许能看到了真实数据里基本就崩了。原因主要有三个。第一是基线漂移传感器预热、温度变化、硬件滤波都会让整段信号的零电平缓慢上下移动固定阈值跟着失效。第二是尖峰幅度差异大真实场景里既有幅度很大的强脉冲也有贴着噪声上沿的弱脉冲一条阈值线只能照顾其中一类。第三是噪声本身不是平稳的不同时间段背景噪声方差差别很大固定的绝对幅度阈值在噪声大的时段触发频繁、噪声小的时段又漏检。我最早吃这个亏是在处理一批生理信号记录时前一分钟数据安静得像小湖泊后一分钟受运动干扰噪声幅度翻了五六倍。固定阈值下前一分钟的尖峰倒是没漏后一分钟的尖峰全被淹没在误检里了。后来我把阈值改成随着局部噪声水平自适应变化问题才真正解决。这也是 Jerry_Spike 切换成自适应阈值的直接起因。1.3 工具选型与可移植性技术选型上我选了 Python搭配 NumPy 和 SciPy。原因很实际SciPy 的 signal 模块里已经有成熟的滤波器函数和峰值搜索函数不需要自己重新实现底层算法省掉一大截调试时间而且 Python 做数据分析、画图、模型验证都非常顺手处理完检测结果直接就能用 matplotlib 可视化复盘。可能有人会问为什么不用深度学习模型来识别尖峰。神经网络确实能学习复杂波形模式但代价是训练数据要足够多、标注成本高、模型在不同采样率下通用性差而且推理过程不像规则方法那样容易解释。尖峰检测本质上是一个先验很明确的问题——尖峰就是局部极大值、持续时间短、和周围背景差异大用信号处理的方式完全可以把这件事讲清楚也方便维护。后期真要叠加更复杂的分类任务直接在 Jerry_Spike 的事件输出基础上加一层分类器就可以了。还有一点是代码可移植性。选 Python 不意味着只能跑 Python核心算法拆出来后MAD 阈值计算、滤波、峰值搜索这几步用的都是非常标准的信号处理逻辑将来要移植到 C 或者嵌入式平台也不难。做工具选型的时候我会刻意提醒自己别为了酷炫把简单问题复杂化一个能稳定跑一年的工具远比一个三天两头要重新调参的模型有价值。2. 核心算法与关键参数2.1 信号预处理先过滤再检测尖峰检测的第一步不是直接找峰而是先把原始信号里不想要的成分去掉。原始波形里通常混着两类东西一类是低频基线漂移频率可能低到零点几赫兹另一类是高频噪声可能分布在几千赫兹以上。这两类成分都会严重干扰阈值计算所以预处理环节我固定做两件事高通滤波压掉漂移低通滤波压掉高频毛刺。具体实现上我用了 SciPy 的 Butterworth 滤波器阶数选二阶再通过filtfilt做零相位滤波。为什么要强调零相位因为普通滤波会引入相位延迟尖峰在滤波后的时间位置会被整体平移检测出来的时间戳就不准了。filtfilt先正向再反向滤波相位偏移互相抵消尖峰位置基本不偏。代价是计算量翻倍但对离线处理来说完全值得。滤波截止频率的选择很关键。高通截止频率一般取 20 到 30 Hz能压住呼吸、运动、热漂移等低频干扰低通截止频率要结合采样率和信号固有特征比如神经信号常见 0.1 到 5 kHz 的频带能量采样率 10 kHz 时低通设 1500 Hz 已经够用。需要注意一点滤波不是越狠越好低通设得太低会把尖峰头部削圆宽度特征失真幅度也被压低后面再好的阈值算法也救不回来。我自己习惯在滤波后立刻画一张对比图确认滤波后的尖峰形状还在、位置没动、幅度没有明显衰减再做下一步。2.2 自适应阈值MAD 的底层逻辑自适应阈值部分Jerry_Spike 用的核心统计量是 MAD即中位数绝对偏差。MAD 的定义是median(|x - median(x)|)也就是每个样本减去该区间中位数之后取绝对值再对这些绝对值取中位数。它描述的是信号波动的典型幅度但和标准差不同的是中位数本身对离群值非常不敏感。为什么这点重要因为尖峰本身就是离群值。假如用均值加标准差来做阈值几个大幅度尖峰会把标准差拉得很大阈值抬升之后反而把后续的尖峰漏掉了这是典型的被强者带偏。中位数就不一样即使区间里混进了几个非常大的尖峰中位数几乎不受影响得到的噪声水平估计稳健很多。计算最终阈值时通常把 MAD 乘以一个归一化常数 1.4826再乘以阈值因子。1.4826 来源于正态分布下 MAD 与标准差的关系乘上它之后阈值就相当于约 k 倍标准差。实际使用中 k 取 3 到 5 比较合理取 3 时灵敏度高但误检多取 5 时更保守。我在 Jerry_Spike 里默认取 4大部分数据的表现都还可以。阈值计算需要一个滑动窗口窗口太短估计出的噪声水平波动大窗口太长遇到信号状态变化时跟不上比如运动干扰突变那一段阈值会残留很久的高噪声记忆。我一般把窗口设成 0.5 到 2 秒。算法上用的是分块计算每块内部算一个 MAD 阈值块和块之间阈值会有台阶这在稀疏尖峰场景下问题不大但如果你的信号背景噪声变化非常平滑可以考虑用更细的滑动窗口加线性插值平滑阈值曲线。2.3 不应期和波形约束减少重复与假峰自适应阈值能解决大部分问题但还有一个非常实际的坑同一个实际尖峰检测时可能被标记成两个事件。因为一个尖峰通常先上升后下降上升沿超过阈值一次下降沿回落后再波动一下可能又超一次如果不加约束就会把一个事件误报成两个。解决办法是引入不应期。思想很简单两个合理的尖峰检测点之间必须间隔至少一段时间默认我取 5 到 8 毫秒。在 SciPy 的find_peaks里这对应distance参数意思是两个峰之间至少间隔多少个采样点。这个参数必须在波形形态分析之后设定因为不同信号的尖峰宽度不一样不应期要比单个尖峰的宽度稍大又不能大到把紧挨着的真实连续事件漏掉。除了不应期还要做宽度约束。真实尖峰通常持续几毫秒到几十毫秒太窄的检测点多半是高频噪声毛刺太宽的往往是基线抬升或者电极移动伪迹。在 Jerry_Spike 里我通过find_peaks的width参数限定可接受的峰宽度范围比如宽度范围设为 2 到 50 毫秒。这一步能把大量高而窄的噪声点直接过滤掉效果立竿见影。3. 实操从数据到事件序列3.1 环境准备与数据组织运行 Jerry_Spike 的环境非常朴素Python 3.8 以上就够了核心依赖只有 NumPy 和 SciPy画图建议再加一个 Matplotlib。处理的数据格式我是统一成两种要么只给一个一维数组加一个采样率要么给一个二维数组加一列时间戳。二维数组适合多通道记录每个通道是独立的检测对象不要把所有通道混在一起算阈值否则某一个通道出现大幅度强脉冲会把所有通道的阈值都抬高。数据导入这一步我习惯先把原始波形转换成 float64避免整数截断影响幅度判断。然后确认采样率的单位Jerry_Spike 所有时间相关的参数都以秒为单位内部再换算成样本数这样切换不同采样率的数据时参数不用改来改去。曾经有次我从旧数据里拿到的信号标称采样率写的是 1024实际存储的时候已经降采样到 256检测出来的宽度全部失真后来我加了数据自带采样率校验的日志才算治住这个问题。3.2 核心代码SpikeDetector 实现下面给出 Jerry_Spike 的核心类实现代码不复杂胜在结构清楚import numpy as np from scipy import signal class SpikeDetector: def __init__(self, fs10000, hp_cutoff20.0, lp_cutoff1500.0, mad_window1.0, threshold_factor4.0, refractory5e-3, width_range(2e-3, 50e-3)): self.fs fs self.hp_cutoff hp_cutoff self.lp_cutoff lp_cutoff self.mad_window int(mad_window * fs) self.threshold_factor threshold_factor self.refractory refractory self.width_range width_range self._init_filters() def _init_filters(self): nyq self.fs / 2 if self.hp_cutoff and self.hp_cutoff nyq: self.b_hp, self.a_hp signal.butter(2, self.hp_cutoff / nyq, high) else: self.b_hp self.a_hp None if self.lp_cutoff and self.lp_cutoff nyq: self.b_lp, self.a_lp signal.butter(2, self.lp_cutoff / nyq, low) else: self.b_lp self.a_lp None def preprocess(self, x): if self.b_hp is not None: x signal.filtfilt(self.b_hp, self.a_hp, x) if self.b_lp is not None: x signal.filtfilt(self.b_lp, self.a_lp, x) return x def detect(self, x): x self.preprocess(x) abs_x np.abs(x) n len(abs_x) threshold np.zeros(n) last_thr np.finfo(float).max win self.mad_window for i in range(0, n, win): seg abs_x[i:i win] if len(seg) 3: mad np.median(np.abs(seg - np.median(seg))) if mad 0: last_thr self.threshold_factor * 1.4826 * mad threshold[i:i len(seg)] last_thr min_w int(self.width_range[0] * self.fs) max_w int(self.width_range[1] * self.fs) peaks, _ signal.find_peaks( abs_x, heightthreshold, distanceint(self.refractory * self.fs), width(min_w, max_w) ) return peaks / self.fs, threshold调用方式很简单detector SpikeDetector(fs10000, hp_cutoff30, lp_cutoff1500, threshold_factor4.0, refractory5e-3) event_seconds, threshold_curve detector.detect(raw_signal)代码里find_peaks的height参数传的是数组含义是每个样本点的动态阈值峰只有在abs_x[i]大于threshold[i]时才会被保留。distance按样本数换算不应期width则用采样率把秒转换成样本区间。整体流程就是预处理、MAD 分块阈值、峰值搜索、宽度约束四步全部耦合在类里新增数据直接砸进来就能跑。3.3 用模拟信号验证检测效果写算法最怕在真实数据上翻车所以我做任何改动之前都会先用人工构造的模拟信号测一遍。模拟数据的优势是真实事件位置完全已知可以准确计算召回率和精确率。构造的思路很简单先造一个低频正弦波模拟真实信号背景叠加高斯白噪声模拟干扰再在固定的随机位置注入已知宽度和幅度的尖峰。下面是验证过程的核心代码fs 10000 t np.arange(0, 10, 1 / fs) x np.sin(2 * np.pi * 3 * t) 0.1 * np.random.randn(len(t)) true_peaks np.array([20000, 35000, 52000, 73000]) for idx in true_peaks: x[idx:idx 50] np.exp(-np.arange(50) / 8) * 2.0 detector SpikeDetector(fs, hp_cutoff30, lp_cutoff1500, threshold_factor4.0, refractory5e-3) peak_times, _ detector.detect(x)检测结果里四个注入的尖峰全部被找到没有多出额外的假事件。我后来又跑了多组不同噪声水平的模拟观察到一个规律在信噪比大于 5 的情况下MAD 自适应阈值的检测效果非常稳定信噪比降到 2 附近时弱尖峰开始出现漏检这时优先降低阈值因子到 3或者放宽宽度约束可以找回一部分漏检但误检会相应上升。建议实际使用时用模拟信号先标定一遍参数再把同一组参数套到真实数据上这样心里才有底。4. 常见问题与排查实录4.1 误检率居高不下误检多通常会集中在两种场景一是噪声不是高斯白噪声而是时不时出现的机械抖动、工频干扰或者电极伪迹这些噪声的形状和尖峰非常像二是阈值因子设得太低把噪声的局部极值也当成事件了。排查顺序我一般是先画图把原始波形、滤波后波形、阈值曲线和最终检测标记叠在一起看一眼误检点全都落在哪里。如果误检点都是又窄又高的纯随机噪声优先调高threshold_factor比如从 4 调到 5如果误检点集中在某个固定频率的振荡干扰那就是滤波没滤干净应该加一个陷波滤波器或者收紧低通截止频率。还有一个小技巧可以统计检测出来的事件宽度分布大部分真实尖峰的宽度是相对集中的从宽度直方图里能很直观看出异常点。4.2 尖峰被吞漏检问题漏检最典型的原因是两个阈值因子设太大或者预处理阶段把尖峰削得太狠。前者容易理解后者我会单独检查高通滤波器如果有明显超过 50 Hz 的截止频率很容易把宽一点的慢尖峰当成基线漂移直接滤掉低通滤波器截止太低则会把尖峰幅度压小。如果确认滤波没问题那就要考虑幅度差异大的情况。某些真实数据里强尖峰和弱尖峰的幅度可以差一个数量级这时单一阈值根本不可能同时抓住两类事件。我的处理方法是把检测分成两轮第一轮用高阈值抓强尖峰第二轮把第一轮检测出的强尖峰时间点附近的数据挖掉再用低阈值重新检测弱尖峰最后合并结果。这个方法在 Jerry_Spike 里作为可选模式实现了处理神经信号时效果很明显。4.3 基线漂移把阈值带偏分块 MAD 对基线漂移其实有不错的抗性因为高通滤波已经先把绝大部分漂移滤掉了。但有些数据源的低频漂移特别剧烈滤波后仍然留下残差这时分块阈值会在漂移上升段整体抬高导致弱尖峰被漏掉。遇到这种情况我建议分两步。第一步把高通截止频率适当提高一点比如从 20 Hz 提到 35 Hz压住漂移第二步减小 MAD 窗口比如从 1 秒缩短到 0.5 秒让阈值跟着局部噪声水平快速变化。注意窗口缩短会让阈值曲线波动更剧烈如果误检增加就配合提高阈值因子来平衡。4.4 长数据跑得慢MAD 分块计算加上filtfilt处理几十 MB 级别的长时间记录时会感觉明显变慢。就 Jerry_Spike 而言优化最有效的办法是把分块 MAD 循环向量化。因为每块之间的阈值只有块边界处跳变实际可以用reshape切块、批量计算中位数再repeat展开成完整阈值曲线速度能提升几十倍。如果数据大到单线程处理都很吃力还可以考虑用多进程按通道并行因为各通道检测完全独立。我在实际项目中处理 16 通道、每通道 2 小时的数据单通道串行大约要几分钟并行后压到十几秒完全可接受。还有一点容易被忽略filtfilt的计算量比一阶段滤波大一倍如果时间戳精度、定位的准确度要求没那么苛刻改成lfilter也能省不少时间。下面把排查问题整理成一个速查表方便实际调参时快速定位现象可能原因推荐处理误检多阈值因子过低提高 threshold_factor 到 5~6误检多噪声有周期性干扰加陷波滤波或收紧低通误检点宽度偏窄高频毛刺被当成峰调高宽度范围下限漏检弱尖峰阈值因子过高降低阈值因子到 3漏检弱尖峰高通滤波太狠降低高通截止频率阈值整体被抬高剧烈基线漂移提高高通截止频率、缩小 MAD 窗口一个尖峰两个事件不应期太短增大 refractory时间戳偏移滤波器引入相位使用零相位 filtfilt数据量大跑不动串行处理按通道并行、向量化 MAD 计算最后说一个我自己的调参习惯每次拿到新数据不会直接跑完整检测而是先取三十秒数据把原始波形、滤波后波形、阈值曲线、检测标记四条线叠在一起画一张图肉眼过一遍再动参数。尖峰检测这个活儿算法再花哨最后看的还是波形本身。Jerry_Spike 能做的其实就是把一个熟练工的眼睛功夫变成可重复的代码。把这套链路跑通之后后面换数据源、加通道、做事件分类都会轻松很多。
Jerry_Spike 是我给一个尖峰信号检测工具起的项目代号。Jerry 是当时开发机上的用户名节选Spike 就是信号里那些短促有力的尖峰脉冲起名的时候没想太多。结果这个工具后来陪着我处理了神经电生理峰电位、工业振动冲击波形、音频瞬态噪声三类完全不同源的数据从最初一个几百行的小脚本慢慢长成了带滤波、自适应阈值、波形筛选、事件导出的完整流程。这篇文章直接拿 Jerry_Spike 当例子把尖峰检测从原理到落地完整拆一遍。如果你手里也有一堆波形数据要在杂讯里找脉冲这里应该有一份能直接抄作业的方案。1. 项目概述与设计思路1.1 尖峰检测要解决的真实问题尖峰检测不是超过一条线就算数它指的是在一段连续采样的信号里找出持续时间很短、幅度明显突出的瞬态事件。最常见的场景就是神经电生理数据里的峰电位神经元放电表现为几十到几百微秒的快速脉冲幅度不大却携带着关键的生理信息。换到工业侧机械设备撞击、齿轮啮合冲击、轴承早期故障的信号也往往以尖峰的形式呈现。再比如音频里的爆音、瞬态噪声本质上也是尖峰。这类信号最让人头疼的地方在于目标形态不稳定。同样是放电尖峰不同的神经元放电波形差别很大同样是轴承故障不同转速下的冲击脉冲宽度也不一样。再加上背景噪声、基线漂移、传感器本身的低频晃动单纯靠肉眼在波形上框选尖峰短数据还好数据一长就完全不可行了。我在做这个工具之前也试过直接调商用软件里的固定阈值检测器效果非常勉强。阈值调高了小幅值的真实尖峰全部漏掉调低了噪声上下一抖就触发一堆假事件。真正让我下决心自己写一套检测链路的原因是我需要同时控制误检和漏检并且希望整个检测流程可以被批量运行、可以重复复现这样后续处理不同批次的数据时标准才能统一。Jerry_Spike 的定位就是一套可配置的批量尖峰检测链路输入原始波形和采样率输出一个事件时间戳列表以及每个事件对应的幅度、宽度等特征。1.2 固定阈值为什么不可靠很多刚接触尖峰检测的人第一反应是设一条固定阈值线超过阈值就标记成尖峰。这个方法在实验室理想信噪比条件下也许能看到了真实数据里基本就崩了。原因主要有三个。第一是基线漂移传感器预热、温度变化、硬件滤波都会让整段信号的零电平缓慢上下移动固定阈值跟着失效。第二是尖峰幅度差异大真实场景里既有幅度很大的强脉冲也有贴着噪声上沿的弱脉冲一条阈值线只能照顾其中一类。第三是噪声本身不是平稳的不同时间段背景噪声方差差别很大固定的绝对幅度阈值在噪声大的时段触发频繁、噪声小的时段又漏检。我最早吃这个亏是在处理一批生理信号记录时前一分钟数据安静得像小湖泊后一分钟受运动干扰噪声幅度翻了五六倍。固定阈值下前一分钟的尖峰倒是没漏后一分钟的尖峰全被淹没在误检里了。后来我把阈值改成随着局部噪声水平自适应变化问题才真正解决。这也是 Jerry_Spike 切换成自适应阈值的直接起因。1.3 工具选型与可移植性技术选型上我选了 Python搭配 NumPy 和 SciPy。原因很实际SciPy 的 signal 模块里已经有成熟的滤波器函数和峰值搜索函数不需要自己重新实现底层算法省掉一大截调试时间而且 Python 做数据分析、画图、模型验证都非常顺手处理完检测结果直接就能用 matplotlib 可视化复盘。可能有人会问为什么不用深度学习模型来识别尖峰。神经网络确实能学习复杂波形模式但代价是训练数据要足够多、标注成本高、模型在不同采样率下通用性差而且推理过程不像规则方法那样容易解释。尖峰检测本质上是一个先验很明确的问题——尖峰就是局部极大值、持续时间短、和周围背景差异大用信号处理的方式完全可以把这件事讲清楚也方便维护。后期真要叠加更复杂的分类任务直接在 Jerry_Spike 的事件输出基础上加一层分类器就可以了。还有一点是代码可移植性。选 Python 不意味着只能跑 Python核心算法拆出来后MAD 阈值计算、滤波、峰值搜索这几步用的都是非常标准的信号处理逻辑将来要移植到 C 或者嵌入式平台也不难。做工具选型的时候我会刻意提醒自己别为了酷炫把简单问题复杂化一个能稳定跑一年的工具远比一个三天两头要重新调参的模型有价值。2. 核心算法与关键参数2.1 信号预处理先过滤再检测尖峰检测的第一步不是直接找峰而是先把原始信号里不想要的成分去掉。原始波形里通常混着两类东西一类是低频基线漂移频率可能低到零点几赫兹另一类是高频噪声可能分布在几千赫兹以上。这两类成分都会严重干扰阈值计算所以预处理环节我固定做两件事高通滤波压掉漂移低通滤波压掉高频毛刺。具体实现上我用了 SciPy 的 Butterworth 滤波器阶数选二阶再通过filtfilt做零相位滤波。为什么要强调零相位因为普通滤波会引入相位延迟尖峰在滤波后的时间位置会被整体平移检测出来的时间戳就不准了。filtfilt先正向再反向滤波相位偏移互相抵消尖峰位置基本不偏。代价是计算量翻倍但对离线处理来说完全值得。滤波截止频率的选择很关键。高通截止频率一般取 20 到 30 Hz能压住呼吸、运动、热漂移等低频干扰低通截止频率要结合采样率和信号固有特征比如神经信号常见 0.1 到 5 kHz 的频带能量采样率 10 kHz 时低通设 1500 Hz 已经够用。需要注意一点滤波不是越狠越好低通设得太低会把尖峰头部削圆宽度特征失真幅度也被压低后面再好的阈值算法也救不回来。我自己习惯在滤波后立刻画一张对比图确认滤波后的尖峰形状还在、位置没动、幅度没有明显衰减再做下一步。2.2 自适应阈值MAD 的底层逻辑自适应阈值部分Jerry_Spike 用的核心统计量是 MAD即中位数绝对偏差。MAD 的定义是median(|x - median(x)|)也就是每个样本减去该区间中位数之后取绝对值再对这些绝对值取中位数。它描述的是信号波动的典型幅度但和标准差不同的是中位数本身对离群值非常不敏感。为什么这点重要因为尖峰本身就是离群值。假如用均值加标准差来做阈值几个大幅度尖峰会把标准差拉得很大阈值抬升之后反而把后续的尖峰漏掉了这是典型的被强者带偏。中位数就不一样即使区间里混进了几个非常大的尖峰中位数几乎不受影响得到的噪声水平估计稳健很多。计算最终阈值时通常把 MAD 乘以一个归一化常数 1.4826再乘以阈值因子。1.4826 来源于正态分布下 MAD 与标准差的关系乘上它之后阈值就相当于约 k 倍标准差。实际使用中 k 取 3 到 5 比较合理取 3 时灵敏度高但误检多取 5 时更保守。我在 Jerry_Spike 里默认取 4大部分数据的表现都还可以。阈值计算需要一个滑动窗口窗口太短估计出的噪声水平波动大窗口太长遇到信号状态变化时跟不上比如运动干扰突变那一段阈值会残留很久的高噪声记忆。我一般把窗口设成 0.5 到 2 秒。算法上用的是分块计算每块内部算一个 MAD 阈值块和块之间阈值会有台阶这在稀疏尖峰场景下问题不大但如果你的信号背景噪声变化非常平滑可以考虑用更细的滑动窗口加线性插值平滑阈值曲线。2.3 不应期和波形约束减少重复与假峰自适应阈值能解决大部分问题但还有一个非常实际的坑同一个实际尖峰检测时可能被标记成两个事件。因为一个尖峰通常先上升后下降上升沿超过阈值一次下降沿回落后再波动一下可能又超一次如果不加约束就会把一个事件误报成两个。解决办法是引入不应期。思想很简单两个合理的尖峰检测点之间必须间隔至少一段时间默认我取 5 到 8 毫秒。在 SciPy 的find_peaks里这对应distance参数意思是两个峰之间至少间隔多少个采样点。这个参数必须在波形形态分析之后设定因为不同信号的尖峰宽度不一样不应期要比单个尖峰的宽度稍大又不能大到把紧挨着的真实连续事件漏掉。除了不应期还要做宽度约束。真实尖峰通常持续几毫秒到几十毫秒太窄的检测点多半是高频噪声毛刺太宽的往往是基线抬升或者电极移动伪迹。在 Jerry_Spike 里我通过find_peaks的width参数限定可接受的峰宽度范围比如宽度范围设为 2 到 50 毫秒。这一步能把大量高而窄的噪声点直接过滤掉效果立竿见影。3. 实操从数据到事件序列3.1 环境准备与数据组织运行 Jerry_Spike 的环境非常朴素Python 3.8 以上就够了核心依赖只有 NumPy 和 SciPy画图建议再加一个 Matplotlib。处理的数据格式我是统一成两种要么只给一个一维数组加一个采样率要么给一个二维数组加一列时间戳。二维数组适合多通道记录每个通道是独立的检测对象不要把所有通道混在一起算阈值否则某一个通道出现大幅度强脉冲会把所有通道的阈值都抬高。数据导入这一步我习惯先把原始波形转换成 float64避免整数截断影响幅度判断。然后确认采样率的单位Jerry_Spike 所有时间相关的参数都以秒为单位内部再换算成样本数这样切换不同采样率的数据时参数不用改来改去。曾经有次我从旧数据里拿到的信号标称采样率写的是 1024实际存储的时候已经降采样到 256检测出来的宽度全部失真后来我加了数据自带采样率校验的日志才算治住这个问题。3.2 核心代码SpikeDetector 实现下面给出 Jerry_Spike 的核心类实现代码不复杂胜在结构清楚import numpy as np from scipy import signal class SpikeDetector: def __init__(self, fs10000, hp_cutoff20.0, lp_cutoff1500.0, mad_window1.0, threshold_factor4.0, refractory5e-3, width_range(2e-3, 50e-3)): self.fs fs self.hp_cutoff hp_cutoff self.lp_cutoff lp_cutoff self.mad_window int(mad_window * fs) self.threshold_factor threshold_factor self.refractory refractory self.width_range width_range self._init_filters() def _init_filters(self): nyq self.fs / 2 if self.hp_cutoff and self.hp_cutoff nyq: self.b_hp, self.a_hp signal.butter(2, self.hp_cutoff / nyq, high) else: self.b_hp self.a_hp None if self.lp_cutoff and self.lp_cutoff nyq: self.b_lp, self.a_lp signal.butter(2, self.lp_cutoff / nyq, low) else: self.b_lp self.a_lp None def preprocess(self, x): if self.b_hp is not None: x signal.filtfilt(self.b_hp, self.a_hp, x) if self.b_lp is not None: x signal.filtfilt(self.b_lp, self.a_lp, x) return x def detect(self, x): x self.preprocess(x) abs_x np.abs(x) n len(abs_x) threshold np.zeros(n) last_thr np.finfo(float).max win self.mad_window for i in range(0, n, win): seg abs_x[i:i win] if len(seg) 3: mad np.median(np.abs(seg - np.median(seg))) if mad 0: last_thr self.threshold_factor * 1.4826 * mad threshold[i:i len(seg)] last_thr min_w int(self.width_range[0] * self.fs) max_w int(self.width_range[1] * self.fs) peaks, _ signal.find_peaks( abs_x, heightthreshold, distanceint(self.refractory * self.fs), width(min_w, max_w) ) return peaks / self.fs, threshold调用方式很简单detector SpikeDetector(fs10000, hp_cutoff30, lp_cutoff1500, threshold_factor4.0, refractory5e-3) event_seconds, threshold_curve detector.detect(raw_signal)代码里find_peaks的height参数传的是数组含义是每个样本点的动态阈值峰只有在abs_x[i]大于threshold[i]时才会被保留。distance按样本数换算不应期width则用采样率把秒转换成样本区间。整体流程就是预处理、MAD 分块阈值、峰值搜索、宽度约束四步全部耦合在类里新增数据直接砸进来就能跑。3.3 用模拟信号验证检测效果写算法最怕在真实数据上翻车所以我做任何改动之前都会先用人工构造的模拟信号测一遍。模拟数据的优势是真实事件位置完全已知可以准确计算召回率和精确率。构造的思路很简单先造一个低频正弦波模拟真实信号背景叠加高斯白噪声模拟干扰再在固定的随机位置注入已知宽度和幅度的尖峰。下面是验证过程的核心代码fs 10000 t np.arange(0, 10, 1 / fs) x np.sin(2 * np.pi * 3 * t) 0.1 * np.random.randn(len(t)) true_peaks np.array([20000, 35000, 52000, 73000]) for idx in true_peaks: x[idx:idx 50] np.exp(-np.arange(50) / 8) * 2.0 detector SpikeDetector(fs, hp_cutoff30, lp_cutoff1500, threshold_factor4.0, refractory5e-3) peak_times, _ detector.detect(x)检测结果里四个注入的尖峰全部被找到没有多出额外的假事件。我后来又跑了多组不同噪声水平的模拟观察到一个规律在信噪比大于 5 的情况下MAD 自适应阈值的检测效果非常稳定信噪比降到 2 附近时弱尖峰开始出现漏检这时优先降低阈值因子到 3或者放宽宽度约束可以找回一部分漏检但误检会相应上升。建议实际使用时用模拟信号先标定一遍参数再把同一组参数套到真实数据上这样心里才有底。4. 常见问题与排查实录4.1 误检率居高不下误检多通常会集中在两种场景一是噪声不是高斯白噪声而是时不时出现的机械抖动、工频干扰或者电极伪迹这些噪声的形状和尖峰非常像二是阈值因子设得太低把噪声的局部极值也当成事件了。排查顺序我一般是先画图把原始波形、滤波后波形、阈值曲线和最终检测标记叠在一起看一眼误检点全都落在哪里。如果误检点都是又窄又高的纯随机噪声优先调高threshold_factor比如从 4 调到 5如果误检点集中在某个固定频率的振荡干扰那就是滤波没滤干净应该加一个陷波滤波器或者收紧低通截止频率。还有一个小技巧可以统计检测出来的事件宽度分布大部分真实尖峰的宽度是相对集中的从宽度直方图里能很直观看出异常点。4.2 尖峰被吞漏检问题漏检最典型的原因是两个阈值因子设太大或者预处理阶段把尖峰削得太狠。前者容易理解后者我会单独检查高通滤波器如果有明显超过 50 Hz 的截止频率很容易把宽一点的慢尖峰当成基线漂移直接滤掉低通滤波器截止太低则会把尖峰幅度压小。如果确认滤波没问题那就要考虑幅度差异大的情况。某些真实数据里强尖峰和弱尖峰的幅度可以差一个数量级这时单一阈值根本不可能同时抓住两类事件。我的处理方法是把检测分成两轮第一轮用高阈值抓强尖峰第二轮把第一轮检测出的强尖峰时间点附近的数据挖掉再用低阈值重新检测弱尖峰最后合并结果。这个方法在 Jerry_Spike 里作为可选模式实现了处理神经信号时效果很明显。4.3 基线漂移把阈值带偏分块 MAD 对基线漂移其实有不错的抗性因为高通滤波已经先把绝大部分漂移滤掉了。但有些数据源的低频漂移特别剧烈滤波后仍然留下残差这时分块阈值会在漂移上升段整体抬高导致弱尖峰被漏掉。遇到这种情况我建议分两步。第一步把高通截止频率适当提高一点比如从 20 Hz 提到 35 Hz压住漂移第二步减小 MAD 窗口比如从 1 秒缩短到 0.5 秒让阈值跟着局部噪声水平快速变化。注意窗口缩短会让阈值曲线波动更剧烈如果误检增加就配合提高阈值因子来平衡。4.4 长数据跑得慢MAD 分块计算加上filtfilt处理几十 MB 级别的长时间记录时会感觉明显变慢。就 Jerry_Spike 而言优化最有效的办法是把分块 MAD 循环向量化。因为每块之间的阈值只有块边界处跳变实际可以用reshape切块、批量计算中位数再repeat展开成完整阈值曲线速度能提升几十倍。如果数据大到单线程处理都很吃力还可以考虑用多进程按通道并行因为各通道检测完全独立。我在实际项目中处理 16 通道、每通道 2 小时的数据单通道串行大约要几分钟并行后压到十几秒完全可接受。还有一点容易被忽略filtfilt的计算量比一阶段滤波大一倍如果时间戳精度、定位的准确度要求没那么苛刻改成lfilter也能省不少时间。下面把排查问题整理成一个速查表方便实际调参时快速定位现象可能原因推荐处理误检多阈值因子过低提高 threshold_factor 到 5~6误检多噪声有周期性干扰加陷波滤波或收紧低通误检点宽度偏窄高频毛刺被当成峰调高宽度范围下限漏检弱尖峰阈值因子过高降低阈值因子到 3漏检弱尖峰高通滤波太狠降低高通截止频率阈值整体被抬高剧烈基线漂移提高高通截止频率、缩小 MAD 窗口一个尖峰两个事件不应期太短增大 refractory时间戳偏移滤波器引入相位使用零相位 filtfilt数据量大跑不动串行处理按通道并行、向量化 MAD 计算最后说一个我自己的调参习惯每次拿到新数据不会直接跑完整检测而是先取三十秒数据把原始波形、滤波后波形、阈值曲线、检测标记四条线叠在一起画一张图肉眼过一遍再动参数。尖峰检测这个活儿算法再花哨最后看的还是波形本身。Jerry_Spike 能做的其实就是把一个熟练工的眼睛功夫变成可重复的代码。把这套链路跑通之后后面换数据源、加通道、做事件分类都会轻松很多。
