简介面向电机故障诊断中振动和电流信号的时频分析需求这份压缩包给出了希尔伯特黄变换HHT的完整Matlab实现适合从事电机故障诊断、设备健康监测的研究人员及工程师参考。包内共6个文件全部为Matlab的m脚本按功能分为经验模态分解、希尔伯特谱分析、主流程调用、峰值查找、结果绘图和FFT辅助分析模块覆盖从信号分解、瞬时频率提取到故障特征呈现的完整链路可直接用于处理电机内圈、外圈故障所诱发的非线性非平稳信号识别谐波失真或新增频率分量。压缩包仅3KB代码精简、无额外依赖运行主脚本即可对采集信号完成分解与谱分析便于按需修改后嵌入自己的故障诊断流程。已有231人学习下载适合具备一定信号处理基础、希望快速上手HHT并落地到电机故障识别场景的读者。1. 电机故障诊断为什么要选小波包变换傅立叶变换的失速与故障信号的本质电机故障诊断的常规路子是先抓信号再做变换。下面这套方案里小波包变换是性价比最高的一个它不像傅立叶变换那样只告诉你“有没有某个频率”而是把故障信号按频带切开告诉你“故障能量落在哪个频段”。对于轴承磨损、转子断条这类早期故障故障特征往往藏在高频段里傅立叶变换看全局频谱很难发现小波包把频带细分之后特征一下就暴露了。这篇笔记从一个可运行的“程序.zip”出发把信号变换、特征提取到故障分类的完整链路讲清楚。新手能照着跑通熟手可以拿参数边界去调自己的数据。2. 小波包变换的原理与选型为什么它能抓住傅立叶丢掉的故障特征2.1 傅立叶变换在电机故障信号上的三个失速点傅立叶变换从来没说它不能用于电机故障诊断而是它在两类最常见的故障信号上表现疲软。第一类是冲击型信号比如轴承外圈点蚀滚珠经过缺陷点时产生一个瞬态冲击这个冲击在时域是一个宽度不到一个毫秒的尖峰在频域却覆盖几百赫兹到几千赫兹。傅立叶变换得到的是整段信号的平均频谱冲击能量只有相对很小的占比平均之后被基频和背景噪声淹没。早期故障里冲击强度本来就弱再用全局平均一稀释频谱上几乎看不到特征。第二类是调制型信号典型的是转子断条。断条故障会在定子电流中产生边频带特征频率与供电基频只差 2 倍的转差频率也就是几赫兹的量级。要在频谱里把这两个相距几赫兹的谱峰分开频率分辨率至少要达到 1Hz 甚至更细而频率分辨率等于采样时长的倒数——想要 0.5Hz 的分辨率就要采样至少 2 秒。更麻烦的是基频的谱泄漏会把边频带整个糊掉window 加得再认真泄漏都在。第三类是非平稳工况。电机的负载波动、转速爬升过程中信号本身就在变傅立叶变换假设信号在观测窗口内是平稳的这个前提在变工况下直接不成立。上面这三点就是傅立叶变换在电机故障信号上的失速点。需要的是一个既能看频带、又能保留时间位置的变换这正是小波包变换的核心价值。2.2 小波包变换的分频逻辑从全频带到等宽子带连续小波变换的思路是给信号戴上一个可变镜头的放大镜高频段用窄时间窗去卡瞬态低频段用宽时间窗去看周期。离散小波变换DWT把这条思路工程化了但它只对低频部分不断二分高频部分分解到第一层就扔在那里不管。轴承故障的冲击成分恰恰落在高频段DWT 等于先把最有用的信息给跳过了。小波包变换WPT纠正了这个问题它对每一层的低频和高频同时做二分。采样率 8192Hz 的信号奈奎斯特频率是 4096Hz第 1 层分成 02048Hz 和 20484096Hz 两段第 2 层再把这两段各自对半分成四段到第 4 层就有 16 个等宽频带每个带宽 256Hz。用程序包里的话说这就是“把故障信号按频带切碎了看”。原始信号 (fs8192, 0~4096 Hz) ├── 第1层: [0~2048] [2048~4096] ├── 第2层: [0~1024] [1024~2048] [2048~3072] [3072~4096] └── 第4层: 16 个频带每个 256Hz正常电机运行时能量集中在基频和少量低次谐波上频带能量分布比较“靠前”轴承外圈故障出现时高频段能量占比会明显抬升转子断条则是基频附近的边频带能量上升。每个故障类型在小波包频带上都有不同的能量分布形态这就是后续分类的特征基础。提示小波包分出的频带宽度 fs / 2^(level1)。例fs8192第4层每个频带宽 8192/32256Hz。这个公式在调参时反复用到建议记下来。2.3 小波基、分解层数、采样频率怎么定三个参数里小波基的选择最像“玄学”但它有章可循。我处理电机振动和电流信号默认从 db4 和 sym5 两个基函数入手。db4 是工程界最常用的 Daubechies 小波时域支撑短能比较好地卡住冲击特征sym5 是近似对称小波相位畸变小适合电流信号的调制分析。Haardb1在冲击信号上会把瞬态“方块化”特征方差大别用它做轴承故障。分解层数取决于你要看哪段频带。目标频带中心频率 f_target分解层数 level floor(log2(fs / (2 * f_target))) 附近。以 8192Hz 采样为例要看 2000Hz 附近的轴承高频谐振第 3 层就够要分辨 3Hz 的转子断条边频带第 4 层也不够需要第 8 层以上但层数越高每层的小波系数越短末尾几个系数全是边界效应特征就“糊了”。我的习惯是分两步先用第 4 层做整体能量分布诊断锁定异常频段再对目标频段单独放大分析。采样频率的上限更实际。轴承故障的高频冲击常伴有几 kHz 的固有谐振采样率低于 8kHz 时这些谐振被混叠吃掉。但采样率高也有代价同样的频带宽度需要更多层数数据量成倍上涨。我做便携式采集设备时用 8192Hz 采样做精密实验室分析用 25.6kHz没有“万能采样率”这回事。三个参数之间是联动关系改采样率就必须重新算层数不能只调一个。3. 把故障信号变成特征向量程序包里的数据生成与特征提取代码3.1 第一步生成带噪声的仿真故障信号压缩包里的data_gen.py负责生成用于验证的特征信号。真实故障数据不一定人人都有仿真信号的目的是先把整条链路跑通——信号里面有明确的故障频率成分特征提取和分类代码只要能识别出这些频率逻辑上就成立。仿真模型分三部分基频正弦、谐波成分、随机冲击序列。import numpy as np def simulate_motor_signal(fs8192, duration4.0, modebearing): 生成电机故障仿真信号。 包含基频、低次谐波、随机冲击和背景噪声。 mode: normal 正常工况 / bearing 轴承外圈故障 / rotor_bar 转子断条 t np.arange(0, duration, 1 / fs) n len(t) f_supply 50.0 # 供电基频 slider 0.03 # 转差率 s # 基波 3、5 次谐波对应电机正常运行的主要频谱成分 signal (1.0 * np.sin(2 * np.pi * f_supply * t) 0.15 * np.sin(2 * np.pi * 3 * f_supply * t) 0.08 * np.sin(2 * np.pi * 5 * f_supply * t)) if mode bearing: # 轴承故障外圈特征频率约 93Hz每转一次发生一次冲击 f_bpfo 93.0 impact_time np.arange(0, duration, 1 / f_bpfo) impulse_indices (impact_time * fs).astype(int) impulse np.zeros(n) decay np.exp(-t[: int(0.008 * fs)] / 0.0016) # 指数衰减冲击 for idx in impulse_indices: if idx n - len(decay): impulse[idx : idx len(decay)] 0.5 * decay signal impulse # 冲击叠加到振动信号上 elif mode rotor_bar: # 转子断条基频附近出现 ±2*s*f_supply 的边频带 f_side 2 * slider * f_supply # 约 3Hz sideband 0.25 * np.sin(2 * np.pi * (f_supply f_side) * t) sideband 0.20 * np.sin(2 * np.pi * (f_supply - f_side) * t) signal sideband # 边频带叠加到基频上 # 背景噪声高斯白噪声模拟传感器与电气干扰 rng np.random.default_rng(42) signal rng.normal(0, 0.03, n) return t, signal这个仿真脚本有三个用途。第一是验证特征提取代码有没有写错如果仿真信号里明明加了 93Hz 的轴承冲击提取结果却在别的频带里占了主导说明频带切分或重构逻辑有问题。第二是为分类器提供有标签的样本后面做模板库时正好用上。第三是帮助你理解“故障特征落到哪个频带”把参数调成自己现场工况的样子比直接套真实数据容易定位问题。提示仿真只是链路验证不能替代现场数据。现场振动信号里的轴系不对中、齿轮啮合、负载波动仿真里没有务必用实际采集数据做最终标定。3.2 第二步小波包分解与频带能量提取拿到信号后的核心动作是小波包分解。PyWavelets 库的WaveletPacket可以一次性完成多层分解关键是要把分解结果按频率顺序取出来。这一步有两个细节最容易出错一是获取第 level 层所有节点二是节点顺序按频率而不是按树的自然顺序排列。import numpy as np import pywt def wpt_energy_ratio(data, fs8192, waveletdb4, maxlevel4): 对单通道信号做小波包分解返回最深层各频带的能量占比。 参数: data : 一维信号数组长度至少 2^maxlevel fs : 采样率用于计算频带边界 wavelet : 小波基推荐 db4 / sym5 maxlevel : 分解层数第 L 层频带宽度 fs / 2^(L1) 返回: (ratio, band_edges) ratio : 归一化能量占比shape 为 (2^maxlevel,) band_edges : 每个频带的起止频率用于画图和核对 wp pywt.WaveletPacket(data, wavelet, maxlevelmaxlevel) # 关键点orderfreq 按频率排序取节点而不是按树的自然顺序 nodes wp.get_level(maxlevel, orderfreq) energies [] for node in nodes: rec node.reconstruct(updateFalse) # 重构该频带对应的时间序列 energies.append(np.sum(rec**2)) # 该频带的总能量 energies np.array(energies) total energies.sum() if total 1e-12: # 信号全零时直接返回均匀占比避免除零 return np.ones_like(energies) / len(energies), None ratio energies / total # 归一化为能量占比 # 频带边界第 L 层每个频带宽度为 fs / 2^(L1) band_width fs / (2**(maxlevel 1)) band_edges [(i * band_width, (i 1) * band_width) for i in range(len(energies))] return ratio, band_edges代码里有几个点要说明。node.reconstruct(updateFalse)是把第 L 层某个频带的小波系数重建成与原始信号等长的时间序列updateFalse表示重构时不修改父节点的系数避免相互覆盖。能量用的是重构信号的平方和没有除以时间窗长度因为后面要做归一化绝对数值大小无所谓。归一化成占比之后信号幅值刻度、传感器灵敏度差异都被消掉了同一特征向量才能在不同设备间对比——这一条是现场诊断能落地的前提。3.3 第三步把多段样本组织成特征矩阵单段信号的ratio是 16 维向量第 4 层一次诊断至少要采集几十秒信号切成多段窗口取平均才能压住随机噪声带来的抖动。下面这段代码把 N 段窗口的特征堆成一个矩阵每行是一个样本的特征向量。def build_feature_matrix(signal, fs8192, window_len2, overlap0.5, maxlevel4, waveletdb4): 将长信号滑窗切片逐片提取小波包能量占比返回特征矩阵。 窗口 2s50% 重叠fs8192 时每窗 16384 点。 win_points int(window_len * fs) step int(win_points * (1 - overlap)) if len(signal) win_points: raise ValueError(信号太短至少需要一个完整窗口) features [] start 0 while start win_points len(signal): segment signal[start : start win_points] ratio, _ wpt_energy_ratio(segment, fs, wavelet, maxlevel) features.append(ratio) start step return np.array(features)窗口长度是一个需要现场权衡的参数。2 秒窗口在 8192Hz 采样下是 16384 个点能满足第 4 层分解后每个子频带 256Hz 分辨率的要求。窗口太短频带能量估计的方差大同一个故障样本反复测出来的特征向量的波动很大窗口太长信号里的非平稳成分被平均掉冲击性故障的瞬态特征反而被弱化。50% 重叠算出来的多个窗口特征不是独立的但用于求均值方差没有问题。特征矩阵的 shape 是(窗口数, 2^maxlevel)。程序包里的约定是每类工况采样 60 秒信号切成 59 个窗口2 秒窗、50% 重叠把 59 行的特征向量取平均得到该类工况的模板向量。这个模板向量就是后面分类判据的基础。4. 故障类型怎么分出来特征模板匹配与分类结果的可信度4.1 模板匹配分类不依赖深度学习的快速方案拿到各类工况的特征模板之后分类可以做得非常朴素新样本的特征向量与哪个模板的欧氏距离最近就判为哪类故障。这个方法不必训练神经网络也不需要人工标注大量数据对冷启动的诊断项目尤其合适。硬件条件有限的现场跑一个小波包提取加距离计算在单片机上都能完成。import numpy as np def classify_by_template(feature, templates, labels, threshold_ratio0.15): 模板匹配分类。 参数: feature : 当前样本的特征向量shape (n_bands,) templates : 二维数组每一行是一类工况的模板特征 labels : 模板对应的工况标签如 [正常, 轴承外圈, 转子断条] threshold_ratio: 距离异常阈值最小距离超过该比例的判定为未知故障 返回: label, min_dist dists np.linalg.norm(templates - feature, axis1) min_idx np.argmin(dists) min_dist dists[min_idx] # 模板向量的平均模长用于把绝对距离归一化避免受信号幅值影响 avg_norm np.mean(np.linalg.norm(templates, axis1)) if avg_norm 0 and min_dist / avg_norm threshold_ratio: return 未知故障, min_dist return labels[min_idx], min_dist这里引入了一个“未知故障”的通道是我在做诊断项目时踩出来的经验。模板库里只有正常、轴承外圈、转子断条三类现场来了一个不对中故障它离所有模板都很远这时强行归到最近的一类会造成误报。与其给一个错误答案不如报“未知故障”让工程师去人工复核。threshold_ratio取 0.15 是比较保守的初值实际使用时要拿正常工况的历史数据做校准。注意模板匹配只适合故障特征相对固定的场景。如果电机负载变化大同一故障的特征向量会漂移这时需要按负载区间分别建模板而不是用一个全局模板硬扛。4.2 特征向量怎么组织才能看出门道特征向量不是拿到 16 维能量占比就直接丢给分类器。我的习惯是同时算两种视图一是单个频带的相对变化二是相邻频带的比值。反映在代码里就是除了ratio原始向量之外还派生一个相邻频带差向量np.diff(ratio)。原因很实在正常电机的能量分布是“低频很高、高频很低”原始向量第一维永远是最大分类器很容易被这个绝对主导特征带偏用相邻频带差就能凸显出高频段“某一段意外突起”的细节这对早期轴承故障特别敏感。组织特征时还要注意样本量。模板匹配对每类工况只有十几个样本都能工作但 T 检验或后续要接机器学习模型时每类至少要有 30 个以上独立窗口。独立的意思是窗口之间不重叠或者重叠率不超过 25%不然算出来的方差是虚低的模型评估会虚高。这个细节在文献里经常被跳过但实际做项目时它决定了你调参的方向对不对。4.3 分类效果怎么评估混淆矩阵与单类准确率的陷阱分类代码写完了自然要问一句效果怎么样最简单的做法是把每类样本各取一半一部分做模板另一部分做测试算出四分类的准确率。但只看总准确率会掩盖问题如果正常样本占了总数的一半一个“什么都判正常”的分类器也能拿到 50% 以上的准确率。正确做法是看混淆矩阵的每一行确认每个故障类别的召回率都达到预期尤其是轴承故障这是最容易和正常状态混淆的一类。# 假设 test_features: 测试特征矩阵test_labels: 真实标签 from sklearn.metrics import confusion_matrix import numpy as np pred_labels [classify_by_template(f, templates, labels)[0] for f in test_features] cm confusion_matrix(test_labels, pred_labels, labelslabels) print(混淆矩阵 (行真实, 列预测):) print(cm) # 单类召回率 对角线元数 / 该行总样本数逐行打印查看哪类被漏报 recall np.diag(cm) / cm.sum(axis1) for i, lab in enumerate(labels): print(f{lab}: 召回率 {recall[i]:.2%})单类召回率有一个常见的坑如果轴承故障样本里的冲击幅值偏小它的特征向量会落在“正常”和“轴承外圈”两个模板的中间地带距离上更靠近正常模板于是被错判。这时不要急着调分类器阈值先回头检查特征提取层数是不是太低或者窗口长度是不是把冲击能量平均掉了。大多数分类效果不佳的问题根子在前面特征提取的参数上不在分类器本身。5. 电机故障诊断代码里的 5 个常见坑现象、原因与解决办法5.1 频带顺序错位小波包节点编号不是频率顺序现象画出来的特征向量里明明是高频段的位置出了一个大峰值但对照仿真的故障频率一看对不上号。原因小波包树的自然顺序是按“低-高-低-高”交替排列的不是从低到高平铺。直接按树的顺序取节点并当成频率顺序用拿到的是错位的频带。PyWavelets 的get_level(level, orderfreq)参数可以解决但如果库版本较老或换用了别的语言实现这个参数不一定可靠。解决最稳妥的办法是拿一个已知频率的单频正弦信号做全流程测试。生成一个 1000Hz 正弦做完分解之后看哪个节点能量最高再和代码里算出的频带边界对照。我每次新装环境都会先跑这个“标定脚本”确认band_edges算对了再继续。5.2 采样率不足导致故障特征被混叠吃掉现象轴承故障信号在 4kHz 采样率下做出来的特征向量与正常工况几乎无差别换成 16kHz 采样差异立刻出现。原因轴承冲击的固有谐振频率常在 2kHz 以上4kHz 采样的奈奎斯特频率只有 2kHz谐振成分被混叠到低频段能量散落在多个频带里指纹被冲淡。解决先查传感器的可用频响再决定采样率。加速度传感器量程内最高频响是 10kHz就用 25.6kHz 采样只有 2kHz 频响的工业变送器再怎么加密采样也没用高频成分在传感器端就已经被滤掉了。这是硬件层面的边界软件无法弥补。5.3 分解层数设太高特征反而“糊了”现象从 4 层改成 6 层分解之后各频带能量占比变得差别很小分类准确率不升反降。原因层数一旦超过 log2(信号长度) 的上限高频子带的小波系数就没有足够的数据点重构时大量依赖边界填充值。这些填充值是“假数据”做出来的能量里面有相当一部分是边界效应特征自然失真。另外层数越高单个频带越窄频谱泄漏让相邻频带之间的能量串扰也越严重。解决控制分解层数 L 满足 2^L 小于等于信号采样点数。更实际的经验是选好目标频带宽度之后用前面那个公式倒推需要的层数够用就不要再往上加。想分辨 3Hz 边频带时正确做法是先降采样而不是单纯增加分解层数。5.4 边界效应让特征向量两端失真现象同一段长信号切成窗口逐窗提特征前几个窗口和最后几个窗口的特征向量明显偏离中间窗口的均值且偏离方向一致。原因小波分解默认在信号两端做边界延拓补零或周期延拓当窗口长度有限时两端延拓区占了真实系数的一大部分重构信号两端几十个点是“编造”出来的。这些点能量计入总能量特征就偏了。解决两种处理方式。一是在重构能量计算时每端舍弃约 2^L / 2 个采样点再求平方和二是尽量用symmetric对称延拓模式它比补零在大多数电机信号上更接近真实趋势。程序包里默认用modesymmetric但如果你用别的库别忘了显式指定。5.5 小波基选错Haar 在冲击信号上的“方块化”现象用 Haar 小波db1做轴承故障同一个工况反复测的特征向量方差极大分类器很难稳定。原因Haar 小波是分段常数形状对轴承冲击这类低脉宽瞬态信号它用矩形的台阶去逼近尖锐脉冲逼近误差大且对脉冲的相位位置极其敏感。脉冲位置在窗口内偏移几个采样点系数分布就变化明显特征自然不稳定。解决换 db4 或 sym5 这类有平滑振荡形态的小波基。选择标准不是某本书上的固定答案而是取一段真实轴承故障数据分别用 db4、db5、sym5 跑特征提取看那个基函数下同类样本的方差最小。这个实验做一次比抄参数管用得多。6. 让诊断结果更可信重构残差检查与阈值标定的一个实操技巧6.1 用重构残差检查分解是否可信小波包分解有一个天然的验证手段把最深层所有节点重构的信号相加应当等于原始信号。这个等式几乎不成立因为边界延拓和滤波器组的重构误差会留下残差但残差的量级可以告诉你分解质量。残差超过原始信号 RMS 的 1%先去查采样点数是不是太少、边界模式是不是设错了而不是急着看特征。def check_reconstruction(data, waveletdb4, maxlevel4): 重构残差校验各节点重构信号之和减去原始信号RMS 应低于原始信号 RMS 的 1%。 wp pywt.WaveletPacket(data, wavelet, maxlevelmaxlevel) nodes wp.get_level(maxlevel, orderfreq) recon_sum np.zeros_like(data) for node in nodes: recon_sum node.reconstruct(updateFalse) residual data - recon_sum rms_signal np.sqrt(np.mean(data**2)) rms_residual np.sqrt(np.mean(residual**2)) ratio rms_residual / rms_signal if rms_signal 0 else float(inf) return ratio, residual我现在的习惯是每换一组参数都会先跑一次这个检查。它不给出故障结论而是给结论一个置信基础。曾经有一次我在现场发现特征向量异常波动排查了很久最后是残差检查暴露了问题——我用的信号截取函数把数组长度截成了奇数小波包分解时最后一点被延拓填充残差到了 5% 的量级。修掉截取逻辑特征立刻稳定了。6.2 阈值标定找“最安静”工况做基准而不是用仿真均值程序包里的模板库直接用仿真数据做了均值这是验证链路用的。现场部署时阈值和模板必须重新标定。具体做法是选一个电机状态最干净的时段——比如新保养完、空载稳定运行、轴承温度稳定之后连续采集 10 分钟信号滑窗提取特征向量。这批特征向量的均值和标准差就是这台电机自己的“健康基线”。判定故障的门槛取 3 倍标准差即当前窗口的特征向量与健康基线的马氏距离超过 3σ 时才算异常。为什么不用仿真或别人项目的阈值因为传感器安装位置、电机功率、负载率、现场环境振动这些因素都会改变特征的绝对尺度同一个轴承故障在 5kW 电机和 200kW 电机上表现出的能量占比差异可以到几倍。基准必须来自同一台设备的健康状态这是现场诊断最朴素也最可靠的标定方式。这个标定过程花的时间不多但价值很大。我见过太多诊断系统上线后第一周疯狂误报原因就是拿实验室的阈值去套现场数据。反过来做过“最安静工况”标定的系统后续几个月的误报都能控制在个位数。踩过这个坑之后我每次给新产线部署这套程序流程第一件事不是看故障数据而是先把健康基线采足、采准。希望这些参数边界和调试习惯能帮你少走几段弯路。本文还有配套的精品资源点击获取