Documentation of trfsssd

Global Index (short | long) | Local contents | Local Index (short | long)

Function Synopsis

[G,PV] = trfsssd(th,nnu,nny,w)

Help text

TRFSSSD        Auxiliary routine to TH2FF

   [G,PHI]=trfsssd(th)

Cross-Reference Information

Listing of function trfsssd

function [G,PV] = trfsssd(th,nnu,nny,w)

%   L. Ljung 10-2-90,11-11-93
%   Copyright (c) 1986-98 by The MathWorks, Inc.
%   $Revision: 2.3 $  $Date: 1997/12/02 03:42:18 $

[a,b,c,d,k]=eta2ss(th); [nx,nu]=size(b);[ny,nx]=size(c);
if any(nnu>nu),error('There are not that many inputs in the model!'),end
if any(nny>ny),error('There are not that many outputs in the model!'),end
T=gett(th);Ncap=getncap(th);if isempty(Ncap),Ncap=inf;end
if isempty(nnu),nnu=1:nu;end
if isempty(nny),nny=1:ny;end
if length(nnu)==1,if nnu<0,nnu=1:nu;end,end
if length(nny)==1,if nny<0,nny=1:ny;end,end
lny=length(nny);
[dum,dum,lambda]=th2par(th);
G=[];PV=[];
for ky=nny
    thbb=thss2th(th,ky,'noises');
    if nargout==1&nu>0
       gy=th2ff(thbb,nnu,w);
    else
       gy=th2ff(thbb,[nnu nu+[1:ny]],w);
       pp=zeros(length(w),1);psd=pp;
       for jj=nu+[1:ny]
         [w,amp,fas,sdamp]=getff(gy,jj,1);
         dum=zeros(length(w),1);dumsd=dum;
         for kk=nu+[1:ny]
           [w,ampkk,faskk,sdampkk]=getff(gy,kk,1);
           term=ampkk.*cos((faskk-fas)/180*pi)*lambda(jj-nu,kk-nu);
           dum=dum+term;
           if kk==jj,factor=4;else factor=2;end
           dumsd=dumsd+factor*term.^2;
         end
         pp=pp+dum.*amp;
         psd=psd+dumsd.*sdamp.^2;
       end
       gy=gy(:,1:5*length(nnu));
       pvy(1,:)=[100 0 50]+(ky-1)*1000;
       pvy(2:length(w)+1,:)=[w abs(T)*pp ...
                             sqrt(2/Ncap)*abs(T)*pp+abs(T)*sqrt(psd)];
       PV=[PV pvy];
    end
    if ~isempty(gy),    gy(1,:)=gy(1,:)+(ky-1)*1000;G=[G gy]; end
end

if length(nnu)>0
     for ku=nnu
        [g,p]=trfsaux(a,b,c(nny,:),d(nny,:),ku,w,T);
        scount=1;
        for ky=nny
          G(2:length(w)+1,find(G(1,:)==1000*(ky-1)+ku))=g(:,scount);
          G(2:length(w)+1,find(G(1,:)==1000*(ky-1)+20+ku))=p(:,scount);
          scount=scount+1;
        end
     end
end

if nargout>1 | nu==0
    lam=sqrtm(lambda);g=[];
    for ke=1:ny
       g=[g,trfsaux(a,k*lam,c,lam,ke,w,T).^2];
    end
    for ky=nny
       index=[ky:ny:ny*ny]; li=length(index);
       if li==1,col=g(:,index);else col=sum(g(:,index)')';end
       PV(2:length(w)+1,find(PV(1,:)==1000*(ky-1)))=...
          abs(T)*col;
    end
end
if nu==0 G=PV;end