getarma.m

matlab算法集 matlab算法集

function theta = getarma (u,y,n,m)
%-----------------------------------------------------------------------
% Usage:       theta = arma (u,y,n,m)
%
% Description: Identify the parameters an auto-regressive moving-
%              average (ARMA) model having the following transfer
%              function:
%
%              Y(z)      theta(n+1) + ... + theta(n+m+1)*z^(-m)
%              ---- = -------------------------------------------
%              U(z)   1 + theta(1)*z^(-1) + ... + theta(n)*z^(-n)   
%
% Inputs:      u = p by 1 vector containing input samples
%              y = p by 1 vector containing output samples
%              n = number of past outputs used (n >= 0)
%              m = number of past inputs used (m >= 0)
%
% Outputs      theta = (n+m+1) by 1 vector containing ARMA model
%                      parameters.
%            
% Notes:       1. The input u must be rich in frequency content:
%                 the magnitude spectrum must be nonzero at at
%                 least p frequencies 
%              2. Use the function arma to compute the output of
%                 ARMA model. */
%-----------------------------------------------------------------------
   
% Initialize

   chkvec (u,1,'getarma');
   chkvec (y,2,'getarma');
   n = args (n,0,n,3,'getarma');
   m = args (m,0,m,4,'getarma');
   e = 0;
   q = n + m + 1;
   p = length(u);
   if length(y) ~= p
      disp ('Arguments 1 and 2 of getarma must be of the same length.')
      return
   end
   if p < q
      disp ('Insufficient number of samples in getarma.')
      return
   end 
   X = zeros (p,q);
  
% Compute state matrix X

   for i = 1 : p
      for j = 1 : n
         if j < i
            X(i,j) = -y(i-j);
         end
      end 
      for j = 0 : m
         if j < i
            X(i,n+1+j) = u(i-j);
         end
      end
   end

% Compute the least-squares estimate, theta  

   if rank(X) < q
      disp ('Insufficient frequency content in input signal for getarma.')
      return
   end
   theta = gauss (X,y);

% Compute error

   for i = 1 : p
      r = y(i);	
      for j = 1 : q
         r = r - X(i,j)*theta(j);
      end
      e = e + r*r;
   end
%-----------------------------------------------------------------------