latest

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([])



#copied EMG signal filtering based on Empirical Mode Decomposition

'''
(a) x (an auxiliary variable) is set to the signal S(t) to be
analysed, and a variable k, which is the number of estimated
IMFs, is set to zero.
(b) Splines are fitted to the upper extrema and the lower
extrema. This will define the lower (LE) and upper
envelopes (UE).
(c) The average envelope, m, is calculated as the arithmetic
mean between UE and LE.
(d) A candidate IMF, h, is estimated as the difference between x
and m.
(e) If h does not fulfill the criteria defining an IMF, it is assigned
to the variable x and the steps (b)–(d) are repeated.
Otherwise, if h is an IMF then the procedure moves to step
(f).
(f) If h is an IMF it is saved as ck, where k is the kth component.
(g) The mean squared error, mse, between two consecutive IMFs
Ck
i and ck, is calculated, and this value is compared to a
stopping condition (usually a very small value, i.e. 10
5
).
(h) If the stopping condition is not reached, the partial residue,
rk, is estimated as the difference between a previous partial
residue rk
1 and ck, and its content is assigned to the dummy
variable x and the steps of (b)–(d) are repeated.
(i) If the stopping condition is reached then the sifting process
is finalized and the final residue rfinal can be estimated as the
difference between S(t) and the sum of all IMFs.

'''


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(110*np.pi*x)+np.sin(np.pi*440*x)
N = 3 #decimals to round mean (higher number means more precision)

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