home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / FBODE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  5.4 KB  |  179 lines
  1. function [magout,phase,w] = fbode(a,b,c,d,iu,w)
  2. %FBODE  Bode frequency response for continuous-time linear systems.
  3. %    FBODE(A,B,C,D,IU) produces a Bode plot from the single input IU to
  4. %    all the outputs of the continuous 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.  The frequency range and 
  7. %    number of points are chosen automatically.  FBODE is a faster, but
  8. %    less accurate, version of BODE.
  9. %
  10. %    FBODE(NUM,DEN) produces the Bode plot for the polynomial transfer
  11. %    function G(s) = NUM(s)/DEN(s) where NUM and DEN contain the 
  12. %    polynomial coefficients in descending powers of s. 
  13. %
  14. %    FBODE(A,B,C,D,IU,W) or FBODE(NUM,DEN,W) uses the user-supplied 
  15. %    frequency vector W which must contain the frequencies, in 
  16. %    radians/sec, at which the Bode response is to be evaluated.  See 
  17. %    LOGSPACE to generate log. spaced frequency vectors.  When invoked
  18. %    with left hand arguments,
  19. %        [MAG,PHASE,W] = FBODE(A,B,C,D,...)
  20. %        [MAG,PHASE,W] = FBODE(NUM,DEN,...) 
  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 and BODE.
  25.  
  26. %     J.N. Little 12-5-88
  27. %    Revised CMT 8-2-90, ACWG 6-21-92
  28. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  29.  
  30.  
  31. nargs = nargin;
  32. if nargs==0, eval('exresp(''fbode'')'), return, end
  33.  
  34. error(nargchk(2,6,nargs));
  35.  
  36. % --- Determine which syntax is being used ---
  37. if (nargs==1),
  38.     error('Wrong number of input arguments.');
  39.  
  40. elseif (nargs==2),    % Transfer function form without frequency vector
  41.     num = a; den = b;
  42.     w = freqint(num,den,20);
  43.     [ny,nn] = size(num); nu = 1;
  44.  
  45. elseif (nargs==3),    % Transfer function form with frequency vector
  46.     num = a; den = b;
  47.     w = c;
  48.     [ny,nn] = size(num); nu = 1;
  49.  
  50. elseif (nargs==4),    % State space system, w/o iu or frequency vector
  51.     error(abcdchk(a,b,c,d));
  52.     w = freqint(a,b,c,d,30);
  53.     [iu,nargs,mag,phase]=mulresp('bode',a,b,c,d,w,nargout,1);
  54.     if ~iu, if nargout, magout = mag; end, return, end
  55.     [ny,nu] = size(d);
  56.  
  57. elseif (nargs==5),    % State space system, with iu but w/o freq. vector
  58.     error(abcdchk(a,b,c,d));
  59.     w = freqint(a,b,c,d,30);
  60.     [ny,nu] = size(d);
  61.  
  62. else            % State space system, with iu and freq. vector
  63.     error(abcdchk(a,b,c,d));
  64.     [ny,nu] = size(d);
  65.     
  66. end
  67.  
  68. if nu*ny==0, phase=[]; w=[]; if nargout~=0, magout=[]; end, return, end
  69.  
  70. % --- Evaluate the frequency response ---
  71. if (nargs == 2)|(nargs == 3)
  72.     % It is in transfer function form.  Do directly, using Horner's method
  73.     % of polynomial evaluation at the frequency points, for each row in
  74.     % the numerator.  Then divide by the denominator.
  75.     [ma,na] = size(num);
  76.     s = sqrt(-1)*w(:);
  77.     for i=1:ma
  78.         g(:,i) = polyval(num(i,:),s);
  79.     end
  80.     g = polyval(den,s)*ones(1,ma).\g;
  81. else
  82.     [no,nu] = size(d);
  83.     [ns,na] = size(a);
  84.     nw = max(size(w));
  85.  
  86.     [p,a] = eig(a);       % Reduce A to diagonal form
  87.     if rcond(p)<sqrt(eps), 
  88.         disp('Warning: Diagonalization matrix singular.  Use BODE instead.')
  89.     end
  90.     % Apply similarity transformations to B and C
  91.     if ns>0,
  92.         b = p \ b(:,iu);  
  93.         c = c * p;
  94.         d = d(:,iu);
  95.         s = (w*sqrt(-1)).';
  96.         s = s(ones(ns,1),:);
  97.         a2 = diag(a)*ones(1,nw);
  98.         b = b*ones(1,nw);
  99.         g = (c*((1 ./(s-a2)).*b) + diag(d)*ones(no,nw)).';
  100.     else
  101.         d = d(:,iu);
  102.         g = (diag(d)*ones(no,nw)).';
  103.     end                
  104. end
  105.  
  106. mag = abs(g);
  107. phase = (180./pi)*unwrap(atan2(imag(g),real(g)));
  108.  
  109. % Uncomment out the following statement if you don't want the phase to  
  110. % be unwrapped.  Note that phase unwrapping will not always work; it is
  111. % only a "guess" as to whether +-360 should be added to the phase 
  112. % to make it more aesthetically pleasing.  (See UNWRAP.M)
  113.  
  114. %phase = (180./pi)*atan2(imag(g),real(g));
  115.  
  116. %Try to correct phase anomaly for plants with integrators
  117. %by adding multiples of -360 degrees.
  118. if (nargs == 2 | nargs == 3) 
  119.     nd = length(den);
  120.     f = find(fliplr(den) == zeros(1,nd));
  121.     nintgs = sum(f == 1:length(f));    
  122.     if phase(1) > 0 & nintgs > 0
  123.         phase = phase - 360;
  124.     end
  125. else 
  126.     if abs(det(a)) < eps
  127.         if phase(1) > 0 & nintgs > 0
  128.             phase = phase - 360;
  129.         end
  130.     end
  131. end
  132.  
  133. % If no left hand arguments then plot graph.
  134. if nargout==0
  135.     holdon = ishold;
  136.     subplot(211) 
  137.     if holdon
  138.         hold on
  139.     end
  140.     semilogx(w,20*log10(mag),w,zeros(1,length(w)),'w:')
  141.     % If hold is set to on on the current axis then set it on the first axis too.
  142.     % This enables two bode response to be superimposed on each other with
  143.     % the following commands:
  144.     %    bode(num, den); hold on; bode(num2, den2)
  145.     grid
  146.     xlabel('Frequency (rad/sec)'), ylabel('Gain dB')
  147.  
  148.     subplot(212), 
  149.     semilogx(w,phase)
  150.     xlabel('Frequency (rad/sec)'), ylabel('Phase deg')
  151.  
  152.  
  153.     % Set tick marks up to be in multiples of 30, 90, 180, 360 ... degrees.
  154.     ytick = get(gca, 'ytick');
  155.     ylim = get(gca, 'ylim');
  156.     yrange = ylim(2) - ylim(1);
  157.     n = round(log(yrange/(length(ytick)*90))/log(2));
  158.  
  159.     set(gca, 'ylimmode', 'manual')
  160.  
  161.     if n >=0 
  162.         ytick = [-90*2^n:-(90*2^n):ylim(1), 0:(90*2^n):ylim(2)];
  163.         set(gca,'ytick',ytick);
  164.     elseif n >= -2 
  165.         ytick = [-30:-30:ylim(1), 0:30:ylim(2)];
  166.         set(gca,'ytick',ytick);
  167.     end
  168.     grid
  169.     % Reset the graph: subplot(111)
  170.     subplot(111)
  171.     return % Suppress output 
  172. end
  173.  
  174. % Uncomment the following line for decibels, but be warned that the
  175. % MARGIN function will not work with decibel data.
  176. % mag = 20*log10(mag);
  177.  
  178. magout = mag; 
  179.