function [GH,PV]=spa(z,M,W,maxsize,T)
%SPA	Performs spectral analysis.
%
%	G = spa(Z)  or	[G,NSP] = spa(Z)
%
%	Z : The output-input data Z=[y u], with y and u being column vectors.
%	For multi-variable systems Z=[y1 y2 .. u1 u2 ..]. For time-series Z=y.
%	If an input is present G is returned as the transfer function estimate,
%	and NSP (if specified) is the estimated noise spectrum. 
%	Estimated standard deviations are also supplied. The functions are of
%	the standard frequency function format (see HELP FREQFUNC)
%	If Z is a time series G is returned as the estimated spectrum, along
%	with its estimated standard deviation.	A Hamming window of lag size 
%	min(length(Z)/10,30) is used, which can be changed to M by
%	G = spa(Z,M)  or [G,NSP] = spa(Z,M)
%       For multi-output systems M is entered as a row vector of the same 
%	length as the number of outputs. For default windows use M=[-1 -1 ..].
%	The functions are computed at 128 equally spaced frequency-values
%	between 0 (excluded) and pi. The functions can be computed at arbitrary
%	frequencies w (a row vector, generated e.g. by LOGSPACE) by
%	G = SPA(Z,M,w)  or  [G,NSP] = SPA(Z,M,w)
%	With G = spa(Z,M,w,maxsize,T) also the memory trade-off and the
%	sampling interval can be changed. See HELP AUXVAR and SETT.

%	L. Ljung 10-1-86, 29-3-92 (mv-case)
%	Copyright (c) 1986-92 by the MathWorks, Inc.
%	All Rights Reserved.

[Ncap,nz]=size(z);  i=sqrt(-1); dw=0;
if nz>Ncap, error(' The data should be entered in column vectors'),return,end
% *** Set up default values ***
maxsdef = idmsize(Ncap); %Mod 89-10-31
if nargin<5, T=[]; end
if isempty(T), T=1; end,if T<0, T=1;end
if nargin<2, M=[]; end %Mod 89-10-31
Mdef=min(30,floor(Ncap/10));
if isempty(M), M=Mdef;end       %Mod 89-10-31
if any(M-M(1)~=0), 
disp('The same window size has to be applied to all outputs')
disp('The first entry will by used by spa')
end
if any(M<0), M=Mdef*ones(length(M));end       %Mod 89-10-31
if nargin<4, maxsize=[]; end
if isempty(maxsize), maxsize=maxsdef;end
if maxsize<0, maxsize=maxsdef;end
if nargin<3, W=[];end
if isempty(W), W=[1:128]*pi/128/T; dw=1;end
if any(W<0), W=[1:128]*pi/128/T; dw=1;end
ny=length(M);nu=nz-ny;M=M(1);
if nu<0,error('You have specified too many outputs in M!'),end
if dw==1 & M>128 & ny==1 & nu<2,
   disp('Suggestion: For large M, use etfe(z,M) instead!')
end %mod 91-11-11
[nrw,rcw]=size(W);if nrw>rcw,W=W.';end 
n=length(W);

% *** Compute the covariance functions ***

R=covf(z,M,maxsize);

% *** Form the Hamming lag window ***

pd=pi/(M-1);
window=[1 0.5*ones(1,M-1)+0.5*cos(pd*[1:M-1])];
if dw~=1 | nz>2, window(1)=window(1)/2;end

if nz<3
% *** For a SISO system with equally spaced frequency points
%     we compute the spectra using FFT for highest speed ***

noll=zeros(1,2*(n-M)+1);
order=2:n+1;
nfft = 2.^nextpow2(2*n);
if nz==2
    if dw==1
       rtem=R(4,:).*window;
       FIU=real(fft([rtem noll rtem(M:-1:2)],nfft));FIU=FIU(order);
       rtem=R(3,:).*window;
       FIYU=fft([R(2,1:M).*window noll rtem(M:-1:2)],nfft);
       FIYU=FIYU(order);
    end

 % *** For arbitrary frequency points we use POLYVAL for the
%      computation of spectra ****

    if dw~=1
        rtem=R(4,:).*window;
        FIU=2*real(polyval(rtem(M:-1:1),exp(-i*W*T)));
        rtem=R(2,:).*window; rtemp=R(3,:).*window;
        FIYU=polyval(rtem(M:-1:1),exp(-i*W*T)) + polyval(rtemp(M:-1:1),exp(i*W*T));
    end

% *** Now compute the transfer function estimate ***

       G=FIYU./FIU;
       GH(1,1:3)=[101,1,21];
       GH(2:n+1,1)=W'; GH(2:n+1,2)=abs(G'); 
       GH(2:n+1,3)=phase(G)'*180/pi;
end

% *** Compute the noise spectrum (= the output spectrum for
%     a time series) ***


rtem=R(1,:).*window;
   if dw==1
FIY=real(fft([rtem noll rtem(M:-1:2)],nfft)); FIY=FIY(order);
   end
   if dw~=1
            FIY=2*real(polyval(rtem(M:-1:1),exp(-i*W*T)));
   end
if nz==2, FFV=(FIYU.*conj(FIYU))./FIU;,else FFV=zeros(1,n);end
PV(1,1:3)=[100 0 50];
PV(2:n+1,1)=W'; PV(2:n+1,2)=T*abs(FIY-FFV)';%mod 9009 
PV(2:n+1,3)=sqrt(0.75*M/Ncap*2)*PV(2:n+1,2);
if nz==2,GH(:,4)=GH(:,3);GH(1,1:5)=[101,1,51,21,71];
	GH(2:n+1,3)=sqrt(0.75*M/Ncap/2)*sqrt((PV(2:n+1,2)/T)./FIU');
	 GH(2:n+1,5)=(180/pi)*GH(2:n+1,3)./GH(2:n+1,2);
end
if nz==1 , GH=PV;end

return
end  % end case nz<3
%
% 
% *** Now we have to deal with the multi-variable case ***
%
if ny==1, 
disp(['Your data is interpreted as single output with ',int2str(nu) ,' inputs'])
disp('If you have multi-output data, let the argument M have as many elements as the number of outputs!')
end
FI=zeros(nz,n);
for k=1:nz^2
    rtem=R(k,:).*window;
    FI(k,:)=polyval(rtem(M:-1:1),exp(-i*W*T));
end
outputs=[1:ny];inputs=[ny+1:nz];
for k=1:n
    FIZ=zeros(nz); FIZ(:)=FI(:,k); FIZ=FIZ+FIZ';
    if nu>0, GCk=(FIZ(outputs,inputs)/FIZ(inputs,inputs))';
       PHIUINV=inv(FIZ(inputs,inputs));
       for kd=1:nu,diagCOV(k,kd)=PHIUINV(kd,kd);end
       PVCk=(FIZ(outputs,outputs) -FIZ(outputs,inputs)*GCk)';
       GC(k,:)=GCk(:)';PVC(k,:)=diag(PVCk)';
    else PVC(k,:)=diag(FIZ)';
    end
end

PA=abs(PVC);
if nu>0
GH=zeros(n+1,5*nu*ny);
for ky=1:ny
for k=1:5:nu*5
    ku=(k+4)/5;kky=(ky-1)*nu;kkz=(ky-1)*nu*5;
    GH(1,kkz+[k:k+4])=(ky-1)*1000+[100+ku,ku,50+ku,20+ku,70+ku];
    GH(2:n+1,kkz+k)=W';
    GH(2:n+1,kkz+k+1)=abs(GC(:,ku+kky));
    GH(2:n+1,kkz+k+2)=sqrt(0.75*M/Ncap/2*(PA(:,ky).*diagCOV(:,(k+4)/5)));
    GH(2:n+1,kkz+k+3)=phase(GC(:,ku+kky)')'*180/pi;
    GH(2:n+1,kkz+k+4)=(180/pi)*GH(2:n+1,kkz+k+2)./GH(2:n+1,kkz+k+1);
end
end    
  if nargout==1, return,end
end
for ky=1:ny
PV(1,(ky-1)*3+[1:3])=1000*(ky-1)+[100,0,50];
PV(2:n+1,1+(ky-1)*3)=W';
PV(2:n+1,2+(ky-1)*3)=T*PA(:,ky);%mod 9009
PV(2:n+1,3+(ky-1)*3)=sqrt(1.5*M/Ncap)*PA(:,ky)*T;
end
if nu==0, GH=PV;end

