Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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.1
- 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):
- 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=1)))
- 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
- Cki and ck, is calculated, and this value is compared to a
- stopping condition (usually a very small value, i.e. 105
- ).
- (h) If the stopping condition is not reached, the partial residue,
- rk, is estimated as the difference between a previous partial
- residue rk1 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)
- pl.plot(x,y,'b-',label = "Original signal")
- pl.legend()
- for i in range(2,k+1):
- pl.subplot(k+1,1,i)
- pl.plot(x,c[i-1],'r-',label = "IMF"+str(i-1))
- pl.legend()
- pl.subplot(k+1,1,k+1)
- pl.plot(x,residue,'g-', label = "Residue")
- pl.legend()
- x= np.arange(0,1,0.001)
- y = np.cos(4*np.pi*x)+np.sin(4*np.pi*x)
- emd(y,y)
- drawIMF(x,y)
- pl.tight_layout()
- pl.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement