SHARE
TWEET

latest

a guest Nov 14th, 2019 135 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import numpy as np
  2. import matplotlib.pyplot as pl
  3. from scipy.signal import argrelextrema
  4. from scipy.interpolate import CubicSpline
  5. import sys
  6. sys.setrecursionlimit(10000)
  7.  
  8.  
  9. def findExtremum(sgnl):
  10.    
  11.     extPoints = np.array([(x,y) for x in argrelextrema(sgnl, np.less) for y in argrelextrema(sgnl, np.greater) ])
  12.     return extPoints;
  13.  
  14. def getEnvelopeMean(sgnl):
  15.     extremePoints = findExtremum(sgnl)
  16.     #define lower envelope Xs
  17.     lx = extremePoints[0,0]
  18.     lx = np.insert(lx,0,0)
  19.     lx = np.append(lx,(len(sgnl)-1))
  20.    
  21.     #define lower envelope Ys
  22.     ly = np.array([(sgnl[x]) for x in lx])
  23.  
  24.     #define upper envelope Xs
  25.     ux = extremePoints[0,1]
  26.     ux = np.insert(ux,0,0)  
  27.     ux = np.append(ux,(len(sgnl)-1))
  28.    
  29.     #define upper envelope Ys
  30.     uy = np.array([(sgnl[x]) for x in ux])
  31.     #FIND LOWER ENVELOPE using cubic spline interpolation
  32.     y1 = CubicSpline(lx,ly)
  33.     x1 = np.linspace(lx.min(), lx.max(), sgnl.size)
  34.    
  35.     #FIND UPPER ENVELOPE using cubic spline interpolation
  36.     y2 = CubicSpline(ux,uy)
  37.     x2 = np.linspace(ux.min(), ux.max(), sgnl.size)
  38.    
  39.     #returns mean envelope signal
  40.     return np.array(y1(x1)+y2(x2)/2.0)
  41.  
  42. def sdCondition(sgnl):
  43.     #Returns true if stopping condition of standart deviation is true
  44.     return np.std(sgnl)<0.3
  45.  
  46. def isMonotonic(sgnl):
  47.     dx = np.diff(sgnl) #dx is an array. np.diff(sgnl) calculates discrete difference along given axis
  48.     return np.all(dx <= 0) or np.all(dx >= 0)   #returns true if array is either increasing or decreasing
  49.    
  50. def isSignalZero(sgnl):
  51.     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.
  52.  
  53.  
  54. def isImf(sgnl):
  55.     global N
  56.     extremas = findExtremum(sgnl)
  57.    
  58.     #counts how many times signal crosses 0 by counting sign changes
  59.     zero_crossings = (np.diff(np.sign(sgnl)) != 0).sum()
  60.      
  61.     #count number of optimas (minima count + maxima count)
  62.     extremaCount = extremas[0,0].size + extremas[0,1].size    
  63.    
  64.     #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
  65.     return (((zero_crossings==extremaCount)or(np.abs(zero_crossings-extremaCount))==1) and isSignalZero(np.around(getEnvelopeMean(sgnl),decimals=N)))
  66.    
  67.  
  68.  
  69. sumc = 0
  70. k= 0
  71. c = np.zeros(100000,dtype=object)
  72. r = np.zeros(100000,dtype=object)
  73. residue = np.array([])
  74.  
  75.  
  76.  
  77. #copied EMG signal filtering based on Empirical Mode Decomposition
  78.  
  79. '''
  80. (a) x (an auxiliary variable) is set to the signal S(t) to be
  81. analysed, and a variable k, which is the number of estimated
  82. IMFs, is set to zero.
  83. (b) Splines are fitted to the upper extrema and the lower
  84. extrema. This will define the lower (LE) and upper
  85. envelopes (UE).
  86. (c) The average envelope, m, is calculated as the arithmetic
  87. mean between UE and LE.
  88. (d) A candidate IMF, h, is estimated as the difference between x
  89. and m.
  90. (e) If h does not fulfill the criteria defining an IMF, it is assigned
  91. to the variable x and the steps (b)–(d) are repeated.
  92. Otherwise, if h is an IMF then the procedure moves to step
  93. (f).
  94. (f) If h is an IMF it is saved as ck, where k is the kth component.
  95. (g) The mean squared error, mse, between two consecutive IMFs
  96. Cki and ck, is calculated, and this value is compared to a
  97. stopping condition (usually a very small value, i.e. 105
  98. ).
  99. (h) If the stopping condition is not reached, the partial residue,
  100. rk, is estimated as the difference between a previous partial
  101. residue rk1 and ck, and its content is assigned to the dummy
  102. variable x and the steps of (b)–(d) are repeated.
  103. (i) If the stopping condition is reached then the sifting process
  104. is finalized and the final residue rfinal can be estimated as the
  105. difference between S(t) and the sum of all IMFs.
  106.  
  107. '''
  108.  
  109.  
  110. def emd(x,orgs):
  111.     global k, c, r, residue, sumc
  112.     m = getEnvelopeMean(x)
  113.     h= x-m
  114.     if(not isImf(h)):
  115.         x = h
  116.         emd(x,orgs)    
  117.     else:
  118.         k+=1
  119.         c[k] = h
  120.         c[0] = c[1]
  121.         if(not sdCondition(c[k])):
  122.             r[0] = orgs
  123.             r[k] = r[k-1] - c[k]
  124.             x = r[k]
  125.             emd(x,orgs)
  126.         else:
  127.             for i in range(1,k+1):
  128.                 sumc += c[i]
  129.             residue = orgs - sumc
  130.            
  131.  #Draws Original signal, IMFS, and Residue  
  132. def drawIMF(x,y):
  133.     global k , c, residue
  134.     pl.subplot(k+1,1,1).set_title("Original signal")
  135.     pl.plot(x,y,'b-')
  136.     for i in range(2,k+1):
  137.         pl.subplot(k+1,1,i).set_title("IMF"+str(i-1))
  138.         pl.plot(x,c[i-1],'r-')
  139.        
  140.     pl.subplot(k+1,1,k+1).set_title("Residue")
  141.     pl.plot(x,residue,'g-')
  142.     print(k)
  143.  
  144.  
  145. x= np.arange(0,1,0.01)
  146. y = np.sin(110*np.pi*x)+np.sin(np.pi*440*x)
  147. N = 3 #decimals to round mean (higher number means more precision)
  148.  
  149. emd(y,y)
  150. drawIMF(x,y)
  151. pl.tight_layout()
  152. pl.show()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top