ngrid.m

数字通信第四版原书的例程

function [] = ngrid(s)
%NGRID	NGRID generates grid lines for a Nichols chart. The chart relates
%	the complex number H/(1+H) to H, where H is any complex number.
%
%	NGRID generates a grid over the region -40 db to 40 db in 
%	magnitude and -360 degrees to 0 degrees in phase with lines of 
%	constant mag(H/(1+H)) and angle(H/(1+H)) drawn in.  NGRID plots 
%	the Nichols chart grid over an existing Nichols plot such as one
%	generated with NICHOLS or DNICHOLS.  
%
%	NGRID('new') clears the graphics screen before plotting the grid
%	and sets HOLD ON so that the Nichols response can be plotted over
%	the grid using
%		ngrid('new')
%		nichols(num,den); or nichols(a,b,c,d,iu);
%
%	See also: NICHOLS and DNICHOLS.

%	J.N. Little 2-23-88
%	Revised: CMT 7-12-90, ACWG 6-21-92, Wes W 8-17-92.
%	Copyright (c) 1986-93 by the MathWorks, Inc.

status = ishold;
limit=axis;
if nargin == 1
  hold off
else
  axis(axis)
  hold on
  if limit(3) >= 0,
    disp('This is not a Nichols Chart. Use ngrid(''new'') for new Nichols Chart grid');
    hold off
    return
  end;
end
% Define desired gain and phase curves
mm = [6 3 1 .5 .25 0 -1 -3 -6 -12 -20 -40 ];
mx = [6 3 2 1 .75 .5 .4 .3 .25 .2 .15 .1 .05 0 -.05 -.1 -.15 -.2 -.25 -.3 -.4 -.5 -.75 -1 -2 -3 -4 -5 -6 -9 -12 -16 -20 -30 -40 ];
mm = 10 .^(mm/20);
mx = 10 .^(mx/20);

pp = [-359 -355 -350 -340 -330 -300 -270 -240 -210 -180 -150 -120 -90 -60 -30 -20 -10 -5 -1];
px = [-359 -358 -357 -356 -355 -352.5 -350 -345 -340 -335 -330 -315 -300 -285 -270 -255 -240 -225 -210 -195 -180 -165 -150 -135 -120 -105 -90 -75 -60 -45 -30 -25 -20 -15 -10 -5 -4 -3 -2 -1];

% --- Perform mapping and plot Nichols grid lines ---

i = sqrt(-1);

% Plot phase lines
[p,m] = meshgrid(pp,mx);
z = m .*exp(i*p/180*pi);
g = z./(1-z);
gain = 20*log10(abs(g));
phase = rem(angle(g)/pi*180+360,360)-360;
plot(phase,gain,'w:');
hold on

% Plot Magnitude lines
[p,m] = meshgrid(px,mm);
z = m .*exp(i*p/180*pi);
g = z./(1-z);
gain = 20*log10(abs(g));
phase = rem(angle(g)/pi*180+360,360)-360;
plot(phase',gain','w:')

% Now put magnitude labels (in db) on the curve.
[mp,np]=size(p);
mm = 20*log10(mm);
nm=length(mm)+1;
for i=1:nm-1
%  text(phase(i,np),gain(i,np),sprintf('%.3g',mm(nm-i)))
  text(phase(i,np),gain(i,np),sprintf('%.3g db',mm(i)))
end
xlabel('Open-Loop Phase (deg)')
ylabel('Open-Loop Gain (db)')

if nargin == 1 %	ngrid('new')
	set(gca,'xlim',[-360, 0]);
	set(gca,'ylim',[-40, 40]);
end

if (nargin==0) & (~status), hold off, end	    % Return hold to previous status