EMD分解及预测

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
import pandas as pd

#判定当前的时间序列是否是单调序列
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 False
    else:
        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 False
    else:
        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):
        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)
    
#经验模态分解方法
def emd(x):
    imf=[]
    while not ismonotonic(x):
        x1=x
        sd=np.inf
        while sd>0.1 or  (not isImf(x1)):
#             print isImf(x1)
            s1=getspline(x1)
            s2=-getspline(-1*x1)
            x2=x1-(s1+s2)/2
            sd=np.sum((x1-x2)**2)/np.sum(x1**2)
            x1=x2
        
        imf.append(x1)
        x=x-x1
    imf.append(x)
    return imf

def get_line(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):
        t=interpolate.splrep(peaks,y=x[peaks], w=None, xb=None, xe=None,k=len(peaks)-1)
        return t
    t=interpolate.splrep(peaks,y=x[peaks])
    return t

def emd_line(imf):
    imf_line = []
    for i in range(len(imf)):
        imf_line.append(get_line(imf[i]))
    return imf_line
        
def emd_predict(x,imf_line,k):
    predict_imf = []
    for i in range(len(x)):
        predict_x = 0
        for j in range(k):
            predict_x = predict_x + interpolate.splev(x[i],imf_line[j])
        predict_imf.append(predict_x)
    return predict_imf

imf1=emd(train_y)
emd_lines = emd_line(imf1)
predict_imf = emd_predict(train_y,emd_lines,10)

plt.plot(train_x,imf1[5])
plt.show()