def levinson_1d(r, order):
r = np.atleast_1d(r)<EOL>if r.ndim > <NUM_LIT:1>:<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>n = r.size<EOL>if n < <NUM_LIT:1>:<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>elif order > n - <NUM_LIT:1>:<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>if not np.isreal(r[<NUM_LIT:0>]):<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>elif not np.isfinite(<NUM_LIT:1> / r[<NUM_LIT:0>]):<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>a = np.empty(order + <NUM_LIT:1>, '<STR_LIT>')<EOL>t = np.empty(order + <NUM_LIT:1>, '<STR_LIT>')<EOL>k = np.empty(order, '<STR_LIT>')<EOL>a[<NUM_LIT:0>] = <NUM_LIT:1.><EOL>e = r[<NUM_LIT:0>]<EOL>for i in range(<NUM_LIT:1>, order + <NUM_LIT:1>):<EOL><INDENT>acc = r[i]<EOL>for j in range(<NUM_LIT:1>, i):<EOL><INDENT>acc += a[j] * r[i - j]<EOL><DEDENT>k[i - <NUM_LIT:1>] = -acc / e<EOL>a[i] = k[i - <NUM_LIT:1>]<EOL>for j in range(order):<EOL><INDENT>t[j] = a[j]<EOL><DEDENT>for j in range(<NUM_LIT:1>, i):<EOL><INDENT>a[j] += k[i - <NUM_LIT:1>] * np.conj(t[i - j])<EOL><DEDENT>e *= <NUM_LIT:1> - k[i - <NUM_LIT:1>] * np.conj(k[i - <NUM_LIT:1>])<EOL><DEDENT>return a, e, k<EOL>
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).