Scalar Levinson Algorithm and Autocorrelation in MATLAB

Hello,

I am writitng a script to solve the levinson recursion equations and am having trouble with formatting the indices.

These are the equations I am using typed into word:

Here is my code:

data = csvread('RV_sequence2.csv');
dataT = data(:,1);
dataD = data(:,2);
[n,dat] = size(data);
numlag = 10;

% Autocorrelation Sequence - this is correct
for ii = 1:numlag
    lag = ii - 1;
        for iii = 1:n-lag
            seq1(iii,ii) = dataD(iii)* dataD(iii+lag);         
            numer = sum(seq1);
        end
        autocorr = numer/n;
end

% initializing
a(1,1) = -autocorr(2)/autocorr(1);
P(1,1) = autocorr(1)*(-(autocorr(2))^2/autocorr(1));

% Recursion Sequence
for N = 1:4
%     lag = N -1;
for i = 1:n
    alpha(N) = autocorr(N+1) + sum(a(N)*autocorr(N+1-1));
    a(N+1) = a(N,i) - (alpha(N)/P(N))*a(N,N+1-i);
end
    P(N+1) = P(N) - (alpha(N))^2/P(N);
    a(N+1) = - alpha(N)/P(N);
end

Answer 1

[John Stanfield]

I was able to solve the indices problem and correctly solve the levinson recursion equations. Posting just in case someone needs this in the future!

% Initialization
P = zeros(1,3);
alpha = zeros(1,3);
a(1,1) = -autocorr(2)/autocorr(1);
P(1) = autocorr(1);
P(2) = autocorr(1) - ((autocorr(2)^2)/autocorr(1));



% Start Algorithm
for N = 1:numlag-2
    sum = 0;
    for j = 1:N
        sum = sum + a(N,j)*autocorr(N+2-j);
    end
    alpha(N,1) = autocorr(N+2) + sum;
for i = 1:N
    a(N+1,i) = a(N,i) - (alpha(N)/P(N+1))*a(N,N+1-i);
end
    a(N+1,i) = a(N,i) - (alpha(N)/P(N+1))*a(N,N+1-i);
    a(N+1,N+1) = - alpha(N)/P(N+1);
    P(N+2) = P(N+1) - (alpha(N))^2/P(N+1);
    

% Tolerance Checks    
    zerochecker = abs(P(2) - P(1));
    if zerochecker < 0.5 
        N = 0;
        break
    end
    check = abs(abs(P(N+2)) - abs(P(N+1)));
    if  check < 1
        break
    end
end

As a BONUS I am posting my autocorrelation sequence below:

function [autocorr] = corr(numlag,dataD,n)
% autocorrelation sequence
for ii = 1:numlag
    lag = ii - 1;
        for iii = 1:n-lag
            seq1(iii,ii) = dataD(iii).* dataD(iii+lag);         
            numer = sum(seq1);
        end
        autocorr_unrounded = round(numer,2)';
end
autocorr = autocorr_unrounded/n;
end