home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / NICHOLS.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.3 KB  |  137 lines
  1. function [magout,phase,w] = nichols(a,b,c,d,iu,w)
  2. %NICHOLS Nichols frequency response for continuous-time linear systems.
  3. %    NICHOLS(A,B,C,D,IU) produces a Nichols plot from the single
  4. %    input IU to all the outputs of the continuous state-space system 
  5. %    (A,B,C,D).  IU is an index into the inputs of the system and 
  6. %    specifies which input to use for the Nichols response.  The 
  7. %    frequency range and number of points are chosen automatically.
  8. %
  9. %    NICHOLS(NUM,DEN) produces the Nichols plot for the polynomial 
  10. %    transfer function G(s) = NUM(s)/DEN(s) where NUM and DEN contain
  11. %    the polynomial coefficients in descending powers of s. 
  12. %
  13. %    NICHOLS(A,B,C,D,IU,W) or NICHOLS(NUM,DEN,W) uses the user-supplied
  14. %    freq. vector W which must contain the frequencies, in radians/sec,
  15. %    at which the Nichols response is to be evaluated.  When invoked 
  16. %    with left hand arguments,
  17. %        [MAG,PHASE,W] = NICHOLS(A,B,C,D,...)
  18. %        [MAG,PHASE,W] = NICHOLS(NUM,DEN,...) 
  19. %    returns the frequency vector W and matrices MAG and PHASE (in 
  20. %    degrees) with as many columns as outputs and length(W) rows.  No 
  21. %    plot is drawn on the screen.  Draw the Nichols grid with NGRID.
  22. %
  23. %    See also: LOGSPACE,NGRID,MARGIN,BODE, and NYQUIST.
  24.  
  25. %    Clay M. Thompson  7-6-90
  26. %    Revised ACWG 6-21-92
  27. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  28.  
  29. if nargin==0, eval('exresp(''nichols'')'), return, end
  30.  
  31. error(nargchk(2,6,nargin));
  32.  
  33. % --- Determine which syntax is being used ---
  34. if (nargin==1),
  35.     error('Wrong number of input arguments.');
  36.  
  37. elseif (nargin==2),    % Transfer function form without frequency vector
  38.     num = a; den = b;
  39.     w = freqint2(num,den,30);
  40.     [ny,nn] = size(num); nu=1;
  41.  
  42. elseif (nargin==3),    % Transfer function form with frequency vector
  43.     num = a; den = b;
  44.     w = c;
  45.     [ny,nn] = size(num); nu=1;
  46.  
  47. elseif (nargin==4),    % State space system, w/o iu or frequency vector
  48.     error(abcdchk(a,b,c,d));
  49.     w = freqint2(a,b,c,d,30);
  50.     [iu,nargin,mag,phase]=mulresp('nichols',a,b,c,d,w,nargout,1);
  51.     if ~iu, if nargout, magout = mag; end, return, end
  52.     [ny,nu] = size(d);
  53.  
  54. elseif (nargin==5),    % State space system, with iu but w/o freq. vector
  55.     error(abcdchk(a,b,c,d));
  56.     w = freqint2(a,b,c,d,30);
  57.     [ny,nu] = size(d);
  58.  
  59. else
  60.     error(abcdchk(a,b,c,d));
  61.     [ny,nu] = size(d);
  62.  
  63. end
  64.  
  65. if nu*ny==0, phase=[]; w=[]; if nargout~=0, magout=[]; end, return, end
  66.  
  67. % --- Evaluate the frequency response ---
  68. if (nargin==2)|(nargin==3),
  69.     g = freqresp(num,den,sqrt(-1)*w);
  70. else
  71.     g = freqresp(a,b,c,d,iu,sqrt(-1)*w);
  72. end
  73.  
  74. mag = abs(g);
  75. phase = -180+180./pi*atan2(-imag(g),-real(g));
  76.  
  77. % If no left hand arguments then plot graph.
  78. if nargout==0
  79.     status = ishold;
  80.     magdb = 20*log10(mag);
  81.  
  82.     % Find the points where the phase wraps.  Plot each branch separately.
  83.     % The following code is based on the algorithm in the unwrap function.
  84.  
  85.     cutoff = 200;
  86.     [m, n] = size(phase);                % Assume column orientation.
  87.  
  88.     pmin = min(phase); pmin = pmin(ones(m, 1), :);% To force REM to behave.
  89.     p = rem(phase - pmin, 360) + pmin;            % Phases modulo 360.
  90.  
  91.     dp = [p(1,:);diff(p)];            % Differentiate phases.
  92.     jumps = (dp > cutoff) + (dp < -cutoff);    % Array of jumps locations
  93.     if n==1, jvec = jumps; else jvec = (sum(jumps')~=0)'; end
  94.     seg = cumsum(jvec);                % Used to identify segments
  95.  
  96.     jloc = find(jvec~=0);
  97.  
  98.     % Strip off each continuous piece and plot separately
  99.     pj=[]; mj=[]; jj=[];
  100.     for k=0:seg(length(seg)),
  101.         mseg = magdb(seg==k,:); pseg = p(seg==k,:);
  102.  
  103.         % Add last point of previous segments so plots are continuous
  104.         if ~isempty(pj), pj(jj)=pseg(1,jj); mj(jj)=mseg(1,jj); end
  105.         pseg=[pj;pseg]; mseg=[mj;mseg];
  106.  
  107.         % Remember last point
  108.         [np,mp] = size(pseg);
  109.         if np~=0,
  110.             pj = pseg(np,:); mj = mseg(np,:); 
  111.             if k<length(jloc), jj=(jumps(jloc(k+1),:)~=0); end
  112.         end
  113.  
  114.         if np>1,
  115.             plot(pseg,mseg),         % Plot segment
  116.             hold on
  117.         end
  118.     end
  119.  
  120.     axstate = axis('state');
  121.     if axstate(1)=='a'  % axis is set to auto range
  122.         set(gca,'xlim',[-360 0]);
  123.     end
  124.     set(gca,'xtick',[-360, -270, -180, -90 0]);
  125.     xlabel('Open-Loop Phase (deg)')
  126.     ylabel('Open-Loop Gain (db)')
  127.  
  128.     if ~status, hold off, end    % Return hold to previous status
  129.     return % Suppress output
  130. end
  131.  
  132. % Uncomment the following line for decibels, but be warned that the 
  133. % MARGIN function will not work with decibel data.
  134.  
  135. % mag = 20*log10(mag);
  136. magout = mag; 
  137.