dtimvec.m

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

function npts=dtimvec(a,b,c,x0,st)
%DTIMVEC Returns a time vector for use with discrete time response
%        functions.
%
%	NPTS=DTIMVEC(a,b,c,X0,ST,PRECISION) returns a suggestion  
%	for the number of samples, NPTS. that should be used 
%	with discrete state space systems (a,b,c) that have known 
%	initial conditions X0. NPTS can be used,
%	for instance, for step and impulse responses. 
%	
%	DTIMVEC attempts to produce a number of points that ensures the 
%	output responses have decayed to approx. ST % of their initial 
%	or peak values. 

%	A.C.W.Grace 9-21-89
%	Copyright (c) 1986-93 by the MathWorks, Inc.

 	[r,n]=size(c);
% Work out scale factors by looking at d.c. gain, 
% eigenvectors and eigenvalues of the system.
	[m,r]=eig(a);
	r=diag(r);
% Equate the  responses to the sum of a set of 
% exponentials(r) multiplied by a vector of magnitude terms(vec).
	[r1,n]=size(c);
	if r1>1, c1=max(abs(c)); else c1=c; end
% Cater for the case when m is singular
	if rcond(m)<eps 
		vec = (c1*m).'.*(1e4*ones(n,1));
	else
		vec=(c1*m).'.*(m\x0);
	end
% Cater for the case when poles and zeros cancel each other out:
	vec=vec+1e-20*(vec==0);
% d.c. gain (not including d matrix)
	dcgain =c1*x0; % or dcgain=c\(eye(n)-a)*b;
	ind=find(imag(r)>0);
% If the d.c gain is small then base the scaling
% factor on the estimated maximum peak in the response.
	if abs(dcgain)<1e-8;
% Estimation of the maximum peak in the response.
		peak=2*imag(vec(ind)).*exp(-abs(0.5*real(r(ind)).*pi./imag(r(ind))));
		pk=max([abs(dcgain);abs(peak);eps]);
	else
		pk=dcgain;
	end
% Work out the st% settling time for the responses.
% (or st% peak value settling time for low d.c. gain systems).
	lt=(st*pk./abs(vec));
	n1=log(lt)./log(abs(r)+(abs(r)==1));
	n=max(real(n1));
% For unstable problems n will be negative
	if n<=3, n=max(abs(n1))+2; end
	n=min([n,1000]);
% Round the maximum time to an appropriate value for plotting.
	nn=chop(n,1,5);
	if abs((n-nn)/n)>0.2, 
		nn=chop(n,1,1);
		if abs((n-nn)/n)>0.2
			nn=chop(n,2,2);
		end
	end
	npts=nn+1;