arma_last.m

针对ARMA的仿真源程序 对P, Q阶的求解

%--------------------------------------------------------------------------
%
%   ARMA_LAST.M
%
%   This function provides an ARMA spectral estimate which is maximum
% entropy satisfying correlation constraint (number of poles)and cepstrum 
% constrains (number of ceros)
%
%   [sigma,outn,outd]=arma_last(xinf,np,nq)
%
% The function requires:
%
%   xinf(nsam)  Input signal (real valued) of nsam samples
%   nq          Order of the denominator (# of poles)
%   np          Order of the numerator (# of ceros)
%
% The function provides:
%
%   outd(nq)    The numerator coefficients outn(1)=1
%   outn(np)    The denominator coefficients outd(1)=1
%   sigma       The square-root of input noise power
%
% The function requires spa_corc.m to compute the autocorrealtion matrix
% estimate
%
%   Miguel Angel LAGUNAS and Petre STOICA                June 2007
%
%--------------------------------------------------------------------------
function [sigma,outn,outd]=arma_last(xinf,np,nq);
nsam=length(xinf);ra1=spa_corc(xinf,nq);nle=(nsam/2);
ra2=spa_corc(xinf,nle);[u d]=eig(ra2);
d1=diag(log(diag(d)));ca1=u*d1*u';clear u d1 d ra2 nsam;
for i=1:np
  cuso(i)=0;
  for j=1:nle-i;
    cuso(i)=cuso(i)+ca1(j,i+j);
  end;cuso(i)=cuso(i)/(nle);
end;
cep=(1:1:np).*cuso(1:np);clear ca1 nle cuso;
i=2;h(1)=1;
while i<np+0.5;
  h(i)=0;
  for k=1:i-1;
    h(i)=h(i)+h(k)*cep(i-k);
  end;
  h(i)=h(i)/(i-1);i=i+1;
end;
hh=flipud(hankel(fliplr(h)));clear h cep;
if nq>np;
   hh=[hh zeros(np,nq-np)];
else;
   hh=hh(1:np,1:nq);
end;
[u d]=eig(ra1,hh'*hh);sas=diag(d);ua=[];sasn=[];
for j=1:nq;
  [auso luso]=min(sas);sasn=[sasn;sas(luso)];
  sas(luso)=10e9;ua=[ua u(:,luso)];
end;
outd=real(ua(:,1)');outd=outd/outd(1);sigma=1/sqrt(sasn(1));
hha=hh/sigma;outn=outd*hha';sigma=sigma*outn(1);outn=outn/outn(1);
clear hh auso luso sas sasn u d ua ra1;
%---------------------------------------------------------------------