Mercurial > hg > python-allan
diff _allan.pyx @ 0:3dcd4cfd8f82
Import
author | Daniele Nicolodi <daniele.nicolodi@obspm.fr> |
---|---|
date | Fri, 18 Apr 2014 18:44:54 +0200 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_allan.pyx Fri Apr 18 18:44:54 2014 +0200 @@ -0,0 +1,86 @@ +import numpy as np +cimport numpy as np +cimport cython + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +def _xmvar(np.ndarray[np.double_t, ndim=1] x, int m): + + cdef double c + cdef double v + cdef int n = x.shape[0] + cdef unsigned int j + cdef unsigned int i + + c = 0 + for j in range(n - 3 * m + 1): + v = 0 + for i in range(j, j + m): + v += (x[i + 2 * m] - 2 * x[i + m] + x[i]) / m + c += v * v + return c / (2.0 * m*m * (n - 3 * m + 1)) + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +def _xmvar2(np.ndarray[np.double_t, ndim=1] x, np.ndarray[np.int_t, ndim=1] mm): + + cdef double z + cdef double zi + cdef unsigned int n = x.shape[0] + cdef unsigned int i + cdef unsigned int k + cdef unsigned int m + cdef np.ndarray[np.double_t, ndim=1] w = np.empty([n + 1, ], np.double) + cdef np.ndarray[np.double_t, ndim=1] mvar = np.empty([mm.shape[0], ], np.double) + + mvar.fill(np.nan) + + w[0] = 0.0 + for i in range(n): + w[i + 1] = w[i] + x[i] + + for k, m in enumerate(mm): + if 3 * m > n: + continue + z = 0 + for i in range(n - 3 * m + 1): + zi = (w[i + 3 * m] - 3 * w[i + 2 * m] + 3 * w[i + m] - w[i]) / m + z += zi * zi + mvar[k] = z / (2.0 * m*m * (n - 3 * m + 1)) + + return mvar + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +def _xavar(double[:] x, int[:] mm): + + cdef double z + cdef double zi + cdef unsigned int n = x.shape[0] + cdef unsigned int i + cdef unsigned int k + cdef unsigned int m + cdef int N + + r = np.empty([mm.shape[0], ], np.double) + r.fill(np.nan) + cdef double[::1] avar = r + + for k, m in enumerate(mm): + if m < 1: + continue + N = (n // m) - 2 + if N < 1: + continue + z = 0 + for i in range(N): + zi = x[i * m] - 2 * x[(i + 1) * m] + x[(i + 2) * m] + z += (zi * zi) + avar[k] = z / (2.0 * m*m * N) + + return r