www.gusucode.com > signal 工具箱matlab源码程序 > signal/private/arparest.m

    function [a,e,msg,msgobj] = arparest( x, p, method)
%ARPAREST   AR parameter estimation via a specified method.
%   A = ARPAREST(X,ORDER,METHOD) returns the polynomial A corresponding to 
%   the AR parametric signal model estimate of vector X using the specified
%   METHOD.  ORDER is the model order of the AR system.
%
%   Supported methods are: 'covariance' and 'modified' although all of the
%   methods of CORRMTX will work. In particular if 'autocorrelation' is
%   used, the results should be the same as those of ARYULE (but slower).
%
%   [A,E] = ARPAREST(...) returns the variance estimate E of the white noise
%   input to the AR model.

%   Ref: S. Kay, MODERN SPECTRAL ESTIMATION,
%              Prentice-Hall, 1988, Chapter 7
%        S. Marple, DIGITAL SPECTRAL ANALYSIS WITH APPLICATION,
%              Prentice-Hall, 1987, Chapter 8.
%        P. Stoica and R. Moses, INTRODUCTION TO SPECTRAL ANALYSIS,
%              Prentice-Hall, 1997, Chapter 3

%   Author(s): R. Losada and P. Pacheco
%   Copyright 1988-2014 The MathWorks, Inc.

narginchk(3,3)

% Initialize in case we return early
a = []; e = [];

% Assign msg in case there are no errors
msg ='';
msgobj = [];

% enforce backwards compatibility with row vectors
if isvector(x)
  x = x(:);
end

validateattributes(x,{'numeric'},{'nonempty','finite','2d'},'arpest','X');

% Set up necessary but not sufficient conditions for the correlation
% matrix to be nonsingular. From (Marple)
switch method,
case 'covariance',
   minlength_x = 2*p;
case 'modified',
   minlength_x = 3*p/2;
otherwise
   msgobj = message('signal:arparest:UnknMethod');
   msg = getString(msgobj);
   return
end

% Do some data sanity testing
if size(x,1) < minlength_x
   if strcmp(method, 'modified')
       multiplier = '3/2';
    else
       multiplier = '2';
    end
    if isvector(x)
       msgID = 'signal:arparest:VectorTooSmallForModel';
    else
       msgID = 'signal:arparest:MatrixTooSmallForModel';
    end
    msgobj = message(msgID,'X',multiplier);
    msg = getString(msgobj);
    return
end
if issparse(x),
   msgobj = message('signal:arparest:InputSignalCannotBeSparse');
   msg = getString(msgobj);
   return
end
if isempty(p) || p ~= round(p),
   msgobj = message('signal:arparest:ModelOrderMustBeInteger');
   msg = getString(msgobj);
   return
end

% 'like' will also copy over the complexity we only want the class
a = zeros(p+1,size(x,2),class(x)); %#ok<ZEROLIKE>
e = zeros(1,size(x,2),class(x)); %#ok<ZEROLIKE>

for chan=1:size(x,2)
   % Generate the appropriate data matrix
   XM = corrmtx(x(:,chan),p,method);
   Xc = XM(:,2:end);
   X1 = XM(:,1);

   % Coefficients estimated via the covariance method
   a(:,chan) = [1; -Xc\X1];

   % Estimate the input white noise variance
   Cz = X1'*Xc;
   e(:,chan) = X1'*X1 + Cz*a(2:end,chan);

   % Ignore the possible imaginary part due to numerical errors and force
   % the variance estimate of the white noise to be positive
   e(:,chan) = abs(real(e(:,chan)));  
end

a = a.'; % By convention all polynomials are row vectors

% [EOF] arparest.m