dnyquist.m

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

function [reout,im,w] = dnyquist(a,b,c,d,Ts,iu,w)
%DNYQUIST Nyquist frequency response for discrete-time linear systems.
%	DNYQUIST(A,B,C,D,Ts,IU) produces a Nyquist plot from the single
%	input IU to all the outputs of the system:      -1
%         x[n+1] = Ax[n] + Bu[n]    G(w) = C(exp(jwT)I-A) B + D  
%         y[n]   = Cx[n] + Du[n]    RE(w) = real(G(w)), IM(w) = imag(G(w))
%	The frequency range and number of points are chosen automatically.
%
%	DNYQUIST(NUM,DEN,Ts) produces the Nyquist plot for the polynomial
%	transfer function G(z) = NUM(z)/DEN(z) where NUM and DEN contain
%	the polynomial coefficients in descending powers of z. 
%
%	DNYQUIST(A,B,C,D,Ts,IU,W) or DNYQUIST(NUM,DEN,Ts,W) uses the user-
%	supplied freq. vector W which must contain the frequencies, in 
%	rad/sec, at which the Nyquist response is to be evaluated.  
%	Aliasing will occur at frequencies greater than the Nyquist 
%	frequency (pi/Ts). With left hand arguments,
%		[RE,IM,W] = DNYQUIST(A,B,C,D,Ts,...)
%		[RE,IM,W] = DNYQUIST(NUM,DEN,Ts,...) 
%	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,DBODE, and DNICHOLS.

% 	J.N. Little 10-11-85
%	Revised A.C.W.Grace 8-15-89, 2-12-91, 6-21-92
%	Revised Clay M. Thompson 7-9-90
%	Copyright (c) 1986-93 by the MathWorks, Inc.

if nargin==0, eval('dexresp(''dnyquist'')'), return, end

error(nargchk(3,7,nargin));

% --- Determine which syntax is being used ---
if (nargin==3),		% Transfer function without frequency vector
	num = a; den = b;
	Ts=c;
	w=dfrqint2(num,den,Ts,30);
	[ny,nn] = size(num); nu = 1;

elseif (nargin==4), 	% Transfer function form with frequency vector
	num = a; den = b;
	Ts=c; w=d;
	[ny,nn] = size(num); nu = 1;

elseif (nargin==5),	% State space system without iu or freq. vector
	error(abcdchk(a,b,c,d));
	w=dfrqint2(a,b,c,d,Ts,30);
	[iu,nargin,re,im]=dmulresp('dnyquist',a,b,c,d,Ts,w,nargout,0);
	if ~iu, if nargout, reout = re; end, return, end
	[ny,nu] = size(d);

elseif (nargin==6),	% State space system, with iu but w/o freq. vector
	error(abcdchk(a,b,c,d));
	w=dfrqint2(a,b,c,d,Ts,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

% Compute frequency response
if (nargin==3)|(nargin==4)
	g=freqresp(num,den,exp(sqrt(-1)*w*Ts));
else
	g=freqresp(a,b,c,d,iu,exp(sqrt(-1)*w*Ts));
end
re = real(g);
im = imag(g);

% If no left hand arguments then plot graph.
if nargout==0,
	status = ishold;

	% Make arrows
	sx=max(re)-min(re);   [sy,sample]=max(abs(2*im));
	arrow=[-1;0;-1]*sx/50+sqrt(-1)*[1;0;-1]*sy/50;
	sample=sample+(sample==1);
	reim=diag(g(sample,:));
	d=diag(g(sample,:)-g(sample-1,:)); 
	rot = sign(d);  % Arrow is always horizontal i.e. -1 or +1
	% rot=d./abs(d)  Use this when arrow is not horizontal
	xy =ones(3,1)*reim' + ones(3,1)*rot'.*arrow;
	xy2=ones(3,1)*reim' - ones(3,1)*rot'.*arrow;
	[m,n]=size(g);
	if (n==1)
		plot(re,im,'r-',re,-im,'g--',real(xy),-imag(xy),'r-',real(xy2),imag(xy2),'g-')
	else
		plot(re,im,re,-im,'--');
		hold on
		plot(real(xy),-imag(xy),'-',real(xy2),imag(xy2),'-')
	end
	hold on

	xlabel('Real Axis'), ylabel('Imag Axis')

	% Make cross at s = -1 + j0, i.e the -1 point
	limits = axis;
	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-')

% Axis
	limits = axis;
	plot([limits(1:2);0,0]',[0,0;limits(3:4)]','w:');

	if ~status, hold off, end	% Return hold to previous status
	return % Suppress output
end
reout = re;