Mercurial > hg > python-allan

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