home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / AR.M next >
Encoding:
Text File  |  1993-03-11  |  4.3 KB  |  117 lines
  1. function [th,ref]=ar(y,n,approach,win,maxsize,T)
  2. %AR    Computes AR-models of signals using various approaches
  3. %
  4. %    TH = ar(y,n)  or  TH = ar(y,n,approach)  or TH = ar(y,n,approach,win)
  5. %    
  6. %    TH: returned as the estimated parameters of the AR-model,see HELP THETA
  7. %    y: The time series to be modelled, a column vector
  8. %    n: The order of the AR-model
  9. %    approach: The method used, one of the following ones:
  10. %       'fb' : The forward-backward approach (default)
  11. %       'ls' : The Least Squares method
  12. %       'yw' : The Yule-Walker method
  13. %       'burg': Burg's method
  14. %       'gl' : A geometric lattice method
  15. %       For the two latter ones, reflection coefficients and loss functions
  16. %       are returned in REFL by [TH,REFL]=ar(y,n,approach)
  17. %       If any of these arguments end with a zero (like 'burg0'), the 
  18. %       computation of the covariance is suppressed.
  19. %    win : Windows employed, one of the following ones:
  20. %       'now' : No windowing (default, except when approach='yw')
  21. %       'prw' : Prewindowing
  22. %       'pow' : Postwindowing
  23. %       'ppw' : pre- and post-windowing
  24. %    for TH = ar(y,n,approach,win,maxsize,T) see HELP AUXVAR
  25.  
  26. %    L. Ljung 10-7-87
  27. %    Copyright (c) 1987-90 by the MathWorks, Inc.
  28. %    All Rights Reserved.
  29.  
  30. % Some initial tests on the input arguments
  31. [Ncap,ny]=size(y);Ncap0=Ncap;
  32. maxsdef=idmsize(Ncap,n);
  33. if nargin<6,T=[];end
  34. if nargin<5,maxsize=[];end
  35. if nargin<4,win=[];end
  36. if nargin<3,approach=[];end
  37. if isempty(T),T=1;end,if T<0,T=1;end
  38. if isempty(maxsize), maxsize=maxsdef;end,if maxsize<0,maxsize=maxsdef;end
  39. if isempty(win);win='now';end ,if win<0,win='now';end
  40. if isempty(approach),approach='fb';end,if approach<0,approach='fb';end
  41. pt=1; %pt=1 means that covariances should be computed
  42. if approach(length(approach))=='0',pt=0;end
  43. if length(approach)>1,appr=approach(1:2);else appr=[approach ' '];end
  44. ap=[appr;appr;appr;appr;appr;appr;appr;appr;appr;appr];
  45. ss=(ap==['fb';'FB';'yw';'YW';'ls';'LS';'bu';'BU';'gl';'GL']);
  46. if sum(ss(:,1).*ss(:,2))==0,error('The input argument APPROACH must be one of the following (within quotes): ''fb'', ''ls'',''yw'', ''burg'', ''gl'''),end
  47. ap=[win;win;win;win]; ss=(ap==['prw';'pow';'ppw';'now']);
  48. if sum((ss(:,1).*ss(:,2)).*ss(:,3))==0,error('The input argument WIN must be one of the following (within quotes): ''now'',     ''prw'', ''pow'', ''ppw'''),end
  49. [Ncap,ny]=size(y); if ny>Ncap, y=y';end
  50. if min(Ncap,ny)>1, error('Only single time series can be handled'),end
  51. nstart=1;nend=Ncap;
  52. if appr=='yw' | appr=='YW',win='ppw';end
  53. if win=='prw' | win=='ppw',y=[zeros(n,1);y];nstart=n+1;nend=length(y);end
  54. if win=='pow' | win=='ppw',y=[y;zeros(n,1)];nend=length(y)-n;end
  55.  
  56. [Ncap,ny]=size(y);
  57.  
  58. %build up the theta-matrix
  59. % 
  60. I=[0 T 0 n 0 0];
  61. if pt~=0,th=zeros(n+3,max(7,n));else th=zeros(3,max(7,n));end
  62. th(1,1:6)=I;
  63. ti=fix(clock);ti(1)=ti(1)/100;
  64. th(2,2:6)=ti(1:5);th(2,7)=12;
  65.  
  66. % First the lattice based algorithms
  67.  
  68. if appr=='bu' | appr=='BU' | appr=='gl' | appr=='GL'
  69.    ef=y;eb=y;
  70.    rho(1)=y'*y/Ncap;
  71.    for p=1:n
  72.    nef=ef(p+1:Ncap)'*ef(p+1:Ncap); neb=eb(p:Ncap-1)'*eb(p:Ncap-1);
  73.    if appr=='gl' | appr=='GL',den=sqrt(nef*neb);end
  74.    if appr=='bu' | appr=='BU',den=(nef+neb)/2;end
  75.    r(p)=(-eb(p:Ncap-1)'*ef(p+1:Ncap))/den;
  76.    efold=ef;
  77.    ef(2:Ncap)=ef(2:Ncap)+r(p)*eb(1:Ncap-1);
  78.    A(p)=r(p);
  79.    eb(2:Ncap)=eb(1:Ncap-1)+conj(r(p))*efold(2:Ncap);
  80.    A(1:p-1)=A(1:p-1)+r(p)*conj(A(p-1:-1:1));
  81.    rho(p+1)=rho(p)*(1-r(p)*r(p));
  82.    end
  83. th(3,1:n)=A; 
  84. e=pe(y(nstart:nend),th);lam=e'*e/(nend-nstart+1-n);
  85. th(1,1)=lam; th(2,1)=lam*(1+n/Ncap0)/(1-n/Ncap0);
  86. ref=[0 r;rho];
  87. if pt==0, return,end
  88. end
  89. % Now compute the regression matrix
  90.  
  91. nmax=n;
  92. M=floor(maxsize/n);
  93. R=zeros(n);F=zeros(n,1);
  94. fb=((appr=='fb') | (appr=='FB'));
  95. if fb,RB=zeros(n);FBc=zeros(n,1);yb=conj(y(Ncap:-1:1));end
  96. for k=nmax:M:Ncap-1
  97.     jj=(k+1:min(Ncap,k+M));
  98.     phi=zeros(length(jj),n);if fb,phib=zeros(length(jj),n);end
  99.     for k=1:n,phi(:,k)=-y(jj-k);end
  100.     if fb,for k=1:n,phib(:,k)=-yb(jj-k);end,end
  101.     R=R+phi'*phi;F=F+phi'*y(jj);
  102.     if fb,RB=RB+phib'*phib;FBc=FBc+phib'*yb(jj);end
  103. end
  104. P=inv(R);
  105. if appr~='bu' & appr~='gl'
  106. if ~fb
  107.     if Ncap>M , t=R\F;else t=phi\y(jj);end
  108. end
  109. if fb
  110.     if Ncap>M/2 , t=(R+RB)\(F+FBc);else t=[phi;phib]\[y(jj);yb(jj)];end
  111. end
  112. th(3,1:n)=t.';
  113. e=pe(y(nstart:nend),th);lam=e'*e/(nend-nstart+1-n);
  114. th(1,1)=lam; th(2,1)=lam*(1+n/Ncap0)/(1-n/Ncap0);
  115. end
  116. if pt~=0,P=lam*P;th(4:n+3,1:n)=P;end
  117.