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