home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / control / nyquist.m < prev    next >
Text File  |  1999-03-05  |  7KB  |  187 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 [realp,imagp,w] = nyquist(sys,w,outputs,inputs,atol)
  20. # [realp,imagp,w] = nyquist(sys[,w,outputs,inputs,atol])
  21. # Produce Nyquist plots of a system
  22. #
  23. # Compute the frequency response of a system.
  24. # inputs: (pass as empty to get default values
  25. #   sys: system data structure (must be either purely continuous or discrete;
  26. #        see is_digital)
  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 
  30. #         T=sysgettsam(sys) (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. #   atol: for interactive nyquist plots: atol is a change-in-angle tolerance 
  44. #     (in degrees) for the of asymptotes (default = 0; 1e-2 is a good choice).
  45. #     Consecutive points along the asymptotes whose angle is within atol of
  46. #     the angle between the largest two points are omitted for "zooming in"
  47. #
  48. # outputs:
  49. #    realp, imagp: the real and imaginary parts of the frequency response
  50. #       G(jw) or G(exp(jwT)) at the selected frequency values.
  51. #    w: the vector of frequency values used
  52. #
  53. # If no output arguments are given, nyquist plots the results to the screen.
  54. # If atol != 0 and asymptotes are detected then the user is asked 
  55. #    interactively if they wish to zoom in (remove asymptotes)
  56. # Descriptive labels are automatically placed.  See xlabel, ylable, title,
  57. # and replot.
  58. #
  59. # Note: if the requested plot is for an MIMO system, a warning message is
  60. # presented; the returned information is of the magnitude 
  61. # ||G(jw)|| or ||G(exp(jwT))|| only; phase information is not computed.
  62.    
  63.   # By R. Bruce Tenison, July 13, 1994
  64.   # A. S. Hodel July 1995 (adaptive frequency spacing, 
  65.   #     remove acura parameter, etc.)
  66.   # Revised by John Ingram July 1996 for system format
  67.   #
  68.  
  69.   # Both bode and nyquist share the same introduction, so the common parts are 
  70.   # in a file called bodquist.m.  It contains the part that finds the 
  71.   # number of arguments, determines whether or not the system is SISO, and 
  72.   # computes the frequency response.  Only the way the response is plotted is
  73.   # different between the two functions.
  74.  
  75.   save_val = implicit_str_to_num_ok;    # save for later
  76.   implicit_str_to_num_ok = 1;
  77.  
  78.   # check number of input arguments given
  79.   if (nargin < 1 | nargin > 5)
  80.     usage("[realp,imagp,w] = nyquist(sys[,w,outputs,inputs,atol])");
  81.   endif
  82.   if(nargin < 2)
  83.     w = [];
  84.   endif
  85.   if(nargin < 3)
  86.     outputs = [];
  87.   endif
  88.   if(nargin < 4)
  89.     inputs = [];
  90.   endif
  91.   if(nargin < 5)
  92.     atol = 0;
  93.   elseif(!(is_sample(atol) | atol == 0))
  94.     error("atol must be a nonnegative scalar.")
  95.   endif
  96.  
  97.   # signal to bodquist who's calling
  98.    
  99.   [f,w] = bodquist(sys,w,outputs,inputs,"nyquist");
  100.  
  101.   # Get the real and imaginary part of f.
  102.   realp = real(f);
  103.   imagp = imag(f);
  104.  
  105.   # No output arguments, then display plot, otherwise return data.
  106.   if (nargout == 0)
  107.     dnplot = 0;
  108.     while(!dnplot)
  109.       if(gnuplot_has_multiplot)
  110.         oneplot();
  111.         gset key;
  112.       endif
  113.       clearplot();
  114.       grid ("on");
  115.       gset data style lines;
  116.   
  117.       if(is_digital(sys))
  118.         tstr = " G(e^{jw}) ";
  119.       else
  120.         tstr = " G(jw) ";
  121.       endif
  122.       xlabel(["Re(",tstr,")"]);
  123.       ylabel(["Im(",tstr,")"]);
  124.   
  125.       [stn, inn, outn] = sysgetsignals(sys);
  126.       if(is_siso(sys))
  127.         title(sprintf("Nyquist plot from %s to %s, w (rad/s) in [%e, %e]", ...
  128.       nth(inn,1), nth(outn,1), w(1), w(length(w))) )
  129.       endif
  130.   
  131.       gset nologscale xy;
  132.  
  133.       axis(axis2dlim([[vec(realp),vec(imagp)];[vec(realp),-vec(imagp)]]));
  134.       plot(realp,imagp,"- ;+w;",realp,-imagp,"-@ ;-w;");
  135.  
  136.       # check for interactive plots
  137.       dnplot = 1; # assume done; will change later if atol is satisfied
  138.       if(atol > 0 & length(f) > 2)
  139.  
  140.         # check for asymptotes
  141.         fmax = max(abs(f));
  142.         fi = max(find(abs(f) == fmax));
  143.         
  144.         # compute angles from point to point
  145.         df = diff(f);
  146.         th = atan2(real(df),imag(df))*180/pi;
  147.  
  148.         # get angle at fmax
  149.         if(fi == length(f)) fi = fi-1; endif
  150.         thm = th(fi);
  151.     
  152.         # now locate consecutive angles within atol of thm
  153.         ith_same = find(abs(th - thm) < atol);
  154.         ichk = union(fi,find(diff(ith_same) == 1));
  155.  
  156.         #locate max, min consecutive indices in ichk
  157.         loval = max(complement(ichk,1:fi));
  158.         if(isempty(loval)) loval = fi;
  159.         else               loval = loval + 1;   endif
  160.  
  161.         hival = min(complement(ichk,fi:length(th)));
  162.         if(isempty(hival))  hival = fi+1;      endif
  163.  
  164.         keep_idx = complement(loval:hival,1:length(w));
  165.  
  166.         if(length(keep_idx))
  167.           resp = input("Remove asymptotes and zoom in (y or n): ",1);
  168.           if(resp(1) == "y")
  169.             dnplot = 0;                 # plot again
  170.             w = w(keep_idx);
  171.             f = f(keep_idx);
  172.             realp = real(f);
  173.             imagp = imag(f);
  174.           endif
  175.         endif
  176.  
  177.      endif
  178.     endwhile
  179.     w = [];
  180.     realp=[];
  181.     imagp=[];
  182.   endif
  183.  
  184.   implicit_str_to_num_ok = save_val;    # restore value
  185.  
  186. endfunction
  187.