home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / TH2FF.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.3 KB  |  125 lines
  1. function [G,PV]=th2ff(th,nnu,w,nny)
  2. %TH2FF    Computes a model's frequency function,along with its standard deviation
  3. %
  4. %    G = th2ff(TH)  or  [G,NSP] = th2ff(TH)
  5. %
  6. %    TH: A matrix defining a model, as described in HELP THETA
  7. %    G is returned as the transfer function estimate, and NSP (if specified)
  8. %    as the noise spectrum, corresponding to the model TH. These matrices
  9. %    contain also estimated standard deviations, calculated from the
  10. %    covariance matrix in TH, and are  of the standard frequency function
  11. %    format (see HELP FREQFUNC).If TH describes a time series, G is returned
  12. %    as its spectrum. Both discrete and continuous time models are handled.
  13. %    
  14. %    If the model TH has several inputs, G will be returned as the transfer
  15. %    functions of selected inputs # j1 j2 .. jk by 
  16. %    G = th2ff(TH,[j1 j2 ... jk])  [default is all inputs]. The functions 
  17. %    are computed at 128 equally spaced frequency-values between 0(excluded)
  18. %    and pi/T, where T is the sampling interval specified by TH. The func-
  19. %    tions can be computed at arbitrary frequencies w (a row vector, gene-
  20. %    rated e.g. by LOGSPACE) by G = th2ff(TH,ku,w). The transfer function
  21. %    can be plotted by BODEPLOT. bodeplot(th2ff(TH),sd) is a possible 
  22. %    construction. If the model TH has several outputs, G will be returned
  23. %    as the frequency function at selected outputs ky (a row vector) by
  24. %    G=th2ff(TH,ku,w,ky); (Default is all outputs). See also TRF.
  25.  
  26. %    L. Ljung 7-7-87,1-25-92
  27. %    Copyright (c) 1987-92 by the MathWorks, Inc.
  28. %    All Rights Reserved.
  29.  
  30. T=gett(th);
  31.  
  32. % *** Set up default values ***
  33. if nargin<4, nny=[];end
  34. if nargin<3, w=[];end
  35. if nargin<2 , nnu=[];end
  36. if T>0
  37.     wdef=pi*[1:128]/128/T;
  38. else    wdef=logspace(log10(pi/abs(T)/100),log10(10*pi/abs(T)),128);
  39. end
  40.  
  41. if isempty(w),w=wdef;end
  42. if length(w)==1, if w<0, w=wdef;end,end
  43. [nrw,ncw]=size(w);if nrw>ncw, w=w.';end
  44. if T>0 & any(w>pi/T),disp('WARNING: Frequencies in w exceed the Nyquist frequency!'),end
  45. if isempty(nnu),nnu=1:th(1,3);end
  46. if isthss(th)
  47.    if nargout==1 
  48.       eval('G=trfsssd(th,nnu,nny,w);')
  49.    else eval('[G,PV]=trfsssd(th,nnu,nny,w);')
  50.    end
  51.    return
  52. end
  53. if length(nnu)==1, if nnu==0,nnu=[];end,if nnu<0, nnu=1:th(1,3);end,end
  54.  
  55. % *** Compute the model polynomials and form the basic
  56. %      matrix of frequencies ***
  57. nu=th(1,3);
  58.  
  59. [a,b,c,d,f]=th2poly(th);[par,P]=th2par(th);if isempty(P),isP=0;else isP=1;end
  60. na=th(1,4);nb=th(1,5:4+nu);nc=th(1,5+nu);nd=th(1,6+nu);
  61. nf=th(1,7+nu:6+2*nu);nk=th(1,7+2*nu:6+3*nu);nnd=na+sum(nb)+nc+nd+sum(nf);
  62. nm=max([length(a)+length(f)-1 length(b) length(c) length(d)+length(a)-1]);
  63. i=sqrt(-1);
  64. if T>0,OM=exp(-i*[0:nm-1]'*w*T);end
  65. if T<=0,
  66.    OM=ones(1,length(w));
  67.    for kom=1:nm-1
  68.        OM=[OM;(i*w).^kom];
  69.    end
  70.    [nrth,ncth]=size(th);
  71.    if nrth>3+nnd
  72.        delays=th(nrth,1:length(nk));else delays=zeros(1,length(nk));
  73.    end
  74. end
  75.  
  76. % *** Compute the transfer function(s) GC=B/AF ***
  77.  
  78. sc=1;kf=0;ks=0;nrc=length(w)+1;
  79. for k=nnu
  80.     G(1,sc:sc+4)=[100+k,k,k+50,k+20,k+70];    
  81.     G(2:nrc,sc)=w';
  82.     gn=conv(a,f(k,:));
  83.     if T>0
  84.         indb=1:length(b(k,:));indg=1:length(gn);
  85.     else indb=length(b(k,:)):-1:1; indg=length(gn):-1:1;
  86.     end
  87.     GC=(b(k,:)*OM(indb,:))./(gn*OM(indg,:));
  88.     G(2:nrc,sc+1)=abs(GC');
  89.     G(2:nrc,sc+3)=phase(GC)'*180/pi;
  90.     if T<0, G(2:nrc,sc+3)=G(2:nrc,sc+3)-w'*delays(k)*180/pi;end
  91.     ll=[1:na,na+ks+1:na+ks+nb(k),na+ks+nb(k)+nc+nd+kf+1:na+ks+nb(k)+nc+nd+kf+nf(k)];
  92.     ks=ks+nb(k);kf=kf+nf(k);
  93.     if isP,
  94.         P=th(3+ll,ll);
  95.         [dga,dgp]=ffsdcal(a,b(k,1:nb(k)+nk(k)),f(k,1:nf(k)+1),na,nb(k),nf(k),nk(k),GC,OM,P,T); 
  96.         G(2:nrc,[sc+2 sc+4])=[dga,dgp]; 
  97.     end
  98.     sc=sc+5;
  99. end
  100. if nargout==1 & nu>0, return, end
  101.  
  102. % *** Compute the transfer function H=C/AD *** 
  103.  
  104. hn=conv(a,d);
  105. if T>0,indc=1:length(c);indh=1:length(hn);
  106. else   indc=length(c):-1:1; indh=length(hn):-1:1;
  107. end
  108. H=(c*OM(indc,:))./(hn*OM(indh,:));
  109.  
  110. PV(1,1:3)=[100,0,50];
  111. PV(2:nrc,1)=w';
  112. PV(2:nrc,2)=(abs(H).^2)'*th(1,1)*abs(T);
  113. %
  114. % *** Compute the standard deviation of the spectrum ***
  115. if isP,Ncap=getncap(th);
  116.    if ~isempty(Ncap)
  117.    ll=[1:na,na+sum(nb)+1:na+sum(nb)+nc+nd];
  118.    P=th(3+ll,ll);
  119.    dga=ffsdcal(a,c,d,na,nc,nd,0,H,OM,P,T);
  120.    PV(2:nrc,3)=sqrt(2/Ncap)*PV(2:nrc,2)+2*abs(T)*th(1,1)*(dga.*abs(H)');
  121.    end
  122. end
  123. if isempty(nnu), G=PV;end
  124.  
  125.