def levinson_1d(r, order):
    """Levinson-Durbin recursion, to efficiently solve symmetric linear systems
    with toeplitz structure.

    Parameters
    ---------
    r : array-like
        input array to invert (since the matrix is symmetric Toeplitz, the
        corresponding pxp matrix is defined by p items only). Generally the
        autocorrelation of the signal for linear prediction coefficients
        estimation. The first item must be a non zero real.

    Notes
    ----
    This implementation is in python, hence unsuitable for any serious
    computation. Use it as educational and reference purpose only.

    Levinson is a well-known algorithm to solve the Hermitian toeplitz
    equation:

                       _          _
        -R[1] = R[0]   R[1]   ... R[p-1]    a[1]
         :      :      :          :      *  :
         :      :      :          _      *  :
        -R[p] = R[p-1] R[p-2] ... R[0]      a[p]
                       _
    with respect to a (  is the complex conjugate). Using the special symmetry
    in the matrix, the inversion can be done in O(p^2) instead of O(p^3).
    """
    r = np.atleast_1d(r)
    if r.ndim > 1:
        raise ValueError("Only rank 1 are supported for now.")

    n = r.size
    if n < 1:
        raise ValueError("Cannot operate on empty array !")
    elif order > n - 1:
        raise ValueError("Order should be <= size-1")

    if not np.isreal(r[0]):
        raise ValueError("First item of input must be real.")
    elif not np.isfinite(1 / r[0]):
        raise ValueError("First item should be != 0")

    # Estimated coefficients
    a = np.empty(order + 1, 'float32')
    # temporary array
    t = np.empty(order + 1, 'float32')
    # Reflection coefficients
    k = np.empty(order, 'float32')

    a[0] = 1.
    e = r[0]

    for i in range(1, order + 1):
        acc = r[i]
        for j in range(1, i):
            acc += a[j] * r[i - j]
        k[i - 1] = -acc / e
        a[i] = k[i - 1]

        for j in range(order):
            t[j] = a[j]

        for j in range(1, i):
            a[j] += k[i - 1] * np.conj(t[i - j])

        e *= 1 - k[i - 1] * np.conj(k[i - 1])

    return a, e, k
def levinson_1d ( r , order ) : r = np . atleast_1d ( r ) if r . ndim > 1 : raise ValueError ( "Only rank 1 are supported for now." ) n = r . size if n < 1 : raise ValueError ( "Cannot operate on empty array !" ) elif order > n - 1 : raise ValueError ( "Order should be <= size-1" ) if not np . isreal ( r [ 0 ] ) : raise ValueError ( "First item of input must be real." ) elif not np . isfinite ( 1 / r [ 0 ] ) : raise ValueError ( "First item should be != 0" ) a = np . empty ( order + 1 , 'float32' ) t = np . empty ( order + 1 , 'float32' ) k = np . empty ( order , 'float32' ) a [ 0 ] = 1. e = r [ 0 ] for i in range ( 1 , order + 1 ) : acc = r [ i ] for j in range ( 1 , i ) : acc += a [ j ] * r [ i - j ] k [ i - 1 ] = - acc / e a [ i ] = k [ i - 1 ] for j in range ( order ) : t [ j ] = a [ j ] for j in range ( 1 , i ) : a [ j ] += k [ i - 1 ] * np . conj ( t [ i - j ] ) e *= 1 - k [ i - 1 ] * np . conj ( k [ i - 1 ] ) return a , e , k
Levinson-Durbin recursion, to efficiently solve symmetric linear systems
    with toeplitz structure.

    Parameters
    ---------
    r : array-like
        input array to invert (since the matrix is symmetric Toeplitz, the
        corresponding pxp matrix is defined by p items only). Generally the
        autocorrelation of the signal for linear prediction coefficients
        estimation. The first item must be a non zero real.

    Notes
    ----
    This implementation is in python, hence unsuitable for any serious
    computation. Use it as educational and reference purpose only.

    Levinson is a well-known algorithm to solve the Hermitian toeplitz
    equation:

                       _          _
        -R[1] = R[0]   R[1]   ... R[p-1]    a[1]
         :      :      :          :      *  :
         :      :      :          _      *  :
        -R[p] = R[p-1] R[p-2] ... R[0]      a[p]
                       _
    with respect to a (  is the complex conjugate). Using the special symmetry
    in the matrix, the inversion can be done in O(p^2) instead of O(p^3).