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