這篇文章講解了小波包（WP）與離散小波變換 (DWT) 的不同之處，說明了小波包變換如何執行信号的等帶寬精細濾波，而不是DWT中的較粗糙的倍頻程濾波。由于小波包的等帶寬精細濾波，使得小波包在許多應用中比DWT表現的更好。
離散小波和離散小波包變換
下圖展示了一維信号的DWT分解樹
小波包分解樹如下
在小波包變換中，濾波操作也應用于小波/細節系數，結果就是小波包将輸入信号的子帶分解為更精細的等寬間隔。第j級的子帶帶寬近似為
然而子帶帶寬的近似值還取決于濾波器的頻率局部化程度，對于像“fk18”（18 個系數）這樣的 Fejér-Korovkin 濾波器，近似值非常好，然而對于像 Haar ('haar') 這樣的濾波器來說，近似值是不準确的。下面是一個小波包時頻分析的例子。
小波包時頻分析
由于小波包将頻率軸劃分為比 DWT 更精細的子帶間隔，是以小波包在時頻分析方面更具優勢。考慮加性噪聲中頻率為 150 和 200 Hz 的兩個間歇正弦波，采樣頻率為1 kHz，為減緩小波包變換中的時間分辨率損失，使用未下采樣的小波包變換。
dt = 0.001;
t = 0:dt:1-dt;
x = ...
cos(2*pi*150*t).*(t>=0.2 & t<0.4)+sin(2*pi*200*t).*(t>0.6 & t<0.9);
y = x+0.05*randn(size(t));
[wpt,~,F] = my_udwpt(x,'TimeAlign',true);
contour(t,F.*(1/dt),abs(wpt).^2)
grid on
xlabel('Time (secs)')
ylabel('Hz')
title('Time-Frequency Analysis -- Undecimated Wavelet Packet Transform')
小波包變換能夠分離150和200Hz的分量，但DWT卻不能，因為150 和 200 Hz 屬于同一倍頻程。4級 DWT 的倍頻程為（Hz）
下面采用DWT進行分析，這2個頻率成分在頻率上沒有有效分離
當然，連續小波變換(CWT) 比小波包變換具有更高的時頻分辨率，但連續小波變換計算量過高，且當使用正交小波時，小波包變換還具有能量保持的優勢。
小波包變換中的能量保持特性
首先，導入一段地震信号資料
load kobe
plot(kobe)
grid on
xlabel('Seconds')
title('Kobe Earthquake Data')
axis tight
擷取資料的下采樣和未下采樣的小波包變換，分解到第 3 級。為確定下采樣小波包變換的結果一緻，将邊界擴充模式設定為“periodic”
wptreeDecimated = dwpt(kobe,'Level',3,'Boundary','periodic');
wptUndecimated = my_udwpt (kobe,3);
計算下采樣和未下采樣小波包3級系數的總能量，并與原始信号的能量進行比較
decimatedEnergy = sum(cell2mat(cellfun(@(x) sum(abs(x).^2),wptreeDecimated,'UniformOutput',false)))
decimatedEnergy = 2.0189e+11
undecimatedEnergy = sum(sum(abs(wptUndecimated).^2,2))
undecimatedEnergy = 2.0189e+11
signalEnergy = norm(kobe,2)^2
signalEnergy = 2.0189e+11
得證
軸承故障診斷
首先參考如下幾篇文章
基于包絡譜的軸承故障診斷方法-第1篇 - 哥廷根數學學派的文章 - 知乎 https://zhuanlan.zhihu.com/p/533579665
基于包絡譜的軸承故障診斷方法-第2篇 - 哥廷根數學學派的文章 - 知乎 https://zhuanlan.zhihu.com/p/533984966
基于離散小波變換的滾動軸承故障診斷 - 哥廷根數學學派的文章 - 知乎 https://zhuanlan.zhihu.com/p/534179963
以第3個軸承Z軸70Hz第8組資料為例，其包絡譜如下
x1 = load('#8');x1 = x1.Z;x1 = x1(1:5120*2);x1 = x1-mean(x1);
fs = 10240;
N = length(x1);
t = 0:1/fs:(N-1)/fs;
v=70;
[pEnvInner, fEnvInner, xEnvInner, tEnvInner] = envspectrum(x1, fs);plot(fEnvInner, pEnvInner)
xlim([0 800]);hold on;ncomb = 20;helperPlotCombs(ncomb,6.587)*v);
包絡譜中軸承故障特征頻率及倍頻出現，但相對幅值較低。接下來看一下3層小波包分解及小波包時頻譜圖
[wpt,~,F] = my_udwpt (x1,'TimeAlign',true,3);
contour(t,F.*fs,abs(wpt).^2)
grid on
xlabel('Time (secs)')
ylabel('Hz')
title('Time-Frequency Analysis -- Undecimated Wavelet Packet Transform')
哈哈，到此居然想将小波包時頻譜輸入到CNN進行故障診斷。接下來看一下8個子頻帶波形
subplot(4,2,1)
stem(t,wpt(1,:),'Marker','none','ShowBaseline','off')
ylabel('mV')
subplot(4,2,2)
stem(t,wpt(2,:),'Marker','none','ShowBaseline','off')
ylabel('mV')
subplot(4,2,3)
stem(t,wpt(3,:),'Marker','none','ShowBaseline','off')
ylabel('mV')
subplot(4,2,4)
stem(t,wpt(4,:),'Marker','none','ShowBaseline','off')
ylabel('mV')
subplot(4,2,5)
stem(t,wpt(5,:),'Marker','none','ShowBaseline','off')
ylabel('mV')
subplot(4,2,6)
stem(t,wpt(6,:),'Marker','none','ShowBaseline','off')
ylabel('mV')
subplot(4,2,7)
stem(t,wpt(7,:),'Marker','none','ShowBaseline','off')
ylabel('mV')
subplot(4,2,8)
stem(t,wpt(8,:),'Marker','none','ShowBaseline','off')
ylabel('mV')
看下8個子帶對應的包絡譜
可見第3個子帶和第8個子帶的包絡譜故障特征頻率較為明顯，不妨重新放大看一下
這下顯而易見了，
碼字不易，且行且珍惜，詳細代碼及資料見如下連結
https://mianbaoduo.com/o/bread/YpyXk5Zr