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