Global Index (short | long) | Local contents | Local Index (short | long)
[GH,PV,ZSPE]=spa(z,M,W,maxsize,T,inhib)
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 also 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 also AUXVAR and SETT. A third output argument Z_SPE as in [G,PV,Z_SPE] = SPA(Z) gives directly the spectrum matrix of Z as follows: RESHAPE(Z_SPE(:,k),nz,nz)) = The spectrum S at frequency W(k), where S = sum_over_m [E(z(t+m)*z(t)') * exp(-i*W(k)*m*T)] z' is complex conjugate, transpose. nz is the number of signals in z. See also BODEPLOT, ETFE, FFPLOT and FREQFUNC.
| This function calls | This function is called by |
|---|---|
function [GH,PV,ZSPE]=spa(z,M,W,maxsize,T,inhib)
% L. Ljung 10-1-86, 29-3-92 (mv-case),11-11-94 (complex data case)
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 2.4 $ $Date: 1998/10/23 17:04:23 $
if nargin<1
disp('Usage: G = SPA(Z)')
disp(' [G,NSP,Z_SPECT] = SPA(Z,RES_PAR,FREQUENCIES,MAXSIZE,T)')
return
end
[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<6,inhib=0;end
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<0), M=Mdef*ones(1,length(M));end
M=M(:)';
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 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
W=W(:)';
ny=length(M);nu=nz-ny;M=M(1);
if any(any(imag(z)~=0))|nargout==3,spec_case=1;else spec_case=0;end
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
if M>=Ncap, M=Ncap-1;
disp(['WARNING: Resolution parameter M changed to ',int2str(M)])
end
[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 | spec_case, window(1)=window(1)/2;end
if nz<3&ny==1 & ~spec_case,
% *** 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 %if nz==2
% *** 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)=real(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)=real(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&~inhib&nz>2,
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];ZSPE=zeros(nz^2,n);
for k=1:n
FIZ=zeros(nz); FIZ(:)=FI(:,k); FIZ=FIZ+FIZ';ZSPE(:,k)=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)=real(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)=real(sqrt(1.5*M/Ncap)*PA(:,ky)*T);
end
if nu==0, GH=PV;end