简介：面向电机故障诊断中振动和电流信号的时频分析需求，这份压缩包给出了希尔伯特黄变换（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 层分成 0～2048Hz 和 2048～4096Hz 两段，第 2 层再把这两段各自对半分成四段，到第 4 层就有 16 个等宽频带，每个带宽 256Hz。用程序包里的话说，这就是“把故障信号按频带切碎了看”。

原始信号 (fs=8192, 0~4096 Hz) ├── 第1层: [0~2048] [2048~4096] ├── 第2层: [0~1024] [1024~2048] [2048~3072] [3072~4096] └── 第4层: 16 个频带，每个 256Hz

正常电机运行时，能量集中在基频和少量低次谐波上，频带能量分布比较“靠前”；轴承外圈故障出现时，高频段能量占比会明显抬升；转子断条则是基频附近的边频带能量上升。每个故障类型在小波包频带上都有不同的能量分布形态，这就是后续分类的特征基础。

提示：小波包分出的频带宽度 = fs / 2^(level+1)。例：fs=8192，第4层每个频带宽 8192/32=256Hz。这个公式在调参时反复用到，建议记下来。

2.3 小波基、分解层数、采样频率怎么定

三个参数里，小波基的选择最像“玄学”，但它有章可循。我处理电机振动和电流信号，默认从 db4 和 sym5 两个基函数入手。db4 是工程界最常用的 Daubechies 小波，时域支撑短，能比较好地卡住冲击特征；sym5 是近似对称小波，相位畸变小，适合电流信号的调制分析。Haar（db1）在冲击信号上会把瞬态“方块化”，特征方差大，别用它做轴承故障。

分解层数取决于你要看哪段频带。目标频带中心频率 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(fs=8192, duration=4.0, mode="bearing"): """ 生成电机故障仿真信号。 包含基频、低次谐波、随机冲击和背景噪声。 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, fs=8192, wavelet="db4", maxlevel=4): """ 对单通道信号做小波包分解，返回最深层各频带的能量占比。 参数: data : 一维信号数组，长度至少 2^maxlevel fs : 采样率，用于计算频带边界 wavelet : 小波基，推荐 db4 / sym5 maxlevel : 分解层数，第 L 层频带宽度 = fs / 2^(L+1) 返回: (ratio, band_edges) ratio : 归一化能量占比，shape 为 (2^maxlevel,) band_edges : 每个频带的起止频率，用于画图和核对 """ wp = pywt.WaveletPacket(data, wavelet, maxlevel=maxlevel) # 关键点：order="freq" 按频率排序取节点，而不是按树的自然顺序 nodes = wp.get_level(maxlevel, order="freq") energies = [] for node in nodes: rec = node.reconstruct(update=False) # 重构该频带对应的时间序列 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^(L+1) 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(update=False)是把第 L 层某个频带的小波系数重建成与原始信号等长的时间序列，update=False表示重构时不修改父节点的系数，避免相互覆盖。能量用的是重构信号的平方和，没有除以时间窗长度，因为后面要做归一化，绝对数值大小无所谓。归一化成占比之后，信号幅值刻度、传感器灵敏度差异都被消掉了，同一特征向量才能在不同设备间对比——这一条是现场诊断能落地的前提。

3.3 第三步：把多段样本组织成特征矩阵

单段信号的ratio是 16 维向量（第 4 层），一次诊断至少要采集几十秒信号，切成多段窗口取平均，才能压住随机噪声带来的抖动。下面这段代码把 N 段窗口的特征堆成一个矩阵，每行是一个样本的特征向量。

def build_feature_matrix(signal, fs=8192, window_len=2, overlap=0.5, maxlevel=4, wavelet="db4"): """ 将长信号滑窗切片，逐片提取小波包能量占比，返回特征矩阵。 窗口 2s，50% 重叠，fs=8192 时每窗 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_ratio=0.15): """ 模板匹配分类。 参数: feature : 当前样本的特征向量，shape (n_bands,) templates : 二维数组，每一行是一类工况的模板特征 labels : 模板对应的工况标签，如 ["正常", "轴承外圈", "转子断条"] threshold_ratio: 距离异常阈值，最小距离超过该比例的判定为"未知故障" 返回: label, min_dist """ dists = np.linalg.norm(templates - feature, axis=1) min_idx = np.argmin(dists) min_dist = dists[min_idx] # 模板向量的平均模长用于把绝对距离归一化，避免受信号幅值影响 avg_norm = np.mean(np.linalg.norm(templates, axis=1)) 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, labels=labels) print("混淆矩阵 (行=真实, 列=预测):") print(cm) # 单类召回率 = 对角线元数 / 该行总样本数，逐行打印查看哪类被漏报 recall = np.diag(cm) / cm.sum(axis=1) for i, lab in enumerate(labels): print(f"{lab}: 召回率 {recall[i]:.2%}")

单类召回率有一个常见的坑：如果轴承故障样本里的冲击幅值偏小，它的特征向量会落在“正常”和“轴承外圈”两个模板的中间地带，距离上更靠近正常模板，于是被错判。这时不要急着调分类器阈值，先回头检查特征提取层数是不是太低，或者窗口长度是不是把冲击能量平均掉了。大多数分类效果不佳的问题，根子在前面特征提取的参数上，不在分类器本身。

5. 电机故障诊断代码里的 5 个常见坑：现象、原因与解决办法

5.1 频带顺序错位：小波包节点编号不是频率顺序

现象：画出来的特征向量里，明明是高频段的位置出了一个大峰值，但对照仿真的故障频率一看，对不上号。

原因：小波包树的自然顺序是按“低-高-低-高”交替排列的，不是从低到高平铺。直接按树的顺序取节点并当成频率顺序用，拿到的是错位的频带。PyWavelets 的get_level(level, order="freq")参数可以解决，但如果库版本较老或换用了别的语言实现，这个参数不一定可靠。

解决：最稳妥的办法是拿一个已知频率的单频正弦信号做全流程测试。生成一个 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对称延拓模式，它比补零在大多数电机信号上更接近真实趋势。程序包里默认用mode="symmetric"，但如果你用别的库，别忘了显式指定。

5.5 小波基选错：Haar 在冲击信号上的“方块化”

现象：用 Haar 小波（db1）做轴承故障，同一个工况反复测的特征向量方差极大，分类器很难稳定。

原因：Haar 小波是分段常数形状，对轴承冲击这类低脉宽瞬态信号，它用矩形的台阶去逼近尖锐脉冲，逼近误差大且对脉冲的相位位置极其敏感。脉冲位置在窗口内偏移几个采样点，系数分布就变化明显，特征自然不稳定。

解决：换 db4 或 sym5 这类有平滑振荡形态的小波基。选择标准不是某本书上的固定答案，而是取一段真实轴承故障数据，分别用 db4、db5、sym5 跑特征提取，看那个基函数下同类样本的方差最小。这个实验做一次，比抄参数管用得多。

6. 让诊断结果更可信：重构残差检查与阈值标定的一个实操技巧

6.1 用重构残差检查分解是否可信

小波包分解有一个天然的验证手段：把最深层所有节点重构的信号相加，应当等于原始信号。这个等式几乎不成立，因为边界延拓和滤波器组的重构误差会留下残差，但残差的量级可以告诉你分解质量。残差超过原始信号 RMS 的 1%，先去查采样点数是不是太少、边界模式是不是设错了，而不是急着看特征。

def check_reconstruction(data, wavelet="db4", maxlevel=4): """ 重构残差校验：各节点重构信号之和减去原始信号，RMS 应低于原始信号 RMS 的 1%。 """ wp = pywt.WaveletPacket(data, wavelet, maxlevel=maxlevel) nodes = wp.get_level(maxlevel, order="freq") recon_sum = np.zeros_like(data) for node in nodes: recon_sum += node.reconstruct(update=False) 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 电机上表现出的能量占比差异可以到几倍。基准必须来自同一台设备的健康状态，这是现场诊断最朴素也最可靠的标定方式。

这个标定过程花的时间不多，但价值很大。我见过太多诊断系统上线后第一周疯狂误报，原因就是拿实验室的阈值去套现场数据。反过来，做过“最安静工况”标定的系统，后续几个月的误报都能控制在个位数。踩过这个坑之后，我每次给新产线部署这套程序流程，第一件事不是看故障数据，而是先把健康基线采足、采准。希望这些参数边界和调试习惯，能帮你少走几段弯路。

本文还有配套的精品资源，点击获取