home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / nyquist.m < prev    next >
Text File  |  1999-12-24  |  7KB  |  207 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{realp}, @var{imagp}, @var{w}] =} nyquist (@var{sys}@{, @var{w}, @var{out_idx}, @var{in_idx}, @var{atol}@})
  21. ## @deftypefnx {Function File } {} nyquist (@var{sys}@{, @var{w}, @var{out_idx}, @var{in_idx}, @var{atol}@})
  22. ## Produce Nyquist plots of a system; if no output arguments are given, Nyquist
  23. ## plot is printed to the screen.
  24. ## 
  25. ## Compute the frequency response of a system.
  26. ## @strong{Inputs} (pass as empty to get default values)
  27. ## @table @var
  28. ## @item sys
  29. ##        system data structure (must be either purely continuous or discrete;
  30. ##        see is_digit)
  31. ## @item w 
  32. ##      frequency values for evaluation.
  33. ##      if sys is continuous, then bode evaluates @math{G(jw)}
  34. ##      if sys is discrete, then bode evaluates @math{G(exp(jwT))}, where 
  35. ##         @math{@var{T}=sysgetts(@var{sys})} (the system sampling time)
  36. ## @item default
  37. ##      the default frequency range is selected as follows: (These
  38. ##      steps are NOT performed if @var{w} is specified)
  39. ## @end table
  40. ## @enumerate 
  41. ## @item via routine bodquist, isolate all poles and zeros away from 
  42. ## @var{w}=0 (@var{jw}=0 or @math{exp(@var{jwT})=1}) and select the frequency
  43. ## range based on the breakpoint locations of the frequencies.
  44. ## @item if @var{sys} is discrete time, the frequency range is limited
  45. ## to @var{jwT} in 
  46. ## @ifinfo
  47. ## [0,2p*pi]
  48. ## @end ifinfo
  49. ## @iftex
  50. ## $[0,2p*\pi]$
  51. ## @end iftex
  52. ## @item A "smoothing" routine is used to ensure that the plot phase does
  53. ##              not change excessively from point to point and that singular
  54. ##              points (e.g., crossovers from +/- 180) are accurately shown.
  55. ## @end enumerate
  56. ##   outputs, inputs: the indices of the output(s) and input(s) to be used in
  57. ##     the frequency response; see sysprune.
  58. ## 
  59. ## @strong{Inputs} (pass as empty to get default values)
  60. ## @table @var
  61. ## @item   atol
  62. ## for interactive nyquist plots: atol is a change-in-slope tolerance 
  63. ## for the of asymptotes (default = 0; 1e-2 is a good choice).  This allows
  64. ## the user to ``zoom in'' on portions of the Nyquist plot too small to be
  65. ## seen with large asymptotes.
  66. ## @end table
  67. ## @strong{Outputs}
  68. ## @table @var
  69. ## @item    realp, imagp
  70. ## the real and imaginary parts of the frequency response
  71. ##        @math{G(jw)} or @math{G(exp(jwT))} at the selected frequency values.
  72. ## @item    w
  73. ##  the vector of frequency values used
  74. ## @end table
  75. ## 
  76. ##  If no output arguments are given, nyquist plots the results to the screen.
  77. ##  If @var{atol} != 0 and asymptotes are detected then the user is asked 
  78. ##     interactively if they wish to zoom in (remove asymptotes)
  79. ##  Descriptive labels are automatically placed.
  80. ## 
  81. ##  Note: if the requested plot is for an MIMO system, a warning message is
  82. ##  presented; the returned information is of the magnitude 
  83. ##  ||G(jw)|| or ||G(exp(jwT))|| only; phase information is not computed.
  84. ## 
  85. ## @end deftypefn
  86.  
  87. function [realp, imagp, w] = nyquist (sys, w, outputs, inputs, atol)
  88.    
  89.   ## By R. Bruce Tenison, July 13, 1994
  90.   ## A. S. Hodel July 1995 (adaptive frequency spacing, 
  91.   ##     remove acura parameter, etc.)
  92.   ## Revised by John Ingram July 1996 for system format
  93.  
  94.   ## Both bode and nyquist share the same introduction, so the common
  95.   ## parts are in a file called bodquist.m.  It contains the part that
  96.   ## finds the number of arguments, determines whether or not the system
  97.   ## is SISO, andd computes the frequency response.  Only the way the
  98.   ## response is plotted is different between the two functions.
  99.  
  100.   ## check number of input arguments given
  101.   if (nargin < 1 | nargin > 5)
  102.     usage("[realp,imagp,w] = nyquist(sys[,w,outputs,inputs,atol])");
  103.   endif
  104.   if(nargin < 2)
  105.     w = [];
  106.   endif
  107.   if(nargin < 3)
  108.     outputs = [];
  109.   endif
  110.   if(nargin < 4)
  111.     inputs = [];
  112.   endif
  113.   if(nargin < 5)
  114.     atol = 0;
  115.   elseif(!(is_sampl(atol) | atol == 0))
  116.     error("atol must be a nonnegative scalar.")
  117.   endif
  118.  
  119.   ## signal to bodquist who's calling
  120.    
  121.   [f,w] = bodquist(sys,w,outputs,inputs,"nyquist");
  122.  
  123.   ## Get the real and imaginary part of f.
  124.   realp = real(f);
  125.   imagp = imag(f);
  126.  
  127.   ## No output arguments, then display plot, otherwise return data.
  128.   if (nargout == 0)
  129.     dnplot = 0;
  130.     while(!dnplot)
  131.       if(gnuplot_has_multiplt)
  132.         oneplot();
  133.         gset key;
  134.       endif
  135.       clearplot();
  136.       grid ("on");
  137.       gset data style lines;
  138.   
  139.       if(is_digit(sys))
  140.         tstr = " G(e^{jw}) ";
  141.       else
  142.         tstr = " G(jw) ";
  143.       endif
  144.       xlabel(["Re(",tstr,")"]);
  145.       ylabel(["Im(",tstr,")"]);
  146.   
  147.       [stn, inn, outn] = sysgetsg(sys);
  148.       if(is_siso(sys))
  149.         title(sprintf("Nyquist plot from %s to %s, w (rad/s) in [%e, %e]", ...
  150.       nth(inn,1), nth(outn,1), w(1), w(length(w))) )
  151.       endif
  152.   
  153.       gset nologscale xy;
  154.  
  155.       axis(ax2dlim([[vec(realp),vec(imagp)];[vec(realp),-vec(imagp)]]));
  156.       plot(realp,imagp,"- ;+w;",realp,-imagp,"-@ ;-w;");
  157.  
  158.       ## check for interactive plots
  159.       dnplot = 1; # assume done; will change later if atol is satisfied
  160.       if(atol > 0 & length(f) > 2)
  161.  
  162.         ## check for asymptotes
  163.         fmax = max(abs(f));
  164.         fi = max(find(abs(f) == fmax));
  165.         
  166.         ## compute angles from point to point
  167.         df = diff(f);
  168.         th = atan2(real(df),imag(df))*180/pi;
  169.  
  170.         ## get angle at fmax
  171.         if(fi == length(f)) fi = fi-1; endif
  172.         thm = th(fi);
  173.     
  174.         ## now locate consecutive angles within atol of thm
  175.         ith_same = find(abs(th - thm) < atol);
  176.         ichk = union(fi,find(diff(ith_same) == 1));
  177.  
  178.         ## locate max, min consecutive indices in ichk
  179.         loval = max(complmnt(ichk,1:fi));
  180.         if(isempty(loval)) loval = fi;
  181.         else               loval = loval + 1;   endif
  182.  
  183.         hival = min(complmnt(ichk,fi:length(th)));
  184.         if(isempty(hival))  hival = fi+1;      endif
  185.  
  186.         keep_idx = complmnt(loval:hival,1:length(w));
  187.  
  188.         if(length(keep_idx))
  189.           resp = input("Remove asymptotes and zoom in (y or n): ",1);
  190.           if(resp(1) == "y")
  191.             dnplot = 0;                 # plot again
  192.             w = w(keep_idx);
  193.             f = f(keep_idx);
  194.             realp = real(f);
  195.             imagp = imag(f);
  196.           endif
  197.         endif
  198.  
  199.      endif
  200.     endwhile
  201.     w = [];
  202.     realp=[];
  203.     imagp=[];
  204.   endif
  205.  
  206. endfunction
  207.