home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / bode.m < prev    next >
Text File  |  1999-12-24  |  7KB  |  211 lines
  1. ## Copyright (C) 1996,1998 Auburn University.  All Rights Reserved
  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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. 
  18.  
  19. ## -*- texinfo -*-
  20. ## @deftypefn {Function File } {[@var{mag}, @var{phase}, @var{w}] =} bode(@var{sys}@{,@var{w}, @var{out_idx}, @var{in_idx}@})
  21. ## If no output arguments are given: produce Bode plots of a system; otherwise,
  22. ## compute the frequency response of a system data structure
  23. ## 
  24. ## @strong{Inputs}
  25. ## @table @var
  26. ## @item   sys
  27. ##  a system data structure (must be either purely continuous or discrete;
  28. ##      see is_digit)
  29. ## @item   w
  30. ##  frequency values for evaluation.
  31. ## 
  32. ## if @var{sys} is continuous, then bode evaluates @math{G(jw)} where
  33. ## @math{G(s)} is the system transfer function.
  34. ## 
  35. ## if @var{sys} is discrete, then bode evaluates G(@code{exp}(jwT)), where 
  36. ## @itemize @bullet
  37. ## @item @var{T}=@code{sysgetts(@var{sys})} (the system sampling time) and
  38. ## @item @math{G(z)} is the system transfer function.
  39. ## @end itemize
  40. ## 
  41. ## @strong{ Default} the default frequency range is selected as follows: (These
  42. ##         steps are NOT performed if @var{w} is specified)
  43. ## @enumerate
  44. ## @item via routine bodquist, isolate all poles and zeros away from
  45. ## @var{w}=0 (@var{jw}=0 or @math{@code{exp}(jwT)}=1) and select the frequency
  46. ## range based on the breakpoint locations of the frequencies.
  47. ## @item if @var{sys} is discrete time, the frequency range is limited
  48. ##               to @math{jwT} in 
  49. ## @ifinfo
  50. ## [0,2 pi /T]
  51. ## @end ifinfo
  52. ## @iftex
  53. ## @tex 
  54. ## $[0,2\pi/T]$
  55. ## @end tex
  56. ## @end iftex
  57. ## @item A "smoothing" routine is used to ensure that the plot phase does
  58. ##               not change excessively from point to point and that singular
  59. ##               points (e.g., crossovers from +/- 180) are accurately shown.
  60. ## 
  61. ## @end enumerate
  62. ## @item out_idx, in_idx
  63. ##  the indices of the output(s) and input(s) to be used in
  64. ##      the frequency response; see @code{sysprune}.
  65. ## @end table
  66. ## @strong{Outputs}
  67. ## @table @var
  68. ## @item  mag, phase
  69. ##  the magnitude and phase of the frequency response
  70. ##        @math{G(jw)} or @math{G(@code{exp}(jwT))} at the selected frequency values.
  71. ## @item w
  72. ## the vector of frequency values used
  73. ## @end table
  74. ## 
  75. ## @strong{Notes}
  76. ## @enumerate
  77. ## @item If no output arguments are given, e.g.,
  78. ## @example
  79. ## bode(sys);
  80. ## @end example
  81. ## bode plots the results to the 
  82. ## screen.  Descriptive labels are automatically placed. 
  83. ## 
  84. ## Failure to include a concluding semicoln will yield some garbage
  85. ## being printed to the screen (@code{ans = []}).
  86. ## 
  87. ## @item If the requested plot is for an MIMO system, mag is set to
  88. ##  @math{||G(jw)||} or @math{||G(@code{exp}(jwT))||}
  89. ## and phase information is not computed.
  90. ## @end enumerate
  91. ## @end deftypefn 
  92.  
  93. function [mag_r, phase_r, w_r] = bode (sys, w, outputs, inputs, plot_style)
  94.  
  95.   ## Written by John Ingram  July 10th, 1996
  96.   ## Based on previous code
  97.   ## By R. Bruce Tenison, July 13, 1994
  98.   ## Modified by David Clem November 13, 1994
  99.   ## again by A. S. Hodel July 1995 (smart plot range, etc.)
  100.   ## Modified by Kai P. Mueller September 28, 1997 (multiplt mode)
  101.  
  102.   ## check number of input arguments given
  103.   if (nargin < 1 | nargin > 5)
  104.     usage("[mag,phase,w] = bode(sys[,w,outputs,inputs,plot_style])");
  105.   endif
  106.   if(nargin < 2)
  107.     w = [];
  108.   endif
  109.   if(nargin < 3)
  110.     outputs = [];
  111.   endif
  112.   if(nargin < 4)
  113.     inputs = [];
  114.   endif
  115.   if(nargin < 5)
  116.     plot_style = "dB";
  117.   endif
  118.  
  119.   if (strcmp (plot_style, "dB"))
  120.     do_db_plot = 1;
  121.   elseif (strcmp (plot_style, "mag"))
  122.     do_db_plot = 0;
  123.   else
  124.     error ("bode: invalid value of plot_style specified");
  125.   endif
  126.  
  127.   [f, w] = bodquist(sys,w,outputs,inputs,"bode");
  128.  
  129.   [stname,inname,outname] = sysgetsg(sys);
  130.   systsam = sysgetts(sys);
  131.  
  132.   ## Get the magnitude and phase of f.
  133.   mag = abs(f);
  134.   phase = arg(f)*180.0/pi;
  135.  
  136.   if (nargout < 1),
  137.     ## Plot the information
  138.     if(gnuplot_has_multiplt)
  139.       oneplot();
  140.     endif
  141.     gset autoscale;
  142.     if(gnuplot_has_multiplt)
  143.       gset nokey;
  144.     endif
  145.     clearplot();
  146.     gset data style lines;
  147.     if(is_digit(sys))
  148.       xlstr = ["Digital frequency w=rad/sec.  pi/T=",num2str(pi/systsam)];
  149.       tistr = "(exp(jwT)) ";
  150.     else
  151.       xlstr = "Frequency in rad/sec";
  152.       tistr = "(jw)";
  153.     endif
  154.     xlabel(xlstr);
  155.     if(is_siso(sys))
  156.       if (gnuplot_has_multiplt)
  157.         subplot(2,1,1);
  158.       endif
  159.       title(["|[Y/U]",tistr,"|, u=", nth(inname,1),", y=",nth(outname,1)]);
  160.     else
  161.       title([ "||Y(", tistr, ")/U(", tistr, ")||"]);
  162.       disp("MIMO plot from")
  163.       disp(outlist(inname,"    "));
  164.       disp("to")
  165.       disp(outlist(outname,"    "));
  166.     endif
  167.     wv = [min(w), max(w)];
  168.     if(do_db_plot && max(mag) > 0)
  169.       ylabel("Gain in dB");
  170.       md = 20*log10(mag);
  171.       axvec = ax2dlim([vec(w),vec(md)]);
  172.       axvec(1:2) = wv;
  173.       axis(axvec);
  174.     else
  175.       ylabel("Gain |Y/U|")
  176.       md = mag;
  177.     endif
  178.  
  179.     grid("on");
  180.     if (do_db_plot)
  181.       semilogx(w,md);
  182.     else
  183.       loglog(w,md);
  184.     endif
  185.     if (is_siso(sys))
  186.       if (gnuplot_has_multiplt)
  187.         subplot(2,1,2);
  188.       else
  189.         prompt('Press any key for phase plot');
  190.       endif
  191.       axvec = ax2dlim([vec(w),vec(phase)]);
  192.       axvec(1:2) = wv;
  193.       axis(axvec);
  194.       xlabel(xlstr);
  195.       ylabel("Phase in deg");
  196.       title([ "phase([Y/U]", tistr, ...
  197.      "), u=", nth(inname,1),", y=",nth(outname,1)]);
  198.       grid("on");
  199.       semilogx(w,phase);
  200.       ## This should be the default for subsequent plot commands.
  201.       if(gnuplot_has_multiplt)
  202.         oneplot();
  203.       endif
  204.     endif
  205.   else
  206.     mag_r = mag;
  207.     phase_r = phase;
  208.     w_r = w;
  209.   endif
  210. endfunction
  211.