1

I'm trying wrap my head around linear prediction and figured I'd code up a basic example in Python to test my understanding. The idea behind linear predictive coding is to estimate future samples of a signal based on linear combinations of past samples.

I'm using the lpc module in scikits.talkbox so I don't have to write any of the algorithm myself. Here's my code:

import math
import numpy as np
from scikits.talkbox.linpred.levinson_lpc import levinson, acorr_lpc, lpc

x = np.linspace(0,11,12)

order = 5
"""
a = solution of the inversion
e = prediction error
k = reflection coefficients
"""

(a,e,k) = lpc(x,order,axis=-1)
recon = []

for i in range(order,len(x)):
    sum = 0
    for j in range(order):
        sum += -k[j]*x[i-j-1]
    sum += math.sqrt(e)
    recon.append(sum)

print(recon) 
print(x[order:len(x)])

which gives an output of

[5.618790615323507, 6.316875690307965, 7.0149607652924235, 
7.713045840276882, 8.411130915261339, 9.109215990245799, 9.807301065230257, 
10.505386140214716]
[ 4.  5.  6.  7.  8.  9. 10. 11.]

My concern is that I'm implementing this incorrectly somehow because I figured that if my input array is a linear signal, it should have no issue predicting future values based on past values. However, it does seem to have a particularly high error, especially for the first few values. Would anyone be able to tell me if I'm implementing this correctly or point me to a few examples where this is done in Python? Any help is greatly appreciated, thanks!

1 Answer 1

2

Linear prediction algorithm extends the original sequence with infinite amount of zeros in both directions. So, unless your input signal is constant zero, the extended sequence is not linear and you should expect a nonzero error. Here is my Python implementation:

def lpc(y, m):
    "Return m linear predictive coefficients for sequence y using Levinson-Durbin prediction algorithm"
    #step 1: compute autoregression coefficients R_0, ..., R_m
    R = [y.dot(y)] 
    if R[0] == 0:
        return [1] + [0] * (m-2) + [-1]
    else:
        for i in range(1, m + 1):
            r = y[i:].dot(y[:-i])
            R.append(r)
        R = np.array(R)
    #step 2: 
        A = np.array([1, -R[1] / R[0]])
        E = R[0] + R[1] * A[1]
        for k in range(1, m):
            if (E == 0):
                E = 10e-17
            alpha = - A[:k+1].dot(R[k+1:0:-1]) / E
            A = np.hstack([A,0])
            A = A + alpha * A[::-1]
            E *= (1 - alpha**2)
        return A

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.