home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / DSIGMA.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.9 KB  |  143 lines
  1. function [svout,w] = dsigma(Z1,Z2,Z3,Z4,Z5,Z6,Z7)
  2. %function [svout,w] = dsigma(a,b,c,d,Ts,w,invflag)
  3. %DSIGMA Singular value frequency response of discrete-time linear systems.
  4. %    DSIGMA(A,B,C,D,Ts) (or optional SIGMA(SS_,Ts) in RCT) produces a 
  5. %       singular value plot of matrix:   -1
  6. %                    G(w) = C(exp(jwT)I-A) B + D 
  7. %    as a function of frequency.  The singular values are an extension
  8. %    of Bode magnitude response for MIMO system.  The frequency range 
  9. %       and number of points are chosen automatically. For square systems, 
  10. %       DSIGMA(A,B,C,D,Ts,'inv') produces the singular values of the inverse 
  11. %       matrix     -1                    -1      -1
  12. %              G (w) = [ C(exp(jwT)I-A) B + D ]
  13. %    DSIGMA(A,B,C,D,Ts,W) or DSIGMA(A,B,C,D,Ts,W,'inv') use the user-
  14. %    supplied frequency vector W which must contain the frequencies, in
  15. %    radians/sec, at which the singular value response is to be 
  16. %    evaluated. Aliasing will occur at frequencies greater than the 
  17. %    Nyquist frequency (pi/Ts).  When invoked with left hand arguments,
  18. %           [SV,W] = DSIGMA(A,B,C,D,Ts,...)
  19. %    or  [SV,W] = SIGMA(SS_,Ts,...)    (for Robust Control Toolbox user)
  20. %    returns the frequency vector W and the matrix SV with MIN(NU,NY)
  21. %    columns and length(W) rows, where NU is the number of inputs and
  22. %    NY is the number of outputs.  No plot is drawn on the screen. 
  23. %    The singular values are returned in descending order.
  24. %
  25. %    See also: LOGSPACE, SEMILOGX, DNICHOLS, DNYQUIST and DBODE.
  26.  
  27. %    Clay M. Thompson  7-10-90
  28. %    Revised A.Grace 2-12-91, 6-21-92
  29. %    Revised W.Wang 7/24/92
  30. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  31.  
  32. if nargin==0, eval('dexresp(''dsigma'',1)'), return, end
  33.  
  34. if nargin <= 4  %if not enough parameters, check rct data structure
  35.   cmd=[];
  36.   %The following if..end is add to match up RCT format change
  37.   if exist('mkargs') == 2, %If RCT installed
  38.     inargs='(a,b,c,d,Ts,w,invflag)';
  39.     eval('cmd=mkargs(inargs,nargin,''ss'');')
  40.     eval(cmd);
  41.   end;
  42. else                       % If RCT not installed
  43.   % assign a=Z1;...,w=Z6;invflag=Z7;
  44.   a=Z1; b=Z2; c=Z3; d=Z4;
  45.   if nargin>6,
  46.     invflag=Z7; w=Z6; Ts=Z5;
  47.   elseif nargin>5,
  48.     w=Z6; Ts=Z5;
  49.   elseif nargin>4,
  50.     Ts=Z5;
  51.   end
  52. end;
  53.  
  54. if nargin==7,        % Trap call to RCT function
  55.     if ~isstr(invflag),
  56.         eval('svout = dsigma2(a,b,c,d,Ts,w,invflag);')
  57.         return
  58.     end
  59. end
  60.  
  61. error(nargchk(5,7,nargin));
  62. error(abcdchk(a,b,c,d));
  63. if ~(length(d) | (length(b) &  length(c)))
  64.     return;
  65. end
  66.  
  67. % Determine status of invflag
  68. if nargin==5, 
  69.     invflag = [];
  70.     w = [];
  71. elseif (nargin==6)
  72.     if (isstr(w)),
  73.         invflag = w;
  74.         w = [];
  75.         [ny,nu] = size(d);
  76.         if (ny~=nu), error('The state space system must be square when using ''inv''.'); end
  77.     else
  78.         invflag = [];
  79.     end
  80.  
  81. else
  82.     [ny,nu] = size(d);
  83.     if (ny~=nu), error('The state space system must be square when using ''inv''.'); end
  84. end
  85.  
  86. % Generate frequency range if one is not specified.
  87.  
  88. % If frequency vector supplied then use Auto-selection algorithm
  89. % Fifth argument determines precision of the plot.
  90. if ~length(w)  
  91.     w=dfrqint(a,b,c,d,Ts,30);
  92. end
  93.  
  94. [nx,na] = size(a);
  95. [no,ns] = size(c);
  96. nw = max(size(w));
  97.  
  98. % Balance A
  99. [t,a] = balance(a);
  100. b = t \ b;
  101. c = c * t;
  102.  
  103. % Reduce A to Hessenberg form:
  104. [p,a] = hess(a);
  105.  
  106. % Apply similarity transformations from Hessenberg
  107. % reduction to B and C:
  108. b = p' * b;
  109. c = c * p;
  110.  
  111. s = exp(sqrt(-1)*w*Ts);
  112. I=eye(length(a));
  113. if nx > 0,
  114.     for i=1:length(s)
  115.         if isempty(invflag),
  116.             sv(:,i)=svd(c*((s(i)*I-a)\b) + d);
  117.         else
  118.             sv(:,i)=svd(inv(c*((s(i)*I-a)\b) + d));
  119.         end
  120.     end
  121. else
  122.     for i=1:length(s)
  123.         if isempty(invflag),
  124.             sv(:,i)=svd(d);
  125.         else
  126.             sv(:,i)=svd(inv(d));
  127.         end
  128.     end
  129. end
  130.  
  131. % If no left hand arguments then plot graph.
  132. if nargout==0
  133.     semilogx(w, 20*log10(sv), [w(1), w(length(w))], [0,0], 'w:')
  134.     xlabel('Frequency (rad/sec)')
  135.     ylabel('Singular Values dB')
  136.     return % Suppress output
  137. end
  138. svout =sv; 
  139.  
  140.  
  141.  
  142.  
  143.