%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 線形予測分析
%
% 離散入力からの線形予測係数のパラメータ推定
%
%                                                      2003.12.17 Katto Lab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

% 変数設定
%% 信号数
N = 1000;
%% 信号ベクトル
y = zeros(N,1);
%% 線形予測係数 (推定目標)
a = [-0.9 0.6 -0.5 0.45 0.32 -0.24 -0.15 0.1 0.08 0.02];
% a = [-0.8 0.6 -0.4];
p = length(a);


% p次のARモデル入力生成
for n=1:N
    tmp = 0;
    for i=1:length(a)
        if((n-i) <= 0) continue;
        else
            tmp = tmp -a(i) * y(n-i);
        end;
    end;
    y(n) = tmp + randn;
end;
figure(1);
plot(1:N,y);

% 共分散法によるパラメータ推定
t1 = cputime;
V = xcov(y,y);
for i=1:p
    v(i) = V(N+i-1);
    r1(i) = V(N+i);
end;
for i=1:p
    for j=1:p
        R1(i,j) = v(abs(i-j)+1);
    end;
end;
A1 = -inv(R1)*r1';
t1 = cputime - t1;

% コレスキー分解
t2 = cputime;
[T0,T1] = lu(R1);
c = -inv(T0)*r1';
A2 = inv(T1)*c;
t2 = cputime - t2;

% 自己相関法によるパラメータ推定
t3 = cputime;
[X,RR] = corrmtx(y,p);
for i=1:p
    for j=1:p
        R3(i,j) = RR(i,j);
    end;
end;
for i=1:p
    r3(i) = RR(i+1,1);
end;
A3 = -inv(R3)*r3';
t3 = cputime - t3;

% Levinson-Durbinアルゴリズム
t4 = cputime;
%LD = lpc(y,p);
%for i=1:p
%    A4(i) = real(LD(i+1));
%end;
for i=1:p+1
    s(i) = RR(1,i);
end;
sigma = s(1);
A4 = zeros(p,1);
A4tmp = zeros(p,1);
for m=1:p
    tmp = 0.0;
    for k=1:m-1
        tmp = tmp + A4(k)*s(m-k+1);
    end;
    d = s(m+1) + tmp;
    rho = - d/sigma;
    A4_tmp(m) = rho;
    sigma = sigma + rho * d;
    for k=1:m-1
        A4_tmp(k) = A4(k) + rho*A4(m-k);
    end;
    for k=1:m
        A4(k) = A4_tmp(k);
    end;
end;
t4 = cputime - t4;


% パラメータ推定結果の比較
figure(2);
plot(1:length(a), a,  'ko:'); hold on;
plot(1:length(a), A1, 'rx:'); hold on;
plot(1:length(a), A2, 'g+:'); hold on;
plot(1:length(a), A3, 'b*:'); hold on;
plot(1:length(a), A4, 'cs:'); hold on;

% 計算時間の比較
figure(3);
T = [t1, t2, t3, t4];
plot(1:4, T, 'o-'); 
