function [magout,phase,w] = dbode(a,b,c,d,Ts,iu,w)
%DBODE	Bode frequency response for discrete-time linear systems.
%	DBODE(A,B,C,D,Ts,IU) produces a Bode plot from the single input IU
%	to all the outputs of the discrete state-space system (A,B,C,D).
%	IU is an index into the inputs of the system and specifies which
%	input to use for the Bode response.  Ts is the sample period.  
%	The frequency range and number of points are chosen automatically.
%
%	DBODE(NUM,DEN,Ts) produces the Bode 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. 
%
%	DBODE(A,B,C,D,Ts,IU,W) or DBODE(NUM,DEN,Ts,W) uses the user-
%	supplied freq. vector W which must contain the frequencies, in 
%	radians/sec, at which the Bode response is to be evaluated.  
%	Aliasing will occur at frequencies greater than the Nyquist 
%	frequency (pi/Ts).  See LOGSPACE to generate log. spaced frequency
%	vectors.  With left hand arguments,
%		[MAG,PHASE,W] = DBODE(A,B,C,D,Ts,...)
%		[MAG,PHASE,W] = DBODE(NUM,DEN,Ts,...) 
%	returns the frequency vector W and matrices MAG and PHASE (in 
%	degrees) with as many columns as outputs and length(W) rows.  No 
%	plot is drawn on the screen. 
%	See also: LOGSPACE,SEMILOGX,MARGIN,DNICHOLS, and DNYQUIST.

%	J.N. Little 10-11-85
%	Revised Andy Grace 8-15-89 2-4-91 6-20-92
%	Revised Clay M. Thompson 7-10-90
%	Copyright (c) 1986-93 by the MathWorks, Inc.


nargs = nargin;
if nargs==0, eval('dexresp(''dbode'')'), return, end

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

% --- Determine which syntax is being used ---
if (nargs==3),
	num = a; den = b;
	if length(c)==1,	% Transfer function without frequency vector
		Ts=c;
		w=dfrqint(num,den,Ts,30);
	else  		% Assume this is the old syntax DBODE(NUM,DEN,W)
		disp('Warning: You are using the old syntax.  Use DBODE(NUM,DEN,Ts,W) instead.')
		Ts = 1; w=c;
	end
	[ny,nn] = size(num); nu = 1;

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

elseif (nargs==5),	% State space system without iu or freq. vector
	error(abcdchk(a,b,c,d));
	w=dfrqint(a,b,c,d,Ts,30);
	[iu,nargs,mag,phase]=dmulresp('dbode',a,b,c,d,Ts,w,nargout,1);
	if ~iu, if nargout, magout = mag; end, return, end
	[ny,nu] = size(d);

elseif (nargs==6),
	error(abcdchk(a,b,c,d));
	if length(iu)==1,	% State space system, with iu but w/o freq. vector
		w=dfrqint(a,b,c,d,Ts,30);
	else			% Assume this is the old syntax DBODE(A,B,C,D,IU,W)
		disp('Warning: You are using the old syntax.  Use DBODE(A,B,C,D,Ts,IU,W) instead.')
		w=iu; iu=Ts; Ts=1;
	end
	[ny,nu] = size(d);

else
	error(abcdchk(a,b,c,d));
	[ny,nu] = size(d);

end

if nu*ny==0, phase=[]; w=[]; if nargout~=0, magout=[]; end, return, end

% Compute frequency response
if (nargs==3)|(nargs==4)
	g=freqresp(num,den,exp(sqrt(-1)*w*Ts));
else
	g=freqresp(a,b,c,d,iu,exp(sqrt(-1)*w*Ts));
end
mag = abs(g);
phase = (180/pi)*unwrap(atan2(imag(g),real(g)));

% Uncomment the following statement if you don't want the phase to  
% be unwrapped.  Note that phase unwrapping may not always work; it
% is only a "guess" as to whether +-360 should be added to the phase
% to make it more aestheticly pleasing.
% phase = (180/pi)*imag(log(g));

% If no left hand arguments then plot graph.
if nargout==0
	holdon = ishold;
	if ~holdon
		newplot;
	end
	subplot(211)
	if holdon
		hold on
	end

	% Magnitude plot with 0db line
	semilogx(w,20*log10(mag),[min(w(:)),max(w(:))],[0 0],'w:')
	grid
	xlabel('Frequency (rad/sec)')
	ylabel('Gain dB')

	subplot(212)
	semilogx(w,phase)

	xlabel('Frequency (rad/sec)')
	ylabel('Phase deg')

	subplot(212)
	semilogx(w,phase)


	% Set tick marks up to be in multiples of 30, 90, 180, 360 ... degrees.
	ytick = get(gca, 'ytick');
	ylim = get(gca, 'ylim');
	yrange = ylim(2) - ylim(1);
	n = round(log(yrange/(length(ytick)*90))/log(2));

	set(gca, 'ylimmode', 'manual')

	if n >=0
		ytick = [-90*2^n:-(90*2^n):ylim(1), 0:(90*2^n):ylim(2)];
		set(gca,'ytick',ytick);
	elseif n >= -2
		ytick = [-30:-30:ylim(1), 0:30:ylim(2)];
		set(gca,'ytick',ytick);
	end

	grid
	xlabel('Frequency (rad/sec)'), ylabel('Phase deg')
	subplot(111)
	return % Suppress output
end

% Uncomment the following line for decibels, but be warned that the
% MARGIN function will not work with decibel data.
% mag = 20*log10(mag);

magout = mag;