home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / control / freqresp.m < prev    next >
Text File  |  1999-03-05  |  4KB  |  118 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 [ff,w] = freqresp(sys,USEW,w);
  20.   # function [ff,w] = freqresp(sys,USEW{,w});
  21.   # Frequency response function - used internally by bode, nyquist.
  22.   # minimal argument checking; "do not attempt to do this at home"
  23.   # USEW returned by freqchkw 
  24.   # w: optional, must be present if USEW is given
  25.   #
  26.   # returns: ff = vector of finite G(j*w) entries (or || G(j*w) || for MIMO)
  27.   #          w = vector of frequencies used
  28.   #      ff and w are both returned as row vectors
  29.  
  30.   #  Written by: R. Bruce Tenison July 11, 1994
  31.   # SYS_INTERNAL accesses members of system data structure
  32.  
  33.   save_val = empty_list_elements_ok;
  34.   empty_list_elements_ok = 1;
  35.  
  36.   # Check Args
  37.   if( (nargin < 2) || (nargin > 4) )
  38.     usage ("[ff,w] = freqresp(sys,USEW{,w})");
  39.   elseif( USEW & (nargin < 3) )
  40.     error("USEW=1 but w was not passed.");
  41.   elseif( USEW & isempty(w))
  42.     warning("USEW=1 but w is empty; setting USEW=0");
  43.     USEW=0;
  44.   endif
  45.  
  46.   DIGITAL = is_digital(sys);
  47.  
  48.   # compute default w if needed
  49.   if(!USEW)
  50.     if(is_siso(sys))
  51.       sys = sysupdate(sys,"zp");
  52.       [zer,pol] = sys2zp(sys);
  53.     else
  54.       zer = tzero(sys);
  55.       pol = eig(sys2ss(sys));
  56.     endif
  57.  
  58.     # get default frequency range
  59.     [wmin,wmax] = bode_bounds(zer,pol,DIGITAL,sysgettsam(sys));
  60.     w = logspace(wmin,wmax,50);
  61.   else
  62.     w = reshape(w,1,length(w));     # make sure it's a row vector
  63.   endif
  64.  
  65.   # now get complex values of s or z
  66.   if(DIGITAL)
  67.     jw = exp(i*w*sysgettsam(sys));
  68.   else
  69.     jw = i*w;
  70.   endif
  71.  
  72.   [nn,nz,mm,pp] = sysdimensions(sys);
  73.  
  74.   # now compute the frequency response - divide by zero yields a warning
  75.   if (strcmp(sysgettype(sys),"zp"))
  76.     # zero-pole form (preferred)
  77.     [zer,pol,sysk] = sys2zp(sys);
  78.     ff = ones(size(jw));
  79.     l1 = min(length(zer)*(1-isempty(zer)),length(pol)*(1-isempty(pol)));
  80.     for ii=1:l1
  81.       ff = ff .* (jw - zer(ii)) ./ (jw - pol(ii));
  82.     endfor
  83.  
  84.     # require proper  transfer function, so now just get poles.
  85.     for ii=(l1+1):length(pol)
  86.       ff = ff ./ (jw - pol(ii));
  87.     endfor
  88.     ff = ff*sysk;
  89.  
  90.   elseif (strcmp(sysgettype(sys),"tf"))
  91.     # transfer function form 
  92.     [num,den] = sys2tf(sys);
  93.     ff = polyval(num,jw)./polyval(den,jw);
  94.   elseif (mm==pp)
  95.     # The system is square; do state-space form bode plot
  96.     [sysa,sysb,sysc,sysd,tsam,sysn,sysnz] = sys2ss(sys);
  97.     n = sysn + sysnz;
  98.     for ii=1:length(jw);
  99.       ff(ii) = det(sysc*((jw(ii).*eye(n)-sysa)\sysb)+sysd);
  100.     endfor;
  101.   else
  102.     # Must be state space... bode                            
  103.     [sysa,sysb,sysc,sysd,tsam,sysn,sysnz] = sys2ss(sys);
  104.     n = sysn + sysnz;
  105.     for ii=1:length(jw);
  106.       ff(ii) = norm(sysc*((jw(ii)*eye(n)-sysa)\sysb)+sysd);
  107.     endfor
  108.     
  109.   endif
  110.  
  111.   w = reshape(w,1,length(w));
  112.   ff = reshape(ff,1,length(ff));
  113.  
  114.   # restore global variable
  115.   empty_list_elements_ok = save_val;
  116. endfunction
  117.  
  118.