Mercurial > hg > python-allan

annotate allan.py @ 2:ff2192f47448

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Minor cleanup
author Daniele Nicolodi <daniele@grinta.net>
date Tue, 24 Mar 2015 18:10:40 +0100
parents ebc4e1d7c32f
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
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',
1
ebc4e1d7c32f Add convenience functions to generate an integration times vector
Daniele Nicolodi <daniele@grinta.net>
parents: 0
diff changeset
7 'adev', 'xadev', 'xmdev', 'xtdev', 'xavar', 'xmvar', 'tau' ]
0
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
1
ebc4e1d7c32f Add convenience functions to generate an integration times vector
Daniele Nicolodi <daniele@grinta.net>
parents: 0
diff changeset
53 def tau(expmin=0, expmax=4):
ebc4e1d7c32f Add convenience functions to generate an integration times vector
Daniele Nicolodi <daniele@grinta.net>
parents: 0
diff changeset
54 return np.ravel(np.power(10, np.arange(expmin, expmax))[:,np.newaxis] * np.arange(0, 10))
ebc4e1d7c32f Add convenience functions to generate an integration times vector
Daniele Nicolodi <daniele@grinta.net>
parents: 0
diff changeset
55
ebc4e1d7c32f Add convenience functions to generate an integration times vector
Daniele Nicolodi <daniele@grinta.net>
parents: 0
diff changeset
56
ebc4e1d7c32f Add convenience functions to generate an integration times vector
Daniele Nicolodi <daniele@grinta.net>
parents: 0
diff changeset
57 def tau124(expmin=0, expmax=4):
ebc4e1d7c32f Add convenience functions to generate an integration times vector
Daniele Nicolodi <daniele@grinta.net>
parents: 0
diff changeset
58 return np.ravel(np.power(10.0, np.arange(expmin, expmax))[:,np.newaxis] * np.array([1, 2, 4]))
ebc4e1d7c32f Add convenience functions to generate an integration times vector
Daniele Nicolodi <daniele@grinta.net>
parents: 0
diff changeset
59
ebc4e1d7c32f Add convenience functions to generate an integration times vector
Daniele Nicolodi <daniele@grinta.net>
parents: 0
diff changeset
60
0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
61 def adev(x, tau, sampl=1.0):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
62 """Allan deviation"""
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
63
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
64 x = np.asarray(x)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
65 tau = np.asarray(tau)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
66
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
67 # allocate output vectors
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
68 adev = np.zeros(tau.size)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
69 dadev = np.zeros(tau.size)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
70
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
71 # samples
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
72 n = x.size
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
73 # partitioning
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
74 p = np.floor(tau * sampl).astype(int)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
75
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
76 for i, m in enumerate(p):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
77 d = x[0:n - n % m].reshape(-1, m)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
78 y = np.mean(d, axis=1)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
79 adev[i] = np.sqrt(0.5 * np.mean(np.diff(y)**2))
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
80 dadev[i] = adev[i] / np.sqrt(y.size)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
81
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
82 return adev, dadev
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
83
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
84
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
85 from _allan import _xmvar, _xmvar2, _xavar
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
86
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
87 # Functions to compute Allan deviation and modified Allan deviation
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
88 # from time error data. For an oscillator described by the equation
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
89 # v(t) = v0 * cos(2*pi * f0 * t + phi(t)) where f0 is the nominal
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
90 # oscillator frequency the time deviation is computed as: x(t) =
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
91 # phi(t) / (2*pi * f0)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
92 #
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
93 # The Symmetricom 5115A phase analyzer outputs data in units of cycles
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
94 # of the input signal. To obtain time deviation from the data stream
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
95 # 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
96 # nominal frequency of the oscillator under test
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
97
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
98 def xmvar0(x, tau, tau0=1.0):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
99
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
100 x = np.asarray(x)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
101 tau = np.asarray(tau)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
102
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
103 # allocate output vectors
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
104 mvar = np.zeros(tau.size)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
105
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
106 # partitioning
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
107 p = (tau // tau0).astype(int)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
108
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
109 for k, m in enumerate(p):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
110 mvar[k] = _xmvar(x, m)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
111 # c = 0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
112 # for j in xrange(x.size - 3 * m + 1):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
113 # v = 0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
114 # for i in xrange(j, j + m):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
115 # v += (x[i + 2 * m] - 2 * x[i + m] + x[i]) / m
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
116 # c += v*v
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
117 # mavar[k] = c / (2.0 * m*m * (n - 3 * m + 1))
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
118
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
119 return np.sqrt(mvar) / tau0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
120
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
121
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
122 def xmvar(x, tau, tau0=1.0):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
123
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
124 x = np.asarray(x)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
125 tau = np.asarray(tau)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
126
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
127 # partitioning
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
128 p = (tau // tau0).astype(int)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
129
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
130 return _xmvar2(x, p) / (tau0 * tau0)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
131
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
132
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
133 def xmdev(x, tau, tau0=1.0):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
134
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
135 return np.sqrt(xmvar(x, tau, tau0))
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
136
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
137
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
138 def xtdev(x, tau, tau0=1.0):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
139
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
140 return np.sqrt(xmvar(x, tau, tau0)) * tau / np.sqrt(3.0)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
141
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
142
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
143 def xavar(x, tau, tau0=1.0):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
144
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
145 x = np.asarray(x)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
146 tau = np.asarray(tau)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
147
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
148 # partitioning
2
ff2192f47448 Minor cleanup
Daniele Nicolodi <daniele@grinta.net>
parents: 1
diff changeset
149 p = (tau // tau0).astype(np.int32)
0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
150
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
151 return _xavar(x, p) / (tau0 * tau0)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
152
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
153
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
154 def xadev(x, tau, tau0=1.0):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
155 return np.sqrt(xavar(x, tau, tau0))