-
Notifications
You must be signed in to change notification settings - Fork 1
/
empirical_mode_decomposition.py
70 lines (55 loc) · 1.26 KB
/
empirical_mode_decomposition.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
from numpy import *
from matplotlib.pyplot import *
from scipy.signal import *
from scipy.interpolate import splrep, splev
f = 1
f2 = 1
t = linspace(0,10,10000)
x = zeros(len(t))
xmin = zeros(len(t))
xmax = zeros(len(t))
k = 4 # number of mod
newsi = zeros([k+1,len(x)])
mod = zeros([k,len(x)])
noise = random.uniform(-0.05,0.05,10000)
print(noise)
signal = sin (2*pi*f*t) + noise
x = signal
newsi[0,:] = x
for i in range(k):
x = newsi[i,:]
data1 = argrelmax(x)[0]
xmax = x.take(data1)
tmax = t.take(data1)
data2 = argrelmin(x)[0]
xmin = x.take(data2)
tmin = t.take(data2)
coef = splrep(tmax,xmax)
resmax = splev(t,coef)
coef = splrep(tmin,xmin)
resmin = splev(t,coef)
mu = (resmax + resmin)/2
mod[i,:] = x - mu
newsi[i+1,:] = x - mod[i,:]
#draw graph
figure('Test Signal')
xlabel('t')
ylabel('Amplitude')
plot(t, newsi[0,:])
grid()
figure('Empirical Mode Decomposition')
xlabel('t')
ylabel('Amplitude')
for i in range(k):
plot(t, mod[i,:])
grid()
figure('Result')
xlabel('t')
ylabel('Amplitude')
for i in range(k):
plot(t, newsi[i,:])
# plot(t, newsi[0,:], 'y', label = 'all modes')
# plot(t, newsi[3,:], 'b', label = 'signal without first 3 modes')
# legend(loc = 'best')
grid()
show()