使用EMD【经验模态分解】对一维波形信号进行滤波去噪以及Python实现代码
2)在任意时刻,由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的平均值为零,即上、下包络线相对于时间轴局部对称。
相比于其他分解方法, EMD不需预设基函数,可直接对原始信号进行分解,分解得到IMF分量与最终余项。
这里给出python环境下的实现方法:
两种方式导入:
- git clone 下载到本地环境, python setup.py安装
git clone https://github.com/laszukdawid/PyEMD python setup.py install
- 直接使用Pip install 安装emd-signal数据包到指定环境
pip install EMD-signal
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)*t2) + 6*t2 #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()
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') '''
发布者:全栈程序员-站长,转载请注明出处:https://javaforall.net/225696.html原文链接:https://javaforall.net
