fbode.m

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

function [magout,phase,w] = fbode(a,b,c,d,iu,w)
%FBODE  Bode frequency response for continuous-time linear systems.
%	FBODE(A,B,C,D,IU) produces a Bode plot from the single input IU to
%	all the outputs of the continuous 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.  The frequency range and 
%	number of points are chosen automatically.  FBODE is a faster, but
%	less accurate, version of BODE.
%
%	FBODE(NUM,DEN) produces the Bode 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. 
%
%	FBODE(A,B,C,D,IU,W) or FBODE(NUM,DEN,W) uses the user-supplied 
%	frequency vector W which must contain the frequencies, in 
%	radians/sec, at which the Bode response is to be evaluated.  See 
%	LOGSPACE to generate log. spaced frequency vectors.  When invoked
%	with left hand arguments,
%		[MAG,PHASE,W] = FBODE(A,B,C,D,...)
%		[MAG,PHASE,W] = FBODE(NUM,DEN,...) 
%	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 and BODE.

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


nargs = nargin;
if nargs==0, eval('exresp(''fbode'')'), return, end

error(nargchk(2,6,nargs));

% --- Determine which syntax is being used ---
if (nargs==1),
	error('Wrong number of input arguments.');

elseif (nargs==2),	% Transfer function form without frequency vector
	num = a; den = b;
	w = freqint(num,den,20);
	[ny,nn] = size(num); nu = 1;

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

elseif (nargs==4),	% State space system, w/o iu or frequency vector
	error(abcdchk(a,b,c,d));
	w = freqint(a,b,c,d,30);
	[iu,nargs,mag,phase]=mulresp('bode',a,b,c,d,w,nargout,1);
	if ~iu, if nargout, magout = mag; end, return, end
	[ny,nu] = size(d);

elseif (nargs==5),	% State space system, with iu but w/o freq. vector
	error(abcdchk(a,b,c,d));
	w = freqint(a,b,c,d,30);
	[ny,nu] = size(d);

else			% State space system, with iu and freq. vector
	error(abcdchk(a,b,c,d));
	[ny,nu] = size(d);
	
end

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

% --- Evaluate the frequency response ---
if (nargs == 2)|(nargs == 3)
	% It is in transfer function form.  Do directly, using Horner's method
	% of polynomial evaluation at the frequency points, for each row in
	% the numerator.  Then divide by the denominator.
	[ma,na] = size(num);
	s = sqrt(-1)*w(:);
	for i=1:ma
		g(:,i) = polyval(num(i,:),s);
	end
	g = polyval(den,s)*ones(1,ma).\g;
else
	[no,nu] = size(d);
	[ns,na] = size(a);
	nw = max(size(w));

	[p,a] = eig(a);	   % Reduce A to diagonal form
	if rcond(p)<sqrt(eps), 
		disp('Warning: Diagonalization matrix singular.  Use BODE instead.')
	end
	% Apply similarity transformations to B and C
	if ns>0,
		b = p \ b(:,iu);  
		c = c * p;
		d = d(:,iu);
		s = (w*sqrt(-1)).';
		s = s(ones(ns,1),:);
		a2 = diag(a)*ones(1,nw);
		b = b*ones(1,nw);
		g = (c*((1 ./(s-a2)).*b) + diag(d)*ones(no,nw)).';
	else
		d = d(:,iu);
		g = (diag(d)*ones(no,nw)).';
	end				
end

mag = abs(g);
phase = (180./pi)*unwrap(atan2(imag(g),real(g)));

% Uncomment out the following statement if you don't want the phase to  
% be unwrapped.  Note that phase unwrapping will not always work; it is
% only a "guess" as to whether +-360 should be added to the phase 
% to make it more aesthetically pleasing.  (See UNWRAP.M)

%phase = (180./pi)*atan2(imag(g),real(g));

%Try to correct phase anomaly for plants with integrators
%by adding multiples of -360 degrees.
if (nargs == 2 | nargs == 3) 
	nd = length(den);
	f = find(fliplr(den) == zeros(1,nd));
	nintgs = sum(f == 1:length(f));	
	if phase(1) > 0 & nintgs > 0
		phase = phase - 360;
	end
else 
	if abs(det(a)) < eps
		if phase(1) > 0 & nintgs > 0
		    phase = phase - 360;
		end
	end
end

% If no left hand arguments then plot graph.
if nargout==0
	holdon = ishold;
	subplot(211) 
	if holdon
		hold on
	end
	semilogx(w,20*log10(mag),w,zeros(1,length(w)),'w:')
	% If hold is set to on on the current axis then set it on the first axis too.
	% This enables two bode response to be superimposed on each other with
	% the following commands:
	%	bode(num, den); hold on; bode(num2, den2)
	grid
	xlabel('Frequency (rad/sec)'), ylabel('Gain dB')

	subplot(212), 
	semilogx(w,phase)
	xlabel('Frequency (rad/sec)'), ylabel('Phase deg')


	% 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
	% Reset the graph: subplot(111)
	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;