基于未下采样小波包分析的轴承故障诊断研究

这篇文章讲解了小波包（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