home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / ETFE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.6 KB  |  73 lines
  1. function g = etfe(z,M,N,T)
  2. %ETFE    Computes the Empirical Transfer Function Estimate and Periodogram
  3. %
  4. %    G = etfe(Z)   or   G = etfe(Z,M)
  5. %
  6. %    Z contains the input-output data [y u] or a time series y. Only
  7. %    single (or no) input systems can be handled. If an input is present
  8. %    G is returned as the ETFE (the ratio of the output Fourier transform
  9. %    to the input Fourier transform) for the data. For a time series G
  10. %    is returned as the periodogram (the normed absolute square of the
  11. %    Fourier transform) of the data. G is returned in the standard frequency
  12. %    function format. See HELP FREQFUNC.
  13. %
  14. %    With M specified,a smoothing operation is performed on the raw spectral
  15. %    estimates. The effect of M is then analogous to the same argument in
  16. %    SPA.
  17. %
  18. %    The transfer function is estimated at 128 equally spaced frequencies
  19. %    between 0 (excluded) and pi. This number can be changed to N (a power
  20. %    of two) and the sampling interval can be changed from 1 (default) to T
  21. %    by    G = etfe(Z,M,N,T).
  22. %    The transfer function can be plotted by BODEPLOT. bodeplot(etfe(z))
  23. %    is a possible construction. See also SPA.
  24.  
  25. %    L. Ljung 7-7-87
  26. %    Revised 3-8-89, 4-20-91
  27. %    Copyright (c) 1987-90 The MathWorks, Inc.
  28. %    All Rights Reserved
  29.  
  30. %    Some initial checks on the arguments
  31.  
  32. if nargin<4,T=1;end, if T<0,T=1;end,if isempty(T),T=1;end
  33. if nargin<3,N=128;end, if N<0, N=128;end,if isempty(N),N=128;end
  34. if nargin<2,M=[];end
  35. [Ncap,nz]=size(z);
  36. if nz>2, error('This routine works only for single-input systems!'),return,end
  37.  
  38. % ** Add zeros to data so that it can support N frequency points **
  39.  
  40. if Ncap<2*N,z=[z;zeros(2*N-Ncap,nz)];end 
  41.  
  42. % ** Compute the Fourier transforms by FFT **
  43. nfft = 2^nextpow2(max(Ncap,2*N)); % Corr 891007
  44. if nz==1,Y=abs(fft(z(:,1),nfft)).^2;else Y=fft(z(:,1),nfft);end,l=length(Y);
  45. if M<0,M=l;end,if isempty(M),M=l;end
  46. M1=fix(l/M);sc=l/(2*N);
  47. if fix(sc)~=sc, error('N must be a power of two!'),return,end
  48. if M1>1,
  49. ha = .54 - .46*cos(2*pi*(0:M1)'/M1);
  50. ha=ha/(norm(ha)^2);Y=[Y(l-M1+2:l);Y];
  51.     if nz==1,Y=filter(ha,1,Y);end
  52. end
  53. Yd=Y(M1+fix(M1/2)+sc:sc:M1+fix(M1/2)+l/2);
  54. if nz==1,
  55.     g(1,1:2)=[100 0];
  56.     g(2:N+1,1)=(1:N)'*pi/T/N;%%%%%%Correction made here
  57.     g(2:N+1,2)=T*Yd/Ncap;%Corr 910420
  58.     return
  59. end
  60. if M1==1,clear Y,end
  61. U=fft(z(:,2),nfft);
  62. if M1>1,
  63.     U=[U(l-M1+2:l);U];
  64.     Y=filter(ha,1,Y.*conj(U));
  65.     U=filter(ha,1,abs(U).^2);
  66.     Yd=Y(M1+fix(M1/2)+sc:sc:M1+fix(M1/2)+l/2);
  67. end
  68. Ud=U(M1+fix(M1/2)+sc:sc:l/2+M1+fix(M1/2));clear U
  69. g(1,1:3)=[101 1 21];
  70. g(2:N+1,1)=(1:N)'*pi/N/T; %%%%%Correction made here
  71. g(2:N+1,2)=abs(Yd./Ud);
  72. g(2:N+1,3)=-180*phase((Yd./Ud)')'/pi;
  73.