home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts / control / bodquist.m < prev    next >
Encoding:
Text File  |  1999-12-15  |  5.0 KB  |  157 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{f}, @var{w}] =} bodquist (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx})
  21. ##  used internally by bode, nyquist; compute system frequency response.
  22. ## 
  23. ## @strong{Inputs}
  24. ## @table @var
  25. ## @item sys
  26. ## input system structure
  27. ## @item w
  28. ## range of frequencies; empty if user wants default
  29. ## @item out_idx
  30. ## list of outputs; empty if user wants all
  31. ## @item in_idx
  32. ## list of inputs; empty if user wants all
  33. ## @item rname
  34. ## name of routine that called bodquist ("bode" or "nyquist")
  35. ## @end table
  36. ## @strong{Outputs}
  37. ## @table @var
  38. ## @item w
  39. ##  list of frequencies 
  40. ## @item f
  41. ##  frequency response of sys; @math{f(ii) = f(omega(ii))}
  42. ## @end table
  43. ## @strong{Note} bodquist could easily be incorporated into a Nichols
  44. ## plot function; this is in a "to do" list.
  45. ##
  46. ## Both bode and nyquist share the same introduction, so the common parts are 
  47. ## in bodquist.  It contains the part that finds the number of arguments, 
  48. ## determines whether or not the system is SISO, and computes the frequency 
  49. ## response.  Only the way the response is plotted is different between the 
  50. ## two functions.
  51. ## @end deftypefn
  52.  
  53. function [f, w] = bodquist (sys, w, outputs, inputs, rname)
  54.  
  55.   ## check number of input arguments given
  56.   if (nargin != 5)
  57.     usage("[f,w] = bodquist(sys,w,outputs,inputs,rname)");
  58.   endif
  59.  
  60.   ## check each argument to see if it's in the correct form
  61.   if (!is_struct(sys))
  62.     error("sys must be a system data structure");
  63.   endif
  64.     
  65.   ## let freqresp determine w if it's not already given
  66.   USEW = freqchkw(w);
  67.  
  68.   ## get initial dimensions (revised below if sysprune is called)
  69.   [nn,nz,mm,pp ] = sysdimensions(sys);
  70.  
  71.   ## check for an output vector and to see whether it`s correct
  72.   if (!isempty(outputs))
  73.     if (isempty(inputs))
  74.       inputs = 1:mm;            # use all inputs
  75.       warning([rname,": outputs specified but not inputs"]);
  76.     endif
  77.     sys = sysprune(sys,outputs,inputs);
  78.     [nn,nz,mm,pp ] = sysdimensions(sys);
  79.   endif
  80.  
  81.   ## for speed in computation, convert local copy of 
  82.   ## SISO state space systems to zero-pole  form
  83.   if( is_siso(sys) & strcmp( sysgettype(sys), "ss") )
  84.     [zer,pol,k,tsam,inname,outname] = sys2zp(sys);
  85.     sys = zp2sys(zer,pol,k,tsam,inname,outname);
  86.   endif
  87.  
  88.   ## get system frequency response
  89.   [f,w] = freqresp(sys,USEW,w);   
  90.  
  91.   phase = arg(f)*180.0/pi;
  92.  
  93.   if(!USEW)
  94.     ## smooth plots
  95.     pcnt = 5;        # max number of refinement steps
  96.     dphase = 5;        # desired max change in phase
  97.     dmag = 0.2;        # desired max change in magnitude
  98.     while(pcnt)
  99.       pd = abs(diff(phase));            # phase variation
  100.       pdbig = vec(find(pd > dphase));
  101.  
  102.       lp = length(f);  lp1 = lp-1;        # relative variation
  103.       fd = abs(diff(f));
  104.       fm = max(abs([f(1:lp1); f(2:lp)]));
  105.       fdbig = vec(find(fd > fm/10));
  106.  
  107.       bigpts = union(fdbig, pdbig);
  108.  
  109.       if(isempty(bigpts) )
  110.         pcnt = 0;
  111.       else
  112.         pcnt = pcnt - 1;
  113.         wnew = [];
  114.         crossover_points = find ( phase(1:lp1).*phase(2:lp) < 0);
  115.         pd(crossover_points) = abs(359.99+dphase - pd(crossover_points));
  116.         np_pts = max(3,ceil(pd/dphase)+2);        # phase points
  117.         nm_pts = max(3,ceil(log(fd./fm)/log(dmag))+2);     # magnitude points
  118.         npts = min(5,max(np_pts, nm_pts));
  119.  
  120.         w1 = log10(w(1:lp1));
  121.         w2 = log10(w(2:lp));
  122.         for ii=bigpts
  123.           if(npts(ii))
  124.             wtmp = logspace(w1(ii),w2(ii),npts(ii));
  125.             wseg(ii,1:(npts(ii)-2)) = wtmp(2:(npts(ii)-1));
  126.           endif
  127.         endfor
  128.         wnew = vec(wseg)'; # make a row vector
  129.         wnew = wnew(find(wnew != 0));
  130.         wnew = sort(wnew);
  131.         wnew = create_set(wnew);
  132.         if(isempty(wnew))   # all small crossovers
  133.           pcnt = 0;
  134.         else
  135.           [fnew,wnew] = freqresp(sys,1,wnew);    # get new freq resp points
  136.           w = [w,wnew];            # combine with old freq resp
  137.           f = [f,fnew];
  138.           [w,idx] = sort(w);        # sort into order
  139.           f = f(idx);
  140.           phase = arg(f)*180.0/pi;
  141.         endif
  142.       endif
  143.     endwhile
  144.   endif
  145.  
  146.   ## ensure unique frequency values
  147.   [w,idx] = sort(w);
  148.   f = f(idx);
  149.  
  150.   w_diff = diff(w);
  151.   w_dup = find(w_diff == 0);
  152.   w_idx = complement(w_dup,1:length(w));
  153.   w = w(w_idx);
  154.   f = f(w_idx);
  155.     
  156. endfunction
  157.