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