home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / COVF.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.6 KB  |  54 lines
  1. function R=covf(z,M,maxsize)
  2. %COVF  Computes the covariance function estimate for a data matrix.
  3. %
  4. %    R = covf(Z,M)
  5. %
  6. %    Z : An   N x nz data matrix, typically Z=[y u]
  7. %
  8. %    M: The maximum delay - 1, for which the covariance function is estimated
  9. %    R: The covariance function of Z, returned so that the entry
  10. %    R((i+(j-1)*nz,k+1) is  the estimate of E Zi(t) * Zj(t+k)
  11. %    The size of R is thus nz^2 x M
  12. %
  13. %    For nz<3, an FFT algorithm is used, memory size permitting.
  14. %    For nz>2, straightforward summation is used (in COVF2)
  15. %
  16. %    The memory trade-off is affected by
  17. %    R = COVF(Z,M,maxsize)
  18. %    See HELP AUXVAR  for how to use this option.
  19.  
  20. %    L. Ljung 10-1-86
  21. %    Copyright (c) 1986-90 by the MathWorks, Inc.
  22. %    All Rights Reserved.
  23. [Ncap,nz]=size(z);
  24. maxsdef=idmsize(Ncap);
  25. if nargin<3,maxsize=maxsdef;end
  26. if isempty(maxsize),maxsize=maxsdef;end
  27. if nz>Ncap, error('The data should be arranged in columns'),return,end
  28.  
  29. if nz>2  |  Ncap>maxsize/8, R=covf2(z,M); return,end
  30.  
  31. nfft = 2.^ceil(log(2*Ncap)/log(2));
  32. Ycap=fft(([z(:,1)' zeros(1,Ncap)]),nfft);
  33. if nz==2, Ucap=fft(([z(:,2)'  zeros(1,Ncap)]),nfft);
  34.            YUcap=Ycap.*conj(Ucap);
  35.            UAcap=abs(Ucap).^2;
  36.            clear Ucap
  37. end
  38. YAcap=abs(Ycap).^2;
  39. clear Ycap
  40. RYcap=fft(YAcap,nfft);
  41. n=length(RYcap);
  42. R=real(RYcap(1:M))/Ncap/n;
  43. clear RYcap
  44. if nz==1,return,end
  45.          RUcap=fft(UAcap,nfft);
  46.          ru=real(RUcap(1:M))/Ncap/n;
  47.          clear RUcap
  48.          RYUcap=fft(YUcap,nfft);
  49.          ryu=real(RYUcap(1:M))/Ncap/n;
  50.          ruy=real(RYUcap(n:-1:n-M+2))/Ncap/n;
  51.          clear RYUcap
  52. R=[R;[ryu(1) ruy];ryu;ru];
  53.  
  54.