home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / SPA.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  5.7 KB  |  173 lines
  1. function [GH,PV]=spa(z,M,W,maxsize,T)
  2. %SPA    Performs spectral analysis.
  3. %
  4. %    G = spa(Z)  or    [G,NSP] = spa(Z)
  5. %
  6. %    Z : The output-input data Z=[y u], with y and u being column vectors.
  7. %    For multi-variable systems Z=[y1 y2 .. u1 u2 ..]. For time-series Z=y.
  8. %    If an input is present G is returned as the transfer function estimate,
  9. %    and NSP (if specified) is the estimated noise spectrum. 
  10. %    Estimated standard deviations are also supplied. The functions are of
  11. %    the standard frequency function format (see HELP FREQFUNC)
  12. %    If Z is a time series G is returned as the estimated spectrum, along
  13. %    with its estimated standard deviation.    A Hamming window of lag size 
  14. %    min(length(Z)/10,30) is used, which can be changed to M by
  15. %    G = spa(Z,M)  or [G,NSP] = spa(Z,M)
  16. %       For multi-output systems M is entered as a row vector of the same 
  17. %    length as the number of outputs. For default windows use M=[-1 -1 ..].
  18. %    The functions are computed at 128 equally spaced frequency-values
  19. %    between 0 (excluded) and pi. The functions can be computed at arbitrary
  20. %    frequencies w (a row vector, generated e.g. by LOGSPACE) by
  21. %    G = SPA(Z,M,w)  or  [G,NSP] = SPA(Z,M,w)
  22. %    With G = spa(Z,M,w,maxsize,T) also the memory trade-off and the
  23. %    sampling interval can be changed. See HELP AUXVAR and SETT.
  24.  
  25. %    L. Ljung 10-1-86, 29-3-92 (mv-case)
  26. %    Copyright (c) 1986-92 by the MathWorks, Inc.
  27. %    All Rights Reserved.
  28.  
  29. [Ncap,nz]=size(z);  i=sqrt(-1); dw=0;
  30. if nz>Ncap, error(' The data should be entered in column vectors'),return,end
  31. % *** Set up default values ***
  32. maxsdef = idmsize(Ncap); %Mod 89-10-31
  33. if nargin<5, T=[]; end
  34. if isempty(T), T=1; end,if T<0, T=1;end
  35. if nargin<2, M=[]; end %Mod 89-10-31
  36. Mdef=min(30,floor(Ncap/10));
  37. if isempty(M), M=Mdef;end       %Mod 89-10-31
  38. if any(M-M(1)~=0), 
  39. disp('The same window size has to be applied to all outputs')
  40. disp('The first entry will by used by spa')
  41. end
  42. if any(M<0), M=Mdef*ones(length(M));end       %Mod 89-10-31
  43. if nargin<4, maxsize=[]; end
  44. if isempty(maxsize), maxsize=maxsdef;end
  45. if maxsize<0, maxsize=maxsdef;end
  46. if nargin<3, W=[];end
  47. if isempty(W), W=[1:128]*pi/128/T; dw=1;end
  48. if any(W<0), W=[1:128]*pi/128/T; dw=1;end
  49. ny=length(M);nu=nz-ny;M=M(1);
  50. if nu<0,error('You have specified too many outputs in M!'),end
  51. if dw==1 & M>128 & ny==1 & nu<2,
  52.    disp('Suggestion: For large M, use etfe(z,M) instead!')
  53. end %mod 91-11-11
  54. [nrw,rcw]=size(W);if nrw>rcw,W=W.';end 
  55. n=length(W);
  56.  
  57. % *** Compute the covariance functions ***
  58.  
  59. R=covf(z,M,maxsize);
  60.  
  61. % *** Form the Hamming lag window ***
  62.  
  63. pd=pi/(M-1);
  64. window=[1 0.5*ones(1,M-1)+0.5*cos(pd*[1:M-1])];
  65. if dw~=1 | nz>2, window(1)=window(1)/2;end
  66.  
  67. if nz<3
  68. % *** For a SISO system with equally spaced frequency points
  69. %     we compute the spectra using FFT for highest speed ***
  70.  
  71. noll=zeros(1,2*(n-M)+1);
  72. order=2:n+1;
  73. nfft = 2.^nextpow2(2*n);
  74. if nz==2
  75.     if dw==1
  76.        rtem=R(4,:).*window;
  77.        FIU=real(fft([rtem noll rtem(M:-1:2)],nfft));FIU=FIU(order);
  78.        rtem=R(3,:).*window;
  79.        FIYU=fft([R(2,1:M).*window noll rtem(M:-1:2)],nfft);
  80.        FIYU=FIYU(order);
  81.     end
  82.  
  83.  % *** For arbitrary frequency points we use POLYVAL for the
  84. %      computation of spectra ****
  85.  
  86.     if dw~=1
  87.         rtem=R(4,:).*window;
  88.         FIU=2*real(polyval(rtem(M:-1:1),exp(-i*W*T)));
  89.         rtem=R(2,:).*window; rtemp=R(3,:).*window;
  90.         FIYU=polyval(rtem(M:-1:1),exp(-i*W*T)) + polyval(rtemp(M:-1:1),exp(i*W*T));
  91.     end
  92.  
  93. % *** Now compute the transfer function estimate ***
  94.  
  95.        G=FIYU./FIU;
  96.        GH(1,1:3)=[101,1,21];
  97.        GH(2:n+1,1)=W'; GH(2:n+1,2)=abs(G'); 
  98.        GH(2:n+1,3)=phase(G)'*180/pi;
  99. end
  100.  
  101. % *** Compute the noise spectrum (= the output spectrum for
  102. %     a time series) ***
  103.  
  104.  
  105. rtem=R(1,:).*window;
  106.    if dw==1
  107. FIY=real(fft([rtem noll rtem(M:-1:2)],nfft)); FIY=FIY(order);
  108.    end
  109.    if dw~=1
  110.             FIY=2*real(polyval(rtem(M:-1:1),exp(-i*W*T)));
  111.    end
  112. if nz==2, FFV=(FIYU.*conj(FIYU))./FIU;,else FFV=zeros(1,n);end
  113. PV(1,1:3)=[100 0 50];
  114. PV(2:n+1,1)=W'; PV(2:n+1,2)=T*abs(FIY-FFV)';%mod 9009 
  115. PV(2:n+1,3)=sqrt(0.75*M/Ncap*2)*PV(2:n+1,2);
  116. if nz==2,GH(:,4)=GH(:,3);GH(1,1:5)=[101,1,51,21,71];
  117.     GH(2:n+1,3)=sqrt(0.75*M/Ncap/2)*sqrt((PV(2:n+1,2)/T)./FIU');
  118.      GH(2:n+1,5)=(180/pi)*GH(2:n+1,3)./GH(2:n+1,2);
  119. end
  120. if nz==1 , GH=PV;end
  121.  
  122. return
  123. end  % end case nz<3
  124. %
  125. % 
  126. % *** Now we have to deal with the multi-variable case ***
  127. %
  128. if ny==1, 
  129. disp(['Your data is interpreted as single output with ',int2str(nu) ,' inputs'])
  130. disp('If you have multi-output data, let the argument M have as many elements as the number of outputs!')
  131. end
  132. FI=zeros(nz,n);
  133. for k=1:nz^2
  134.     rtem=R(k,:).*window;
  135.     FI(k,:)=polyval(rtem(M:-1:1),exp(-i*W*T));
  136. end
  137. outputs=[1:ny];inputs=[ny+1:nz];
  138. for k=1:n
  139.     FIZ=zeros(nz); FIZ(:)=FI(:,k); FIZ=FIZ+FIZ';
  140.     if nu>0, GCk=(FIZ(outputs,inputs)/FIZ(inputs,inputs))';
  141.        PHIUINV=inv(FIZ(inputs,inputs));
  142.        for kd=1:nu,diagCOV(k,kd)=PHIUINV(kd,kd);end
  143.        PVCk=(FIZ(outputs,outputs) -FIZ(outputs,inputs)*GCk)';
  144.        GC(k,:)=GCk(:)';PVC(k,:)=diag(PVCk)';
  145.     else PVC(k,:)=diag(FIZ)';
  146.     end
  147. end
  148.  
  149. PA=abs(PVC);
  150. if nu>0
  151. GH=zeros(n+1,5*nu*ny);
  152. for ky=1:ny
  153. for k=1:5:nu*5
  154.     ku=(k+4)/5;kky=(ky-1)*nu;kkz=(ky-1)*nu*5;
  155.     GH(1,kkz+[k:k+4])=(ky-1)*1000+[100+ku,ku,50+ku,20+ku,70+ku];
  156.     GH(2:n+1,kkz+k)=W';
  157.     GH(2:n+1,kkz+k+1)=abs(GC(:,ku+kky));
  158.     GH(2:n+1,kkz+k+2)=sqrt(0.75*M/Ncap/2*(PA(:,ky).*diagCOV(:,(k+4)/5)));
  159.     GH(2:n+1,kkz+k+3)=phase(GC(:,ku+kky)')'*180/pi;
  160.     GH(2:n+1,kkz+k+4)=(180/pi)*GH(2:n+1,kkz+k+2)./GH(2:n+1,kkz+k+1);
  161. end
  162. end    
  163.   if nargout==1, return,end
  164. end
  165. for ky=1:ny
  166. PV(1,(ky-1)*3+[1:3])=1000*(ky-1)+[100,0,50];
  167. PV(2:n+1,1+(ky-1)*3)=W';
  168. PV(2:n+1,2+(ky-1)*3)=T*PA(:,ky);%mod 9009
  169. PV(2:n+1,3+(ky-1)*3)=sqrt(1.5*M/Ncap)*PA(:,ky)*T;
  170. end
  171. if nu==0, GH=PV;end
  172.  
  173.