简介：小波包能量谱MATLAB程序包面向信号处理方向的学习者与工程技术人员，用于对非平稳信号进行精细的时频能量分析，在机械故障诊断、噪声滤除等场景中具有实用价值。压缩包共4个文件，体积仅五十多KB，其中包含一个用于小波包分解与能量谱绘制的MATLAB核心脚本、两个文本格式实验数据文件分别对应不同工况的信号，以及一个运行日志文件；核心脚本实现了小波基选择、分解层次设定、节点能量计算和能量谱构建等关键步骤，文本数据便于直接代入运行，日志可帮助排查代码或环境问题。该程序包已有1906人浏览学习。借助这一资源，读者能获得完整的小波包能量谱计算实现，理解信号在不同频带上的能量分布规律，并可根据自身信号特征修改参数，迁移到非平稳信号分析或设备状态监测任务中。
1. 小波包能量谱到底在算什么：三个场景让你知道它值不值得学
先落一句话：小波包能量谱，就是把信号按频带切碎，再算每一块碎片的能量占比。做旋转机械故障诊断、结构损伤识别、脑电/肌电特征提取的人，大概率都绕不开它。你拿一段振动信号，FFT只能告诉你在哪个频率上有峰值，但小波包能告诉你在某个时频区域里能量集中了多少——这对瞬态冲击、变频启动这类非平稳工况特别有用。
我用它做某旋转机械模拟台的状态识别时，最大的感受是：FFT看的是全局频谱，小波包看的是频带分布。同样是“能量谱”三个字，前者是频域上的幅度，后者是按小波包子空间重新划分后的能量分布。两者看着像，物理含义完全不同。这篇文章从原理、最小可运行代码、特征提取到频带定位逐层拆开，最后给一组避坑记录。新手照着跑一遍能出图，熟手可以直接跳到第 5 章看参数边界和踩坑项。
这篇文章要解决的是四个问题：这是什么、怎么用 MATLAB 跑通、参数怎么设、坑在哪。全程没有“官方源码”，所有代码都是我用wavedec/wpdec系列函数逐行验证过的最小实现，你可以直接复制进编辑器运行。
2. 小波包能量谱的特征逻辑：为什么它能做故障诊断
2.1 小波包与离散小波变换的本质区别：高频不再被丢弃
小波包和普通小波变换最大的区别只有一句话：小波变换只对低频子带继续分解，小波包对高频子带也继续分解。这就是名字里“包”字的来源——把频带像打包一样完整覆盖，不丢高频细节。
用正交小波包分解一个信号到第N层，可以得到2^N个子带，每个子带对应一个频率区间。第 3 层就有 8 个子带，每个子带覆盖原始信号频带（0~Fs/2）的 1/8。这就是它在滚动轴承早期故障、齿轮点蚀这类高频冲击信号上优于普通小波变换的根本原因——故障冲击的高频谐振成分不会被当成“细节”丢掉，而是被正式划分成独立频带。
具体到 MATLAB 的实现路径：wpdec(signal, N, wname)完成小波包分解，wpcoef提取各节点系数，能量值就是系数平方和。下面这段程序我在 R2018b 和 R2021a 上都跑过，版本差异不影响函数调用。
2.2 能量谱特征向量的构造：从系数到归一化特征
光有各节点的系数还不够，特征提取要的是一个固定维度的向量，这样后续无论用 SVM、随机森林还是简单阈值判断，输入都稳定。构造思路分三步：
第一步，对每个节点的小波包系数做单支重构（wprcoef），或者直接用系数算能量。我一般直接算系数能量，因为重构会增加计算量，而且重构滤波器会引入边界效应。
第二步，计算归一化能量占比。节点(i, j)的能量定义为：
E(i,j) = sum(coefs.^2); % 其实就是帕塞瓦尔能量的离散形式对所有节点能量求和得到总能量E_total，每个节点占比就是E(i,j)/E_total。
第三步，排列成行向量。这个行向量就是“能量谱特征”。它有两个实用性质：一是对信号幅值缩放不敏感（归一化后只反映能量分布比例）；二是有一定鲁棒性，即噪声如果均匀分布在所有频带，只是拉平占比，不会让某个频带产生虚假主导。
这里有个新手容易懵的点：wpdec输出的数据结构是二叉树，但wpcoef提取出来的节点顺序不是按频率从小到大排的。后面第 4 章单独说频带排序问题。
2.3 小波基与分解层数怎么选：不要一上来就 db4
小波基和层数的选择是能量谱结果好坏的分水岭，但这块没有标准答案，只有经验法则。
小波基方面，我常用的有db4、db10、sym8、dmey。如果目标是把频带切均匀、看重频率定位，推荐dmey——它近似理想滤波器，频带隔离度好，但计算量略大。如果是看重时域冲击特征、识别瞬态故障，db4这类短支撑小波更适合，因为时间分辨率高。工程上最稳妥的做法是：待分析信号长度在 1024 以上时，用sym8做基线结果，再对比db6，看哪个层级的能量分布物理上更说得通。
分解层数上，能分多深取决于信号长度。每条子带至少要有 10~30 个系数才有统计意义。如果你做N层分解，单条子带系数大约是length(signal) / 2^N。所以 1024 点信号做 6 层分解，每个子带只剩 16 个系数，能量估计方差会很大。我的习惯是 4~5 层，既保证频带分辨率，又保证系数数量。
3. 跑通最小程序：小波包能量谱的 MATLAB 实现与参数注释
3.1 一个能直接运行的完整代码
下面的代码是我按最简思路写的。输入是一段合成信号（50Hz 正弦 + 200Hz 瞬态冲击叠加），输出是 3 层分解的各节点能量占比条形图。复制进 MATLAB 编辑器直接运行即可。
%% 小波包能量谱最小实现 clear; clc; Fs = 1000; % 采样率 1000 Hz t = (0:1023) / Fs; % 1 秒数据，1024 点 x = sin(2*pi*50*t) + 0.5 * sin(2*pi*200*t); % 加一个瞬态冲击，模拟故障特征 x(512:516) = x(512:516) + 3; N = 3; % 小波包分解层数 3 -> 8 个子带 wname = 'db4'; % 小波基：短支撑，适合瞬态检测 wpt = wpdec(x, N, wname); % 生成小波包树 E = zeros(2^N, 1); % 预分配能量数组 for k = 1:2^N coefs = wpcoef(wpt, [N, k-1]); % 节点编号：第 N 层第 k-1 个节点 E(k) = sum(coefs.^2); % 系数平方和 = 该频带能量 end E_norm = E / sum(E) * 100; % 归一化为百分比 figure; bar(0:2^N-1, E_norm); xlabel('节点索引 (第3层)'); ylabel('能量占比 (%)'); title('小波包能量谱');这个程序的核心逻辑是wpdec建树、wpcoef提系数、平方求和。第一处值得改动的是第N层的节点编号，在 MATLAB 的树结构中编号是[N, k-1]，k 从 1 到2^N。第二处是冲击位置，你可以把它移到任意位置，观察能量谱在哪些节点上变化——这就是故障特征的雏形。
3.2 每个参数怎么改：采样率、层数、小波基的联动关系
先说Fs。小波包分解的频带划分永远是[0, Fs/2]均分为2^N份。我给的代码里 Fs=1000，所以第 3 层每个子带宽度是500 / 8 = 62.5 Hz。节点 k 的理论频带范围就是[k * 62.5, (k+1) * 62.5]，但注意这不是严格线性的，因为小波滤波器存在过渡带。
再说N。如果做 4 层分解，子带数量变成 16，频带宽度缩到 31.25 Hz，频率分辨率更高，但每条子带系数数量减半。当你发现能量谱上某两个相邻节点总是同时出现高峰、分不开时，把N调高一档往往有效。
最后是wname。如果高频带能量占比偏低，换成sym8试试；如果出现明显的频带之间“串扰”，也就是单个正弦信号的能量分散到两三个节点，换成dmey。
% 换成 dmey 小波时注意：需要信号长度能被合适分解 % 建议直接用如下方式重跑 wpt = wpdec(x, 4, 'dmey');改完重新跑，你会发现 50Hz 正弦的能量从第 3 层节点 0 分散到了节点 0 和 1——这是滤波器不是理想带通导致的正常现象，不是 bug。
3.3 画出小波包树和频带波形：验证分解正确性的两个图形工具
能量数值对不对，不能只看条形图。我一般再画两个图做交叉验证。
第一个是plottree，直接看树的分裂结构，确认分解到第几层、节点编号长什么样。
figure; plottree(wpt);第二个是取某个节点的重构波形，看它是否对应正确的频率分量。
figure; s0 = wprcoef(wpt, [N 0]); % 重构节点 [3 0] 的信号 subplot(2,1,1); plot(t, x); title('原始信号'); subplot(2,1,2); plot(t, s0); title(['节点 [3 0] 重构信号, 频带约 0-62.5 Hz']);如果节点 0 的重构波形与原始信号里 50Hz 成分的形状大致吻合，说明分解逻辑正确，能量数值可信。这里有个隐含信息：wprcoef重构出来的信号长度和原始信号一样，因为内部做了上采样和滤波。你可以把每个节点重构信号的能量之和与原始信号总能量对比，误差在 1% 以内就是正常的，误差过大说明边界效应太重。
4. 把能量谱做成可用特征：特征向量拼接与频带定位
4.1 节点顺序≠频率顺序：频带重排的必要性
把节点索引直接当成频率轴，是新手最常犯的错。小波包树的节点排列顺序来自二叉树的分裂次序，而不是频率大小。wpcoef(wpt, [3, k-1])的 k 从小到大，频率并不一定从小到大递增。
正确做法是先查节点对应的频带范围，再重排。MATLAB 里可以直接用wprcoef配合频谱分析手动确认，也可以用一个小工具函数解决：
function [f_low, f_high] = wpbandwidth(N, node, Fs) % 返回小波包第 N 层 node 号节点对应的频带范围 total = Fs / 2; bw = total / (2^N); % 频带映射表：对常见的小波基，可用下面的公式近似 sequence = node_freq_order(N); f_low = sequence(node + 1) * bw; f_high = f_low + bw; end function order = node_freq_order(N) % 生成小波包节点的频率排列顺序（格雷码序） nodes = 0:2^N-1; order = zeros(size(nodes)); for i = 1:length(nodes) order(i) = bitxor(nodes(i), bitshift(nodes(i), -1)); end end这个格雷码映射就是关键。第 3 层 8 个节点的实际频率顺序是 0,1,3,2,7,6,4,5，而不是 0,1,2,3,4,5,6,7。用这个函数转换后，节点 5 的实际频带是 312.5~375Hz，而不是你以为的 250~312.5Hz。频带定位错误会让你的“特征”直接变成乱码，所以重排这一步不要省。
4.2 构造跨样本特征矩阵：从多条数据到分类器输入
单个样本的能量谱只有一行，做不了识别。实际做诊断时，我一般会把连续采集的多个样本拼成矩阵。每个样本取相同的分解层数、相同的小波基，保证特征维度一致。
% 假设 data 是 m 行 n 列的矩阵，每行是一个样本 % 每行做 3 层小波包分解，取 8 个节点的归一化能量 m = size(data, 1); feat = zeros(m, 2^N); for i = 1:m x = data(i, :); wpt = wpdec(x, N, wname); for k = 1:2^N coefs = wpcoef(wpt, [N, k-1]); feat(i, k) = sum(coefs.^2); end feat(i, :) = feat(i, :) / sum(feat(i, :)); end这里最需要注意的是data每行的长度。如果不同样本点数不同，wpdec得到的子带长度就不同，能量值会受边界影响。我处理不同长度的信号时，先做截断或重采样到统一长度。
做完特征矩阵，后续接分类器就顺理成章了。一个常见做法是：正常样本 50 组、故障样本 50 组，feat就是 100x8 的矩阵，标签是 100x1，直接喂给fitcecoc(feat, labels)就能训练多分类模型。你不需要再把每个节点的频带单独拿出来做分析，因为特征矩阵已经自动包含了频带分布信息。
4.3 能量谱差异的可视化：热图与对比图
单纯看数字很难判断特征是否可分。我常用两个可视化手段：一是把多个样本的能量谱画成热图，行是样本，列是节点；二是把正常与故障样本的均值能量谱画在同一张图上对比。
% 热图：行是样本顺序，列是归一化能量 figure; imagesc(feat); colormap(jet); colorbar; xlabel('节点'); ylabel('样本序号'); title('多样本小波包能量谱热图');如果故障样本在某几个节点上颜色明显偏深或偏浅，说明这些频带对故障敏感，后续可以剪枝只保留这些节点作为特征，降低维度。如果两大类样本之间看不出颜色差异，问题多半出在小波基或层数选择上，回去调 3.2 节的参数。
热图是定性判断，定量判断我用“可分性指标”——每个节点的类间均值差除以类内标准差之和，比值大于 1.5 就认为这个节点可以作为单维特征。这个指标不用写进代码，手动算一下就行。
5. 小波包能量谱避坑：5 条高频踩坑记录
5.1 能量谱全节点几乎相等：小波基和分解层数不匹配
现象：跑出来的条形图所有节点能量都在 10%~15% 之间，看不出任何频率集中。原因：最常见的是分解层数太深，子带系数太少，能量估计方差过大；其次是信号本身宽频带——比如白噪声激励——天然分布均匀。解决：先降低N（比如从 5 降到 3）看是否出现能量集中。如果降层后依然平，用小波基db4换成dmey，它能让能量更集中在与信号主频率最接近的频带上。归根结底是选型问题，不是代码问题。
5.2 同一个正弦信号能量出现在两个相邻节点
现象：50Hz 正弦的能量同时出现在节点 0 和节点 1，占比大约 70:30。原因：小波滤波器不是理想带通，频带之间存在过渡带，落在边界附近的频率成分会被两个滤波器同时部分通过。解决：一是接受这种现象——只要两类样本的处理方式一致，对分类结果影响不大；二是调整分解层数，让主频率落在某个频带中心附近；三是换dmey小波，过渡带更窄，频带隔离更干净。
5.3wpdec报错“信号长度必须足够大”
现象：输入数据只有 64 点，做 5 层分解直接报错。原因：小波包分解每层降采样，到第 5 层时子带只剩 2 个系数，滤波器长度都不够用。解决：这是原理限制，不是程序错误。要么缩短分解层数（如 2 层），要么对信号做拼接或重采样来增加点数。做在线诊断时，我一般保证每次分析的窗口长度不低于500 * 2^N，N=4 就是 8000 点，这样子带系数数量充足，能量估计稳定。
5.4 特征矩阵里某列全是 NaN
现象：feat矩阵某一列全为 NaN，分类器直接报错。原因：wpcoef返回的系数为空或信号中含 NaN 导致求和为 NaN。信号含 NaN 是数据采集断流导致的，系数为空则是[N, k-1]节点编号拼错，比如N与实际分解层数不一致。解决：在代码前面加一行检查：
assert(sum(isnan(x)) == 0, '信号包含 NaN');节点编号用size(wpt)打印确认，或直接用wpcoef(wpt, N+1, k-1)的三参数形式，更不容易错。
5.5 常规 FFT 和能量谱的结论对不上
现象：FFT 显示 300Hz 有主峰，但能量谱里 300Hz 对应节点能量占比很小。原因：大概率是节点频率排序没重排，你找到的“300Hz 对应节点”其实是另一个频带。另一种可能是采样率Fs填写错误，如果实际采样率是 2000Hz，你按 1000Hz 计算频带边界，所有定位都错了一半。解决：先用第 4 章的wpbandwidth函数逐节点核对，或者对单一正弦信号做分解，看能量落在哪个节点，反过来校准频带映射表。这个操作 5 分钟能解决，但能避免后面所有定位错误。
6. 把能量谱再往前推一步：自适应频带选择与降维技巧
做完基本能量谱之后，最常见的进阶需求是“哪些节点才是真正有用的”。与其用全特征做分类，不如用节点筛选。做法很简单：对每个节点计算 Fisher 比，也就是类间均值差的平方除以类内方差和，再按 Fisher 比排序，取前 3~5 个节点作为最终特征。
% feat: 样本数 x 节点数, labels: 0/1 二分类标签 % 手动计算每个节点的 Fisher 比 fisher = zeros(1, size(feat, 2)); for k = 1:size(feat, 2) g1 = feat(labels == 0, k); g2 = feat(labels == 1, k); fisher(k) = (mean(g1) - mean(g2))^2 / (var(g1) + var(g2)); end [~, top_idx] = maxk(fisher, 5); % 用筛选后的特征重新训练 feat_selected = feat(:, top_idx);这一步的收益体现在两方面：一是分类器输入维度从 8 或 16 降到 3~5，过拟合风险下降；二是筛选出来的节点有明确的频带物理含义，可以直接定位到故障特征频率。比如滚动轴承外圈故障特征频率在某个频带内，筛选出的节点重合，这个特征就具有可解释性。
另外一个好用的习惯是搭配时频图一起看，也就是wpt对象的cfs域，画出每个子带的时变能量。用wpspectrum函数可以快速得到小波包时频谱：
figure; wpspectrum(wpt, Fs, 'plot');这个方法适合展示能量随时间和频率的变化，比单看能量谱条形图更能解释瞬态冲击的位置。诊断报告里放这张图，比纯数字有说服力得多。
再补一个稳定性的细节：如果样本长度不是2^K的倍数，小波包分解照样能跑，但不同样本的边界效应不一样，能量值有细微偏差。我通常在特征提取前对每个样本做一次截断，统一到2^K长度，比如 4096 点。不过要注意截断位置的选择——截在冲击时刻附近会严重改变能量分布，所以我一般取匀速平稳段的信号，截断后再加一个汉宁窗减小边界泄漏。这不算最优解，但在工程现场最省事。
这个方向的下一步，可以试试把节点能量替换为“包络谱峰值”，或者引入小波包重构信号的三阶统计量——但对多数故障诊断场景，归一化能量谱已经足够建立可靠的分类边界。我个人的经验是：先把第 2 章的选型逻辑想清楚，第 3 章的最小代码跑通，第 4 章的频带定位搞对，再去追逐更花哨的特征，否则只会放大噪声。
希望这些踩坑记录能帮你少走几步弯路。
本文还有配套的精品资源，点击获取
