Mercurial > hg > python-allan

annotate allan.py @ 0:3dcd4cfd8f82

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import
author Daniele Nicolodi <daniele.nicolodi@obspm.fr>
date Fri, 18 Apr 2014 18:44:54 +0200
parents
children ebc4e1d7c32f
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',
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))