Function lnyquist: Plotting the Nyquist Frequency Response on a Log-base-2 Scale
Below is the function lnyquist.m. This function is a modified version of MATLAB's nyquist command, and has the same attributes as the original, with a few improvements. The function lnyquist.m plots
(log2(1+abs(G(jw))),angle(G(jw))
in polar coordinates by taking the log of the magnitude, the magnitude scale is compressed and the overall shape of the Nyquist plot is easier to see on the screen. We use log base 2 and add one to the magnitude so as to preserve the key attributes of the -1 point for the Nyquist plot.
The lnyquist function also takes poles on the imaginary axis into account when creating the Nyquist plot, and plots around them.
Copy the following text into a file lnyquist.m. Put the file in the same directory as the MATLAB software, or in a directory which is contained in MATLAB's search path.
function [reout,imt,w] = lnyquist(a,b,c,d,iu,w)
%LNYQUIST Nyquist frequency response for continuous-time linear systems.
%
% This Version of the NYQUIST Command takes into account poles at the
% jw-axis and loops around them when creating the frequency vector in order
% to produce the appropriate Nyquist Diagram (The NYQUIST command does
% not do this and therefore produces an incorrect plot when we have poles in the
% jw axis).
%
% NOTE: This version of LNYQUIST does not account for pole-zero
% cancellation. Therefore, the user must simplify the transfer function before using
% this command.
%
% LNYQUIST(A,B,C,D,IU) produces a Nyquist plot from the single input
% IU to all the outputs of the system:
% . -1
% x = Ax + Bu G(s) = C(sI-A) B + D
% y = Cx + Du RE(w) = real(G(jw)), IM(w) = imag(G(jw))
%
% The frequency range and number of points are chosen automatically.
%
% LNYQUIST(NUM,DEN) produces the Nyquist plot for the polynomial
% transfer function G(s) = NUM(s)/DEN(s) where NUM and DEN contain
% the polynomial coefficients in descending powers of s.
%
% LNYQUIST(A,B,C,D,IU,W) or LNYQUIST(NUM,DEN,W) uses the user-supplied
% freq. vector W which must contain the frequencies, in radians/sec,
% at which the Nyquist response is to be evaluated. When invoked
% with left hand arguments,
% [RE,IM,W] = LNYQUIST(A,B,C,D,...)
% [RE,IM,W] = LNYQUIST(NUM,DEN,...)
% returns the frequency vector W and matrices RE and IM with as many
% columns as outputs and length(W) rows. No plot is drawn on the
% screen.
% See also: LOGSPACE,MARGIN,BODE, and NICHOLS.
%
% J.N. Little 10-11-85
% Revised ACWG 8-15-89, CMT 7-9-90, ACWG 2-12-91, 6-21-92,
% AFP 2-23-93
% WCM 8-30-97
%
********************************************************************
Modifications made to the nyquist - takes into account poles on jw axis.
then goes around these to make up frequency vector
%
%
if nargin==0, eval('exresp(''nyquist'')'), return, end
--- Determine which syntax is being used ---
nargin1 = nargin;
nargout1 = nargout;
if (nargin1==1), % System form without frequency vector
[num,den]=tfdata(a,'v');
z = roots(num);
p = roots(den);
zp = [z;p];
wpos = zp(find(abs(zp)>0));
if(min(abs(p)) == 0)
wstart = max(eps, 0.03*min([1;wpos]));
else
wstart = max(eps, 0.03*min(abs(zp)));
end
wstop = max([1000;30*wpos]);
w = logspace(log10(wstart),log10(wstop),max(51,10*max(size(zp))+1));
%w = freqint2(num,den,30);
[ny,nn] = size(num); nu = 1;
%error('Wrong number of input arguments.');
elseif (nargin1==2),
if(isa(a,'ss')|isa(a,'tf')|isa(a,'zpk')) % System with frequency vector
[num,den]=tfdata(a,'v');
w = b;
else % Transfer function form without frequency vector
num = a; den = b;
z = roots(num);
p = roots(den);
zp = [z;p];
wpos = zp(find(abs(zp)>0));
if(min(abs(p)) == 0)
wstart = max(eps, 0.03*min([1;wpos]));
else
wstart = max(eps, 0.03*min(abs(zp)));
end
wstop = max([1000;30*wpos]);
w = logspace(log10(wstart),log10(wstop),max(51,10*max(size(zp))+1));
%w = freqint2(num,den,30);
end
[ny,nn] = size(num); nu = 1;
elseif (nargin1==3), % Transfer function form with frequency vector
num = a; den = b;
w = c;
[ny,nn] = size(num); nu = 1;
elseif (nargin1==4), % State space system, w/o iu or frequency vector
error(abcdchk(a,b,c,d));
[num,den] = ss2tf(a,b,c,d);
[z,p,k]= ss2zp(a,b,c,d);
[num,den] = zp2tf(z,p,k);
zp = [z;p];
wpos = zp(find(abs(zp)>0));
if(min(abs(p)) == 0)
wstart = max(eps, 0.03*min([1;wpos]));
else
wstart = max(eps, 0.03*min(abs(zp)));
end
wstop = max([1000;30*wpos]);
w = logspace(log10(wstart),log10(wstop),max(51,10*max(size(zp))+1));
%w = freqint2(a,b,c,d,30);
nargin1 = 2;%[iu,nargin,re,im]=mulresp('nyquist',a,b,c,d,w,nargout1,0);
%if ~iu, if nargout, reout = re; end, return, end
[ny,nu] = size(d);
elseif (nargin1==5), % State space system, with iu but w/o freq. vector
error(abcdchk(a,b,c,d));
z = tzero(a,b,c,d);
p = eig(a);
zp = [z;p];
wpos = zp(find(abs(zp)>0));
if(min(abs(p)) == 0)
wstart = max(eps, 0.03*min([1;wpos]));
else
wstart = max(eps, 0.03*min(abs(zp)));
end
wstop = max([1000;30*wpos]);
w = logspace(log10(wstart),log10(wstop),max(51,10*max(size(zp))+1));
%w = freqint2(a,b,c,d,30);
[ny,nu] = size(d);
else
error(abcdchk(a,b,c,d));
[ny,nu] = size(d);
end
if nu*ny==0, im=[]; w=[]; if nargout~=0, reout=[]; end, return, end
*********************************************************************
depart from the regular nyquist program here
now we have a frequency vector, a numerator and denominator
now we create code to go around all poles and zeroes here.
if (nargin1==5) | (nargin1 ==4) | (nargin1 == 6)
[num,den]=ss2tf(a,b,c,d);
end
tol = 1e-6; %defined tolerance for finding imaginary poles
z = roots(num);
p = roots(den);
***** If all of the poles are at the origin, just move them a tad to the left***
if norm(p) == 0
if(isempty(z))
tad = 1e-1;
else
tad = min([1e-1; 0.1*abs(z)]);
end
length_p = length(p);
p = -tad*ones(length_p,1);
den = den(1,1)*[1 tad];
for ii = 2:length_p
den = conv(den,[1 tad]);
end
zp = [z;p];
wpos = zp(find(abs(zp)>0));
if(min(abs(p)) == 0)
wstart = max(eps, 0.03*min([1;wpos]));
else
wstart = max(eps, 0.03*min(abs(zp)));
end
wstop = max([1000;30*wpos]);
w = logspace(log10(wstart),log10(wstop),max(51,10*max(size(zp))+1));
%w = freqint2(num,den,30);
end
zp = [z;p]; % combine the zeros and poles of the system
nzp = length(zp); % number of zeros and poles
ones_zp=ones(nzp,1);
%Ipo = find((abs(real(p))<1e-6) & (imag(p)>=0)) %index poles with zero real part + non-neg imag part
Ipo = find((abs(real(p)) < tol) & (imag(p)>=0)); %index poles with zero real part + non-neg imag part
if ~isempty(Ipo) %
**** only if we have such poles do we do the following:*************************
po = p(Ipo); % poles with 0 real part and non-negative imag part
check for distinct poles
[y,ipo] = sort(imag(po)); % sort imaginary parts
po = po(ipo);
dpo = diff(po);
idpo = find(abs(dpo)>tol);
idpo = [1;idpo+1]; % indexes of the distinct poles
po = po(idpo); % only distinct poles are used
nIpo = length(idpo); % # of such poles
originflag = find(imag(po)==0); % locate origin pole
s = []; % s is our frequency response vector
%w = sqrt(-1)*w; % create a jwo vector to evaluate all frequencies with
for ii=1:nIpo % for all Ipo poles
[nrows,ncolumns]=size(w);
if nrows == 1
w = w.'; % if w is a row, make it a column
end;
if nIpo == 1
r(ii) = .1;
else % check distances to other poles and zeroes
pdiff = zp-po(ii)*ones_zp; % find the differences between
% poles you are checking and other poles and zeros
ipdiff = find(abs(pdiff)>tol); % ipdiff is all nonzero differences
r(ii)=0.2*min(abs(pdiff(ipdiff))); % take half this difference
r(ii)=min(r(ii),0.1); % take the minimum of this diff.and .1
end;
t = linspace(-pi/2,pi/2,25);
if (ii == originflag)
t = linspace(0,pi/2,25);
end; % gives us a vector of points around each Ipo
s1 = po(ii)+r(ii)*(cos(t)+sqrt(-1)*sin(t)); % detour here
s1 = s1.'; % make sure it is a column
% Now here I reconstitute s - complex frequency - and
% evaluate again.
bottomvalue = po(ii)-sqrt(-1)*r(ii); % take magnitude of imag part
if (ii == originflag) % if this is an origin point
bottomvalue = 0;
end;
topvalue = po(ii)+sqrt(-1)*r(ii); % the top value where detour stops
nbegin = find(imag(w) < imag(bottomvalue)); %
nnbegin = length(nbegin); % find all the points less than encirclement
if (nnbegin == 0)& (ii ~= originflag) % around jw root
sbegin = 0
else sbegin = w(nbegin);
end;
nend = find(imag(w) > imag(topvalue)); % find all points greater than
nnend = length(nend); % encirclement around jw root
if (nnend == 0)
send = 0
else send = w(nend);
end
w = [sbegin; s1; send]; % reconstituted half of jw axis
end
else
w = sqrt(-1)*w;
end
%end % this ends the loop for imaginary axis poles
*************************************************************
back to the regular nyquist program here
Compute frequency response
if (nargin1==1)|(nargin1==2)|(nargin1==3)
gt = freqresp(num,den,w);
else
gt = freqresp(a,b,c,d,iu,w);
end
***********************************************************
nw = length(gt);
mag = abs(gt); % scaling factor added
ang = angle(gt);
mag = log2(mag+1); % scale by log2(mag) throughout
for n = 1:nw
h(n,1) = mag(n,1)*(cos(ang(n,1))+sqrt(-1)*sin(ang(n,1)));
end; % recalculate G(jw) with scaling factor
gt = h;
***********************************************************
ret=real(gt);
imt=imag(gt);
If no left hand arguments then plot graph.
if nargout==0,
status = ishold;
plot(ret,imt,'r-',ret,-imt,'g--')
plot(real(w),imag(w))
modifications added here
%*******************************************
%
set(gca, 'YLimMode', 'auto')
limits = axis;
% Set axis hold on because next plot command may rescale
set(gca, 'YLimMode', 'auto')
set(gca, 'XLimMode', 'manual')
hold on
% Make arrows
for k=1:size(gt,2),
g = gt(:,k);
re = ret(:,k);
im = imt(:,k);
sx = limits(2) - limits(1);
[sy,sample]=max(abs(2*im));
arrow=[-1;0;-1] + 0.75*sqrt(-1)*[1;0;-1];
sample=sample+(sample==1);
reim=diag(g(sample,:));
d=diag(g(sample+1,:)-g(sample-1,:));
% Rotate arrow taking into account scaling factors sx and sy
d = real(d)*sy + sqrt(-1)*imag(d)*sx;
rot=d./abs(d); % Use this when arrow is not horizontal
arrow = ones(3,1)*rot'.*arrow;
scalex = (max(real(arrow)) - min(real(arrow)))*sx/50;
scaley = (max(imag(arrow)) - min(imag(arrow)))*sy/50;
arrow = real(arrow)*scalex + sqrt(-1)*imag(arrow)*scaley;
xy =ones(3,1)*reim' + arrow;
xy2=ones(3,1)*reim' - arrow;
[m,n]=size(g);
hold on
plot(real(xy),-imag(xy),'r-',real(xy2),imag(xy2),'g-')
end
xlabel('Real Axis'), ylabel('Imag Axis')
limits = axis;
% Make cross at s = -1 + j0, i.e the -1 point
if limits(2) >= -1.5 & limits(1) <= -0.5 % Only plot if -1 point is not far out.
line1 = (limits(2)-limits(1))/50;
line2 = (limits(4)-limits(3))/50;
plot([-1+line1, -1-line1], [0,0], 'w-', [-1, -1], [line2, -line2], 'w-')
end
% Axis
plot([limits(1:2);0,0]',[0,0;limits(3:4)]','w:');
plot(-1,0,'+k');
if ~status, hold off, end % Return hold to previous status
return % Suppress output
end
%reout = ret;
% plot(real(p),imag(p),'x',real(z),imag(z),'o');
