/trunk/octave-forge/extra/control-legacy/inst/nyquist.m
MATLAB | 210 lines | 204 code | 6 blank | 0 comment | 18 complexity | ca8c74a6a522b750416b871a328378c3 MD5 | raw file
Possible License(s): GPL-2.0, BSD-3-Clause, LGPL-2.1, GPL-3.0, LGPL-3.0
- ## Copyright (C) 1996, 1998, 2000, 2003, 2004, 2005, 2006, 2007
- ## Auburn University. All rights reserved.
- ##
- ##
- ## This program is free software; you can redistribute it and/or modify it
- ## under the terms of the GNU General Public License as published by
- ## the Free Software Foundation; either version 3 of the License, or (at
- ## your option) any later version.
- ##
- ## This program is distributed in the hope that it will be useful, but
- ## WITHOUT ANY WARRANTY; without even the implied warranty of
- ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- ## General Public License for more details.
- ##
- ## You should have received a copy of the GNU General Public License
- ## along with this program; see the file COPYING. If not, see
- ## <http://www.gnu.org/licenses/>.
- ## -*- texinfo -*-
- ## @deftypefn {Function File} {[@var{realp}, @var{imagp}, @var{w}] =} nyquist (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx}, @var{atol})
- ## @deftypefnx {Function File} {} nyquist (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx}, @var{atol})
- ## Produce Nyquist plots of a system; if no output arguments are given, Nyquist
- ## plot is printed to the screen.
- ##
- ## Compute the frequency response of a system.
- ##
- ## @strong{Inputs} (pass as empty to get default values)
- ## @table @var
- ## @item sys
- ## system data structure (must be either purely continuous or discrete;
- ## see @code{is_digital})
- ## @item w
- ## frequency values for evaluation.
- ## If sys is continuous, then bode evaluates @math{G(@var{jw})};
- ## if sys is discrete, then bode evaluates @math{G(exp(@var{jwT}))},
- ## where @var{T} is the system sampling time.
- ## @item default
- ## the default frequency range is selected as follows: (These
- ## steps are @strong{not} performed if @var{w} is specified)
- ## @enumerate
- ## @item via routine @command{__bodquist__}, isolate all poles and zeros away from
- ## @var{w}=0 (@var{jw}=0 or @math{exp(@var{jwT})=1}) and select the frequency
- ## range based on the breakpoint locations of the frequencies.
- ## @item if @var{sys} is discrete time, the frequency range is limited
- ## to @var{jwT} in
- ## @ifinfo
- ## [0,2p*pi]
- ## @end ifinfo
- ## @iftex
- ## @tex
- ## $ [ 0,2 p \pi ] $
- ## @end tex
- ## @end iftex
- ## @item A ``smoothing'' routine is used to ensure that the plot phase does
- ## not change excessively from point to point and that singular
- ## points (e.g., crossovers from +/- 180) are accurately shown.
- ## @end enumerate
- ## @item atol
- ## for interactive nyquist plots: atol is a change-in-slope tolerance
- ## for the of asymptotes (default = 0; 1e-2 is a good choice). This allows
- ## the user to ``zoom in'' on portions of the Nyquist plot too small to be
- ## seen with large asymptotes.
- ## @end table
- ## @strong{Outputs}
- ## @table @var
- ## @item realp
- ## @itemx imagp
- ## the real and imaginary parts of the frequency response
- ## @math{G(jw)} or @math{G(exp(jwT))} at the selected frequency values.
- ## @item w
- ## the vector of frequency values used
- ## @end table
- ##
- ## If no output arguments are given, nyquist plots the results to the screen.
- ## If @var{atol} != 0 and asymptotes are detected then the user is asked
- ## interactively if they wish to zoom in (remove asymptotes)
- ## Descriptive labels are automatically placed.
- ##
- ## Note: if the requested plot is for an @acronym{MIMO} system, a warning message is
- ## presented; the returned information is of the magnitude
- ## @iftex
- ## @tex
- ## $ \Vert G(jw) \Vert $ or $ \Vert G( {\rm exp}(jwT) \Vert $
- ## @end tex
- ## @end iftex
- ## @ifinfo
- ## ||G(jw)|| or ||G(exp(jwT))||
- ## @end ifinfo
- ## only; phase information is not computed.
- ## @end deftypefn
- ## Author: R. Bruce Tenison <btenison@eng.auburn.edu>
- ## Created: July 13, 1994
- ## A. S. Hodel July 1995 (adaptive frequency spacing,
- ## remove acura parameter, etc.)
- ## Revised by John Ingram July 1996 for system format
- function [realp, imagp, w] = nyquist (sys, w, outputs, inputs, atol)
- ## Both bode and nyquist share the same introduction, so the common
- ## parts are in a file called __bodquist__.m. It contains the part that
- ## finds the number of arguments, determines whether or not the system
- ## is SISO, andd computes the frequency response. Only the way the
- ## response is plotted is different between the two functions.
- ## check number of input arguments given
- if (nargin < 1 || nargin > 5)
- print_usage ();
- endif
- if (nargin < 2)
- w = [];
- endif
- if (nargin < 3)
- outputs = [];
- endif
- if (nargin < 4)
- inputs = [];
- endif
- if (nargin < 5)
- atol = 0;
- elseif (! (is_sample (atol) || atol == 0))
- error ("nyquist: atol must be a nonnegative scalar")
- endif
- ## signal to __bodquist__ who's calling
- [f, w, sys] = __bodquist__ (sys, w, outputs, inputs, "nyquist");
- ## Get the real and imaginary part of f.
- realp = real (f);
- imagp = imag (f);
- ## No output arguments, then display plot, otherwise return data.
- if (nargout == 0)
- dnplot = 0;
- while (! dnplot)
- plot (realp, imagp, "- ;+w;", realp, -imagp, "-@ ;-w;");
- grid ("on");
- if (is_digital (sys))
- tstr = " G(e^{jw}) ";
- else
- tstr = " G(jw) ";
- endif
- xlabel (sprintf ("Re(%s)", tstr));
- ylabel (sprintf ("Im(%s)", tstr));
- [stn, inn, outn] = sysgetsignals (sys);
- if (is_siso (sys))
- title (sprintf ("Nyquist plot from %s to %s, w (rad/s) in [%e, %e]",
- inn{1}, outn{1}, w(1), w(end)));
- endif
- axis (__axis2dlim__ ([[realp(:), imagp(:)]; [realp(:), -imagp(:)]]));
- ## check for interactive plots
- dnplot = 1; # assume done; will change later if atol is satisfied
- if (atol > 0 && length (f) > 2)
- ## check for asymptotes
- fmax = max (abs (f));
- fi = find (abs (f) == fmax, 1, "last");
- ## compute angles from point to point
- df = diff (f);
- th = atan2 (real (df), imag (df)) * 180 / pi;
- ## get angle at fmax
- if (fi == length (f))
- fi = fi-1;
- endif
- thm = th(fi);
- ## now locate consecutive angles within atol of thm
- ith_same = find (abs (th - thm) < atol);
- ichk = union (fi, find (diff (ith_same) == 1));
- ## locate max, min consecutive indices in ichk
- loval = max (complement (ichk, 1:fi));
- if (isempty (loval))
- loval = fi;
- else
- loval = loval + 1;
- endif
- hival = min (complement (ichk, fi:length(th)));
- if (isempty (hival))
- hival = fi+1;
- endif
- keep_idx = complement (loval:hival, 1:length(w));
- if (length (keep_idx))
- resp = input ("Remove asymptotes and zoom in (y or n): ", 1);
- if (resp(1) == "y")
- dnplot = 0; # plot again
- w = w(keep_idx);
- f = f(keep_idx);
- realp = real (f);
- imagp = imag (f);
- endif
- endif
- endif
- endwhile
- w = realp = imagp = [];
- endif
- endfunction