python时间序列EMD分解预测

本文主要是介绍python时间序列EMD分解预测，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

经验模态分解
经验模态分解的python实现

安装包

pyhht

github地址

pip install pyhht

from pyhht.emd import EMD
from pyhht.visualization import plot_imfs
emd = EMD(data.RUL[:10000])
imfs = emd.decompose()

PyEMD

github地址
注：作者将安装包名写错，安装后将文件名改为 PyEMD即可。

pip install EMD-signal

经验模态分解python初步实现

import math
import numpy as np 
import pylab as pl
import matplotlib.pyplot as plt
import scipy.signal as signal
from scipy import fftpack  
import scipy.signal as signal
from scipy import interpolate#判定当前的时间序列是否是单调序列
def ismonotonic(x):max_peaks=signal.argrelextrema(x,np.greater)[0]min_peaks=signal.argrelextrema(x,np.less)[0]all_num=len(max_peaks)+len(min_peaks)if all_num>0:return Falseelse:return True#寻找当前时间序列的极值点
def findpeaks(x):return signal.argrelextrema(x,np.greater)[0]#判断当前的序列是否为 IMF 序列
def isImf(x):N=np.size(x)pass_zero=np.sum(x[0:N-2]*x[1:N-1]<0)#过零点的个数peaks_num=np.size(findpeaks(x))+np.size(findpeaks(-x))#极值点的个数if abs(pass_zero-peaks_num)>1:return Falseelse:return True#获取当前样条曲线
def getspline(x):N=np.size(x)peaks=findpeaks(x)print '当前极值点个数：',len(peaks)if(len(peaks)<=3):if(len(peaks)<2):peaks=np.concatenate(([0],peaks))peaks=np.concatenate((peaks,[N-1]))#这里是为了防止样条次数不够，无法插值的情况t=interpolate.splrep(peaks,y=x[peaks], w=None, xb=None, xe=None,k=len(peaks)-1)return interpolate.splev(np.arange(N),t)t=interpolate.splrep(peaks,y=x[peaks])return interpolate.splev(np.arange(N),t)
#     f=interp1d(np.concatenate(([0,1],peaks,[N+1])),np.concatenate(([0,1],x[peaks],[0])),kind='cubic')
#     f=interp1d(peaks,x[peaks],kind='cubic')
#     return f(np.linspace(1,N,N))#经验模态分解方法
def emd(x):imf=[]while not ismonotonic(x):x1=xsd=np.infwhile sd>0.1 or  (not isImf(x1)):print isImf(x1)s1=getspline(x1)s2=-getspline(-1*x1)x2=x1-(s1+s2)/2sd=np.sum((x1-x2)**2)/np.sum(x1**2)x1=x2imf.append(x1)x=x-x1imf.append(x)return imf

改进算法

import math
import numpy as np 
import pylab as pl
import matplotlib.pyplot as plt
import scipy.signal as signal
from scipy import fftpack  
import scipy.signal as signal
from scipy import interpolate#判定当前的时间序列是否是单调序列
def ismonotonic(x):max_peaks=signal.argrelextrema(x,np.greater)[0]min_peaks=signal.argrelextrema(x,np.less)[0]all_num=len(max_peaks)+len(min_peaks)if all_num>0:return Falseelse:return True#寻找当前时间序列的极值点
def findpeaks(x):#     df_index=np.nonzero(np.diff((np.diff(x)>=0)+0)<0)#     u_data=np.nonzero((x[df_index[0]+1]>x[df_index[0]]))
#     df_index[0][u_data[0]]+=1#     return df_index[0]return signal.argrelextrema(x,np.greater)[0]#判断当前的序列是否为 IMF 序列
def isImf(x):N=np.size(x)pass_zero=np.sum(x[0:N-2]*x[1:N-1]<0)#过零点的个数peaks_num=np.size(findpeaks(x))+np.size(findpeaks(-x))#极值点的个数if abs(pass_zero-peaks_num)>1:return Falseelse:return True#获取当前样条曲线
def getspline(x):N=np.size(x)peaks=findpeaks(x)
#     print '当前极值点个数：',len(peaks)peaks=np.concatenate(([0],peaks))peaks=np.concatenate((peaks,[N-1]))if(len(peaks)<=3):
#         if(len(peaks)<2):
#             peaks=np.concatenate(([0],peaks))
#             peaks=np.concatenate((peaks,[N-1]))
#             t=interpolate.splrep(peaks,y=x[peaks], w=None, xb=None, xe=None,k=len(peaks)-1)
#             return interpolate.splev(np.arange(N),t)t=interpolate.splrep(peaks,y=x[peaks], w=None, xb=None, xe=None,k=len(peaks)-1)return interpolate.splev(np.arange(N),t)t=interpolate.splrep(peaks,y=x[peaks])return interpolate.splev(np.arange(N),t)
#     f=interp1d(np.concatenate(([0,1],peaks,[N+1])),np.concatenate(([0,1],x[peaks],[0])),kind='cubic')
#     f=interp1d(peaks,x[peaks],kind='cubic')
#     return f(np.linspace(1,N,N))#经验模态分解方法
def emd(x):imf=[]while not ismonotonic(x):x1=xsd=np.infwhile sd>0.1 or  (not isImf(x1)):
#             print isImf(x1)s1=getspline(x1)s2=-getspline(-1*x1)x2=x1-(s1+s2)/2sd=np.sum((x1-x2)**2)/np.sum(x1**2)x1=x2imf.append(x1)x=x-x1imf.append(x)return imf