code: def lpc(blk, order=None):
    
    # Calculate the covariance for each lag pair
    phi = lag_matrix(blk, order)
    order = len(phi) - 1
  
    # Inner product for filters based on above statistics
    def inner(a, b):
      return sum(phi[i][j] * ai * bj
                 for i, ai in enumerate(a.numlist)
                 for j, bj in enumerate(b.numlist)
                )
  
    A = ZFilter(1)
    B = [z ** -1]
    beta = [inner(B[0], B[0])]
  
    m = 1
    while True:
      try:
        k = -inner(A, z ** -m) / beta[m - 1] # Last one is really a PARCOR coeff
      except ZeroDivisionError:
        raise ZeroDivisionError("Can't find next coefficient")
      if k >= 1 or k <= -1:
        raise ValueError("Unstable filter")
      A += k * B[m - 1]
  
      if m >= order:
        A.error = inner(A, A)
        return A
  
      gamma = [inner(z ** -(m + 1), B[q]) / beta[q] for q in xrange(m)]
      B.append(z ** -(m + 1) - sum(gamma[q] * B[q] for q in xrange(m)))
      beta.append(inner(B[m], B[m]))
      m += 1
docs: Find the Linear Predictive Coding (LPC) coefficients as a ZFilter object,
    the analysis whitening filter. This implementation is based on the
    covariance method, assuming a zero-mean stochastic process, finding
    the coefficients iteratively and greedily like the lattice implementation
    in Levinson-Durbin algorithm, although the lag matrix found from the given
    block don't have to be toeplitz. Slow, but this strategy don't need NumPy.