home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / DBODE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.4 KB  |  149 lines
  1. function [magout,phase,w] = dbode(a,b,c,d,Ts,iu,w)
  2. %DBODE    Bode frequency response for discrete-time linear systems.
  3. %    DBODE(A,B,C,D,Ts,IU) produces a Bode plot from the single input IU
  4. %    to all the outputs of the discrete state-space system (A,B,C,D).
  5. %    IU is an index into the inputs of the system and specifies which
  6. %    input to use for the Bode response.  Ts is the sample period.  
  7. %    The frequency range and number of points are chosen automatically.
  8. %
  9. %    DBODE(NUM,DEN,Ts) produces the Bode 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. %    DBODE(A,B,C,D,Ts,IU,W) or DBODE(NUM,DEN,Ts,W) uses the user-
  14. %    supplied freq. vector W which must contain the frequencies, in 
  15. %    radians/sec, at which the Bode response is to be evaluated.  
  16. %    Aliasing will occur at frequencies greater than the Nyquist 
  17. %    frequency (pi/Ts).  See LOGSPACE to generate log. spaced frequency
  18. %    vectors.  With left hand arguments,
  19. %        [MAG,PHASE,W] = DBODE(A,B,C,D,Ts,...)
  20. %        [MAG,PHASE,W] = DBODE(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. 
  24. %    See also: LOGSPACE,SEMILOGX,MARGIN,DNICHOLS, and DNYQUIST.
  25.  
  26. %    J.N. Little 10-11-85
  27. %    Revised Andy Grace 8-15-89 2-4-91 6-20-92
  28. %    Revised Clay M. Thompson 7-10-90
  29. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  30.  
  31.  
  32. nargs = nargin;
  33. if nargs==0, eval('dexresp(''dbode'')'), return, end
  34.  
  35. error(nargchk(3,7,nargs));
  36.  
  37. % --- Determine which syntax is being used ---
  38. if (nargs==3),
  39.     num = a; den = b;
  40.     if length(c)==1,    % Transfer function without frequency vector
  41.         Ts=c;
  42.         w=dfrqint(num,den,Ts,30);
  43.     else          % Assume this is the old syntax DBODE(NUM,DEN,W)
  44.         disp('Warning: You are using the old syntax.  Use DBODE(NUM,DEN,Ts,W) instead.')
  45.         Ts = 1; w=c;
  46.     end
  47.     [ny,nn] = size(num); nu = 1;
  48.  
  49. elseif (nargs==4),     % Transfer function form with frequency vector
  50.     num = a; den = b;
  51.     Ts=c; w=d;
  52.     [ny,nn] = size(num); nu = 1;
  53.  
  54. elseif (nargs==5),    % State space system without iu or freq. vector
  55.     error(abcdchk(a,b,c,d));
  56.     w=dfrqint(a,b,c,d,Ts,30);
  57.     [iu,nargs,mag,phase]=dmulresp('dbode',a,b,c,d,Ts,w,nargout,1);
  58.     if ~iu, if nargout, magout = mag; end, return, end
  59.     [ny,nu] = size(d);
  60.  
  61. elseif (nargs==6),
  62.     error(abcdchk(a,b,c,d));
  63.     if length(iu)==1,    % State space system, with iu but w/o freq. vector
  64.         w=dfrqint(a,b,c,d,Ts,30);
  65.     else            % Assume this is the old syntax DBODE(A,B,C,D,IU,W)
  66.         disp('Warning: You are using the old syntax.  Use DBODE(A,B,C,D,Ts,IU,W) instead.')
  67.         w=iu; iu=Ts; Ts=1;
  68.     end
  69.     [ny,nu] = size(d);
  70.  
  71. else
  72.     error(abcdchk(a,b,c,d));
  73.     [ny,nu] = size(d);
  74.  
  75. end
  76.  
  77. if nu*ny==0, phase=[]; w=[]; if nargout~=0, magout=[]; end, return, end
  78.  
  79. % Compute frequency response
  80. if (nargs==3)|(nargs==4)
  81.     g=freqresp(num,den,exp(sqrt(-1)*w*Ts));
  82. else
  83.     g=freqresp(a,b,c,d,iu,exp(sqrt(-1)*w*Ts));
  84. end
  85. mag = abs(g);
  86. phase = (180/pi)*unwrap(atan2(imag(g),real(g)));
  87.  
  88. % Uncomment the following statement if you don't want the phase to  
  89. % be unwrapped.  Note that phase unwrapping may not always work; it
  90. % is only a "guess" as to whether +-360 should be added to the phase
  91. % to make it more aestheticly pleasing.
  92. % phase = (180/pi)*imag(log(g));
  93.  
  94. % If no left hand arguments then plot graph.
  95. if nargout==0
  96.     holdon = ishold;
  97.     if ~holdon
  98.         newplot;
  99.     end
  100.     subplot(211)
  101.     if holdon
  102.         hold on
  103.     end
  104.  
  105.     % Magnitude plot with 0db line
  106.     semilogx(w,20*log10(mag),[min(w(:)),max(w(:))],[0 0],'w:')
  107.     grid
  108.     xlabel('Frequency (rad/sec)')
  109.     ylabel('Gain dB')
  110.  
  111.     subplot(212)
  112.     semilogx(w,phase)
  113.  
  114.     xlabel('Frequency (rad/sec)')
  115.     ylabel('Phase deg')
  116.  
  117.     subplot(212)
  118.     semilogx(w,phase)
  119.  
  120.  
  121.     % Set tick marks up to be in multiples of 30, 90, 180, 360 ... degrees.
  122.     ytick = get(gca, 'ytick');
  123.     ylim = get(gca, 'ylim');
  124.     yrange = ylim(2) - ylim(1);
  125.     n = round(log(yrange/(length(ytick)*90))/log(2));
  126.  
  127.     set(gca, 'ylimmode', 'manual')
  128.  
  129.     if n >=0
  130.         ytick = [-90*2^n:-(90*2^n):ylim(1), 0:(90*2^n):ylim(2)];
  131.         set(gca,'ytick',ytick);
  132.     elseif n >= -2
  133.         ytick = [-30:-30:ylim(1), 0:30:ylim(2)];
  134.         set(gca,'ytick',ytick);
  135.     end
  136.  
  137.     grid
  138.     xlabel('Frequency (rad/sec)'), ylabel('Phase deg')
  139.     subplot(111)
  140.     return % Suppress output
  141. end
  142.  
  143. % Uncomment the following line for decibels, but be warned that the
  144. % MARGIN function will not work with decibel data.
  145. % mag = 20*log10(mag);
  146.  
  147. magout = mag; 
  148.  
  149.