经验小波变换在信号处理及轴承故障诊断中的应用

经验小波变换(EWT)使用所谓的自适应小波细分方案从而创建信号的多分辨率分析 (MRA)，我在科研中用的较多，但是效果也是看运气的，别人都是“数据驱动”，我靠“运气驱动”。EWT从信号频谱分割开始，将输入信号分解为若干的子频带信号，并提供信号的完美重建。EWT由Gilles教授开发，且Gilles和Heal提出了最原始的基于直方图的方法来分割频谱。

多分辨分析MRA将信号分解为不同尺度的分量（子频带），并通过对每个时间点的分量求和来恢复原始信号。目前许多算法都可以归为MRA框架，比如最大重叠离散小波变换，经验模式分解等算法。之所以EWT等基于小波的方法用的多，是因为该类方法具有严谨的数学推导过程，可解释性较强。EWT方法大体可分为4个步骤：

1.预设参数并选择分割频谱的方式；

2.自适应分割信号频谱，获得一组边界；

3.根据边界序列和经验小波构造滤波器组；

4.滤波并重构，获得一系列具有紧支撑的MRA分量

具体的理论推导可以阅读相关的一系列文章，文章末尾会给出相关文献，不再赘述，主要讲一下EWT在Matlab中是如何应用的，以及改进方向

EWT算法

在MATLAB的高版本中，可以直接使用ewt 函数来获取信号的MRA分量， MATLAB中EWT算法的结构如下：

（1）计算信号的Multitaper功率谱估计，这是功率谱的一个平滑、低方差估计，并将估计值归一化为在[0,1]范围内。默认情况下，严格识别所有峰值大于70%的峰值。关于Multitaper功率谱，可参考如下文章。

频谱分析中如何理解taper？ - xiaotu Peng的回答 - 知乎

https://www.zhihu.com/question/372832396/answer/1838971028

（2）经验小波通带的过渡带在相邻峰的几何平均频率处交叉（这个看一下原文就很清楚了），Meyer小波的构造以及确定参数γ的方式见文献[1]，由此形成一个Parseval 紧框架，Parseval 紧框架在现代信号处理中很常用。

（3）为了确定相邻通带之间的边界，可以使用相邻峰之间的第一个局部最小值，如果没有确定局部最小值，默认为几何平均值。注意，局部最小值方法对噪声极为敏感，实际应用时效果不好。

因为EWT的小波形成 Parseval 紧框架，所以滤波器组是自对偶的，即分析滤波器组等于综合滤波器组。 EWT在频域中对信号进行滤波，然后进行逆变换以获得分析系数。

下面给出一个简单的例子。首先，设置一个测试信号

rng default
fs = 500;
t = 0:1/fs:1-1/fs;
f1 = 1./(6/5+cos(2*pi*t));
f2 = 1./(3/2+sin(2*pi*t));
f3 = cos(32*pi*t+cos(64*pi*t));
sig = f1+f2.*f3;
sig = sig+randn(1,length(sig))/2;
sig = sig-mean(sig);
plot(t,sig)
xlabel('Time (sec)')
ylabel('Amplitude')
title('Test Signal')

然后绘制信号的周期谱图和平滑Multitaper估计值

[Pxx,F] = periodogram(sig,[],[],500);
Pxxmt = pmtm(sig,5,[],500,'Tapers','sine','power');
subplot(2,1,1)
plot(F,Pxx)
title('Periodogram')
subplot(2,1,2)
plot(F,Pxxmt)
title('Smoothed Estimate')
xlabel('Frequency (Hz)')

EWT分解

为了描述方便，使用默认设置的ewt 函数来获取信号的 MRA 和有关滤波器组的信息，当然还有好几种算法可以控制ewt函数并获得信号的MRA。

[mra,~,~,info] = ewt(sig);
size(mra)

可以看到，ewt将sig信号默认分解为了两个MRA分量，看一下分解波形，频谱就不画了

subplot(2,1,1);plot(t,mra(:,1));ylabel('mra1');
subplot(2,1,2);plot(t,mra(:,2));ylabel('mra2');

指定最大峰值数量

默认情况下，ewt 会找到2个MRA成分，观察一下ewt的滤波器组通带（中心频率与带宽），因为频率是归一化的，所以要乘以采样频率

info.FilterBank.Passbands*fs

注意到在 22 Hz处有一个边界，第1段有两个峰，因此将 MaxNumPeaks设置为3，以便ewt 使用3个最大峰确定滤波器通带

[mra,cfs,~,info] = ewt(sig,'MaxNumPeaks',3);
info.FilterBank.Passbands*fs

看一下分解波形

subplot(3,1,1);plot(t,mra(:,1));ylabel('mra1');
subplot(3,1,2);plot(t,mra(:,2));ylabel('mra2');
subplot(3,1,3);plot(t,mra(:,3));ylabel('mra3');

验证一下是否完美重建

max(abs(sig'-sum(mra,2)))

sum(sum(abs(cfs).^2))

norm(sig,2)^2

指定百分比阈值

此外还可以设置用于确定在Multitaper功率谱中保留哪些峰的百分比阈值，而不是指定最大峰数，指定最大峰数这个方法在噪声较大的信号中效果很拉跨。信号的Multitaper功率谱估计中的局部最大值被归一化到[0,1]范围内，将 百分比阈值设置为 2。

[~,~,~,info] = ewt(sig,'PeakThresholdPercent',2);
info.FilterBank.Passbands*fs

此外，还可以设置频谱的分割方法（例如localmin，geomean方法）和频率分辨率等参数，就不一一列举了。

下面开始重头戏：轴承故障诊断，试试“运气驱动”吧。

轴承数据的话先看下前几篇文章

基于包络谱的轴承故障诊断方法-第1篇 - 哥廷根数学学派的文章 - 知乎
https://zhuanlan.zhihu.com/p/533579665

基于包络谱的轴承故障诊断方法-第2篇 - 哥廷根数学学派的文章 - 知乎
https://zhuanlan.zhihu.com/p/533984966

基于离散小波变换的滚动轴承故障诊断 - 哥廷根数学学派的文章 - 知乎
https://zhuanlan.zhihu.com/p/534179963

以第3轴承Z轴70Hz第8组数据为例，其包络谱如下

下面开始ewt分解，首先，指定最大峰值数量MaxNumPeaks=3，求得ewt的滤波器组通带信息为：

子频带的波形及包络谱如下：

看一下最大峰值数量MaxNumPeaks=4时的各子频带的包络谱：

看一下最大峰值数量MaxNumPeaks=5时的各子频带的包络谱：

可见，MaxNumPeaks>5就不用考虑了

调整百分比阈值

首先，将百分比阈值设置为 2，得到10个mra分量，看一下ewt的滤波器组通带信息：

ans =

1.0e+03 *

4.0840 5.1200

3.9510 4.0840

3.9380 3.9510

3.3400 3.9380

2.7670 3.3400

2.4940 2.7670

2.1300 2.4940

1.9060 2.1300

1.0800 1.9060

0 1.0800

其子频带分量的包络谱如下：

可见第5个分量故障特征频率十分突出,放大看一下：

再看一下将百分比阈值设置为5，得到6个mra分量，看一下ewt的滤波器组通带信息：

ans =

1.0e+03 *

3.9510 5.1200

3.9380 3.9510

3.3400 3.9380

2.7670 3.3400

2.4940 2.7670

0 2.4940

其子频带分量的包络谱如下：

重看第4个分量：

频谱分割方法

当MaxNumPeaks=5，频谱分割方法为localmin（局部极小）时，各mra分量的包络谱如下图：

效果拉跨，localmin受噪声影响太大了。

当MaxNumPeaks=5，频谱分割方法为geomean（几何平均）时，各子频带分量的包络谱如下图：

重看第4个分量

此外，还有很多优秀的频谱分割方法，比如自适应法，就不一一列举了。

参考文献：

[1] Gilles, Jér^ome. “Empirical Wavelet Transform.” IEEE Transactions on Signal Processing 61, no. 16 (August 2013): 3999–4010. https://doi.org/10.1109/TSP.2013.2265222.

[2] Gilles, Jér^ome, Giang Tran, and Stanley Osher. “2D Empirical Transforms. Wavelets, Ridgelets, and Curvelets Revisited.” SIAM Journal on Imaging Sciences 7, no. 1 (January 2014): 157–86. https://doi.org/10.1137/130923774.

[3] Gilles, Jér^ome, and Kathryn Heal. “A Parameterless Scale-Space Approach to Find Meaningful Modes in Histograms — Application to Image and Spectrum Segmentation.” International Journal of Wavelets, Multiresolution and Information Processing 12, no. 06 (November 2014): 1450044. https://doi.org/10.1142/S0219691314500441.

重要重要重要

EWT的不足也很明显，比如EWT会消耗大量的时间，并且分割出很多无效分量；新的模态混叠现象出现了，同一个分量也可能会被分割为两部分等。对此很多学者对此进行了研究和改进，为了对应经验小波变换的四个步骤，各种改进方法也被作者分为四类。第一类：优化输入参数。为了解决分割频谱和图像时的参数问题，Gilles将傅里叶谱变为尺度空间表示中的函数，提出了概率、“Otsu's”、“k-Means”3种方式，从尺度空间函数中获取有意义的模式。上述基于尺度空间函数的频谱分割方式增加了经验小波变换的自适应性，但是会导致边界的个数增多，也就意味着无效分量会增多，运算耗时也会增加。第二类：更换分割谱。经验小波变换在频谱中划分边界，但Gilles提出的三种方式极易受到噪声的干扰而使边界的个数增多。因此，抑制谱中的噪声成分、放大有用成分有利于边界的划分，无效分量的个数也会被减少。因此可以用自回归功率谱、贝塞尔级数展开(Fourier-Bessel Series Expansion, FBSE)等代替频谱。第三类：优化边界。比如采用顺序统计滤波器(Order tatistic Filter, OSF)的包络估计方法来分割频谱。第四类：筛选MRA分量，这个研究的就比较多了，什么改进峭度指标之类的，不再赘述。

好了，EWT就到这里了，以后慢慢更新一些改进算法，码字不易，且行且珍惜。