home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / bode.m < prev    next >
Text File  |  1999-04-29  |  5KB  |  152 lines
  1. # Copyright (C) 1996,1998 A. Scottedward Hodel 
  2. #
  3. # This file is part of Octave. 
  4. #
  5. # Octave is free software; you can redistribute it and/or modify it 
  6. # under the terms of the GNU General Public License as published by the 
  7. # Free Software Foundation; either version 2, or (at your option) any 
  8. # later version. 
  9. # 
  10. # Octave is distributed in the hope that it will be useful, but WITHOUT 
  11. # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
  12. # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
  13. # for more details.
  14. # 
  15. # You should have received a copy of the GNU General Public License 
  16. # along with Octave; see the file COPYING.  If not, write to the Free 
  17. # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 
  18.  
  19. function [mag,phase,w] = bode(sys,w,outputs,inputs)
  20. # [mag,phase,w] = bode(sys[,w,outputs,inputs])
  21. # Produce Bode plots of a system
  22. #
  23. # Compute the frequency response of a system.
  24. # inputs:
  25. #   sys: system data structure (must be either purely continuous or discrete;
  26. #     see is_digit)
  27. #   w: frequency values for evaluation.
  28. #      if sys is continuous, then bode evaluates G(jw)
  29. #      if sys is discrete, then bode evaluates G(exp(jwT)), where T=sys.tsam
  30. #         (the system sampling time)
  31. #      default: the default frequency range is selected as follows: (These
  32. #        steps are NOT performed if w is specified)
  33. #          (1) via routine bodquist, isolate all poles and zeros away from
  34. #              w=0 (jw=0 or exp(jwT)=1) and select the frequency
  35. #             range based on the breakpoint locations of the frequencies.
  36. #          (2) if sys is discrete time, the frequency range is limited
  37. #              to jwT in [0,2p*pi]
  38. #          (3) A "smoothing" routine is used to ensure that the plot phase does
  39. #              not change excessively from point to point and that singular
  40. #              points (e.g., crossovers from +/- 180) are accurately shown.
  41. #   outputs, inputs: the indices of the output(s) and input(s) to be used in
  42. #     the frequency response; see sysprune.
  43. # outputs:
  44. #    mag, phase: the magnitude and phase of the frequency response
  45. #       G(jw) or G(exp(jwT)) at the selected frequency values.
  46. #    w: the vector of frequency values used
  47. # If no output arguments are given, bode plots the results to the screen.
  48. # Descriptive labels are automatically placed.  See xlabel, ylable, title,
  49. # and replot.
  50. #
  51. # Note: if the requested plot is for an MIMO system, mag is set to
  52. # ||G(jw)|| or ||G(exp(jwT))|| and phase information is not computed.
  53.  
  54. # Written by John Ingram  July 10th, 1996
  55. # Based on previous code
  56. # By R. Bruce Tenison, July 13, 1994
  57. # Modified by David Clem November 13, 1994
  58. # again by A. S. Hodel July 1995 (smart plot range, etc.)
  59. # Modified by Kai P. Mueller September 28, 1997 (multiplt mode)
  60.  
  61.   # check number of input arguments given
  62.   if (nargin < 1 | nargin > 4)
  63.     usage("[mag,phase,w] = bode(sys[,w,outputs,inputs])");
  64.   endif
  65.   if(nargin < 2)
  66.     w = [];
  67.   endif
  68.   if(nargin < 3)
  69.     outputs = [];
  70.   endif
  71.   if(nargin < 4)
  72.     inputs = [];
  73.   endif
  74.  
  75.   [f, w] = bodquist(sys,w,outputs,inputs,"bode");
  76.  
  77.   [stname,inname,outname] = sysgetsg(sys);
  78.   systsam = sysgetts(sys);
  79.  
  80.   # Get the magnitude and phase of f.
  81.   mag = abs(f);
  82.   phase = arg(f)*180.0/pi;
  83.  
  84.   if (nargout < 1),
  85.     # Plot the information
  86.     if(gnuplot_has_multiplt)
  87.       oneplot();
  88.     endif
  89.     gset autoscale;
  90.     if(gnuplot_has_multiplt)
  91.       gset nokey;
  92.     endif
  93.     clearplot();
  94.     gset data style lines;
  95.     if(is_digit(sys))
  96.       xlstr = ["Digital frequency w=rad/sec.  pi/T=",num2str(pi/systsam)];
  97.       tistr = "(exp(jwT)) ";
  98.     else
  99.       xlstr = "Frequency in rad/sec";
  100.       tistr = "(jw)";
  101.     endif
  102.     xlabel(xlstr);
  103.     if(is_siso(sys))
  104.       if (gnuplot_has_multiplt)
  105.         subplot(2,1,1);
  106.       endif
  107.       title(["|[Y/U]",tistr,"|, u=", nth(inname,1),", y=",nth(outname,1)]);
  108.     else
  109.       title([ "||Y(", tistr, ")/U(", tistr, ")||"]);
  110.       disp("MIMO plot from")
  111.       disp(outlist(inname,"    "));
  112.       disp("to")
  113.       disp(outlist(outname,"    "));
  114.     endif
  115.     wv = [min(w), max(w)];
  116.     if(max(mag) > 0)
  117.       ylabel("Gain in dB");
  118.       md = 20*log10(mag);
  119.     else
  120.       ylabel("Gain |Y/U|")
  121.       md = mag;
  122.     endif
  123.  
  124.     axvec = ax2dlim([vec(w),vec(md)]);
  125.     axvec(1:2) = wv;
  126.     axis(axvec);
  127.     grid("on");
  128.     semilogx(w,md);
  129.     if (is_siso(sys))
  130.       if (gnuplot_has_multiplt)
  131.         subplot(2,1,2);
  132.       else
  133.         prompt('Press any key for phase plot');
  134.       endif
  135.       axvec = ax2dlim([vec(w),vec(phase)]);
  136.       axvec(1:2) = wv;
  137.       axis(axvec);
  138.       xlabel(xlstr);
  139.       ylabel("Phase in deg");
  140.       title([ "phase([Y/U]", tistr, ...
  141.      "), u=", nth(inname,1),", y=",nth(outname,1)]);
  142.       grid("on");
  143.       semilogx(w,phase);
  144.       # This should be the default for subsequent plot commands.
  145.       if(gnuplot_has_multiplt)
  146.         oneplot();
  147.       endif
  148.     endif
  149.     mag = phase = w = [];
  150.   endif
  151. endfunction
  152.