• API
• FAQ
• Tools
• Archive
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. Cki and ck, is calculated, and this value is compared to a
97. stopping condition (usually a very small value, i.e. 105
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 rk1 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.
Not a member of Pastebin yet?