EMD_Code_Signal_Processing

import numpy as np
import matplotlib.pyplot as pl
from scipy.signal import argrelextrema
from scipy.interpolate import CubicSpline
import sys
sys.setrecursionlimit(10000)


def findExtremum(sgnl):

    extPoints = np.array([(x,y) for x in argrelextrema(sgnl, np.less) for y in argrelextrema(sgnl, np.greater) ])
    return extPoints;

def getEnvelopeMean(sgnl):
    extremePoints = findExtremum(sgnl)
    #define lower envelope Xs
    lx = extremePoints[0,0]
    lx = np.insert(lx,0,0)
    lx = np.append(lx,(len(sgnl)-1))

    #define lower envelope Ys
    ly = np.array([(sgnl[x]) for x in lx])

    #define upper envelope Xs
    ux = extremePoints[0,1]
    ux = np.insert(ux,0,0)
    ux = np.append(ux,(len(sgnl)-1))

    #define upper envelope Ys
    uy = np.array([(sgnl[x]) for x in ux])
    #FIND LOWER ENVELOPE using cubic spline interpolation
    y1 = CubicSpline(lx,ly)
    x1 = np.linspace(lx.min(), lx.max(), sgnl.size)

    #FIND UPPER ENVELOPE using cubic spline interpolation
    y2 = CubicSpline(ux,uy)
    x2 = np.linspace(ux.min(), ux.max(), sgnl.size)

    #returns mean envelope signal
    return np.array(y1(x1)+y2(x2)/2.0)

def sdCondition(sgnl):
    #Returns true if stopping condition of standart deviation is true
    return np.std(sgnl)<0.3

def isMonotonic(sgnl):
    dx = np.diff(sgnl) #dx is an array. np.diff(sgnl) calculates discrete difference along given axis
    return np.all(dx <= 0) or np.all(dx >= 0)   #returns true if array is either increasing or decreasing

def isSignalZero(sgnl):
    return (not np.any(sgnl)) #returns true if all elements of array is zero. np.any tests whether any array element along a given axis evaluates to True.


def isImf(sgnl):
    global N
    extremas = findExtremum(sgnl)

    #counts how many times signal crosses 0 by counting sign changes
    zero_crossings = (np.diff(np.sign(sgnl)) != 0).sum()

    #count number of optimas (minima count + maxima count)
    extremaCount = extremas[0,0].size + extremas[0,1].size

    #returns true if conditions of IMF are true. Those are number of zero crossings and extremas must be equal or differ by 1, and mean envelope should be zero
    return (((zero_crossings==extremaCount)or(np.abs(zero_crossings-extremaCount))==1) and isSignalZero(np.around(getEnvelopeMean(sgnl),decimals=N)))



sumc = 0
k= 0
c = np.zeros(100000,dtype=object)
r = np.zeros(100000,dtype=object)
residue = np.array([])



def emd(x,orgs):
    global k, c, r, residue, sumc
    m = getEnvelopeMean(x)
    h= x-m
    if(not isImf(h)):
        x = h
        emd(x,orgs)
    else:
        k+=1
        c[k] = h
        c[0] = c[1]
        if(not sdCondition(c[k])):
            r[0] = orgs
            r[k] = r[k-1] - c[k]
            x = r[k]
            emd(x,orgs)
        else:
            for i in range(1,k+1):
                sumc += c[i]
            residue = orgs - sumc

 #Draws Original signal, IMFS, and Residue
def drawIMF(x,y):
    global k , c, residue
    pl.subplot(k+1,1,1).set_title("Original signal")
    pl.plot(x,y,'b-')
    for i in range(2,k+1):
        pl.subplot(k+1,1,i).set_title("IMF"+str(i-1))
        pl.plot(x,c[i-1],'r-')

    pl.subplot(k+1,1,k+1).set_title("Residue")
    pl.plot(x,residue,'g-')
    print(k)


x= np.arange(0,1,0.01)
y = np.sin(1020*np.pi*x)+np.sin(np.pi*440*x)
N = 2 #decimals to round mean (higher number means more precision)

emd(y,y)
drawIMF(x,y)
pl.tight_layout()
pl.show()