durlev.m

时间序列工具箱

function [ar,rc,PE,d] = durlev(ACF);
% function  [ar,rc,PE] = durlev(ACF);
% estimates AR(p) model parameter by solving the
% Yule-Walker with the Durbin-Levinson recursion
% for multiple channels
%  INPUT:
% ACF	Autocovariance function [0..P]
%
%  OUTPUT
% ar	autoregressive model parameter	
% rc    reflection coefficients
% PE    remaining error variance
%
% Reference:
% Brockwell P.J. and Davis R. A. "Time Series: Theory and Methods", 2nd ed. Springer, 1991

%	Version 2.42
%	17.05.1998
%	Copyright (c) 1996-1998 by  Alois Schloegl
%	[e-mail]	

% This library is free software; you can redistribute it and/or
% modify it under the terms of the GNU Library General Public
% License as published by the Free Software Foundation; either
% Version 2 of the License, or (at your option) any later version.
%
% This library is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% Library General Public License for more details.
%
% You should have received a copy of the GNU Library General Public
% License along with this library; if not, write to the
% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
% Boston, MA  02111-1307, USA.


[lr,lc]=size(ACF);

% Inititialization
ar=zeros(lr,lc-1);
rc=zeros(lr,lc-1);
PE=ACF;	
d=zeros(lr,1);
if lc<2 return;end;

% Durbin-Levinson Algorithm
K=1; % introduced, because (*) line does not work in Octave 2.0.12
        % for L=1:lr, d(L)=ar(L,1:K-1)*ACF(L,K:-1:2)';end;
        ar(:,K)=(ACF(:,K+1)-d)./PE(:,K);
        rc(:,K)=ar(:,K);
        PE(:,K+1) = PE(:,K).*(1-ar(:,K).*ar(:,K));
        %(*) ar(:,1:K-1)=ar(:,1:K-1)-ar(:,K*ones(K-1,1)).*ar(:,K-1:-1:1);
for K=2:lc-1,
        for L=1:lr, d(L)=ar(L,1:K-1)*ACF(L,K:-1:2)';end;
        ar(:,K)=(ACF(:,K+1)-d)./PE(:,K);
        rc(:,K)=ar(:,K);
        PE(:,K+1) = PE(:,K).*(1-ar(:,K).*ar(:,K));
        ar(:,1:K-1)=ar(:,1:K-1)-ar(:,K*ones(K-1,1)).*ar(:,K-1:-1:1);
end;