自建信号处理代码库：从预处理到故障诊断的实战指南

做信号处理的人早晚会走到这一步手里的测试数据越来越多、场景越来越杂你不再满足于在软件里点按钮出图更不想每次分析都临时翻开工具库文档凑一段滤波代码。我大概是两年前开始整理自己的一套信号处理与分析代码库从最底层的去均值、去趋势到带通滤波器设计、FFT频谱分析再到时域/频域特征提取和故障诊断阈值判定全部收敛成能复用的函数。今天这篇文章就是这套信号处理分析代码的功能说明我会把每一段关键代码到底在解决什么问题、为什么这么设计讲清楚也会把实际运行中踩过的那些文档里不会写的坑一并说出来。如果你刚入门信号处理看这篇文章能获得一条完整的代码组织思路如果已经有经验建议直接跳到第五章和第六章那里是我在工程化改造和故障诊断过程中最真实的教训也是花了不少时间换回来的干货。1. 为什么自己搭一套信号处理代码库而不是随手调现成接口1.1 现成滤波函数很强但它不替你思考很多做信号处理的同行包括早年的我喜欢直接调 scipy.signal.butter 然后 filtfilt 一把梭。这样确实能出图但只能出像样的图。真正的问题从来不在接口长什么样而在于你的截止频率从哪里来带通上下限是多少为什么选4阶而不是8阶滤波后信号边界那两段振荡是什么引起的、要不要处理我最近碰到的例子是这样的一台设备的振动传感器数据里混着明显的直流漂移和零点几赫兹的慢变趋势直接做FFT低频段能量奇高根本看不出轴承故障特征频率。很多人这时候会把频谱放大去看高频段但低频漂移会通过频谱泄漏污染整个频段这是一种隐蔽的错误。如果不把数据预处理当回事后面的分析全是空中楼阁。所以第一个结论是现成库函数只是零件自己搭代码库是为了建立一条可控的分析流程。我需要在每一步都能说清楚这一步在干什么、为什么需要、参数怎么定的而不是把一个黑盒丢给数据。1.2 代码库的整体架构与模块划分这套代码库我按数据流向分成了几个模块目录结构大概是这样的signal_toolkit/ ├── io.py # 数据读取统一格式 ├── preprocess.py # 去均值、去趋势、坏值处理、滤波 ├── spectrum.py # 窗函数、FFT、功率谱、峰值检测 ├── features.py # 时域/频域特征提取 ├── diagnosis.py # 滑动窗口特征序列、阈值判定 └── utils.py # 分帧、结果保存等公共工具io统一数据读取把csv、tdms、mat等不同格式归一成numpy数组同时补上采样率等元信息。preprocess去均值、去趋势、坏值处理、带通滤波。spectrum窗函数、FFT、功率谱、峰值检测。features时域特征RMS、峰值因子、峭度、频域特征重心频率、边带能量。diagnosis基于滑动窗口的特征序列阈值判定与报警去抖。utils公共工具比如分帧函数、保存结果。分层的好处有三个第一每一层都能单独测试比如我可以只测滤波模块而不用管诊断逻辑第二遇到新数据场景大多数时候只是调整参数不是重写代码第三给同事接手时接口清晰不需要读懂每个for循环。一个典型调用链路是这样的data, fs io.load_signal(test.csv) clean preprocess.clean_signal(data, fs) freqs, amp spectrum.compute_spectrum(clean, fs) feats features.extract_features(clean, fs) result diagnosis.run_diagnosis(clean, fs)这也是我在新项目里最快验证数据可用性的路径先把整条链路跑通再回头细化算法。2. 预处理模块的代码设计去均值、去趋势与滤波器的边界处理2.1 去均值和去趋势看似简单坑都在信号里先看第一个函数import numpy as np from scipy import signal def remove_offset(x): 去掉直流分量 return x - np.mean(x)去均值看似简单但它影响巨大FFT的0Hz频点对应信号的直流分量不去掉0Hz附近会被堆起一座高塔低频细节全部被掩盖如果后续还要做包络谱分析直流分量也会干扰包络解调。再往下是去趋势def remove_trend(x, kindlinear): 去掉线性或多项式趋势 if kind linear: return signal.detrend(x, typelinear) # 多项式趋势一般不超过3阶 t np.arange(len(x), dtypenp.float64) coeffs np.polyfit(t, x, deg3) trend np.polyval(coeffs, t) return x - trend实际采集的振动数据经常带有温漂或线缆扰动引起的趋势。我平时用线性去趋势最多因为它不会过度扭曲信号形状。什么时候用多项式我的判断标准是先做一次FFT如果低频段还是明显抬升且原始曲线肉眼可见缓慢弯曲才尝试3阶多项式。阶数太高容易把真实低频成分一起拟掉所以一般不超过3阶。别忘了坏值处理。传感器偶尔会给出NaN或者一个明显的毛刺。我的习惯是先检查NaN占比再用中值滤波处理尖峰def clean_spikes(x, kernel3): 中值滤波去除孤立尖峰kernel必须是奇数 x np.where(np.isfinite(x), x, np.median(x)) return signal.medfilt(x, kernel_sizekernel)这里有个实战经验中值滤波窗口不能开大。窗口大了真实的瞬时冲击比如轴承滚珠敲击外圈也会被当成毛刺滤掉。我一般取3到5个点宁可多留一点尖峰也不牺牲真实冲击成分。2.2 带通滤波器的设计与零相位处理振动分析里带通滤波几乎是标配它用来限定分析频带把不关心的低频和高频噪声隔离开。基础的滤波器函数def bandpass_filter(x, fs, low, high, order4): 带通滤波返回滤波后的信号 b, a signal.butter(order, [low, high], btypebandpass, fsfs) y signal.filtfilt(b, a, x, padtypeodd, padlenNone) return y截止频率 low 和 high 怎么定通常依据对象的工作转速和故障特征频率范围。比如一台转速30Hz左右的旋转设备轴承外圈故障特征频率往往在100到200Hz附近那我就会把带通设为80到300Hz尽量保留特征频带去掉转频低频大能量和高频噪声。为什么默认4阶而不是8阶高阶滤波器衰减更陡但相位失真更大数值稳定性也更差。filtfilt是零相位滤波它把信号正向和反向各滤一次抵消了相位偏移代价是不能对实时数据直接用——它要动用未来的数据来修正当前点。对离线分析来说是福音对实时嵌入式分析就得换lfilter。filtfilt 的边界处理需要特别留意。padtypeodd 表示对边界做奇延拓比默认的反射延拓在多数振动信号上更平稳padlen 是延拓长度如果留空SciPy会自动算一个推荐值。但自动值在短信号上不一定靠谱我后面专门有一章讲这里的坑。大量实际项目里我的做法是先在整段长时间信号上滤波再切出需要的区间而不是直接对短片段滤波。3. 频域分析和特征提取从FFT到频谱峰值检测的工程写法3.1 频谱分辨率与窗函数的选择走到频域这一步很多人第一个问题是为什么同样的数据别人的频谱比我清晰答案大概率出在FFT点数和窗函数上。频谱分辨率是 fs / N其中 fs 是采样率N 是FFT点数。比如采样率是10kHz取2048个点做FFT频率分辨率就是 10000 / 2048约等于4.88Hz。这意味着频域里两个相隔小于4.88Hz的频率成分在谱上会糊成一个峰你看不出它们是两个。我常用的频谱函数def compute_spectrum(x, fs, nfftNone, windowhann): 计算单边幅度谱/功率谱 x x - np.mean(x) if nfft is None: nfft len(x) win signal.get_window(window, nfft) xw x[:nfft] * win X np.fft.rfft(xw, nnfft) freqs np.fft.rfftfreq(nfft, d1.0/fs) amp np.abs(X) * 2.0 / (win.sum()) return freqs, amp这里有个细节幅值校正用 win.sum() 而不是 nfft。因为加窗后信号总能量减小乘以这个系数才能让峰值幅度大致等于真实正弦幅值。这个细节在手册里不显眼但在需要精确测量峰值的场景里直接影响结果。窗函数的选择逻辑分析连续振动信号我默认用汉宁窗它的旁瓣衰减快能有效抑制频率泄漏如果要做幅值精确测量用平顶窗如果只是看脉冲响应矩形窗反而更合适。重点不是哪种窗最好而是你的分析目标是什么。还有一个常被忽略的点FFT点数可以大于信号长度吗可以零填充能让频谱曲线更平滑看起来更细腻但不会提高真实物理分辨率——那么细的峰不是你测出来的是插值出来的。真实分辨率只由信号长度决定。3.2 频谱峰值检测的代码实现信号进频域之后常见的任务是在谱上找几个峰测转速基频、找故障特征频率、跟踪某个边带。很多新手直接写 np.max然后取 argmax这在干净信号上没问题工程信号里却很容易翻车——最强峰值可能是转频或者某个共振频率你真正关心的特征峰只是局部小峰。scipy.signal.find_peaks 是更可靠的工具。它的核心参数有三个from scipy.signal import find_peaks def detect_peaks(freqs, amp, min_prominence0.05, min_distance3): 在频谱中检测峰值返回索引和对应频率 peaks_idx, props find_peaks( amp, prominencemin_prominence, distancemin_distance, ) peaks_freq freqs[peaks_idx] peaks_amp amp[peaks_idx] return peaks_freq, peaks_amp, propsprominence突出度是峰值相对两侧最低谷的高度它比单纯用 height绝对高度更抗噪声即使整段底噪都升高了一个明显凸出的峰依然能被识别。distance 则要求两个峰之间至少间隔多少个采样点用来排除谐波干扰——比如你找基频distance 至少要比基频对应的频点距离大。在我自己的识别脚本里流程一般是先在全谱上找前5个显著峰再看目标特征频率附近有没有峰并记录峰值与底噪的比值。比值比绝对幅值更有意义因为谱的总能量随工况变化很大绝对幅值没有可比性。4. 故障诊断场景下的特征组合滑动窗口、峭度和RMS阈值判定4.1 为什么单独看时域特征会误判特征提取是信号处理和分析的分水岭前面都在把数据变干净从这里开始你要让数据说话。很多诊断项目习惯只算一个RMS或者只看峭度实际用下来误报率非常高。RMS反映的是信号整体能量在设备正常运转和严重磨损时都会有明显变化但对早期微弱冲击反应很慢。峭度对冲击非常敏感轴承早期故障时峭度会先抬升可它也容易受单次随机撞击、电磁干扰脉冲的影响。我在测试治具上就见过一个螺丝松了峭度直接从3飙到8但设备运转其实正常。所以我的诊断特征是组合式的RMS看整体能量、峭度看冲击性、频谱特征峰看是否出现特定频率成分。单一特征会误判组合特征互相验证才能把误报率压下来。4.2 一个完整的诊断循环代码示例先分帧再对每一帧计算特征序列。分帧函数def frame_signal(x, win_len1024, hop512): 把一维信号切成帧win_len帧长hop帧移 n_frames 1 (len(x) - win_len) // hop n_frames max(0, n_frames) frames np.zeros((n_frames, win_len)) for i in range(n_frames): start i * hop frames[i] x[start : start win_len] return frames帧长和hop的选择有讲究。帧长决定频率分辨率和特征统计的样本数工程上我一般取信号周期的4到8倍hop取帧长的一半保证相邻帧有50%重叠特征序列更平滑。然后对每帧算RMS和峭度并对一个窄带区间求频谱峰值def extract_frame_features(frame, fs): rms np.sqrt(np.mean(frame**2)) mu np.mean(frame) sigma np.std(frame) kurt np.mean((frame - mu)**4) / (sigma**4 1e-12) # 频域窄带峰值 freqs, amp compute_spectrum(frame, fs, nfftlen(frame)) band (freqs 100) (freqs 200) if np.any(band): band_peak np.max(amp[band]) else: band_peak 0.0 return rms, kurt, band_peak诊断判定不是简单超阈值就报警那样单帧毛刺就能触发误报。我的做法是连续M帧超阈值才判定故障这叫去抖M一般取3到5对应的时间长度大约几秒足以排除随机干扰。阈值怎么定用正常工况下半小时的数据先统计特征均值加3倍标准差作为基线后续再根据现场反馈微调。为什么流程要这么设计因为设备故障是渐进式的早期特征参数大概率在均值附近缓慢抬升算法必须能捕捉趋势而不是突发脉冲。我亲眼见过只用单帧阈值的诊断方案在一个晚上被电火花干扰连续误报十几次把维护人员折腾得不轻。5. 工程化改造numpy向量化、numba加速与实时处理的取舍5.1 向量化写法和for循环的差距一个信号分析脚本从自己的笔记本里能跑到能在产线数据上稳定跑中间还隔着性能问题。数据量一大比如24小时连续振动记录几百兆甚至上GPython的for循环会慢到让人怀疑人生。拿上面的分帧函数来说Python循环每帧做一次切片和计算100万个点数据要切成近2000帧逐帧用Python算RMS和峭度测下来慢得离谱。改进思路是用numpy的滑动窗口视图from numpy.lib.stride_tricks import sliding_window_view def frame_signal_fast(x, win_len1024, hop512): 利用滑动窗口视图生成矩阵向量化计算 windows sliding_window_view(x, win_len)[::hop] return windowssliding_window_view 返回的是原数组的视图不复制数据所以内存开销极小。配合 np.sqrt(np.mean(windows**2, axis1)) 这种向量化写法特征计算直接从循环变成了矩阵运算速度能有数量级提升。不过要小心视图是只读的别试图往里面写东西而且步进切片 [::hop] 在不同numpy版本上行为略有差异建议写完后先检查一下形状。5.2 numba加速的实际效果如果矩阵运算还满足不了实时性要求就上numba。numba 的 njit 装饰器能把纯Python数值计算编译成机器码典型加速比在几十倍到上百倍。from numba import njit njit(cacheTrue) def frame_features_numba(x, win_len1024, hop512): n_frames 1 (len(x) - win_len) // hop n_frames max(0, n_frames) out np.zeros((n_frames, 3)) for i in range(n_frames): frame x[i*hop : i*hop win_len] rms np.sqrt(np.mean(frame**2)) mu np.mean(frame) sigma np.std(frame) kurt np.mean((frame - mu)**4) / (sigma**4 1e-12) out[i, 0] rms out[i, 1] kurt # 频域部分这里用numpy简化处理 fft np.abs(np.fft.rfft(frame)) out[i, 2] np.max(fft[10:20]) return outnumba 有个特点第一次调用要先编译慢一次后面就快了。所以在程序里可以加个预热调用或者把 cacheTrue 打开把编译结果缓存到磁盘下次启动就不用重新编译。什么时候才需要numba我个人的判断标准是处理耗时超过采集时长的三分之一就要优化。比如你采了10秒数据处理却花了20秒那无论如何谈不上实时。先用向量化向量化到不了目标再上numba不要一上来就追求极致优化维护一段numba代码的心智成本明显高于普通numpy。6. 实际运行中我踩过的那些坑假频、边界振铃与频谱泄露6.1 采样率不足导致假频代码根本防不住假频混叠是我见过最隐蔽、也最难用代码救回来的问题。当信号里存在高于采样率一半奈奎斯特频率的频率成分时这些高频成分会被折叠到低频段以幻影峰的形式出现在频谱里。最坑的是这个幻影峰看起来无比真实你甚至会围绕它开展分析。有一次我分析电机振动数据频谱里出现了一个奇怪的窄带峰频段恰好在轴承外圈特征频率附近。我加窗、去噪、调滤波器全部无效。最后查硬件才明白采样率只有2kHz而现场变频器产生了十几kHz的开关噪声它们混叠回了低频段。从算法层面没有任何滤波函数能区分真实特征峰和混叠峰因为它们确实叠在了一个频率上。代码层面的教训有两条第一接数据前先确认采样率是否覆盖了信号可能出现的最高频率如果覆盖不了必须在采集端加抗混叠低通滤波器这是硬件的事第二任何可疑峰出现在奈奎斯特频率附近时先怀疑假频再怀疑故障。6.2 filtfilt的边界振铃和padlenfiltfilt 的零相位特性很好但它也有自己的脾气。滤波本质上是卷积边界处信号突然截断会产生瞬态双向滤波会把这种瞬态放大结果就是信号开头和结尾出现明显的振铃振荡。我踩过这样一个坑一段4秒振动信号为了观察某个事件从中间切出0.5秒短片段然后直接对这个片段做带通滤波。滤波后片段两端振铃幅度大得离谱甚至把真实的冲击特征淹没了。第一次看到波形我以为是设备出了问题后来才发现是滤波造成的人为振荡。解决办法很直接不要在短片段上滤波而是先在整段长信号上做filtfilt再切出关心的区间。如果一定要对短片段滤波那就显式设置 padlen并把它取到滤波器阶数的5倍以上。我一般还会加一个检查逻辑滤波后对比边界段和中间段的能量差异如果边界能量明显抬升就说明振铃没有完全消除。6.3 短窗下的频谱泄露加窗解决不了物理分辨率问题频谱泄露发生在信号截断时。当你用FFT分析一段有限长度信号如果截取的不是整周期能量就会从一个频点泄漏到周围的旁瓣里让真正的峰变矮、变胖底噪抬高。加汉宁窗能大幅抑制旁瓣但代价是主瓣变宽频率分辨率下降。经常有人问我加了窗为什么两个靠得很近的峰还是分不开答案是窗函数改变的是旁瓣高度改变不了主瓣宽度。两个频率间隔小于1/T的成分其中T是信号时长无论加什么窗都会糊成一个峰。比如你只有1秒数据50Hz和50.5Hz在频谱上就是分不开的因为物理分辨率只有1Hz。我在做齿轮箱边带分析时特别有感触边带间隔通常是转频如果转频只有0.5Hz而分析时长只有1秒边带根本出不来。这时候唯一的办法是加长分析数据比如采4秒甚至10秒短数据上再精巧的算法也没戏。还有一个实战细节做FFT前先看信号的RMS有没有突变如果有说明信号不是平稳的分段后的频谱才有意义。平稳性这个前提很多教程一句话带过但实际数据里往往是大问题。平稳信号做平均谱、非平稳信号做时频分析这两条路的代码结构完全不同一开始就要想清楚。最后分享一个我整理这套代码库时形成的习惯。新接到一个信号处理任务我不会急着写特征提取和诊断逻辑而是先把预处理和频谱分析这一段跑通输出几张关键波形图人工看一遍确认没有假频、没有明显振铃、频谱峰的位置都解释得通再继续往下做。这一步多花半小时能省掉后面好几天排查问题的功夫。再有个小技巧所有的预处理函数都把原始信号和参数一并保存下来通常存成npy或者parquet这样无论何时回看分析结果都能精确还原当时数据进去时长什么样。分析代码可以迭代原始数据和处理参数不能丢。这个习惯让我在好几个复盘类项目里省了大麻烦。