armals.m

统计自适应信号处理的matlab程序

function [a,b,V,FPE,C]=armals(x,p,q,maxiter,tol)
% function [a,b,V,FPE,C]=armals(x,p,q,maxiter,tol)
% ARMA modeling using non-linear LS optimization
%  
% Programmed by: Dimitris Manolakis
% Revised: 7/1/94
%
%-----------------------------------------------------------
% Copyright 2000, by Dimitris G. Manolakis, Vinay K. Ingle,
% and Stephen M. Kogon.  For use with the book
% "Statistical and Adaptive Signal Processing"
% McGraw-Hill Higher Education.
%-----------------------------------------------------------


Lx=length(x);
pqmax=max([p q]);

%tol=0.01;maxiter=10;

% Compute initial AR(p+q) estimate, residuals,
% and cost function
 
[apq,ear,Vpq,FPEpq]=arls(x,p+q);

% Compute initial ARMA(p,q) parameters using
% equation error LS with input ear and output x

[a,b,V,FPE]=pzls(x,ear,p,q);

% Convert b into its minimum phase equivalent
% and then compute the prediction error

b=miphapol(b); b=b(:);
e=pefilter(a,b,x,zeros(pqmax,1));
V2=(e'*e)/(Lx-pqmax);
V=V2;


% Start iterative minimization
%====================================================

% Initialization
%----------------------------------------------------

g=ones(p+q,1);
l=0;
st=0;

% Minimization loop
%----------------------------------------------------
while [norm(g)>tol l<maxiter st~=1]

l=l+1;

% Compute gradient
%----------------------------------------------------

yf=filter(-1,b,x);
ef=filter(1,b,e);

R=zeros(p+q,p+q);
F=zeros(p+q,1);

for n=pqmax+1:Lx
  psi=zeros(p+q,1);
  for k=1:p, psi(k)=yf(n-k); end
  for k=1:q, psi(k+p,:)=ef(n-k); end
  R=R+psi*psi'; F=F+psi*e(n);
end

g=R\F;

% Search along gradient direction
%----------------------------------------------------

[a1,b1,e,V1,st]=gstep(x,a,b,g,V,p,q); 

a=a1;
b=b1;
V=V1;

end

% End iterative minimization
%====================================================

V=(e'*e)/(Lx-pqmax);

FPE=V*(1+(p+q)/Lx)/(1-(p+q)/Lx);

% Covariance error matrix

C=V*inv(R);