home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / DTRSP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.7 KB  |  177 lines
  1. %function y = dtrsp(sys,input,T,tfinal,x0)
  2. %
  3. %  Time response of a discrete system:
  4. %    SYS:    packed form A,B,C,D system or constant.
  5. %    INPUT:    varying format input vector or constant.
  6. %    T:    sample time.
  7. %    TFINAL:    final time value
  8. %        (optional, default = input final time).
  9. %    X0:    initial state
  10. %        (optional, default = 0).
  11. %
  12. %   See also: TRSP, SAMHLD, TUSTIN, VINTERP, and VDCMATE.
  13.  
  14. function y = dtrsp(sys,input,T,tfinal,x0)
  15. if nargin == 0,
  16.     disp('usage: y = dtrsp(sys,input,T,tfinal,x0)')
  17.     return
  18.     end
  19.  
  20. if nargin < 3,
  21.     error('insufficient number of input arguments')
  22.     return
  23.     end
  24.  
  25. [type,ny,nu,nx] = minfo(sys);
  26. if type == 'syst',
  27.     if nargin < 5,
  28.         x0 = zeros(nx,1);
  29.     else
  30.         [m,n] = size(x0);
  31.         if m ~= nx | n ~= 1,
  32.             error('incorrect initial state vector dimension')
  33.             return
  34.                end
  35.         end
  36.     [a,b,c,d] = unpck(sys);
  37. elseif type ~= 'cons',
  38.     error('the system is not in packed or constant form')
  39.     return
  40.     end
  41.  
  42. %       check that the input vector is compatible.
  43.  
  44. [typeu,nru,ncu,nupts] = minfo(input);
  45. if typeu == 'vary' & ncu == 1,
  46.     [udat,uptr,ut] = vunpck(input);
  47. elseif typeu == 'cons' & ncu == 1,
  48.     nupts = 1;
  49.     udat = input;
  50.     ut = 0;
  51. else
  52.     error('invalid input time vector')
  53.     return
  54.     end
  55.  
  56. if nru ~= nu, 
  57.     error('incompatible input vector dimension')
  58.     return
  59.     end
  60.  
  61. %     an epsilon value for time interval calculation is created.
  62.  
  63. teps = T*1e-8;
  64.  
  65. %       The time vector is created.  The final time is derived from
  66. %       this in case the user did not specify an integer number of
  67. %       time steps.
  68.  
  69. if nargin < 4,
  70.     yt = [ut(1):T:max(ut)].';
  71. elseif tfinal < min(ut),
  72.     error('final time less than initial time')
  73.     return
  74. else
  75.     yt = [ut(1):T:tfinal].';
  76.     end
  77. tfinal = max(yt);
  78. nypts = length(yt);
  79.  
  80. %       Check if the final time is less than the maximum time in the input.
  81. %       If it is truncate the input.  Doing this now can save a
  82. %       lot of time if the interpolation is required.
  83.  
  84. if max(ut) - tfinal > teps,
  85.     nupts = max(find(ut <= tfinal+teps));
  86.     ut = ut(1:nupts);
  87.     udat = udat(1:uptr(nupts)+nru-1,:);
  88.     end
  89.  
  90. %        If the input and output time vectors do not
  91. %        match interpolate the input.
  92.  
  93. interflg = 0;
  94. if nypts ~= nupts,
  95.     interflg = 1;
  96. elseif max(abs(yt - ut)) > teps,
  97.     interflg = 1;
  98.     end
  99.  
  100. if interflg == 1,
  101.     uint = zeros(nru*nypts,1);
  102.     disp('interpolating input vector (zero order hold)')
  103.     if nypts <= nupts,
  104.         for i = 1:nypts,
  105.             idx = max(find(ut <= yt(i)+teps));
  106.             uint((i-1)*nru+1:i*nru) = udat((idx-1)*nru+1:idx*nru);
  107.             end
  108.     elseif nupts == 1,
  109.         uint = reshape(udat*ones(1,nypts),nru*nypts,1);
  110.     else
  111.         for i = 1:nupts-1,
  112.             idx = find(yt >= ut(i) & yt < ut(i+1));
  113.             uint((min(idx)-1)*nru+1:max(idx)*nru) = ...
  114.               reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(idx)), ...
  115.               nru*length(idx),1);
  116.             end
  117.         idx = find(yt >= ut(nupts));
  118.         uint((min(idx)-1)*nru+1:max(idx)*nru) = ...
  119.           reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(idx)), ...
  120.           nru*length(idx),1);
  121.         end
  122.     udat = uint;
  123.     end
  124.  
  125. %       reshape the input so that the matrix multiplications
  126. %       line up.  At this point the ut should correspond to
  127. %       yt.
  128.  
  129. udat = reshape(udat,nru,nypts);
  130.     
  131. %       the input vector has now been expanded to a matrix: u
  132. %       constant systems are treated by a multiplication call
  133.  
  134. %  Initialize the x & y matrices to speed up the code
  135.  
  136. y = zeros(nypts*ny+1,2);
  137. if type == 'syst',
  138.     x = [x0 zeros(length(x0),nypts-1)];
  139.     end
  140.  
  141. %       First, the single time point case is dealt with.  Here
  142. %       only the d term of a system and x0 affect the output.
  143.  
  144. if nypts == 1,
  145.     if type == 'syst',
  146.         y(1:ny,1) = c*x + d*udat;
  147.     else
  148.         y(1:ny,1) = sys*udat;
  149.         end
  150.  
  151. %       Now the multiple timepoint case is treated.  Here 
  152. %       a check is made to determine if the system has dynamics.
  153. %       The no dynamics case can be handled by a simple multiplication.
  154.  
  155. elseif type == 'syst',
  156.     for i=2:nypts,
  157.         x(:,i) = a*x(:,i-1) + b*udat(:,i-1);
  158.         end
  159.     y(1:ny*nypts,1) = reshape(c*x + d*udat,ny*nypts,1);
  160. else
  161.     y(1:ny*nypts,1) = reshape(sys*udat,ny*nypts,1);
  162.     end
  163.  
  164. %       In the case of only a single
  165. %       time point, a varying matrix is still created.  This is
  166. %       after all, assumed to be a time response.
  167.  
  168. y(1:nypts,2) = yt;
  169. y((ny*nypts)+1,2) = inf;
  170. y((ny*nypts)+1,1) = nypts;
  171.  
  172. %-----------------------------------------------------------------------------
  173.  
  174.  
  175. %
  176. % Copyright MUSYN INC 1991,  All Rights Reserved
  177.