home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / DNYQUIST.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.7 KB  |  116 lines
  1. function [reout,im,w] = dnyquist(a,b,c,d,Ts,iu,w)
  2. %DNYQUIST Nyquist frequency response for discrete-time linear systems.
  3. %    DNYQUIST(A,B,C,D,Ts,IU) produces a Nyquist plot from the single
  4. %    input IU to all the outputs of the system:      -1
  5. %         x[n+1] = Ax[n] + Bu[n]    G(w) = C(exp(jwT)I-A) B + D  
  6. %         y[n]   = Cx[n] + Du[n]    RE(w) = real(G(w)), IM(w) = imag(G(w))
  7. %    The frequency range and number of points are chosen automatically.
  8. %
  9. %    DNYQUIST(NUM,DEN,Ts) produces the Nyquist plot for the polynomial
  10. %    transfer function G(z) = NUM(z)/DEN(z) where NUM and DEN contain
  11. %    the polynomial coefficients in descending powers of z. 
  12. %
  13. %    DNYQUIST(A,B,C,D,Ts,IU,W) or DNYQUIST(NUM,DEN,Ts,W) uses the user-
  14. %    supplied freq. vector W which must contain the frequencies, in 
  15. %    rad/sec, at which the Nyquist response is to be evaluated.  
  16. %    Aliasing will occur at frequencies greater than the Nyquist 
  17. %    frequency (pi/Ts). With left hand arguments,
  18. %        [RE,IM,W] = DNYQUIST(A,B,C,D,Ts,...)
  19. %        [RE,IM,W] = DNYQUIST(NUM,DEN,Ts,...) 
  20. %    returns the frequency vector W and matrices RE and IM with as many
  21. %       columns as outputs and length(W) rows.  No plot is drawn on the 
  22. %    screen.
  23. %
  24. %    See also: LOGSPACE,MARGIN,DBODE, and DNICHOLS.
  25.  
  26. %     J.N. Little 10-11-85
  27. %    Revised A.C.W.Grace 8-15-89, 2-12-91, 6-21-92
  28. %    Revised Clay M. Thompson 7-9-90
  29. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  30.  
  31. if nargin==0, eval('dexresp(''dnyquist'')'), return, end
  32.  
  33. error(nargchk(3,7,nargin));
  34.  
  35. % --- Determine which syntax is being used ---
  36. if (nargin==3),        % Transfer function without frequency vector
  37.     num = a; den = b;
  38.     Ts=c;
  39.     w=dfrqint2(num,den,Ts,30);
  40.     [ny,nn] = size(num); nu = 1;
  41.  
  42. elseif (nargin==4),     % Transfer function form with frequency vector
  43.     num = a; den = b;
  44.     Ts=c; w=d;
  45.     [ny,nn] = size(num); nu = 1;
  46.  
  47. elseif (nargin==5),    % State space system without iu or freq. vector
  48.     error(abcdchk(a,b,c,d));
  49.     w=dfrqint2(a,b,c,d,Ts,30);
  50.     [iu,nargin,re,im]=dmulresp('dnyquist',a,b,c,d,Ts,w,nargout,0);
  51.     if ~iu, if nargout, reout = re; end, return, end
  52.     [ny,nu] = size(d);
  53.  
  54. elseif (nargin==6),    % State space system, with iu but w/o freq. vector
  55.     error(abcdchk(a,b,c,d));
  56.     w=dfrqint2(a,b,c,d,Ts,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, im=[]; w=[]; if nargout~=0, reout=[]; end, return, end
  66.  
  67. % Compute frequency response
  68. if (nargin==3)|(nargin==4)
  69.     g=freqresp(num,den,exp(sqrt(-1)*w*Ts));
  70. else
  71.     g=freqresp(a,b,c,d,iu,exp(sqrt(-1)*w*Ts));
  72. end
  73. re = real(g);
  74. im = imag(g);
  75.  
  76. % If no left hand arguments then plot graph.
  77. if nargout==0,
  78.     status = ishold;
  79.  
  80.     % Make arrows
  81.     sx=max(re)-min(re);   [sy,sample]=max(abs(2*im));
  82.     arrow=[-1;0;-1]*sx/50+sqrt(-1)*[1;0;-1]*sy/50;
  83.     sample=sample+(sample==1);
  84.     reim=diag(g(sample,:));
  85.     d=diag(g(sample,:)-g(sample-1,:)); 
  86.     rot = sign(d);  % Arrow is always horizontal i.e. -1 or +1
  87.     % rot=d./abs(d)  Use this when arrow is not horizontal
  88.     xy =ones(3,1)*reim' + ones(3,1)*rot'.*arrow;
  89.     xy2=ones(3,1)*reim' - ones(3,1)*rot'.*arrow;
  90.     [m,n]=size(g);
  91.     if (n==1)
  92.         plot(re,im,'r-',re,-im,'g--',real(xy),-imag(xy),'r-',real(xy2),imag(xy2),'g-')
  93.     else
  94.         plot(re,im,re,-im,'--');
  95.         hold on
  96.         plot(real(xy),-imag(xy),'-',real(xy2),imag(xy2),'-')
  97.     end
  98.     hold on
  99.  
  100.     xlabel('Real Axis'), ylabel('Imag Axis')
  101.  
  102.     % Make cross at s = -1 + j0, i.e the -1 point
  103.     limits = axis;
  104.     line1 = (limits(2)-limits(1))/50;
  105.     line2 = (limits(4)-limits(3))/50;
  106.     plot([-1+line1, -1-line1], [0,0], 'w-', [-1, -1], [line2, -line2], 'w-')
  107.  
  108. % Axis
  109.     limits = axis;
  110.     plot([limits(1:2);0,0]',[0,0;limits(3:4)]','w:');
  111.  
  112.     if ~status, hold off, end    % Return hold to previous status
  113.     return % Suppress output
  114. end
  115. reout = re;  
  116.