octave_packages/tsa-4.2.4/arcext.m

   1 function  [AR,RC] = arcext(MX,P);
   2 % ARCEXT extracts AR and RC of order P from Matrix MX
   3 % function  [AR,RC] = arcext(MX,P);
   4 %
   5 %  INPUT:
   6 % MX    AR and RC matrix calculated by durlev
   7 % P     model order (default maximum possible)
   8 %
   9 %  OUTPUT
  10 % AR    autoregressive model parameter
  11 % RC    reflection coefficients (= -PARCOR coefficients)
  12 %
  13 % All input and output parameters are organized in rows, one row
  14 % corresponds to the parameters of one channel
  15 %
  16 % see also ACOVF ACORF DURLEV
  17 %
  18 % REFERENCES:
  19 %  P.J. Brockwell and R. A. Davis "Time Series: Theory and Methods", 2nd ed. Springer, 1991.
  20 %  S. Haykin "Adaptive Filter Theory" 3rd ed. Prentice Hall, 1996.
  21 %  M.B. Priestley "Spectral Analysis and Time Series" Academic Press, 1981.
  22 %  W.S. Wei "Time Series Analysis" Addison Wesley, 1990.
  23
  24 %  $Id: arcext.m 9609 2012-02-10 10:18:00Z schloegl $
  25 %  Copyright (C) 1998-2003,2008,2012 by Alois Schloegl <alois.schloegl@ist.ac.at>
  26 %       This is part of the TSA-toolbox. See also
  27 %       http://pub.ist.ac.at/~schloegl/matlab/tsa/
  28 %       http://octave.sourceforge.net/
  29 %       http://biosig.sourceforge.net/
  30 %
  31 %    This program is free software: you can redistribute it and/or modify
  32 %    it under the terms of the GNU General Public License as published by
  33 %    the Free Software Foundation, either version 3 of the License, or
  34 %    (at your option) any later version.
  35 %
  36 %    This program is distributed in the hope that it will be useful,
  37 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  38 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  39 %    GNU General Public License for more details.
  40 %
  41 %    You should have received a copy of the GNU General Public License
  42 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  43
  44
  45 [lr,lc]=size(MX);
  46
  47 if ~mod(sqrt(8*lc+1)-1,2)       % full number of elements
  48         K = (sqrt(8*lc+1)-1)/2;
  49 elseif ~mod(lc,2),   % compressed form of MX
  50         K = lc/2;
  51 else            % invalid number of elements
  52         fprintf(2,'Warning ARCEXT: Number of elements is different than a triangular matrix\n');
  53 end;
  54
  55 if (K~=P) && (lc~=K*(K+1)/2),
  56         [AR,RC,PE]=rc2ar(MX(:,(1:P).*(2:P+1)/2));
  57 else
  58         AR = MX(:,P*(P-1)/2+(1:P));
  59         RC = MX(:,(1:P).*(2:P+1)/2);
  60 end;