octave_packages/m/signal/durbinlevinson.m

   1 ## Copyright (C) 1995-2012 Friedrich Leisch
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {} durbinlevinson (@var{c}, @var{oldphi}, @var{oldv})
  21 ## Perform one step of the Durbin-Levinson algorithm.
  22 ##
  23 ## The vector @var{c} specifies the autocovariances @code{[gamma_0, @dots{},
  24 ## gamma_t]} from lag 0 to @var{t}, @var{oldphi} specifies the
  25 ## coefficients based on @var{c}(@var{t}-1) and @var{oldv} specifies the
  26 ## corresponding error.
  27 ##
  28 ## If @var{oldphi} and @var{oldv} are omitted, all steps from 1 to
  29 ## @var{t} of the algorithm are performed.
  30 ## @end deftypefn
  31
  32 ## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at>
  33 ## Description: Perform one step of the Durbin-Levinson algorithm
  34
  35 function [newphi, newv] = durbinlevinson (c, oldphi, oldv)
  36
  37   if (! ((nargin == 1) || (nargin == 3)))
  38     print_usage ();
  39   endif
  40
  41   if (columns (c) > 1)
  42     c = c';
  43   endif
  44
  45   newphi = 0;
  46   newv = 0;
  47
  48   if (nargin == 3)
  49
  50     t = length (oldphi) + 1;
  51
  52     if (length (c) < t+1)
  53       error ("durbinlevinson: C too small");
  54     endif
  55
  56     if (oldv == 0)
  57       error ("durbinlevinson: OLDV = 0");
  58     endif
  59
  60     if (rows (oldphi) > 1)
  61       oldphi = oldphi';
  62     endif
  63
  64     newphi = zeros (1, t);
  65     newphi(1) = (c(t+1) - oldphi * c(2:t)) / oldv;
  66     for i = 2 : t
  67       newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1);
  68     endfor
  69     newv = (1 - newphi(1)^2) * oldv;
  70
  71   elseif(nargin == 1)
  72
  73     tt = length (c)-1;
  74     oldphi = c(2) / c(1);
  75     oldv = (1 - oldphi^2) * c(1);
  76
  77     for t = 2 : tt
  78
  79       newphi = zeros (1, t);
  80       newphi(1) = (c(t+1) - oldphi * c(2:t)) / oldv;
  81       for i = 2 : t
  82         newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1);
  83       endfor
  84       newv = (1 - newphi(1)^2) * oldv;
  85
  86       oldv = newv;
  87       oldphi = newphi;
  88
  89     endfor
  90
  91   endif
  92
  93 endfunction