☰
自建信号处理代码库:从预处理到故障诊断的实战指南
2026/10/24 10:37:11 网站建设 项目流程

做信号处理的人,早晚会走到这一步:手里的测试数据越来越多、场景越来越杂,你不再满足于在软件里点按钮出图,更不想每次分析都临时翻开工具库文档凑一段滤波代码。我大概是两年前开始整理自己的一套信号处理与分析代码库,从最底层的去均值、去趋势,到带通滤波器设计、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, kind='linear'): """去掉线性或多项式趋势""" if kind == 'linear': return signal.detrend(x, type='linear') # 多项式趋势,一般不超过3阶 t = np.arange(len(x), dtype=np.float64) coeffs = np.polyfit(t, x, deg=3) trend = np.polyval(coeffs, t) return x - trend

实际采集的振动数据经常带有温漂或线缆扰动引起的趋势。我平时用线性去趋势最多,因为它不会过度扭曲信号形状。什么时候用多项式?我的判断标准是:先做一次FFT,如果低频段还是明显抬升,且原始曲线肉眼可见缓慢弯曲,才尝试3阶多项式。阶数太高容易把真实低频成分一起"拟掉",所以一般不超过3阶。

别忘了坏值处理。传感器偶尔会给出NaN或者一个明显的毛刺。我的习惯是先检查NaN占比,再用中值滤波处理尖峰:

def clean_spikes(x, kernel=3): """中值滤波去除孤立尖峰,kernel必须是奇数""" x = np.where(np.isfinite(x), x, np.median(x)) return signal.medfilt(x, kernel_size=kernel)

这里有个实战经验:中值滤波窗口不能开大。窗口大了,真实的瞬时冲击(比如轴承滚珠敲击外圈)也会被当成毛刺滤掉。我一般取3到5个点,宁可多留一点尖峰,也不牺牲真实冲击成分。

2.2 带通滤波器的设计与零相位处理

振动分析里带通滤波几乎是标配,它用来限定分析频带,把不关心的低频和高频噪声隔离开。基础的滤波器函数:

def bandpass_filter(x, fs, low, high, order=4): """带通滤波,返回滤波后的信号""" b, a = signal.butter(order, [low, high], btype='bandpass', fs=fs) y = signal.filtfilt(b, a, x, padtype='odd', padlen=None) return y

截止频率 low 和 high 怎么定?通常依据对象的工作转速和故障特征频率范围。比如一台转速30Hz左右的旋转设备,轴承外圈故障特征频率往往在100到200Hz附近,那我就会把带通设为80到300Hz,尽量保留特征频带,去掉转频低频大能量和高频噪声。

为什么默认4阶而不是8阶?高阶滤波器衰减更陡,但相位失真更大,数值稳定性也更差。filtfilt是零相位滤波,它把信号正向和反向各滤一次,抵消了相位偏移,代价是不能对实时数据直接用——它要动用"未来"的数据来修正当前点。对离线分析来说是福音,对实时嵌入式分析就得换lfilter。

filtfilt 的边界处理需要特别留意。padtype='odd' 表示对边界做奇延拓,比默认的反射延拓在多数振动信号上更平稳;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, nfft=None, window='hann'): """计算单边幅度谱/功率谱""" 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, n=nfft) freqs = np.fft.rfftfreq(nfft, d=1.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_prominence=0.05, min_distance=3): """在频谱中检测峰值,返回索引和对应频率""" peaks_idx, props = find_peaks( amp, prominence=min_prominence, distance=min_distance, ) peaks_freq = freqs[peaks_idx] peaks_amp = amp[peaks_idx] return peaks_freq, peaks_amp, props

prominence(突出度)是峰值相对两侧最低谷的高度,它比单纯用 height(绝对高度)更抗噪声:即使整段底噪都升高了,一个明显凸出的峰依然能被识别。distance 则要求两个峰之间至少间隔多少个采样点,用来排除谐波干扰——比如你找基频,distance 至少要比基频对应的频点距离大。

在我自己的识别脚本里,流程一般是:先在全谱上找前5个显著峰,再看目标特征频率附近有没有峰,并记录峰值与底噪的比值。比值比绝对幅值更有意义,因为谱的总能量随工况变化很大,绝对幅值没有可比性。

4. 故障诊断场景下的特征组合:滑动窗口、峭度和RMS阈值判定

4.1 为什么单独看时域特征会误判

特征提取是信号处理和分析的"分水岭":前面都在把数据变干净,从这里开始,你要让数据"说话"。很多诊断项目习惯只算一个RMS或者只看峭度,实际用下来误报率非常高。

RMS反映的是信号整体能量,在设备正常运转和严重磨损时都会有明显变化,但对早期微弱冲击反应很慢。峭度对冲击非常敏感,轴承早期故障时峭度会先抬升,可它也容易受单次随机撞击、电磁干扰脉冲的影响。我在测试治具上就见过一个螺丝松了,峭度直接从3飙到8,但设备运转其实正常。

所以我的诊断特征是组合式的:RMS看整体能量、峭度看冲击性、频谱特征峰看是否出现特定频率成分。单一特征会误判,组合特征互相验证,才能把误报率压下来。

4.2 一个完整的诊断循环代码示例

先分帧,再对每一帧计算特征序列。分帧函数:

def frame_signal(x, win_len=1024, hop=512): """把一维信号切成帧,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, nfft=len(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小时连续振动记录,几百兆甚至上G,Python的for循环会慢到让人怀疑人生。

拿上面的分帧函数来说,Python循环每帧做一次切片和计算,100万个点数据要切成近2000帧,逐帧用Python算RMS和峭度,测下来慢得离谱。改进思路是用numpy的滑动窗口视图:

from numpy.lib.stride_tricks import sliding_window_view def frame_signal_fast(x, win_len=1024, hop=512): """利用滑动窗口视图生成矩阵,向量化计算""" windows = sliding_window_view(x, win_len)[::hop] return windows

sliding_window_view 返回的是原数组的视图,不复制数据,所以内存开销极小。配合 np.sqrt(np.mean(windows**2, axis=1)) 这种向量化写法,特征计算直接从循环变成了矩阵运算,速度能有数量级提升。

不过要小心:视图是只读的,别试图往里面写东西;而且步进切片 [::hop] 在不同numpy版本上行为略有差异,建议写完后先检查一下形状。

5.2 numba加速的实际效果

如果矩阵运算还满足不了实时性要求,就上numba。numba 的 @njit 装饰器能把纯Python数值计算编译成机器码,典型加速比在几十倍到上百倍。

from numba import njit @njit(cache=True) def frame_features_numba(x, win_len=1024, hop=512): 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 out

numba 有个特点:第一次调用要先编译,慢一次,后面就快了。所以在程序里可以加个"预热"调用,或者把 cache=True 打开,把编译结果缓存到磁盘,下次启动就不用重新编译。

什么时候才需要numba?我个人的判断标准是:处理耗时超过采集时长的三分之一就要优化。比如你采了10秒数据,处理却花了20秒,那无论如何谈不上实时。先用向量化,向量化到不了目标再上numba,不要一上来就追求极致优化,维护一段numba代码的心智成本明显高于普通numpy。

6. 实际运行中我踩过的那些坑:假频、边界振铃与频谱泄露

6.1 采样率不足导致假频,代码根本防不住

假频(混叠)是我见过最隐蔽、也最难用代码救回来的问题。当信号里存在高于采样率一半(奈奎斯特频率)的频率成分时,这些高频成分会被"折叠"到低频段,以幻影峰的形式出现在频谱里。最坑的是,这个幻影峰看起来无比真实,你甚至会围绕它开展分析。

有一次我分析电机振动数据,频谱里出现了一个奇怪的窄带峰,频段恰好在轴承外圈特征频率附近。我加窗、去噪、调滤波器全部无效。最后查硬件才明白:采样率只有2kHz,而现场变频器产生了十几kHz的开关噪声,它们混叠回了低频段。从算法层面没有任何滤波函数能区分真实特征峰和混叠峰,因为它们确实叠在了一个频率上。

代码层面的教训有两条:第一,接数据前先确认采样率是否覆盖了信号可能出现的最高频率,如果覆盖不了,必须在采集端加抗混叠低通滤波器,这是硬件的事;第二,任何可疑峰出现在奈奎斯特频率附近时,先怀疑假频,再怀疑故障。

6.2 filtfilt的边界振铃和padlen

filtfilt 的零相位特性很好,但它也有自己的脾气。滤波本质上是卷积,边界处信号突然截断,会产生瞬态,双向滤波会把这种瞬态放大,结果就是信号开头和结尾出现明显的振铃振荡。

我踩过这样一个坑:一段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,这样无论何时回看分析结果,都能精确还原当时数据进去时长什么样。分析代码可以迭代,原始数据和处理参数不能丢。这个习惯让我在好几个复盘类项目里省了大麻烦。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询