|
使用EMD【经验模态分解】对一维波形信号进行滤波去噪以及Python实现代码
EMD[ Emprical Mode Decomposition]经验模态分解方法此处不再过多用赘述， 该信号处理方法可以把输入信号分解为若干个IMF之和，除最后一个IMF分量为分解余项外，其余分量均满足如下两个约束条件： 1）在整个数据段内，极值点的个数和过零点的个数必须相等或相差最多不能超过一个。
2）在任意时刻，由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的平均值为零，即上、下包络线相对于时间轴局部对称。
相比于其他分解方法， EMD不需预设基函数，可直接对原始信号进行分解，分解得到IMF分量与最终余项。
这里给出python环境下的实现方法：
python 导入EMD-signal 数据包： 数据包官网：https://pypi.org/project/EMD-signal/ GitHub源码地址：https://github.com/laszukdawid/pyemd
两种方式导入：
git clone 下载到本地环境， python setup.py安装
git clone https://github.com/laszukdawid/PyEMD
python setup.py install
直接使用Pip install 安装emd-signal数据包到指定环境
pip install EMD-signal
下面是官网的一个使用案例： 使用EMD经验模态分解对特定函数数据进行处理：
from PyEMD import EMD
import numpy as np
s = np.random.random(100)
emd = EMD()
IMFs = emd(s)
这边导入emd可能会报No module named emd 的错误，记得把PyEMD改为小写，这个报错可以解决，后续报PyEMD的错误时，把相应的PyEMD改为小写pyemd即可，这个错曾经折腾了我很久，后来发现安装的包site-package里边，EMD的文件夹名称是pyemd,与官网案例PyEMD不一样；
下面是对信号
S
(
t
)
=
c
o
s
(
22
π
t
2
)
+
6
t
2
S(t) = cos(22 \pi t^2) + 6t^2
S(t)=cos(22πt2)+6t2的分解结果：
from pyemd import EMD
#import EMD-signal
import numpy as np
import matplotlib.pyplot as plt
s = np.random.random(100)
emd = EMD()
#IMFs = emd(s)
t = np.linspace(0,1,1000)
Signal = np.cos((22*np.pi)*t**2) + 6*t**2
#print(Signal)
IMFs = emd(Signal)
plt.figure(figsize=[40,20])
plt.subplot(311)
plt.plot(t,Signal,'r',label = "Input Signal")
plt.xlabel("Time")
plt.ylabel("signal")
plt.title("Input Signal")
plt.subplot(312)
plt.plot(t,IMFs[0,:],'r',label = "Input Signal")
plt.xlabel("Time")
plt.ylabel("signal")
plt.title("IMF1")
plt.subplot(313)
plt.plot(t,IMFs[1,:],'r',label = "Input Signal")
plt.xlabel("Time")
plt.ylabel("signal")
plt.title("IMF2")
plt.legend()
plt.show()
运行结果：  下面这段代码对于一维波形信号进行经验模态分解，舍弃前p个高频信号达到滤波去噪的目的 【使用的经验模态分解方法有emd eemd ceemdan,具体原理不再赘述了】
from pyemd import EMD
from pyemd import EEMD
from pyemd import CEEMDAN
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
'''
data_list = np.loadtxt('5_H157_C145_HW16H251_C252_HW6.2H332_C384_HW6.5H417_C492_HW7_SNR21_3.txt')
emd = EMD()
eemd = EEMD()
ceemdan = CEEMDAN()
'''
if __name__ == "__main__":
data_list = np.loadtxt('5_H157_C154_HW16H247_C272_HW6.2H329_C3100_HW6.5H424_C4122_HW7_SNR24_1_n.txt')
emd = EMD()
eemd = EEMD()
ceemdan = CEEMDAN()
imfs_emd = emd(data_list)
imfs_eemd = eemd(data_list)
imfs_ceemd = ceemdan(data_list)
x = np.linspace(0,1,len(data_list))
plt.figure(1)
plt.subplot(1 + np.shape(imfs_emd)[0], 1, 1 )
plt.plot(x, data_list, 'r')
plt.title("Signal Input")
for i in range(np.shape(imfs_emd)[0]):
plt.subplot(1 + np.shape(imfs_emd)[0],1,2+i)
plt.plot(x,imfs_emd[i,:],'b')
plt.title("IMF-emd"+str(i))
# plt.show()
plt.figure(2)
plt.subplot(1 + np.shape(imfs_eemd)[0], 1, 1 )
plt.plot(x, data_list, 'r')
plt.title("Signal Input")
for i in range(np.shape(imfs_eemd)[0]):
plt.subplot(1 + np.shape(imfs_eemd)[0],1,2 + i)
plt.plot(x, imfs_eemd[i, :], 'b')
plt.title("IMF-eemd" + str(i))
# plt.show()
plt.figure(3)
plt.subplot(1 + np.shape(imfs_ceemd)[0], 1, 1 )
plt.plot(x, data_list, 'r')
plt.title("Signal Input")
for i in range(np.shape(imfs_ceemd)[0]):
plt.subplot(1 + np.shape(imfs_ceemd)[0],1,2 + i)
plt.plot(x, imfs_ceemd[i, :], 'b')
plt.title("IMF-ceemdan" + str(i))
# plt.show()
###############p=2 show
plt.figure(4)
plt.subplot(3,2,1)
plt.plot(x, data_list, 'r')
plt.title("Signal Input")
plt.xlabel('Time /s')
plt.ylabel('Intensity')
plt.subplot(2,1,2)
Emd_out = np.zeros(1024,)
EEmd_out = np.zeros(1024,)
CEEmd_out = np.zeros(1024,)
for i in range(np.shape(imfs_emd)[0]-2):
Emd_out += imfs_emd[2+i,:]
plt.subplot(3,2,2)
plt.plot(x,Emd_out,'g')
plt.title('EMD Denoising Result P = 2')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
for i in range(np.shape(imfs_eemd)[0]-2):
EEmd_out += imfs_eemd[2+i,:]
plt.subplot(3,2,4)
plt.plot(x,EEmd_out,'g')
plt.title('EEMD Denoising Result P = 2')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
for i in range(np.shape(imfs_ceemd)[0]-2):
CEEmd_out += imfs_ceemd[2+i,:]
plt.subplot(3,2,6)
plt.plot(x,CEEmd_out,'g')
plt.title('CEEMDAN Denoising Result P = 2')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
###########p=3
plt.figure(5)
plt.subplot(3,2,1)
plt.plot(x, data_list, 'r')
plt.title("Signal Input")
plt.xlabel('Time /s')
plt.ylabel('Intensity')
Emd_out3 = np.zeros(1024,)
EEmd_out3 = np.zeros(1024,)
CEEmd_out3 = np.zeros(1024,)
for i in range(np.shape(imfs_emd)[0]-3):
Emd_out3 += imfs_emd[3+i,:]
plt.subplot(3,2,2)
plt.plot(x,Emd_out3,'g')
plt.title('EMD Denoising Result P = 3')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
for i in range(np.shape(imfs_eemd)[0]-3):
EEmd_out3 += imfs_eemd[3+i,:]
plt.subplot(3,2,4)
plt.plot(x,EEmd_out3,'g')
plt.title('EEMD Denoising Result P = 3')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
for i in range(np.shape(imfs_ceemd)[0]-3):
CEEmd_out3 += imfs_ceemd[3+i,:]
plt.subplot(3,2,6)
plt.plot(x,CEEmd_out3,'g')
plt.title('CEEMDAN Denoising Result P = 3')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
plt.show()
'''
def Ga_fit(x, *param):
# std = H/2.355
return param[0] * np.exp(-(x - param[1]) ** 2 / (param[2]) ** 2) + \
param[3] * np.exp(-(x - param[4]) ** 2 / (param[5]) ** 2) + \
param[6] * np.exp(-(x - param[7]) ** 2 / (param[8]) ** 2) + \
param[9] * np.exp(-(x - param[10]) ** 2 / (param[11]) ** 2) + \
param[12] * np.exp(-(x - param[13]) ** 2 / (param[14]) ** 2)
popt_emd,pcov = curve_fit(Ga_fit,x,Emd_out,p0 = [50,0.2,0.05,40,0.4,0.04,30,0.6,0.08,20,0.8,0.04,10,0.9,0.01],maxfev=20000)
popt_eemd, pcov = curve_fit(Ga_fit, x, EEmd_out,
p0=[50, 0.2, 0.05, 40, 0.4, 0.04, 30, 0.6, 0.08, 20, 0.8, 0.04, 10, 0.9, 0.01],maxfev=20000)
popt_ceemd, pcov = curve_fit(Ga_fit, x, CEEmd_out,
p0=[50, 0.2, 0.05, 40, 0.4, 0.04, 30, 0.6, 0.08, 20, 0.8, 0.04, 10, 0.9, 0.01],maxfev=20000)
popt_emd3,pcov = curve_fit(Ga_fit,x,Emd_out3,p0 = [50,0.2,0.05,40,0.4,0.04,30,0.6,0.08,20,0.8,0.04,10,0.9,0.01],maxfev=20000)
popt_eemd3, pcov = curve_fit(Ga_fit, x, EEmd_out3,
p0=[50, 0.2, 0.05, 40, 0.4, 0.04, 30, 0.6, 0.08, 20, 0.8, 0.04, 10, 0.9, 0.01],maxfev=20000)
popt_ceemd3, pcov = curve_fit(Ga_fit, x, CEEmd_out3,
p0=[50, 0.2, 0.05, 40, 0.4, 0.04, 30, 0.6, 0.08, 20, 0.8, 0.04, 10, 0.9, 0.01],maxfev=20000)
plt.figure(6)
plt.subplot(4,2,1)
plt.plot(x, data_list, 'r')
plt.title("Signal Input")
plt.xlabel('Time /s')
plt.ylabel('Intensity')
plt.subplot(4, 2, 3)
plt.plot(x,Ga_fit(x,*popt_emd),'g')
plt.title('EMD_fit p=2')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
plt.subplot(4, 2, 5)
plt.plot(x,Ga_fit(x,*popt_eemd),'g')
plt.title('EEMD_fit p=2')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
plt.subplot(4, 2, 7)
plt.plot(x,Ga_fit(x,*popt_ceemd),'g')
plt.title('CEEMD_fit p=2')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
plt.subplot(4, 2, 4)
plt.plot(x, Ga_fit(x, *popt_emd3), 'g')
plt.title('EMD_fit p=3')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
plt.subplot(4, 2, 6)
plt.plot(x, Ga_fit(x, *popt_eemd3), 'g')
plt.title('EEMD_fit p=3')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
plt.subplot(4, 2, 8)
plt.plot(x, Ga_fit(x, *popt_ceemd3), 'g')
plt.title('CEEMD_fit p=3')
plt.xlabel('Time /s')
plt.ylabel('Intensity')
'''
运行结果：      对于这个一维波形 采用p=2滤波去噪比较合理！ 完结，撒花！
|