octave_packages/tsa-4.2.4/durlev.m

   1 function [MX,res,arg3] = durlev(AutoCov);
   2 % function  [AR,RC,PE] = durlev(ACF);
   3 % function  [MX,PE] = durlev(ACF);
   4 % estimates AR(p) model parameter by solving the
   5 % Yule-Walker with the Durbin-Levinson recursion
   6 % for multiple channels
   7 %  INPUT:
   8 % ACF   Autocorrelation function from lag=[0:p]
   9 %
  10 %  OUTPUT
  11 % AR    autoregressive model parameter
  12 % RC    reflection coefficients (= -PARCOR coefficients)
  13 % PE    remaining error variance
  14 % MX    transformation matrix between ARP and RC (Attention: needs O(p^2) memory)
  15 %        AR(:,K) = MX(:,K*(K-1)/2+(1:K));
  16 %        RC(:,K) = MX(:,(1:K).*(2:K+1)/2);
  17 %
  18 % All input and output parameters are organized in rows, one row
  19 % corresponds to the parameters of one channel
  20 %
  21 % see also ACOVF ACORF AR2RC RC2AR LATTICE
  22 %
  23 % REFERENCES:
  24 %  Levinson N. (1947) "The Wiener RMS(root-mean-square) error criterion in filter design and prediction." J. Math. Phys., 25, pp.261-278.
  25 %  Durbin J. (1960) "The fitting of time series models." Rev. Int. Stat. Inst. vol 28., pp 233-244.
  26 %  P.J. Brockwell and R. A. Davis "Time Series: Theory and Methods", 2nd ed. Springer, 1991.
  27 %  S. Haykin "Adaptive Filter Theory" 3rd ed. Prentice Hall, 1996.
  28 %  M.B. Priestley "Spectral Analysis and Time Series" Academic Press, 1981.
  29 %  W.S. Wei "Time Series Analysis" Addison Wesley, 1990.
  30
  31 %       $Id: durlev.m 5090 2008-06-05 08:12:04Z schloegl $
  32 %       Copyright (C) 1998-2002,2008 by Alois Schloegl <a.schloegl@ieee.org>
  33 %
  34 %    This program is free software: you can redistribute it and/or modify
  35 %    it under the terms of the GNU General Public License as published by
  36 %    the Free Software Foundation, either version 3 of the License, or
  37 %    (at your option) any later version.
  38 %
  39 %    This program is distributed in the hope that it will be useful,
  40 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  41 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  42 %    GNU General Public License for more details.
  43 %
  44 %    You should have received a copy of the GNU General Public License
  45 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  46
  47
  48 % Inititialization
  49 [lr,lc]=size(AutoCov);
  50
  51 res=[AutoCov(:,1), zeros(lr,lc-1)];
  52 d=zeros(lr,1);
  53
  54 if nargout<3         % needs O(p^2) memory
  55         MX=zeros(lr,lc*(lc-1)/2);
  56         idx=0;
  57         idx1=0;
  58         % Durbin-Levinson Algorithm
  59         for K=1:lc-1,
  60                 %idx=K*(K-1)/2;  %see below
  61                 % for L=1:lr, d(L)=arp(L,1:K-1)*transpose(AutoCov(L,K:-1:2));end;  % Matlab 4.x, Octave
  62                 % d=sum(MX(:,idx+(1:K-1)).*AutoCov(:,K:-1:2),2);              % Matlab 5.x
  63                 MX(:,idx+K)=(AutoCov(:,K+1)-sum(MX(:,idx1+(1:K-1)).*AutoCov(:,K:-1:2),2))./res(:,K);
  64                 %rc(:,K)=arp(:,K);
  65                 %if K>1   %for compatibility with OCTAVE 2.0.13
  66                         MX(:,idx+(1:K-1))=MX(:,idx1+(1:K-1))-MX(:,(idx+K)*ones(K-1,1)).*MX(:,idx1+(K-1:-1:1));
  67                 %end;
  68                 % for L=1:lr, d(L)=MX(L,idx+(1:K))*(AutoCov(L,K+1:-1:2).');end; % Matlab 4.x, Octave
  69                 % d=sum(MX(:,idx+(1:K)).*AutoCov(:,K+1:-1:2),2);              % Matlab 5.x
  70                 res(:,K+1) = res(:,K).*(1-abs(MX(:,idx+K)).^2);
  71                 idx1=idx;
  72                 idx=idx+K;
  73         end;
  74         %arp=MX(:,K*(K-1)/2+(1:K));
  75         %rc =MX(:,(1:K).*(2:K+1)/2);
  76
  77 else            % needs O(p) memory
  78
  79         arp=zeros(lr,lc-1);
  80         rc=zeros(lr,lc-1);
  81
  82         % Durbin-Levinson Algorithm
  83         for K=1:lc-1,
  84                 % for L=1:lr, d(L)=arp(L,1:K-1)*transpose(AutoCov(L,K:-1:2));end;  % Matlab 4.x, Octave
  85                 % d=sum(arp(:,1:K-1).*AutoCov(:,K:-1:2),2);              % Matlab 5.x
  86                 arp(:,K) = (AutoCov(:,K+1)-sum(arp(:,1:K-1).*AutoCov(:,K:-1:2),2))./res(:,K); % Yule-Walker
  87                 rc(:,K)  = arp(:,K);
  88                 %if K>1   %for compatibility with OCTAVE 2.0.13
  89                         arp(:,1:K-1)=arp(:,1:K-1)-arp(:,K*ones(K-1,1)).*arp(:,K-1:-1:1);
  90                 %end;
  91                 %for L=1:lr, d(L)=arp(L,1:K)*(AutoCov(L,K+1:-1:2).');end; % Matlab 4.x, Octave
  92                 % d=sum(arp(:,1:K).*AutoCov(:,K+1:-1:2),2);              % Matlab 5.x
  93                 res(:,K+1) = res(:,K).*(1-abs(arp(:,K)).^2);
  94         end;
  95
  96         % assign output arguments
  97         arg3=res;
  98         res=rc;
  99         MX=arp;
 100 end; %if