时间序列信号处理（二）----经验模态分解EMD及其改进的EEMD

库的安装

from PyEMD import EEMD, EMD

我们直接引用库PyEMD就好，安装库需要注意直接pip install pyemd，这时候会出现No module named 'PyEMD',我们只需要取anaconda里面找到安装库的位置，把文件夹pyemd改成PyEMD就好，就可以继续使用了；或者就把关于emd、emd-signal等全部移除，直接pip install EMD-signal的命令安装好PyEMD。

1 经验模态分解EMD

经验模态分解依据数据自身的时间尺度特征来进行信号分解，无需预先设定任何基函数，适用于非线性信号的处理。

分解步骤：

1. 标出局部所有极值点
2. 通过三次样条插值(cubic spline line)连接极大值点构成上包络线(upper envelope)，连接极小值点构成下包络线(lower envelope)
3. 求上下包络线的均值  ​
4. 用输入信号减去上下包络线均值

这样就可得到分解的IMF分量

data1 = np.loadtxt('E:/12k1/0HP/ball18.txt')
data = data1[2048*2:2048*3]
print(data)
index = list(range(len(data)))
N = len(data)
print(N)

Fs = 12000
Ts = 1.0/Fs
t = np.arange(2048)
k = np.arange(N)
T = N/Fs
frq = k/T
frq1 = frq[range(int(N/2))]

emd = EMD()
eIMFs = emd.emd(data)
nIMFs = eIMFs.shape[0]
print(nIMFs)

画出IMF分量时域图和进行傅里叶变换后的图。

plt.figure(figsize=(10, 8))
for n in range(6):
    plt.subplot(7, 1, 1)
    plt.plot(data)
    plt.title('outer')
    plt.subplot(7, 1, n+2)
    plt.plot(eIMFs[n], 'g')
    plt.ylabel("eIMF %i" %(n+1))
    plt.locator_params(axis='y', nbins=5)
plt.xlabel("Time [s]")
plt.tight_layout()
plt.show()

plt.figure(figsize=(10, 8))
for i in range(0, 6):
    plt.subplot(7, 1, 1)
    data_f = abs(np.fft.fft(data)) / N
    data_f1 = data_f[range(int(N/2))]
    plt.plot(frq1, data_f1, 'red')
    plt.title('outer')
    plt.subplot(7, 1, i+2)
    data_f2 = abs(np.fft.fft(eIMFs[i])) / N
    data_f3 = data_f2[range(int(N / 2))]
    plt.plot(frq1, data_f3, 'red')
    plt.xlabel('pinlv(hz)')
    plt.ylabel('IMF{}'.format(i + 1))
plt.tight_layout()
plt.show()

由上图经验模态分解存在严重的模态混叠，就会导致特征提取的不完善，于是就有了改进----集合经验模态分解EEMD

2.EEMD 

EEMD是利用噪声辅助的方式，在整个信号加入噪声。当噪声充斥满整个信号，对信号每做一次测试所得到的信号含噪都不一样，然后进行多次测量，最后取平均值，大概就能消除噪声，最后保留的就是信号本身。

分解步骤：

1.将正态分布的白噪声加到原始信号;

2.将加入白噪声的信号作为一个整体,然后进行EMD分解,得到各IMF分量；

3.重复步骤1和2,每次加入新的正态分布白噪声序列;

4.将每次得到的IMF做集成平均处理后作为最终结果

代码和EMD一样。如果觉得EEMD运行缓慢，可以去源代码调节测试次数，但可能会影响到分解效果，也有可能提高。