home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / FRSP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.5 KB  |  94 lines
  1. % function out = frsp(sys,omega,T,balflg)
  2. %
  3. %   Complex frequency response of a SYSTEM matrix, produces
  4. %   a VARYING matrix which contains its frequency response:
  5. %
  6. %    SYS     - SYSTEM matrix
  7. %    OMEGA   - response calculated at these frequencies, can
  8. %                also be another VARYING matrix, and then these
  9. %                independent variables are used
  10. %    T       - 0 (default) continuous system, if T ~= 0 then
  11. %          discrete evaluation with sample time T is used.
  12. %    BALFLG  - 0 (default) balances A matrix prior to eval, nonzero
  13. %                 value for BALFLG leaves state-space unchanged
  14. %    OUT     - VARYING frequency response matrix
  15. %
  16. %   See also: BALANCE, HESS, SAMHLD, TUSTIN, and VPLOT.
  17.  
  18. function out = frsp(sys,omega,T,balflg)
  19.   if nargin < 2
  20.     disp('usage: out = frsp(sys,omega)')
  21.     return
  22.   elseif nargin == 2
  23.    balflg = 0;
  24.    T = 0;
  25.   elseif nargin == 3
  26.    balflg = 0;
  27.   end
  28.   [mtype,mrows,mcols,mnum] = minfo(sys);
  29.   if mtype == 'vary'
  30.     error('input matrix is already a VARYING matrix')
  31.     return
  32.   elseif mtype == 'empt'
  33.     out = [];
  34.     return
  35.   end
  36.   [omtype,omrows,omcols,omnum] = minfo(omega);
  37.   if omtype == 'cons'
  38.     if omrows == 1
  39.       omega = omega.';
  40.     end
  41.     npts = length(omega);
  42.   elseif omtype == 'vary'
  43.     omega = omega(1:omnum,omcols+1);
  44.     npts = omnum;
  45.   else
  46.     error('omega should be a vector, or VARYING matrix')
  47.     return
  48.   end
  49.   if mtype == 'cons'
  50.     out == [];
  51.     for i=1:npts
  52.       out = [out; sys];
  53.     end
  54.     out = [out [omega;zeros(npts*(mrows-1),1)];zeros(1,mcols-1) npts inf];
  55.     return
  56.   end
  57.   if isempty(omega)
  58.     out = [];
  59.     return
  60.   else
  61.     comega = sqrt(-1) * omega;
  62.     if T ~= 0
  63.       comega = exp(T*comega);
  64.     end
  65.     [a,b,c,d] = unpck(sys);
  66.     if nargin == 2 | balflg == 0
  67.       [t,a] = balance(a);
  68.       b = t\b;
  69.       c = c*t;
  70.     end
  71.     [states dum]=size(a);
  72.     [outputs inputs] = size(d);
  73.     [p,hesa] = hess(a);
  74.     hesb = p' * b;
  75.     hesc = c * p;
  76.     npts = length(comega);
  77.     idn = eye(states);
  78.     nrout = mrows;
  79.     ncout = mcols;
  80.     out = zeros((nrout*npts)+1,ncout+1);
  81.     fout = (npts+1)*nrout;
  82.     pout = 1:nrout:fout;
  83.     poutm1 = pout(2:npts+1) - 1; 
  84.     for i=1:npts
  85.       g = (comega(i) * idn - hesa) \ hesb;
  86.       out(pout(i):poutm1(i),1:ncout) = d + hesc*g;
  87.     end
  88.     out(1:npts,ncout+1) = omega;
  89.     out((nrout*npts)+1,ncout+1) = inf;
  90.     out((nrout*npts)+1,ncout) = npts;
  91.   end
  92. %
  93. % Copyright MUSYN INC 1991,  All Rights Reserved
  94.