From octave-sources-request at bevo dot che dot wisc dot edu  Sat Jun  5 19:19:51 1999
Subject: updated bode function
From: Marc Compere <compere at mail dot utexas dot edu>
To: octave-sources at bevo dot che dot wisc dot edu
Date: Sat, 05 Jun 1999 19:19:50 -0500

This is a multi-part message in MIME format.
--------------8E2FCEBD26A0BC31287BA690
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

John,

	This is an updated verison of bode() that works for systems whos
phase exceeds -pi.  It's mo'bettah than the last one I mailed.

	This assumes the following settings in octave:

  gnuplot_has_frames=1;
  gnuplot_has_multiplot=1;

  PS1 = ">> "
  PS2 = ""
  beep_on_error = 1.0
  default_eval_print_flag = 0.0
  default_save_format = "mat-binary"
  define_all_return_values = 1.0
  do_fortran_indexing = 1.0
  empty_list_elements_ok = 1.0
  fixed_point_format = 1.0
  implicit_num_to_str_ok = 1.0
  implicit_str_to_num_ok = 1.0
  ok_to_lose_imaginary_part = 1.0
  page_screen_output = 0.0
  prefer_column_vectors = 0.0
  prefer_zero_one_indexing = 1.0
  print_empty_dimensions = 0.0
  treat_neg_dim_as_zero = 1.0
  warn_function_name_clash = 0.0
  whitespace_in_literal_matrix = "traditional"


Regards,
Marc Compere
-- 
_________________________________________________
Marc Compere, The University of Texas at Austin
e-mail: compere at mail dot utexas dot edu
www   : http://nerdlab.me.utexas.edu/~compere
office: ETC 4.134a, (512) 471-7347
_________________________________________________
--------------8E2FCEBD26A0BC31287BA690
Content-Type: text/plain; charset=us-ascii;
 name="bode.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="bode.m"

function bode(num,den,omega_lo,omega_hi,N_points,amplitude_units,phase_units)
%
% This is designed to be a replacement for Matlab's BODE() command.
% It can automatically choose a good frequency sweep step-size given high 
% and low frequency bounds.
% It does not automatically choose frequency bounds.
% Bode plots are used for studying the STEADY_STATE amplitude and phase
% responses of linear systems to a sinusoidal input.
%
% Function argument summary:
%   num,den         - polynomial coefficients of decreasing power in s
%   omega_lo        - lower frequency bound
%   omega_hi        - upper frequency bound
%
% Optional arguments:
%   N_points        - Number of points with which you want num and den 
%                     evaluated along the frequency axis (default is 200)
%                     If you want to specifiy either of the two lower arguments
%                     but still want automatic step-size choice, set
%                     N_points = -1
%   amplitude_units - 0 (default) produces amplitude in dB, 1 produces amplitude
%                     in units specific to the system's transfer function
%   phase_units     - 0 (default) produces phase angles in degrees, 1 produces
%                     phase angles in radians
%
% For example, num and den for a standard second-order system transfer function:
%
%            output                  (omega_n^2)
%     G(s) = ------   = ---------------------------------------
%            input       s^2 + (2*zeta*omega_n)*s + (omega_n^2)
%
% Would be:
%     num = omega_n^2;
%     den = [ 1 (2*zeta*omega_n) (omega_n^2) ];
% where omega_n and zeta are variables already defined in the octave workspace.
%
% Examples:
% (the first 2 lines produce identical results)
%	bode(num,den,0.01,100)          <-- omega_lo=0.01, omega_hi=100 (rad/s)
%	bode(num,den,0.01,100,-1,0,0,0)
%	bode(num,den,0.01,100,100)      <-- same as above, but use 100 points
%	bode(num,den,0.01,100,-1,1)     <-- same as above, but use units native to the transfer function in amplitude
%	bode(num,den,0.01,100,-1,0,1)   <-- use radians in phase
%
%
% Marc Compere
% compere at mail dot utexas dot edu
% 5 June 1999

clear j % j must be used as the imaginary number sqrt(-1)
clear delta_omega omega num_ den_ num_amplitude den_amplitude amplitude phase

% If N_points is not supplied, use 200
if (nargin<5)      N_points=-1;        end % if
if (N_points==-1)  N_points=200;       end % if
if (nargin<6),     amplitude_units=0;  end % if
if (nargin<7),     phase_units=0;      end % if

delta_omega = log10(omega_hi/omega_lo)/N_points;

m=max(size(num));
n=max(size(den));

% Do the schmack:
phase(1)=0;
offset=0;
omega(1)=omega_lo;
i=1;
while omega(i) <= omega_hi
 %omega(i+1) = omega(i) + delta_omega;   % linear increase in omega
 omega(i+1) = omega(i)*10^delta_omega;	% exponential increase in omega
 s=j*omega(i+1);     % <-- this is the mapping from s-domain to time-domain

 num_ = 0;
 for k=1:m
  num_ = num(m-k+1)*s^(k-1) + num_;
 end % for

 den_ = 0;
 for k=1:n
  den_ = den(n-k+1)*s^(k-1) + den_;
 end % for

 % The amplitude of a ratio of complex numbers is the amplitude of the
 % numerator divided by the amplitude of the denominator.
 num_amplitude=sqrt(real(num_)^2 + imag(num_)^2);
 den_amplitude=sqrt(real(den_)^2 + imag(den_)^2);
 amplitude(i+1)=num_amplitude/den_amplitude;
 phase(i+1)=atan2(imag(num_),real(num_)) - atan2(imag(den_),real(den_)) + offset;

 % The if-statements below are a fix for the fact that the atan2 function always
 % returns a solution between +pi and -pi.  Many bode plots go way past -pi.
 % Currently the provision for solutions that rise above +pi is untested.
 % It's there, but I'm not sure if it works properly.
 % Basically, this says that if the numerical derivative is large enough,
 % you can assume the atan2() function just switched over to giving a solution
 % in a new quadrant.
 if (phase(i+1)-phase(i))>pi     % this is a jump up, meaning the phase is dropping
  offset = offset - 2*pi;
  phase(i+1)=phase(i+1) - 2*pi;
 end
 if (phase(i+1)-phase(i))<-pi     % this is a jump down, meaning the phase is rising
  offset = offset + 2*pi;
  phase(i+1)=phase(i+1) + 2*pi;
 end


 i=i+1;
end % for

amplitude(1)=amplitude(2);
phase(1)=phase(2);

% Bode  plots are typically presented as magnitude (dB) vs. log(omega)
% and phase (deg) vs. log(omega)
% where the conversion from amplitude to decibels is given as:
%     1dB = 20*log10(amplitude)
%figure(1)
subplot(2,1,1)
 title('Bode Plot') 
 grid("on")
 axis([omega_lo,omega_hi]);
 if (amplitude_units==0),
  ylabel('Magnitude (dB)')
  semilogx(omega,20*log10(amplitude))    % remember: 1dB = 20*log10(amplitude)
  replot
 else,
  ylabel('Magnitude (system units)')
  semilogx(omega,amplitude)
  replot
 end % if

subplot(2,1,2)
 title('') 
 xlabel('Frequency (rad/s)')
 grid("on")
 axis([omega_lo,omega_hi]);
 if (phase_units==0),
  ylabel('Phase (deg)')
  semilogx(omega,phase*180/pi)
  replot
 else,
  ylabel('Phase (rad)')
  semilogx(omega,phase)
  replot
 end % if



% Mathematical development of the omega stepsize:
%	Assume the variable delta_omega is constant.
%	If you want to generate equally spaced points along a linearly
%	increasing axis, simply create a vector with a constant increase
%	like this:
%		omega(i+1) = omega(i) + delta_omega
%	However, if you're plotting log10(omega) then you are not creating
%	an equally spaced vector with a constant increase for each subsequent
%	point.  You need an exponentially increasing step 
%	size to end up with an equally spaced omega vector on a log10(omega)
%	plot like this:
%		omega(i+1) = omega(i)*10^delta_omega
%
%
%           w0      w1      w2      w3     ...      wn
%	...-|-------|-------|-------|-------|-------|----> log_a(omega)
%	     <-d_1-> <-d_2-> <-d_3->   ...   <-d_n->
%
% Q: How do you create equally spaced points along a log_base_a
%    plot of omega?
% A: Find out the relationship between w0, w1, w2, w3, ..., and wp
%    such that all the a's are equal.  Namely, d_1 = d_2 = d_3 = ... = d_n.
%     
%	log_a(w1) - log_a(w0) = d_1
%	log_a(w2) - log_a(w1) = d_2
%	log_a(w3) - log_a(w2) = d_3
%		...
%    So in general,
%
%	log_a(w(n)) - log_a(w(n-1)) = d_n
%    Or...
%	log_a( w(n) / w(n-1) ) = d_n	------------- eqn. 1
%
%    In general the logarithm of x to the base a is the exponent which we
%    raise a to get x.  Mathematically,
%    	log_a( a^x ) = x, and a^(log_a(x)) = x for x > 0
%
%    So, by raising both sides of eqn.1 by a you get:
%
%	w(n) / w(n-1) = a^(d_n)
%
%    or,
%
%	w(n) = w(n-1)*a^(d_n)		------------- eqn.2
%
%    Eqn.2 is the function used to get the next value of omega such that
%    the complete omega vector is equally spaced when plotted on a
%    log(omega) axis.
%    Eqn.2 is implemented above on line 69 as:
%	omega(i+1) = omega(i)*10^delta_omega;
%
%    Because the logarithm used in Bode plots is log_base_10, a=10.
%
% Performance comparison:
%  For the octave commands:
%	num=1;
%	den=[1.00000  0.20000  1];
%	bode(num,den,.1,10,.01)
%  A linearly increasing omega vector generated 992 points whereas
%  an exponentially increasing omega vector only generated 202 points.
%  The plot using the exponential stepsize looked identical to the plot
%  made with the constant stepsize but was faster and required roughly
%  1/5'th the number of points to make.
%

--------------8E2FCEBD26A0BC31287BA690--