Dr. Boris J. Lurie
HOME  |  RESUME  |  PUBLICATIONS  |  CONSULTING  |  LECTURES  |  SITE MAP 

 
LIST OF PUBLICATIONS

Books Patents Journal Publications NASA Tech Briefs Conference Proceedings Publications
Discussions Readers' Contributions Bode Step Toolbox Function Listings

FUNCTION LISTINGS

NYQLOG BOLAGNYQ TFSHIFT BONYQAS BOSTEP BOCLOS
BOINTEGR BOCOMP NDCP BNDCP NDCB BNDCB

See Appendix 14 of Classical Feedback Control for descriptions of these MATLAB funtions. The functions can be copied freely for nonprofit use. The author bears no responsibility for the results of the functions' applications.

Download BodeStepToolbox.tar
(MATLAB and SPICE scripts for the examples and problems in Classical Feedback Control are available here.)

NYQLOG

function  NYQLOG(wb,numl,denl)

%NYQLOG Nyquist diagram on logarithmic plane
%function  NYQLOG(wb,numl,denl)
%     The loop transfer function is numl(s)/denl(s), where numl and
%     denl contain the polynomial coefficients in descending powers of s. 
%     wb is the crossover frequency in rad/sec, i.e. frequency at which
%     loop gain is 0 dB.
%     Default:  wb = 1, numl = [-26  178  654  1269  429],
%     denl = [1  18.8  132  543  1164  860  0  0]
%
%     Copyright (c) B. J. Lurie 4-1-98, La Crescenta.
%     See also BOSTEP, BOCLOS, BONYQAS, BOINTEGR, TFSHIFT, STEPLOG from the
%     Bode toolbox.
%     Further explanations are given in Classical Feedback Control by
%     B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.

   if nargin < 3,   
      denl = [1  18.8  132  543  1164  860  0  0];
   end 
   if nargin < 2,   
      numl = [-26  178  654  1269  429]; 
   end 
   if nargin == 0,   
      wb = 1; 
   end 

   [mag, phase] = bode(numl, denl);  plot(phase, 20*log10(mag),'w')
   set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90])
   grid
   axis([-270 -90 -20 70])
   hold on
   [a,b] = bode (numl,denl,wb); plot(b,20*log10(a),'wx')
     for w = [0.0156 0.03125 0.0625 0.125 0.25 0.5 2 4],
         [a,b] = bode(numl,denl,wb*w); plot(b,20*log10(a),'w+');
     end
   plot(-180, 0, 'wo');
   hold off
   title('Nyquist diagram, x marks w = wb, + marks octaves')
   xlabel('loop phase shift in degrees')
   ylabel('loop gain in dB')
   zoom on

Back

BOLAGNYQ

function  BOLAGNYQ(wb,numl,denl)
%BOLAGNYQ Gain and -(phase+180) (i.e. lag margin) responses, and Nyquist
%    diagram on logarithmic plane 
%function  BOLAGNYQ(wb,numl,denl)
%    The loop transfer function is numl(s)/denl(s), where numl and
%    denl contain the polynomial coefficients in descending powers of s. 
%    Frequency wb in rad/sec specifies the position of cursor x; it is
%    convenient to put it at the frequency where the loop gain is 0 dB.
%
%    Default/demo:  wb = 1, numl = [-26  178  654  1269  429],
%    denl = [1  18.8  132  543  1164  860  0  0]
%
%    Copyright (c) B. J. Lurie 6-16-98, La Crescenta.
%    See also BOSTEP, BOCLOS, BONYQAS, BOINTEGR, TFSHIFT, STEPLOG, NYQLOG from
%    the Bodestep toolbox. 
%    Further explanations are given in Classical Feedback Control by
%    B. J. Lurie and P. J. Enright, Marcel Dekker, 2000. 

   if nargin < 3,   
      denl = [1  18.8  132  543  1164  860  0  0];
   end 
   if nargin < 2,   
      numl = [-26  178  654  1269  429]; 
   end 
   if nargin == 0,   
      wb = 1; 
   end 

   [mag, phase, w] = bode(numl, denl);

   subplot(2,1,1)
   semilogx(w, 20*log10(mag),'w', w,(-180-phase),'w--');
   hold on
   [mmb,ppb] = bode(numl, denl,wb);
   semilogx(wb,mmb,'wx'); 
   wmin = min(w); wmax = max(w);
   axis([wmin wmax -50 100])
   xlabel('rad/sec')
   ylabel('gain, dB; lag margin, degr')
   grid
   hold off
    
   subplot(2,1,2)
   plot(phase, 20*log10(mag),'w')
   set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90])
   grid
   axis([-270 -90 -20 50])
   hold on
   [a,b] = bode (numl,denl,wb); plot(b,20*log10(a),'wx')
     for w = [0.0156 0.03125 0.0625 0.125 0.25 0.5 2 4],
         [a,b] = bode(numl,denl,wb*w); plot(b,20*log10(a),'w+');
     end
   plot(-180, 0, 'wo');
   hold off
   xlabel('phase, degr; x marks wb, + mark octaves from wb')
   ylabel('loop gain, dB')
   zoom on
   
   wb
   gain_at_wb = mmb
   lag_margin = 180 + ppb

Back

TFSHIFT

function [numsh, densh] = tfshift(num,den,q,p)
%TFSHIFT   Transfer function with shifted frequency response.
%     [numsh, densh] = tfshift(num,den,q,p)
%     TFSHIFT produces numerator and denominator of transfer function 
%     NUMSH(s)/DENSH(s) shifted in frequency q times from NUM(s)/DEN(s) 
%     where NUM and DEN contain the polynomial coefficients in 
%     descending powers of s. NUMSH(s)/DENSH(s) has the same value 
%     at frequency qw as the function NUM(s)/DEN(s) at frequency w. 
%     The numerator and denumerator are both multiplied by 10^p for 
%     the coefficients to have convenient values. Default is p = 0.
%
%     Copyright (c) 1998 A. Ahmed and B. J. Lurie 2-27-97.
%     Further explanations are given in Classical Feedback Control by
%     B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.

   ln = length(num)-1;
   ld = length(den)-1;
   pn = [ln: -1: 0];
   pd = [ld: -1: 0];
   numsh = num./(q.^pn);
   densh = den./(q.^pd);
   if nargin == 4,   
      numsh = numsh*10^p; 
      densh = densh*10^p;
   end 

Back

BONYQAS

function [gain1,phase1] = bonyqas(u,g,lslope,hslope,nml);
%     function [gain1,phase1] = bonyqas(u,g,lslope,hslope,bnm);
% BONYQAS  Asymptotic Bode and Nyquist diagrams corresponding in minimum
%     phase manner to a piece-linear gain response in dB specified by the
%     gain g(h) at k corner frequencies u(h) in rad/sec. The lengths of
%     vectors g and u must be equal. Coefficient lslope defines the
%     asymptotic slope at lower frequencies in integrator units, and
%     hslope defines the asymptotic slope at higher frequencies. nml is
%     nonminimum phase lag in degrees at u(h); the lag is proportional
%     to the frequency, default is nml = 0.
%     Smoothed Bode diagram is plotted in white. Phase is plotted in green.
%     Nyquist diagram is plotted on logarithmic plane.
%     The function is applicable for reasonably smooth responses. The phase
%     accuracy for typical feedback loop shaping problem without sharp
%     peaks or notches is better than 3 degr. 
%       
%     Default: bonyqas without arguments plots the band-pass response:
%     bonyqas([0.2 2 8 32],[-5 10 6 -12],-1,3);
%     Low-pass loop Example 1: bonyqas([3 4 20 30],[16 13 -9 -10],2, 3, 0)
%     Loop Example 2: bonyqas([3 3.8 22 45],[16 13 -10 -11.5],2, 3, 50)
%
%     Copyright (c) B. Lurie 6-13-98.
%     See also BOSTEP, BOCLOS, BOCOMP, BOINTEGR, NYQLOG, TFSHIFT from the
%     Bode toolbox.
%     Further explanations are given in Classical Feedback Control by
%     B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.

%  default
   if nargin == 0,  
      u = [0.2  2  8  32];        %  4 corner frequencies
      g = [-5  10  6 -12];        % gain in dB; band-pass system
      lslope = -1;                % -1 integrator, i.e. 6.02 dB/oct up-slope
      hslope = 3;                 %  3 integrators, i.e. -18.06 dB/oct
   end
   if nargin < 5,
      nml = 0;
   end
 
   k = length(u); h = [1:k];

%  frequency range of calculations and plots
   wmin = floor(log10(u(1)) - 0.33);
   wmax = ceil(log10(u(k)) + 0.33); 
   w = logspace(wmin,wmax);
   
%  low-frequency asymptote
   gain1 = (-lslope * (20 * log10(w/u(1))) + g(1))';
   phase1 = -lslope * 90;
 % phase1 = -lslope * 90 * ones(length(w),1); % can be used instead
 % semilogx(w, gain1, 'r--') % can be activated for purpose of teaching                
 % hold on

%  segment slopes; segment no.h ends at u(h)  
   for h = [2:k]
      sslope(h) = (g(h) - g(h-1))/log2(u(h)/u(h-1));
   end                       % [0 15 -2 -9]
 % sslope                    % can be activated for purpose of teaching

%  ray slopes
      rslope(1) = sslope(2) + 6.02 * lslope;
      rslope(k) = - 6.02 * hslope - sslope(k);
   for h = [2:(k-1)]
      rslope(h) = sslope(h + 1) - sslope(h);     
   end
 % rayslope_round = round(rslope)       % [9 -17 -7 -9], positive means up
                 
%  frequency responses for k rays 
   nray = 1; dray = [1 1 1];   % normalized ray with damping coefficient 0.5;  
   for h = [1:k]
      [nrayf, drayf] = lp2lp(nray,dray,u(h));   % specifying corner frequency
      tfr = freqs(nrayf, drayf, w); 
      rgain = 20*real(log10(tfr)); 
      rphase = (180/pi) * angle(tfr);          
      rasis = (-rslope(h)/12.041) * rgain';   % minus makes positive ray go up
      % semilogx(w, rasis,'w--');  % can be activated for purpose of teaching
      % semilogx(w, rphase,'g:');  % can be activated for purpose of teaching
      gain1 = (-rslope(h)/12.041) * rgain' + gain1;    % slope denormalization
      phase1 = phase1 - (rslope(h)/12.041) * rphase';
   end

 % nonminimum phase lag addition
   phase1 = phase1 - nml * (w/u(h))';

 % Bode diagram and phase response
   subplot(1,2,1)
   semilogx(w,gain1,'w', w,phase1,'g--', u,g,'wo')
   grid;
   xlabel('frequency, rad/sec');
   ylabel('gain, dB, and phase shift, degr');
   title('Bode diagram')
   hold off;
   zoom on;

 % Nyquist diagram on log plane
   subplot(1,2,2)
   plot(phase1, gain1,'w', -180,0,'wo')
   hold on
   gain2 = (-lslope * (20 * log10(u/u(1))) + g(1))'; % marks
   phase2 = -lslope * 90;
     for h = [1:k]           
        [nrayf, drayf] = lp2lp(nray,dray,u(h));   
        tfr = freqs(nrayf, drayf, u); 
        rgain = 20*real(log10(tfr)); 
        rgain = real(rgain);
        rphase = (180/pi) * angle(tfr);          
        gain2 = (-rslope(h)/12.041) * rgain' + gain2;
        phase2 = phase2 - (rslope(h)/12.041) * rphase';
     end
   phase2 = phase2 - nml * (u/u(h))';
   plot(phase2,gain2,'wo');
   hold off;
   set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90])
   grid
   axis([-270 -90 -20 70])
   title('Nyquist diagram')
   xlabel('phase shift, degr')
   ylabel('gain in dB')
   zoom on

Back

BOSTEP

function [wd,wc,width,numc,denc,numl,denl] = bostep(x,y,n,bnc,zetaz,zetap,typ)
% BOSTEP   Normalized loop transfer function with Bode step.
% function  [wd,wc,width,numc,denc,numl,denl] = bostep(x,y,n,bnc,zetaz,zetap,typ)
% Defaults: [wd,wc,width,numc,denc,numl,denl] = bostep(10,0.1667,3,0.1,0.6,0.4,2),
%   x is gain stability margin in dB, default x = 10;
%   y is phase stability margin y*pi rad, default y = 0.1667;
%   n defines asymptotic slope -6n dB/oct, default n = 3;
%   bnc (0 to 1) is n.p. shift in rad at the end of Bode step (at wc),
%   default bnc = 1;
%   zetaz is damping of zeros at beginning of Bode step (at wd), 
%   default zetaz = 0.6;
%   zetap is damping of poles at the end of Bode step (at wc), default
%   zetap = 0.4;
%   typ (1 or 2) is the number of poles at the origin, default typ = 2.
%   BOSTEP calculates the following parameters of the asymptotic Bode diagram 
%   withcrossover frequency 1 rad/sec: frequencies of beginning and end of 
%   the step wd, wc, and the width of the step width.
%   BOSTEP also produces numerator and denominator of rational loop transfer 
%   function numl(s)/denl(s) approximating the asymptotic loop Bode diagram,
%   where numl and denl contain the polynomial coefficients in descending 
%   powers of s. The typ is the number, 1 or 2, of the poles at the origin.
%   BOSTEP also calculates the minimum phase compensator transfer function
%   numc(s)/denc(s) for a single integrator plant with n.p. lag 
%   (-s + 2wc/bnc)/(s + 2wc/bnc).
%
%   Copyright (c) B. J. Lurie 3-1-98, La Crescenta.
%   See also BOCLOS, BONYQAS, BOCOMP, NYQLOG, TFSHIFT from the Bodestep toolbox.
%   Further explanations are given in Classical Feedback Control by
%   B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
  
   if nargin < 7,
      typ = 2;
   end
   if nargin < 6,   
      zetap = 0.4; 
   end 
   if nargin < 5,   
      zetaz = 0.6; 
   end 
   if nargin < 4,   
      bnc = 0.1; 
   end 
   if nargin < 3,   
      n = 3; 
   end 
   if nargin < 2,   
      y = 0.167; 
   end 
   if nargin == 0,   
      x = 10; 
   end 

%  **  Asymptotic diagram parameters
   width = 0.9*(n + (pi/2)*bnc)/(2*(1 - y));
   wd = 2^(x/12/(1 - y));
   wc = wd*width;

%  ** Poles and zeros for typ = 1
   z1 = 0.5*(1-y); 
   z23 = [1 2*zetaz*wd wd^2];         % a pair of complex zeros
   k  = 0.35*wc^2*sqrt(0.17/y);
   p1 = 0.07*(1-y);  p2 = 0.7*wd;  p3 = (wd + wc)/2.7;
   p45 = [1 2*zetap*wc 0.81*wc^2];   ;% a pair of complex poles

%  ** Poles and zeros for typ = 2
   if typ == 2,
       p1 = 0;
       p3 = (wd + wc)/2.2;
      z23 = 0.95*[1 2.2*zetaz*wd*sqrt(sqrt(y/0.17)) wd^2];
      p45 = [1 2.2*zetap*wc 0.81*wc^2];
   end

%  ** Minimum phase compensator
   numc = k * p3 * conv([1 z1], z23);
   denc = conv(conv([1 p1],[1 p2]),conv([1 p3],p45));
   
%  ** Loop with 1/s plant and n.p. lag (-s + 2wc/bnc)/(s + 2wc/bnc): 
   numl = conv(numc,[-1 2*wc/bnc]);
   denl = conv(conv(denc,[1 0]),[1 2*wc/bnc]);
   format short e
   numl
   denl

%   numl_mp = numc; ne nada!
%   denl_mp = conv(denc,[1 0]);ne nada!

%  ** L-plane Nyquist diagram
   [mag, phase] = bode(numl, denl);  plot(phase, 20*log10(mag),'w')
   set(gca,'XTick',[-240 -210 -180 -150 -120 -90])
   grid
   axis([-240 -90 -20 70])
   hold on
   nl = numl;
   dl = denl;
   [a,b] = bode (nl,dl,1); plot(b,20*log10(a),'wx')
       for w = [0.0156 0.03125 0.0625 0.125 0.25 0.5 2 4],
         [a,b] = bode (nl,dl,w); plot(b,20*log10(a),'w+');
       end
   plot(-180, 0, 'wo')
   hold off
   title('Nyquist diagram, x marks w = 1, + marks octaves')
   xlabel('loop phase shift in degrees')
   ylabel('loop gain in dB')
   zoom on

Back

BOCLOS

function [dclos, nprf1, dprf1, nprf2, dprf2, numtot, dentot] = boclos(numl,denl,y)

% BOCLOS  Normalized closed loop transfer function without and with 
%    prefilter, and frequency response of 4th order Bessel filter. Also,
%    step response, optional.
%    function [dclos, nprf1, dprf1, nprf2, dprf2, numtot, dentot] = 
%          boclos(numl,denl,y)
%    BOCLOS calculates denominator dclos of the closed loop transfer
%    function without prefilter (the numerator is numl);
%    prefilter1 and prefilter2 transfer functions nprf1(s)/dprf1(s)
%    and nprf2(s)/dprf2(s); and closed loop transfer functions:
%    numl(s)/dclos(s) without prefilter and numtot(s)/dentot(s)
%    with prefilter; the coefficients are in descending powers of s.
%    The defaults follow the default outputs of BOSTEP:
%    numl = [-9.6082  626.39  1947.8  3326.0  1093]
%    denl = [1  75.102  492.80  1716.3  3222.6  2217  0  0]; y = 0.1667;
%
%    Copyright (c) B. J. Lurie 3-15-98.
%    See also BOSTEP, BOCOMP, BONYQAS, NYQLOG, BOLAGNYQ, TFSHIFT from the
%    Bodestep toolbox.
%    Further explanations are given in Classical Feedback Control by
%    B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
  
 % default:
   if nargin < 3
      y = 0.1667;
   end
   if nargin < 2
      denl = [1  75.102  492.80  1716.3  3222.6  2217  0  0];
   end
   if nargin == 0
      numl = [-9.6082  626.39  1947.8  3326.0  1093];
   end
   format short e

% Normalized closed loop function denominator (without prefilter)
   if length(denl) - length(numl) == 2
       dclos = [0 0 numl] + denl;
   end
   if length(denl) - length(numl) == 3
       dclos = [0 0 0 numl] + denl;
   end
   if length(denl) - length(numl) == 4
       dclos = [0 0 0 0 numl] + denl;
   end

% First and second prefilter notches
  nprf1 = [1 y/0.1667 1];      dprf1 = [1 2 1];     % prefilter 1 
  nprf2 = [1 0.47*y/0.1667 1]; dprf2 = [1 0.6 1];   % prefilter 2

% Prefilter 3 provides gain hump
  nprf3 = [1 2.5 3.6];         dprf3 = [1 1.5 3.6]; % prefilter 3 

% Closed loop transfer function with prefilters 2 and 3
  numtot = conv(conv(nprf3,numl),conv(nprf1,nprf2));
  dentot = conv(conv(dprf3,dclos),conv(dprf1,dprf2));

% Closed loop function with extra 3rd order Bessel lowpass 
% prefilter,w = 0.25
  numtot = conv([15],numtot);
  dentot = conv([64 96 60 15],dentot);

% Open- and closed-loop and Bessel filter frequency responses
  w = logspace(-1,1);
  [magop,phaseop] = bode(numl, denl, w)     % open-loop
  [magcl,phasecl] = bode(numl, dclos, w)    % clos-loop without prefilter
  [magtot,phasetot] = bode(numtot, dentot, w)  % clos-loop with prefilter
  [magbf,phasebf] = bode(105, [1 10 45 105 105], w) % Bessel 4th order
  % [magbf,phasebf] = bode(15, [1 6 15 15], w) % Bessel 3th, option

% plotting
  subplot(2,1,1)
   semilogx(w, 20*log10(magop),'w', w,(phaseop),'w--');
   hold on
   semilogx(w, 20*log10(magcl),'b', w,(phasecl),'b--');  
   semilogx(w, 20*log10(magtot),'g', w,(phasetot),'g--');
   semilogx(w, 20*log10(magbf),'r', w,(phasebf),'r--');
   wmin = min(w); wmax = max(w);
   axis([wmin wmax -40 20])
   xlabel('rad/sec')
   ylabel('gain, dB; lag margin, degr')
   grid
   hold off

  subplot(2,1,2)  
   sstep(numtot, dentot) 
   xlabel('time, sec')
   ylabel('output')            
   grid
  zoom on

Back

BOINTEGR

function [ao, p1, numr, denr] = bointegr(numc,denc,wb)
% BOINTEGR  expresses the compensator transfer function as a sum of the
%    term of the partial fraction expansion ao/(s + p1) which is dominant
%    a lower frequencies and the rest of the terms combined into transfer
%    function numr/denr, where the coefficients in the vectors of 
%    coefficients are in the descending powers of s. numc is up to
%    9th order. wb is the crossover frequency or anything close enough, 
%    it is used here only to scale the frequency axis; default is 1.
%    Bode diagrams for the compensator and its two parallel paths are
%    plotted.
%    The default/demo is
%    numc = [9.6082  29.262  48.975  16.016];
%    denc = [1  6.8631  24.466  46.749  32.488  0];  wb = 1;
% function [ao, p1, numr, denr] = bointegr(numc,denc);
%
%    Copyright (c) B. Lurie 6-19-98.
%    See also BOLAGNYQ, BONYQAS, BOSTEP, BOCLOS, BOCOMP, NYQLOG, TFSHIFT from
%    the Bodestep toolbox.
%    Further explanations are given in Classical Feedback Control by
%    B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.

% default
   if nargin < 3
      wb = 1;
   end
   if nargin < 2
      numc = [9.6082  29.262  48.975  16.016];
      denc = [1  6.8631  24.466  46.749  32.488  0];
   end
   format short e

%  Partial fraction expansion
   ld = length(denc);
   [r,p,k] = residue(numc,denc);
   ao = r(ld-1); p1 = p(ld-1);

%  Combining the rest of the terms
   if ld == 2
      rr = 0;
      pr = 1;
   end
   if ld == 3
      rr = [r(ld-2)];
      pr = [p(ld-2)];
   end
   if ld == 4
      rr = [r(ld-3) r(ld-2)];
      pr = [p(ld-3) p(ld-2)];
   end
   if ld == 5
      rr = [r(ld-4) r(ld-3) r(ld-2)];
      pr = [p(ld-4) p(ld-3) p(ld-2)];
   end
   if ld == 6
      rr = [r(ld-5) r(ld-4) r(ld-3) r(ld-2)];
      pr = [p(ld-5) p(ld-4) p(ld-3) p(ld-2)];
   end
   if ld == 7
      rr = [r(ld-6) r(ld-5) r(ld-4) r(ld-3) r(ld-2)];
      pr = [p(ld-6) p(ld-5) p(ld-4) p(ld-3) p(ld-2)];
   end
   if ld == 8
      rr = [r(ld-7) r(ld-6) r(ld-5) r(ld-4) r(ld-3) r(ld-2)];
      pr = [p(ld-7) p(ld-6) p(ld-5) p(ld-4) p(ld-3) p(ld-2)];
   end
   if ld == 9
      rr = [r(ld-8) r(ld-7) r(ld-6) r(ld-5) r(ld-4) r(ld-3) r(ld-2)];
      pr = [p(ld-8) p(ld-7) p(ld-6) p(ld-5) p(ld-4) p(ld-3) p(ld-2)];
   end
   [numr,denr] = residue(rr,pr,k);
   numr = real(numr); denr = real(denr); % elimination of imaginary
                                 % parts caused by numerical errors        

%  Print the answers
   ao
   p1
   numr
   denr 
   % if further breaking into parallel channels needed, rr and pr

   wmin = floor(log10(wb)-3);
   wmax = ceil(log10(wb)+2);
   w = logspace(wmin,wmax);
  [mc,pc] = bode(numc,denc,w);
  [mr,pr] = bode(numr,denr,w); 
  [mi,pi] = bode(ao,[1 p1],w);

%  plotting
   wmin = min(w); wmax = max(w);

   subplot(211)
   semilogx(w, 20*log10(mc),'w', w, 20*log10(mr),'g', w, 20*log10(mi),'r');
   axis([wmin wmax -40 80])
   xlabel('rad/sec')
   ylabel('gain, dB')
   grid

   subplot(212)
   semilogx(w,(pc),'w', w,(pr),'g', w,(pi),'r');
   axis([wmin wmax -200 20])
   ylabel('phase shift, degr')
   grid
   hold off
   zoom on

Back

BOCOMP

function [numc_pl,denc_pl] = bocomp(wb,numc,denc,kmot,resist,mom_inert);
% BOCOMP  Compensator transfer function numc_pl(s)/denc_pl(s)for a dc 
%   motor system with loop response generated by bostep. The vectors
%   NUMC_PL and DENC_PL contain  polynomial coefficients in descending
%   powers of s, 
%   The driver voltage gain is 1. The sum of its output resistance with
%   the motor winding resistance is resist.The motor constant is kmot. 
%   The rigid body load moment of inertia is mom_inert.
%   The loop crossover frequency in rad/sec is wb.
%   The function numc(s)/denc(s) is the normalized compensator obtained
%   with function bostep for the nominal plant (1/s)(-s+wc/bnc)/(s+wc/bnc). 
%   BOCOMP calculates: the compensator polynomial coefficients, 
%                      the gain coefficient, the poles and the zeros,
%                      the biquad coefficients for complex poles and zeros,
%           and plots: Bode diagrams for the plant (in green) and
%                      compensator (in read).
%   Defaults/demo initiated by typing bocomp: 
%   mom_inert = 0.027; resist = 2; kmot = 0.7; wb = 100;
%   denc = [ 1  6.8631  24.466  46.749  32.488  0];
%   numc = [9.6082  29.262  48.975  16.016];
%   [numc_pl,denc_pl] = bocomp(wb,numc,denc,kmot,resist,mom_inert)
%
%   Copyright (c) B. J. Lurie 6-18-98.
%   See also BOSTEP, BOCLOS, NYQLOG, TFSCHIFT, BOLAGNYQ, BONYQAS, BOINTEGR
%   from the Bodestep toolbox.
%   Further explanations are given in Classical Feedback Control by
%   B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
  
   if nargin < 6,   
      mom_inert = 0.027; ; 
   end 
   if nargin < 5,   
      resist = 2; 
   end 
   if nargin < 4,   
      kmot = 0.7;
   end 
   if nargin < 3, 
      denc =  [1  7.9547  33.536  74.518  58.275  0];  
   end
   if nargin < 2,   
      numc =  [13.046  45.855  76.234  22.827];   % ratio of 4th to 6th
   end 
   if nargin == 0,   
      wb = 100;      
   end 

%  plant transfer function kmot/(resist*mom_inert)/[s(s+a)] where
   a = kmot^2/(resist*mom_inert);  	          %real pole of the plant
   npm = kmot/(resist*mom_inert); dpm = [1 a 0];  % minimum phase plant
                                                  % without flex modes

%  denormalized m.p. compensator for n.p. single integrator plant
   [numc_wb, denc_wb] = tfshift(numc,denc,wb); 

%  denormalized m.p. compensator for motor plant
   numc_pl = wb*deconv(conv(numc_wb,dpm),[1 0]);
   denc_pl = conv(denc_wb,npm);
   format short e

   gain_coeff = numc_pl(1)/denc_pl(1);
   roots_num = roots(numc_pl);
   biquad_num = [1 -2*real(roots_num(2)) abs(roots_num(2))^2];
   roots_den = roots(denc_pl);
   biquad_den = [1 -2*real(roots_den(2)) abs(roots_den(2))^2];

%  printout
   gain_coeff 
   roots_num   
   biquad_num
   roots_den   
   biquad_den
   numc_pl
   denc_pl

%  plotting
   wmin = floor(log10(wb)-3);
   wmax = ceil(log10(wb)+1);
   w = logspace(wmin,wmax);
   [magc,phasc] = bode(numc_pl,denc_pl,w);
   [magnpm, phasenpm] = bode(npm,dpm,w);
   semilogx(w,20*log10(magc),'r', w,phasc,'r--');     % compensator
   hold on
   semilogx(w,20*log10(magnpm),'g', w,phasenpm,'g--') % plant
   hold off
   xlabel('rad/sec; phase plotted with dashed lines')
   ylabel('dB, degr')
   title('gain and phase of plant and compensator')
   grid
   zoom on

Back

NDCP

function NDCP(w1,n1,d1,n2,d2,n3,d3)
%NDCP Nyquist diagram on logarithmic plane for describing function of
%   nonlinear dynamic compensator with parallel channels.
%   -->df-->link1--
%   --< >(+)-->link3--
%   ------->link2--
%function NDCP(w1,n1,d1,n2,d2,n3,d3)
%   The first channel incorporates a link with describing function
%   df=0:0.2:1 and a linear link with numerator n1 and denominator d1
%   of the transfer function. The second channel is a linear link with
%   numerator n2 and denominator d2.
%   The third linear link with transfer function numerator n3 and
%   denominator d3 is connected in series to the parallel connection.
%   The transfer functions are of the type n(s)/d(s), where n and d
%   contain polynomial coefficients in descending powers of s.
%   w1 is a frequency in rad/sec to be marked by both x and o.
%   The default is: w1=10; n1=200; d1=[1 0 0]; n2=10; d2=[1 0]; n3=d3=1;
%
%   Copyright (c) B. Lurie 6-30-99.
%   For further explanations see Classical Feedback Control by
%   B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 7
d3 = 1;
end
if nargin < 6
n3 = 1;
end
if nargin < 1
w1=10; n1=200; d1=[1 0 0]; n2=10; d2=[1 0];
end

denl = conv(conv(d1,d2),d3);
for df = 1 : -0.2 : 0,
n1 = df * n1;
v1 = conv(n1,d2);
v2 = conv(d1,n2);
ll = length(v1) - length(v2);
if ll < 0
adz = zeros(1,-ll);
num1 = [adz v1] + v2;
else
adz = zeros(1,ll);
num1 = v1 + [adz v2];
end
numl = conv(num1,n3);
[mag, phase] = bode(numl, denl);
plot(phase, 20*log10(mag),'w',-180, 0, 'ro')
hold on
[a,b] = bode (numl,denl,w1);
plot(b,20*log10(a),'wx',b,20*log10(a),'wo')
for w = [0.0156 0.03125 0.0625 0.125 0.25 0.5 2 4],
[a,b] = bode(numl,denl,w1*w); plot(b,20*log10(a),'w+');
end
end

hold off
grid
axis([-270 0 -20 70])
set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90 -60 -30 0])
title('Nyquist diagram, x marks w = w1, + marks octaves')
xlabel('loop phase shift in degrees')
ylabel('loop gain in dB')
zoom on

Back

BNDCP

function BNDCP(w1,n1,d1,n2,d2,n3,d3)
%BNDCP Bode diagram for describing function of nonlinear dynamic
%   compensator with parallel channels.
%   ->df-->link1-
%   --< >(+)-->link3--
%   ->--->link2--
%function BNDCP(w1,n1,d1,n2,d2,n3,d3)
%   The first channel incorporates a link with describing function
%   df=0:0.2:1 and a linear link with numerator n1 and denominator d1
%   of the transfer function. The second channel is a linear link with
%   numerator n2 and denominator d2.
%   The third linear link with transfer function numerator n3 and
%   denominator d3 is connected in series to the parallel connection.
%   The transfer functions are of the type n(s)/d(s), where n and d
%   contain polynomial coefficients in descending powers of s.
%   w1 is a frequency in rad/sec to be marked by both x and o.
%
%   The default is: w1=10; n1=200; d1=[1 0 0]; n2=10; d2=[1 0]; n3=d3=1;
%   Copyright (c) B. Lurie 6-30-99.
%   For further explanations see Classical % Feedback Control by
%   B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 7
d3 = 1;
end
if nargin < 6
n3 = 1;
end
if nargin < 1
w1=10; n1=200; d1=[1 0 0]; n2=10; d2=[1 0];
end

denl = conv(conv(d1,d2),d3);
for df = 1 : -0.2 : 0,
n1 = df * n1;
v1 = conv(n1,d2);
v2 = conv(d1,n2);
ll = length(v1) - length(v2);
if ll < 0
adz = zeros(1,-ll);
num1 = [adz v1] + v2;
else
adz = zeros(1,ll);
num1 = v1 + [adz v2];
end
numl = conv(num1,n3);
bode(numl, denl);
hold on
end
hold off
zoom on

Back

NDCB

function NDCB(w1,n1,d1,n2,d2,n3,d3)
%NDCB Nyquist diagram on logarithmic plane for describing function of
%   nonlinear dynamic compensator with nonlinear feedback channel.
%   -->(+)----->link1----->link3--
%   (-) |
%   ^--link2<--df<--
%function NDCB(w1,n1,d1,n2,d2,n3,d3)
%   The first link is linear with numerator n1 and denominator d1
%   of the transfer function. The feedback channel incorporates a
%   link with describing function df=0:0.2:1 an a second linear
%   link with numerator n2 and denominator d2.
%   The third linear link with transfer function numerator n3 and
%   denominator d3 is connected in series to the first link.
%   The transfer functions are of the type n(s)/d(s), where n and d
%   contain polynomial coefficients in descending powers of s.
%   w1 is a frequency in rad/sec to be marked by both x and o.
%
%   The default is: w1=10; n1=100; d1=[1 0]; n2=1; d2=[1];
%   n3=[1 4]; d3=[1 15 0];
%
%   Copyright (c) B. Lurie 6-30-99.
%   For further explanations see Classical Feedback Control by
%   B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 7
d3 = [1 14 0];
end
if nargin < 5
n3 = [1 2];
end
if nargin < 1
w1=8; n1=[100 100]; d1=[1 0 0]; n2=1; d2=[1];
end
numl = conv(conv(n1,d2),n3);
v1 = conv(d1,d2);
v2 = conv(n1,n2);
for df = 1 : -0.2 : 0,
v2 = df * v2;

ll = length(v1) - length(v2);
if ll < 0
adz = zeros(1,-ll);
den1 = [adz v1] + v2;
else
adz = zeros(1,ll);
den1 = v1 + [adz v2];
end
denl = conv(den1,d3);
[mag, phase] = bode(numl, denl);
plot(phase, 20*log10(mag),'w',-180, 0, 'ro')
hold on
[a,b] = bode (numl,denl,w1);
plot(b,20*log10(a),'wx',b,20*log10(a),'wo')
for w = [0.0156 0.03125 0.0625 0.125 0.25 0.5 2 4],
[a,b] = bode(numl,denl,w1*w); plot(b,20*log10(a),'w+');
end
end
hold off
grid
axis([-270 0 -20 70])
set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90 -60 -30 0])
title('Nyquist diagram, x marks w = w1, + marks octaves')
xlabel('loop phase shift in degrees')
ylabel('loop gain in dB')
zoom on

Back

BNDCB

function BNDCB(w1,n1,d1,n2,d2,n3,d3)
%BNDCB Bode diagrams for describing functions of nonlinear dynamic
%   compensator with nonlinear feedback channel.
%   -->(+)----->link1----->link3--
%   (-) |
%   ^--link2<--df<--
%function BNDCB(w1,n1,d1,n2,d2,n3,d3)
%   The first link is linear with numerator n1 and denominator d1
%   of the transfer function. The feedback channel incorporates a
%   link with describing function df = 0:0.2:1 and a second linear
%   link with numerator n2 and denominator d2.
%   The third linear link with transfer function numerator n3 and
%   denominator d3 is connected in series to the first link.
%   The transfer functions are of the type n(s)/d(s), where n and d
%   contain polynomial coefficients in descending powers of s.
%   w1 is a frequency in rad/sec to be marked by both x and o.
%
%   The default is: w1=10; n1=100; d1=[1 0]; n2=1; d2=[1];
%   n3=[1 4]; d3=[1 15 0];
%
%   Copyright (c) B. Lurie 6-30-99.
%   For further explanations see Classical Feedback Control by
%   B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 7
d3 = [1 14 0];
end
if nargin < 5
n3 = [1 2];
end
if nargin < 1
w1=8; n1=[100 100]; d1=[1 0 0]; n2=1; d2=[1];
end
numl = conv(conv(n1,d2),n3);
v1 = conv(d1,d2);
v2 = conv(n1,n2);
for df = 1 : -0.2 : 0,
v2 = df * v2;

ll = length(v1) - length(v2);
if ll < 0
adz = zeros(1,-ll);
den1 = [adz v1] + v2;
else
adz = zeros(1,ll);
den1 = v1 + [adz v2];
end
denl = conv(den1,d3);

bode(numl, denl);
hold on
end
hold off
zoom on

Back