home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / FREQINT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.6 KB  |  58 lines
  1. function w=freqint(a,b,c,d,npts)
  2. %FREQINT Auto-ranging algorithm for Bode frequency response
  3. %    
  4. %    W=FREQINT(A,B,C,D,Npts)
  5. %    W=FREQINT(NUM,DEN,Npts)
  6.  
  7. %    Andy Grace  7-6-90
  8. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  9.  
  10. % Generate more points where graph is changing rapidly.
  11. % Calculate points based on eigenvalues and transmission zeros. 
  12.  
  13. [na,ma] = size(a);
  14.  
  15. if (nargin==3)&(na==1),        % Transfer function form.
  16.   npts=c;
  17.   ep=roots(b);
  18.   tz=roots(a);
  19. else                % State space form
  20.   if nargin==3, npts=c; [a,b,c,d] = tf2ss(a,b); end
  21.   ep=[eig(a)];
  22.   tz=tzero(a,b,c,d);
  23. end
  24.  
  25. if isempty(ep), ep=-1000; end
  26.  
  27. % Note: this algorithm does not handle zeros greater than 1e5
  28. ez=[ep(find(imag(ep)>=0));tz(find(abs(tz)<1e5&imag(tz)>=0))];
  29.  
  30. % Round first and last frequencies to nearest decade
  31. integ = abs(ez)<1e-10; % Cater for systems with pure integrators
  32. highfreq=round(log10(max(3*abs(real(ez)+integ)+1.5*imag(ez)))+0.5);
  33. lowfreq=round(log10(0.1*min(abs(real(ez+integ))+2*imag(ez)))-0.5);
  34.  
  35. % Define a base range of frequencies
  36. diffzp=length(ep)-length(tz);
  37. w=logspace(lowfreq,highfreq,npts+diffzp+10*(sum(abs(imag(tz))<abs(real(tz)))>0));
  38. ez=ez(find(imag(ez)>abs(real(ez))));
  39.  
  40. % Oscillatory poles and zeros
  41. if ~isempty(ez)
  42.   f=w;
  43.   npts2=2+8/ceil((diffzp+eps)/10);
  44.   [dum,ind]=sort(-abs(real(ez)));
  45.   z=[];
  46.   for i=ind'
  47.     r1=max([0.8*imag(ez(i))-3*abs(real(ez(i))),10^lowfreq]);
  48.     r2=1.2*imag(ez(i))+4*abs(real(ez(i)));
  49.     z=z(find(z>r2|z<r1));
  50.     indr=find(w<=r2&w>=r1);
  51.     f=f(find(f>r2|f<r1));
  52.     z=[z,logspace(log10(r1),log10(r2),sum(w<=r2&w>=r1)+npts2)];
  53.   end
  54.   w=sort([f,z]);
  55. end
  56. w = w(:);
  57.  
  58.