home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / signal / durbinlevinson.m < prev    next >
Text File  |  1997-02-19  |  2KB  |  90 lines
  1. ## Copyright (C) 1995  Friedrich Leisch
  2. ## 
  3. ## This program is free software; you can redistribute it and/or modify
  4. ## it under the terms of the GNU General Public License as published by
  5. ## the Free Software Foundation; either version 2, or (at your option)
  6. ## any later version.
  7. ## 
  8. ## This program is distributed in the hope that it will be useful, but
  9. ## WITHOUT ANY WARRANTY; without even the implied warranty of
  10. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11. ## General Public License for more details. 
  12. ## 
  13. ## You should have received a copy of the GNU General Public License
  14. ## along with this file.  If not, write to the Free Software Foundation,
  15. ## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  16.  
  17. ## usage:  durbinlevinson (c, [oldphi, oldv])
  18. ##
  19. ## Performs one step of the Durbin-Levinson algorithm.
  20. ##
  21. ## The vector c_t = [gamma_0, ..., gamma_t] contains the autocovariances
  22. ## from lag 0 to t, oldphi the coefficients based on c_(t-1) and oldv
  23. ## the corresponding error.
  24. ##
  25. ## If oldphi is omitted, all steps from 1 to t of the algorithm are
  26. ## performed.
  27.    
  28. ## Author:  FL <Friedrich.Leisch@ci.tuwien.ac.at>
  29. ## Description:  Perform one step of the Durbin-Levinson algorithm
  30.   
  31. function [newphi, newv] = durbinlevinson (c, oldphi, oldv)
  32.   
  33.   if( !((nargin == 1) || (nargin == 3)) )
  34.     usage ("durbinlevinson (c, [oldphi, oldv])");
  35.   endif
  36.   
  37.   if( columns (c) > 1 )
  38.     c=c';
  39.   endif
  40.  
  41.   newphi = 0;
  42.   newv = 0;
  43.   
  44.   if (nargin == 3)
  45.     
  46.     t = length (oldphi) + 1;
  47.   
  48.     if (length (c) < t+1)
  49.       error ("durbilevinson:  c too small");
  50.     endif
  51.   
  52.     if (oldv == 0)
  53.       error ("durbinlevinson: oldv = 0");
  54.     endif
  55.     
  56.     if (rows (oldphi) > 1 )
  57.       oldphi = oldphi';
  58.     endif
  59.     
  60.     newphi = zeros (1, t);
  61.     newphi(1) = ( c(t+1) - oldphi * c(2:t) ) / oldv;
  62.     for i = 2 : t
  63.       newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1);
  64.     endfor
  65.     newv = ( 1 - newphi(1)^2 ) * oldv;
  66.     
  67.   elseif(nargin == 1)
  68.     
  69.     tt = length (c)-1;
  70.     oldphi = c(2) / c(1);
  71.     oldv = ( 1 - oldphi^2 ) * c(1);
  72.     
  73.     for t = 2 : tt
  74.     
  75.       newphi = zeros (1, t);
  76.       newphi(1) = ( c(t+1) - oldphi * c(2:t) ) / oldv;
  77.       for i = 2 : t
  78.     newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1);
  79.       endfor
  80.       newv = ( 1 - newphi(1)^2 ) * oldv;
  81.       
  82.       oldv = newv;
  83.       oldphi = newphi;
  84.     
  85.     endfor
  86.     
  87.   endif
  88.  
  89. endfunction
  90.