code: def lpc_ref(signal, order):
      
      if signal.ndim > 1:
          raise ValueError("Array of rank > 1 not supported yet")
      if order > signal.size:
          raise ValueError("Input signal must have a lenght >= lpc order")
  
      if order > 0:
          p = order + 1
          r = np.zeros(p, 'float32')
          # Number of non zero values in autocorrelation one needs for p LPC
          # coefficients
          nx = np.min([p, signal.size])
          x = np.correlate(signal, signal, 'full')
          r[:nx] = x[signal.size - 1:signal.size + order]
          phi = np.dot(sp.linalg.inv(sp.linalg.toeplitz(r[:-1])), -r[1:])
          return np.concatenate(([1.], phi))
      else:
          return np.ones(1, dtype='float32')
docs: Compute the Linear Prediction Coefficients.
  
      Return the order + 1 LPC coefficients for the signal. c = lpc(x, k) will
      find the k+1 coefficients of a k order linear filter:
  
        xp[n] = -c[1] * x[n-2] - ... - c[k-1] * x[n-k-1]
  
      Such as the sum of the squared-error e[i] = xp[i] - x[i] is minimized.
  
      Parameters
      ----------
      signal: array_like
          input signal
      order : int
          LPC order (the output will have order + 1 items)
  
      Notes
      ----
      This is just for reference, as it is using the direct inversion of the
      toeplitz matrix, which is really slow