def totdev(data, rate=1.0, data_type="phase", taus=None):
    """ Total deviation.
        Better confidence at long averages for Allan.

    .. math::

        \\sigma^2_{TOTDEV}( m\\tau_0 ) = { 1 \\over 2 (m\\tau_0)^2 (N-2) }
            \\sum_{i=2}^{N-1} ( {x}^*_{i-m} - 2x^*_{i} + x^*_{i+m} )^2


    Where :math:`x^*_i` is a new time-series of length :math:`3N-4`
    derived from the original phase time-series :math:`x_n` of
    length :math:`N` by reflection at both ends.

    FIXME: better description of reflection operation.
    the original data x is in the center of x*:
    x*(1-j) = 2x(1) - x(1+j)  for j=1..N-2
    x*(i)   = x(i)            for i=1..N
    x*(N+j) = 2x(N) - x(N-j)  for j=1..N-2
    x* has length 3N-4
    tau = m*tau0


    FIXME: bias correction http://www.wriley.com/CI2.pdf page 5

    Parameters
    ----------
    phase: np.array
        Phase data in seconds. Provide either phase or frequency.
    frequency: np.array
        Fractional frequency data (nondimensional). Provide either
        frequency or phase.
    rate: float
        The sampling rate for phase or frequency, in Hz
    taus: np.array
        Array of tau values for which to compute measurement


    References
    ----------
    David A. Howe,
    *The total deviation approach to long-term characterization
    of frequency stability*,
    IEEE tr. UFFC vol 47 no 5 (2000)

    NIST SP 1065 eqn (25) page 23

    """
    phase = input_to_phase(data, rate, data_type)
    (phase, m, taus_used) = tau_generator(phase, rate, taus)
    N = len(phase)

    # totdev requires a new dataset
    # Begin by adding reflected data before dataset
    x1 = 2.0 * phase[0] * np.ones((N - 2,))
    x1 = x1 - phase[1:-1]
    x1 = x1[::-1]

    # Reflected data at end of dataset
    x2 = 2.0 * phase[-1] * np.ones((N - 2,))
    x2 = x2 - phase[1:-1][::-1]

    # check length of new dataset
    assert len(x1)+len(phase)+len(x2) == 3*N - 4
    # Combine into a single array
    x = np.zeros((3*N - 4))
    x[0:N-2] = x1
    x[N-2:2*(N-2)+2] = phase # original data in the middle
    x[2*(N-2)+2:] = x2

    devs = np.zeros_like(taus_used)
    deverrs = np.zeros_like(taus_used)
    ns = np.zeros_like(taus_used)

    mid = len(x1)

    for idx, mj in enumerate(m):
        mj = int(mj)
        d0 = x[mid + 1:]
        d1 = x[mid  + mj + 1:]
        d1n = x[mid - mj + 1:]
        e = min(len(d0), len(d1), len(d1n))

        v_arr = d1n[:e] - 2.0 * d0[:e] + d1[:e]
        dev = np.sum(v_arr[:mid] * v_arr[:mid])

        dev /= float(2 * pow(mj / rate, 2) * (N - 2))
        dev = np.sqrt(dev)
        devs[idx] = dev
        deverrs[idx] = dev / np.sqrt(mid)
        ns[idx] = mid

    return remove_small_ns(taus_used, devs, deverrs, ns)
def totdev ( data , rate = 1.0 , data_type = "phase" , taus = None ) : phase = input_to_phase ( data , rate , data_type ) ( phase , m , taus_used ) = tau_generator ( phase , rate , taus ) N = len ( phase ) x1 = 2.0 * phase [ 0 ] * np . ones ( ( N - 2 , ) ) x1 = x1 - phase [ 1 : - 1 ] x1 = x1 [ : : - 1 ] x2 = 2.0 * phase [ - 1 ] * np . ones ( ( N - 2 , ) ) x2 = x2 - phase [ 1 : - 1 ] [ : : - 1 ] assert len ( x1 ) + len ( phase ) + len ( x2 ) == 3 * N - 4 x = np . zeros ( ( 3 * N - 4 ) ) x [ 0 : N - 2 ] = x1 x [ N - 2 : 2 * ( N - 2 ) + 2 ] = phase x [ 2 * ( N - 2 ) + 2 : ] = x2 devs = np . zeros_like ( taus_used ) deverrs = np . zeros_like ( taus_used ) ns = np . zeros_like ( taus_used ) mid = len ( x1 ) for idx , mj in enumerate ( m ) : mj = int ( mj ) d0 = x [ mid + 1 : ] d1 = x [ mid + mj + 1 : ] d1n = x [ mid - mj + 1 : ] e = min ( len ( d0 ) , len ( d1 ) , len ( d1n ) ) v_arr = d1n [ : e ] - 2.0 * d0 [ : e ] + d1 [ : e ] dev = np . sum ( v_arr [ : mid ] * v_arr [ : mid ] ) dev /= float ( 2 * pow ( mj / rate , 2 ) * ( N - 2 ) ) dev = np . sqrt ( dev ) devs [ idx ] = dev deverrs [ idx ] = dev / np . sqrt ( mid ) ns [ idx ] = mid return remove_small_ns ( taus_used , devs , deverrs , ns )
Total deviation.
        Better confidence at long averages for Allan.

    .. math::

        \\sigma^2_{TOTDEV}( m\\tau_0 ) = { 1 \\over 2 (m\\tau_0)^2 (N-2) }
            \\sum_{i=2}^{N-1} ( {x}^*_{i-m} - 2x^*_{i} + x^*_{i+m} )^2


    Where :math:`x^*_i` is a new time-series of length :math:`3N-4`
    derived from the original phase time-series :math:`x_n` of
    length :math:`N` by reflection at both ends.

    FIXME: better description of reflection operation.
    the original data x is in the center of x*:
    x*(1-j) = 2x(1) - x(1+j)  for j=1..N-2
    x*(i)   = x(i)            for i=1..N
    x*(N+j) = 2x(N) - x(N-j)  for j=1..N-2
    x* has length 3N-4
    tau = m*tau0


    FIXME: bias correction http://www.wriley.com/CI2.pdf page 5

    Parameters
    ----------
    phase: np.array
        Phase data in seconds. Provide either phase or frequency.
    frequency: np.array
        Fractional frequency data (nondimensional). Provide either
        frequency or phase.
    rate: float
        The sampling rate for phase or frequency, in Hz
    taus: np.array
        Array of tau values for which to compute measurement


    References
    ----------
    David A. Howe,
    *The total deviation approach to long-term characterization
    of frequency stability*,
    IEEE tr. UFFC vol 47 no 5 (2000)

    NIST SP 1065 eqn (25) page 23