% ***** MATLAB Code Starts Here
%
olz = [-1;-2]; % Transfer function zeros
olp = [0;-.3;-.7+j*2;-.7-j*2;-15]; % Transfer function poles
K = 100; % Transfer function gain
%
num = K * real(poly(olz)); % Numerator polynomial
den = real(poly(olp)); % Denominator polynomial
%
nn = length(num); % Length of polynomials
nd = length(den);
%
KB = num(nn) / den(nd-1); % Gain in "Bode" form
%
w = logspace(-3,3,601); % Frequency vector
%
magKB = abs(KB)*ones(1,length(w)); % Mag & Phase of Bode gain
phKB = zeros(1,length(w));
%
magz1 = abs(polyval([-1/olz(1) 1],j*w)); % Mag & Phase of 1st zero
phz1 = angle(polyval([-1/olz(1) 1],j*w))*180/pi;
%
magz2 = abs(polyval([-1/olz(2) 1],j*w)); % Mag & Phase of 2nd zero
phz2 = angle(polyval([-1/olz(2) 1],j*w))*180/pi;
%
magp1 = 1./w; % Mag & Phase of 1st pole
php1 = -90*ones(1,length(w));
%
magp2 = 1./abs(polyval([-1/olp(2) 1],j*w)); % Mag & Phase of 2nd pole
php2 = -angle(polyval([-1/olp(2) 1],j*w))*180/pi;
%
denp3 = conv([1 -olp(3)],[1 -olp(4)]); % Polynomial for complex poles
denp3 = denp3 / denp3(3);
magp34 = 1./abs(polyval(denp3,j*w)); % Mag & Phase of complex poles
php34 = -angle(polyval(denp3,j*w))*180/pi;
%
magp5 = 1./abs(polyval([-1/olp(5) 1],j*w)); % Mag & Phase of 5th pole
php5 = -angle(polyval([-1/olp(5) 1],j*w))*180/pi;
%
[mag,ph] = bode(num,den,w); % Mag & Phase of total G(s)
%
% Magnitudes of the individual poles & zeros
%
figure(1),clf,semilogx(w,20*log10(magKB),w,20*log10(magz1),...
w,20*log10(magz2),w,20*log10(magp1),w,20*log10(magp2),w,20*log10(magp34),...
w,20*log10(magp5)),grid,xlabel('Frequency (r/s)'),ylabel('Magnitude (db)'),...
title('Magnitudes of Individual Terms'),gtext('KB'),gtext('z1'),
gtext('z2'),gtext('p1'),gtext('p2'),gtext('p34'),gtext('p5')
%
% Phases of the individual poles & zeros
%
figure(2),clf,semilogx(w,phKB,w,phz1,w,phz2,w,php1,w,php2,w,php34,...
w,php5),grid,xlabel('Frequency (r/s)'),ylabel('Phase (deg)'),...
title('Phases of Individual Terms'),gtext('KB'),gtext('z1'),
gtext('z2'),gtext('p1'),gtext('p2'),gtext('p34'),gtext('p5')
%
% Demonstration of asymptotic properties of magnitude and phase
%
wmagz1 = [min(w) 1 max(w)];
wphz1 = [min(w) 0.1 10 max(w)];
asymagz1 = [1 1 1000];
asyphz1 = [0 0 90 90];
%
figure(3),clf,semilogx(w,20*log10(magz1),wmagz1,20*log10(asymagz1)),grid,...
xlabel('Frequency (r/s)'),ylabel('Magnitude (db)'),...
title('Actual & Asymptotic Magnitudes for (s + 1)')
%
figure(4),clf,semilogx(w,phz1,wphz1,asyphz1),grid,...
xlabel('Frequency (r/s)'),ylabel('Phase (deg)'),...
title('Actual & Asymptotic Phases for (s + 1)')
%
% Total magnitude and phase computed from the individual terms
%
mag1 = (magKB.*magz1.*magz2) .* (magp1.*magp2.*magp34.*magp5);
ph1 = (phKB+phz1+phz2) + (php1+php2+php34+php5);
%
figure(5),clf,semilogx(w,20*log10(mag),'y',w,20*log10(mag1),'r:'),grid,...
xlabel('Frequency (r/s)'),ylabel('Magnitude (db)'),...
title('Magnitude of the Transfer Function')
%
figure(6),clf,semilogx(w,ph,'y',w,ph1,'r:'),grid,...
xlabel('Frequency (r/s)'),ylabel('Phase (deg)'),...
title('Phase of the Transfer Function')
%
% ***** MATLAB Code Stops Here
Click the
icon to return to the Dr. Beale's home page
Lastest revision on
Friday, May 26, 2006 9:05 PM