0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
1 import numpy as np
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
2 from scipy import polyfit, polyval
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
3 import matplotlib.pyplot as plt
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
4 import numexpr
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
5
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
6 __all__ = ['outliers', 'detrend', 'fastdetrend',
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
7 'adev', 'xadev', 'xmdev', 'xtdev', 'xavar', 'xmvar', ]
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
8
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
9
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
10 def outliers(y, alpha=6, remove=False):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
11 """Detect and optionally remove outliers"""
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
12
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
13 m = np.mean(y)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
14 v = np.var(y)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
15
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
16 outliers = np.abs(y - m) > (alpha * np.sqrt(v))
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
17
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
18 if remove:
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
19 x = lambda z: z.nonzero()[0]
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
20 y[outliers] = np.interp(x(outliers), x(~outliers), y[~outliers])
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
21
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
22 return np.nonzero(outliers)[0]
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
23
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
24
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
25 def detrend(y, order=1):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
26 """Detrend"""
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
27
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
28 x = np.arange(y.size)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
29 p = polyfit(x, y, order)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
30 y = y - polyval(p, x)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
31 return y
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
32
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
33
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
34 def fastdetrend(y, order=1, downsample=1000):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
35 # time base
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
36 x = np.arange(y.size)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
37 # compute detrend polynomial from downsampled data
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
38 p = polyfit(x[::downsample], y[::downsample], order)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
39
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
40 # construct numexpr expression
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
41 expr = []
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
42 args = {'y': y, 'x': x}
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
43 for i in range(order + 1):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
44 expr.append("p%d" % i + "*x" * i)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
45 args["p%d" % i] = p[order - i]
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
46 expr = "y - (%s)" % " + ".join(expr)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
47
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
48 # evaluate expression
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
49 y = numexpr.evaluate(expr, local_dict=args)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
50 return y
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
51
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
52
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
53 def adev(x, tau, sampl=1.0):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
54 """Allan deviation"""
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
55
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
56 x = np.asarray(x)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
57 tau = np.asarray(tau)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
58
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
59 # allocate output vectors
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
60 adev = np.zeros(tau.size)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
61 dadev = np.zeros(tau.size)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
62
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
63 # samples
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
64 n = x.size
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
65 # partitioning
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
66 p = np.floor(tau * sampl).astype(int)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
67
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
68 for i, m in enumerate(p):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
69 d = x[0:n - n % m].reshape(-1, m)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
70 y = np.mean(d, axis=1)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
71 adev[i] = np.sqrt(0.5 * np.mean(np.diff(y)**2))
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
72 dadev[i] = adev[i] / np.sqrt(y.size)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
73
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
74 return adev, dadev
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
75
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
76
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
77 from _allan import _xmvar, _xmvar2, _xavar
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
78
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
79 # Functions to compute Allan deviation and modified Allan deviation
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
80 # from time error data. For an oscillator described by the equation
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
81 # v(t) = v0 * cos(2*pi * f0 * t + phi(t)) where f0 is the nominal
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
82 # oscillator frequency the time deviation is computed as: x(t) =
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
83 # phi(t) / (2*pi * f0)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
84 #
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
85 # The Symmetricom 5115A phase analyzer outputs data in units of cycles
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
86 # of the input signal. To obtain time deviation from the data stream
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
87 # d(t) it is sufficient to compute: x(t) = d(t) / f0 where f0 is the
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
88 # nominal frequency of the oscillator under test
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
89
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
90 def xmvar0(x, tau, tau0=1.0):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
91
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
92 x = np.asarray(x)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
93 tau = np.asarray(tau)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
94
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
95 # allocate output vectors
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
96 mvar = np.zeros(tau.size)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
97
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
98 # partitioning
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
99 p = (tau // tau0).astype(int)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
100
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
101 for k, m in enumerate(p):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
102 mvar[k] = _xmvar(x, m)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
103 # c = 0
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
104 # for j in xrange(x.size - 3 * m + 1):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
105 # v = 0
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
106 # for i in xrange(j, j + m):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
107 # v += (x[i + 2 * m] - 2 * x[i + m] + x[i]) / m
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
108 # c += v*v
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
109 # mavar[k] = c / (2.0 * m*m * (n - 3 * m + 1))
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
110
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
111 return np.sqrt(mvar) / tau0
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
112
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
113
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
114 def xmvar(x, tau, tau0=1.0):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
115
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
116 x = np.asarray(x)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
117 tau = np.asarray(tau)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
118
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
119 # partitioning
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
120 p = (tau // tau0).astype(int)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
121
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
122 return _xmvar2(x, p) / (tau0 * tau0)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
123
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
124
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
125 def xmdev(x, tau, tau0=1.0):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
126
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
127 return np.sqrt(xmvar(x, tau, tau0))
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
128
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
129
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
130 def xtdev(x, tau, tau0=1.0):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
131
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
132 return np.sqrt(xmvar(x, tau, tau0)) * tau / np.sqrt(3.0)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
133
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
134
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
135 def xavar(x, tau, tau0=1.0):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
136
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
137 x = np.asarray(x)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
138 tau = np.asarray(tau)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
139
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
140 # partitioning
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
141 p = (tau // tau0).astype(int)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
142
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
143 return _xavar(x, p) / (tau0 * tau0)
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
144
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
145
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
146 def xadev(x, tau, tau0=1.0):
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
147
|
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff
changeset
|
148 return np.sqrt(xavar(x, tau, tau0))
|