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