home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / DTIMVEC.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.1 KB  |  65 lines
  1. function npts=dtimvec(a,b,c,x0,st)
  2. %DTIMVEC Returns a time vector for use with discrete time response
  3. %        functions.
  4. %
  5. %    NPTS=DTIMVEC(a,b,c,X0,ST,PRECISION) returns a suggestion  
  6. %    for the number of samples, NPTS. that should be used 
  7. %    with discrete state space systems (a,b,c) that have known 
  8. %    initial conditions X0. NPTS can be used,
  9. %    for instance, for step and impulse responses. 
  10. %    
  11. %    DTIMVEC attempts to produce a number of points that ensures the 
  12. %    output responses have decayed to approx. ST % of their initial 
  13. %    or peak values. 
  14.  
  15. %    A.C.W.Grace 9-21-89
  16. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  17.  
  18.      [r,n]=size(c);
  19. % Work out scale factors by looking at d.c. gain, 
  20. % eigenvectors and eigenvalues of the system.
  21.     [m,r]=eig(a);
  22.     r=diag(r);
  23. % Equate the  responses to the sum of a set of 
  24. % exponentials(r) multiplied by a vector of magnitude terms(vec).
  25.     [r1,n]=size(c);
  26.     if r1>1, c1=max(abs(c)); else c1=c; end
  27. % Cater for the case when m is singular
  28.     if rcond(m)<eps 
  29.         vec = (c1*m).'.*(1e4*ones(n,1));
  30.     else
  31.         vec=(c1*m).'.*(m\x0);
  32.     end
  33. % Cater for the case when poles and zeros cancel each other out:
  34.     vec=vec+1e-20*(vec==0);
  35. % d.c. gain (not including d matrix)
  36.     dcgain =c1*x0; % or dcgain=c\(eye(n)-a)*b;
  37.     ind=find(imag(r)>0);
  38. % If the d.c gain is small then base the scaling
  39. % factor on the estimated maximum peak in the response.
  40.     if abs(dcgain)<1e-8;
  41. % Estimation of the maximum peak in the response.
  42.         peak=2*imag(vec(ind)).*exp(-abs(0.5*real(r(ind)).*pi./imag(r(ind))));
  43.         pk=max([abs(dcgain);abs(peak);eps]);
  44.     else
  45.         pk=dcgain;
  46.     end
  47. % Work out the st% settling time for the responses.
  48. % (or st% peak value settling time for low d.c. gain systems).
  49.     lt=(st*pk./abs(vec));
  50.     n1=log(lt)./log(abs(r)+(abs(r)==1));
  51.     n=max(real(n1));
  52. % For unstable problems n will be negative
  53.     if n<=3, n=max(abs(n1))+2; end
  54.     n=min([n,1000]);
  55. % Round the maximum time to an appropriate value for plotting.
  56.     nn=chop(n,1,5);
  57.     if abs((n-nn)/n)>0.2, 
  58.         nn=chop(n,1,1);
  59.         if abs((n-nn)/n)>0.2
  60.             nn=chop(n,2,2);
  61.         end
  62.     end
  63.     npts=nn+1;
  64.  
  65.