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