octave_packages/signal-1.1.3/levinson.m

   1 ## Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
   2 ## Copyright (C) 2006 Peter V. Lanspeary, <peter.lanspeary@.adelaide.edu.au>
   3 ##
   4 ## This program is free software; you can redistribute it and/or modify it under
   5 ## the terms of the GNU General Public License as published by the Free Software
   6 ## Foundation; either version 3 of the License, or (at your option) any later
   7 ## version.
   8 ##
   9 ## This program is distributed in the hope that it will be useful, but WITHOUT
  10 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  11 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  12 ## details.
  13 ##
  14 ## You should have received a copy of the GNU General Public License along with
  15 ## this program; if not, see <http://www.gnu.org/licenses/>.
  16
  17 ## usage:  [a, v, ref] = levinson (acf [, p])
  18 ##
  19 ## Use the Durbin-Levinson algorithm to solve:
  20 ##    toeplitz(acf(1:p)) * x = -acf(2:p+1).
  21 ## The solution [1, x'] is the denominator of an all pole filter
  22 ## approximation to the signal x which generated the autocorrelation
  23 ## function acf.
  24 ##
  25 ## acf is the autocorrelation function for lags 0 to p.
  26 ## p defaults to length(acf)-1.
  27 ## Returns
  28 ##   a=[1, x'] the denominator filter coefficients.
  29 ##   v= variance of the white noise = square of the numerator constant
  30 ##   ref = reflection coefficients = coefficients of the lattice
  31 ##         implementation of the filter
  32 ## Use freqz(sqrt(v),a) to plot the power spectrum.
  33 ##
  34 ## REFERENCE
  35 ## [1] Steven M. Kay and Stanley Lawrence Marple Jr.:
  36 ##   "Spectrum analysis -- a modern perspective",
  37 ##   Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981
  38
  39 ## Based on:
  40 ##    yulewalker.m
  41 ##    Copyright (C) 1995 Friedrich Leisch <Friedrich.Leisch@ci.tuwien.ac.at>
  42 ##    GPL license
  43
  44 ## TODO: Matlab doesn't return reflection coefficients and
  45 ## TODO:    errors in addition to the polynomial a.
  46 ## TODO: What is the difference between aryule, levinson,
  47 ## TODO:    ac2poly, ac2ar, lpc, etc.?
  48
  49 function [a, v, ref] = levinson (acf, p)
  50   if ( nargin<1 )
  51     print_usage;
  52   elseif( ~isvector(acf) || length(acf)<2 )
  53     error( "levinson: arg 1 (acf) must be vector of length >1\n");
  54   elseif ( nargin>1 && ( ~isscalar(p) || fix(p)~=p ) )
  55     error( "levinson: arg 2 (p) must be integer >0\n");
  56   else
  57     if ((nargin == 1)||(p>=length(acf))) p = length(acf) - 1; endif
  58     if( columns(acf)>1 ) acf=acf(:); endif      # force a column vector
  59
  60     if nargout < 3 && p < 100
  61       ## direct solution [O(p^3), but no loops so slightly faster for small p]
  62       ##   Kay & Marple Eqn (2.39)
  63       R = toeplitz(acf(1:p), conj(acf(1:p)));
  64       a = R \ -acf(2:p+1);
  65       a = [ 1, a.' ];
  66       v = real( a*conj(acf(1:p+1)) );
  67     else
  68       ## durbin-levinson [O(p^2), so significantly faster for large p]
  69       ##   Kay & Marple Eqns (2.42-2.46)
  70       ref = zeros(p,1);
  71       g = -acf(2)/acf(1);
  72       a = [ g ];
  73       v = real( ( 1 - g*conj(g)) * acf(1) );
  74       ref(1) = g;
  75       for t = 2 : p
  76         g = -(acf(t+1) + a * acf(t:-1:2)) / v;
  77         a = [ a+g*conj(a(t-1:-1:1)), g ];
  78         v = v * ( 1 - real(g*conj(g)) ) ;
  79         ref(t) = g;
  80       endfor
  81       a = [1, a];
  82     endif
  83   endif
  84 endfunction