Mercurial > hg > python-allan

comparison _allan.pyx @ 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
comparison
equal deleted inserted replaced
-1:000000000000 0:3dcd4cfd8f82
1 import numpy as np
2 cimport numpy as np
3 cimport cython
4
5 @cython.boundscheck(False)
6 @cython.wraparound(False)
7 @cython.cdivision(True)
8 def _xmvar(np.ndarray[np.double_t, ndim=1] x, int m):
9
10 cdef double c
11 cdef double v
12 cdef int n = x.shape[0]
13 cdef unsigned int j
14 cdef unsigned int i
15
16 c = 0
17 for j in range(n - 3 * m + 1):
18 v = 0
19 for i in range(j, j + m):
20 v += (x[i + 2 * m] - 2 * x[i + m] + x[i]) / m
21 c += v * v
22 return c / (2.0 * m*m * (n - 3 * m + 1))
23
24
25 @cython.boundscheck(False)
26 @cython.wraparound(False)
27 @cython.cdivision(True)
28 def _xmvar2(np.ndarray[np.double_t, ndim=1] x, np.ndarray[np.int_t, ndim=1] mm):
29
30 cdef double z
31 cdef double zi
32 cdef unsigned int n = x.shape[0]
33 cdef unsigned int i
34 cdef unsigned int k
35 cdef unsigned int m
36 cdef np.ndarray[np.double_t, ndim=1] w = np.empty([n + 1, ], np.double)
37 cdef np.ndarray[np.double_t, ndim=1] mvar = np.empty([mm.shape[0], ], np.double)
38
39 mvar.fill(np.nan)
40
41 w[0] = 0.0
42 for i in range(n):
43 w[i + 1] = w[i] + x[i]
44
45 for k, m in enumerate(mm):
46 if 3 * m > n:
47 continue
48 z = 0
49 for i in range(n - 3 * m + 1):
50 zi = (w[i + 3 * m] - 3 * w[i + 2 * m] + 3 * w[i + m] - w[i]) / m
51 z += zi * zi
52 mvar[k] = z / (2.0 * m*m * (n - 3 * m + 1))
53
54 return mvar
55
56
57 @cython.boundscheck(False)
58 @cython.wraparound(False)
59 @cython.cdivision(True)
60 def _xavar(double[:] x, int[:] mm):
61
62 cdef double z
63 cdef double zi
64 cdef unsigned int n = x.shape[0]
65 cdef unsigned int i
66 cdef unsigned int k
67 cdef unsigned int m
68 cdef int N
69
70 r = np.empty([mm.shape[0], ], np.double)
71 r.fill(np.nan)
72 cdef double[::1] avar = r
73
74 for k, m in enumerate(mm):
75 if m < 1:
76 continue
77 N = (n // m) - 2
78 if N < 1:
79 continue
80 z = 0
81 for i in range(N):
82 zi = x[i * m] - 2 * x[(i + 1) * m] + x[(i + 2) * m]
83 z += (zi * zi)
84 avar[k] = z / (2.0 * m*m * N)
85
86 return r