home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts / control / nichols.m < prev    next >
Encoding:
Text File  |  1999-12-15  |  4.0 KB  |  118 lines
  1. ## Copyright (C) 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. ## [mag,phase,w] = nichols(sys[,w,outputs,inputs])
  20. ## Produce Nichols plot of a system
  21. ##
  22. ## Compute the frequency response of a system.
  23. ## inputs:
  24. ##   sys: system data structure (must be either purely continuous or discrete;
  25. ##     see is_digital)
  26. ##   w: frequency values for evaluation.
  27. ##      if sys is continuous, then nichols evaluates G(jw)
  28. ##      if sys is discrete, then nichols evaluates G(exp(jwT)), where T=sys.tsam
  29. ##         (the system sampling time)
  30. ##      default: the default frequency range is selected as follows: (These
  31. ##        steps are NOT performed if w is specified)
  32. ##          (1) via routine bodquist, isolate all poles and zeros away from
  33. ##              w=0 (jw=0 or exp(jwT)=1) and select the frequency
  34. ##             range based on the breakpoint locations of the frequencies.
  35. ##          (2) if sys is discrete time, the frequency range is limited
  36. ##              to jwT in [0,2p*pi]
  37. ##          (3) A "smoothing" routine is used to ensure that the plot phase does
  38. ##              not change excessively from point to point and that singular
  39. ##              points (e.g., crossovers from +/- 180) are accurately shown.
  40. ##   outputs, inputs: the indices of the output(s) and input(s) to be used in
  41. ##     the frequency response; see sysprune.
  42. ## outputs:
  43. ##    mag, phase: the magnitude and phase of the frequency response
  44. ##       G(jw) or G(exp(jwT)) at the selected frequency values.
  45. ##    w: the vector of frequency values used
  46. ## If no output arguments are given, nichols plots the results to the screen.
  47. ## Descriptive labels are automatically placed.  See xlabel, ylable, title,
  48. ## and replot.
  49. ##
  50. ## Note: if the requested plot is for an MIMO system, mag is set to
  51. ## ||G(jw)|| or ||G(exp(jwT))|| and phase information is not computed.
  52.  
  53. function [mag, phase, w] = nichols (sys, w, outputs, inputs)
  54.  
  55.   ## check number of input arguments given
  56.   if (nargin < 1 | nargin > 4)
  57.     usage("[mag,phase,w] = nichols(sys[,w,outputs,inputs])");
  58.   endif
  59.   if(nargin < 2)
  60.     w = [];
  61.   endif
  62.   if(nargin < 3)
  63.     outputs = [];
  64.   endif
  65.   if(nargin < 4)
  66.     inputs = [];
  67.   endif
  68.  
  69.   [f, w] = bodquist(sys,w,outputs,inputs,"nichols");
  70.  
  71.   [stname,inname,outname] = sysgetsignals(sys);
  72.   systsam = sysgettsam(sys);
  73.  
  74.   ## Get the magnitude and phase of f.
  75.   mag = abs(f);
  76.   phase = arg(f)*180.0/pi;
  77.  
  78.   if (nargout < 1),
  79.     ## Plot the information
  80.     if(gnuplot_has_multiplot)
  81.       oneplot();
  82.     endif
  83.     gset autoscale;
  84.     if(gnuplot_has_multiplot)
  85.       gset nokey;
  86.     endif
  87.     clearplot();
  88.     grid("on");
  89.     gset data style lines;
  90.     if(is_digital(sys))
  91.       tistr = "(exp(jwT)) ";
  92.     else
  93.       tistr = "(jw)";
  94.     endif
  95.     xlabel("Phase (deg)");
  96.     if(is_siso(sys))
  97.       title(["Nichols plot of |[Y/U]",tistr,"|, u=", ...
  98.     sysgetsignals(sys,"in",1,1), ", y=",sysgetsignals(sys,"out",1,1)]);
  99.     else
  100.       title([ "||Y(", tistr, ")/U(", tistr, ")||"]);
  101.       printf("MIMO plot from\n%s\nto\n%s\n",outlist(inname,"    "), ...
  102.         outlist(outname,"    "));
  103.     endif
  104.     if(max(mag) > 0)
  105.       ylabel("Gain in dB");
  106.       md = 20*log10(mag);
  107.     else
  108.       ylabel("Gain |Y/U|")
  109.       md = mag;
  110.     endif
  111.  
  112.     axvec = axis2dlim([vec(phase),vec(md)]);
  113.     axis(axvec);
  114.     plot(phase,md);
  115.     mag = phase = w = [];
  116.   endif
  117. endfunction
  118.