home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / TRF.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.4 KB  |  96 lines
  1. function [G,PV]=trf(th,nnu,w,nny)
  2. %TRF    Computes a model's frequency function.
  3. %
  4. %    G = trf(TH)  or  [G,NSP] = trf(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.
  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 = trf(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 RAD/S, (a row vector,
  20. %    generated e.g. by LOGSPACE) by G = trf(TH,ku,w). The transfer function
  21. %    can be plotted by BODEPLOT. bodeplot(trf(TH)) 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=trf(TH,ku,w,ky); (Default is all outputs). See also TH2FF.
  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,wdef=pi*[1:128]/128/T;
  37. else wdef=logspace(log10(pi/abs(T)/100),log10(10*pi/abs(T)),128);
  38. end
  39. if isempty(w),w=wdef;end
  40. if length(w)==1, if w<0, w=wdef;end,end
  41. [nrw,ncw]=size(w);if nrw>ncw, w=w.';end
  42. if T>0 & any(w>pi/T),disp('WARNING: Frequencies in w exceed the Nyquist frequency!'),end
  43.  
  44. if isempty(nnu), nnu=1:th(1,3);end
  45. if length(nnu)==1,if nnu==0,nnu=[];end,if nnu<0,nnu=1:th(1,3);end,end
  46. if isthss(th), if nargout==1, eval('G=trfss(th,nnu,nny,w);')
  47.          else eval('[G,PV]=trfss(th,nnu,nny,w);'),end,return,end
  48. % *** Compute the model polynomials and form the basic
  49. %      matrix of frequencies ***
  50. nu=th(1,3);
  51. [a,b,c,d,f]=th2poly(th);
  52. na=th(1,4);nb=th(1,5:4+nu);nc=th(1,5+nu);nd=th(1,6+nu);
  53. 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);
  54. nm=max([length(a)+length(f)-1 length(b) length(c) length(d)+length(a)-1]);
  55.    i=sqrt(-1);
  56. if T>0,OM=exp(-i*[0:nm-1]'*w*T);end
  57. if T<=0,
  58. OM=ones(1,length(w));
  59. for kom=1:nm-1
  60. OM=[OM;(i*w).^kom];
  61. end
  62. [nrth,ncth]=size(th);
  63. if nrth>3+nnd,delays=th(nrth,1:length(nk));else delays=zeros(1,length(nk));end
  64. end
  65.  
  66. % *** Compute the transfer function(s) GC=B/AF ***
  67.  
  68. sc=1;kf=0;ks=0;nrc=length(w)+1;
  69. for k=nnu
  70.     G(1,sc:sc+2)=[100+k,k,k+20];    
  71.     G(2:nrc,sc)=w';
  72.      gn=conv(a,f(k,:));
  73. if T>0,indb=1:length(b(k,:));indg=1:length(gn);
  74. else indb=length(b(k,:)):-1:1; indg=length(gn):-1:1;end
  75.      GC=(b(k,:)*OM(indb,:))./(gn*OM(indg,:));
  76.     G(2:nrc,sc+1)=abs(GC');
  77.     G(2:nrc,sc+2)=phase(GC)'*180/pi;
  78.     if T<0, G(2:nrc,sc+2)=G(2:nrc,sc+2)-w'*delays(k)*180/pi;end
  79. sc=sc+3;
  80. end
  81. if nargout==1 & nu>0, return, end
  82.  
  83. % *** Compute the transfer function H=C/AD *** 
  84.  
  85. hn=conv(a,d);
  86. if T>0,indc=1:length(c);indh=1:length(hn);
  87. else indc=length(c):-1:1; indh=length(hn):-1:1;end
  88. H=(c*OM(indc,:))./(hn*OM(indh,:));
  89. PV(1,1:3)=[100,0,50];
  90. PV(2:nrc,1)=w';
  91. PV(2:nrc,2)=abs(T)*(abs(H).^2)'*th(1,1);
  92.  
  93. if isempty(nnu) G=PV;end
  94.  
  95.  
  96.