Levinson durbin recursion for LPC calculation

levinson.jl

function levinson(R,L)
    a = zeros(L,L)
    P = zeros(1,L)

    # for m = 1
    a[1,1] = -R[2]/R[1]
    P[1] = R[1]*(1-a[1,1]^2)

    # for m = 2,3,4,..L
    for m = 2:L
        a[m,m] = (-(R[m+1] + dot(vec(a[m-1,1:m-1]),R[m:-1:2]))/P[m-1])
        a[m,1:m-1] = a[m-1,1:m-1]+a[m,m]*a[m-1,m-1:-1:1]
        P[m] = P[m-1]*(1-a[m,m]^2)
    end
    [1 a[L,:]], P[L]
end

function lpc(sig,L)
    R = xcorr(sig,sig)[length(sig):]
    levinson(R,L)
end

sig =[1:20]
L = 5
a,P = lpc(sig,L)
print(a)
print(P)