docs: 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