这篇文章讲解了小波包(WP)与离散小波变换 (DWT) 的不同之处,说明了小波包变换如何执行信号的等带宽精细滤波,而不是DWT中的较粗糙的倍频程滤波。由于小波包的等带宽精细滤波,使得小波包在许多应用中比DWT表现的更好。
离散小波和离散小波包变换
下图展示了一维信号的DWT分解树
![](https://i-blog.csdnimg.cn/blog_migrate/e898c67e9df3015cda53f368303dc3fb.jpeg)
小波包分解树如下
![](https://i-blog.csdnimg.cn/blog_migrate/4b8b45afa2f683110ffda3987b9dbc69.jpeg)
在小波包变换中,滤波操作也应用于小波/细节系数,结果就是小波包将输入信号的子带分解为更精细的等宽间隔。第j级的子带带宽近似为
![](https://i-blog.csdnimg.cn/blog_migrate/f0e9bb63dd7616a4a11e5e7c22f426ce.jpeg)
然而子带带宽的近似值还取决于滤波器的频率局部化程度,对于像“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个频率成分在频率上没有有效分离
![](https://i-blog.csdnimg.cn/blog_migrate/de36ec8fa4cb54aabfab1157491da118.jpeg)
当然,连续小波变换(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);
![](https://i-blog.csdnimg.cn/blog_migrate/7e5b4311bb701dba4df137d713568967.jpeg)
包络谱中轴承故障特征频率及倍频出现,但相对幅值较低。接下来看一下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')
![](https://i-blog.csdnimg.cn/blog_migrate/f15fcc162f19d0fe01bfbc9b856ed067.jpeg)
哈哈,到此居然想将小波包时频谱输入到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')
![](https://i-blog.csdnimg.cn/blog_migrate/43b895c77a27692fbe626abd953fb4b0.jpeg)
看下8个子带对应的包络谱
可见第3个子带和第8个子带的包络谱故障特征频率较为明显,不妨重新放大看一下
![](https://i-blog.csdnimg.cn/blog_migrate/e450189f9caf5ffccf8c6de507eb04f0.jpeg)
这下显而易见了,
码字不易,且行且珍惜,详细代码及数据见如下链接