MATLAB: Scalar Levinson Algorithm and Autocorrelation in MATLAB

levinsonrecursion

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

Best Answer

  • 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