initial.m

经典通信系统仿真书籍《通信系统与 MATLAB (Proakis)》源代码

function [yout,x,t] = initial(a,b,c,d,x0,t)
%INITIAL Initial condition response of continuous-time linear systems.
%	INITIAL(A,B,C,D,X0) plots the time response of the linear system
%		.
%		x = Ax + Bu
%		y = Cx + Du
%
%	due to an initial condition on the states.  The time vector is 
%	automatically determined based on the system poles and zeros.  
%
%	INITIAL(A,B,C,D,X0,T) plots the initial condition response for the
%	times specified in the vector T.  The time vector must be 
%	regularly spaced.  When invoked with left hand arguments:
%	
%		[Y,X,T] = INITIAL(A,B,C,D,X0,...)
%
%	returns the output and state responses (Y and X), and the time 
%	vector (T).  No plot is drawn on the screen.  The matrix Y has as
%	many columns as outputs and one row for element in T.  Similarly,
%	the matrix X has as many columns as states and length(T) rows.
%	
%	See also: IMPULSE,STEP,LSIM, and DINITIAL.

%	Clay M. Thompson  7-6-90
%	Revised ACWG 6-21-92
%	Copyright (c) 1986-93 by the MathWorks, Inc.

if nargin==0, eval('exresp(''initial'',1)'), return, end

error(nargchk(4,6,nargin));
error(abcdchk(a,b,c,d));

[ny,nu] = size(d);
if (nu*ny==0)|isempty(a),
  x = []; t = []; if nargout~=0, yout=[]; end, return, 
end

[nx,na] = size(a);
if nargin==4		% No x0 specified, use random initial condition.
  x0 = randn(nx,1);

else
  x0 = x0(:);		% Make sure x0 is a column vector
  if length(x0)~=nx,
    error('The length of X0 must be equal to the number of states.');
  end
end

if (nargin==4)|(nargin==5)  % Workout time vector if not supplied.
  % The next two constants control the precision of the plot
  % and the time interval of the plot.
  st=0.005; % Set settling time bound  = 0.5%
  precision=30; % Show approx 30 points for simple graph
  t=timvec(a,b,c,x0,st,precision);
end

% --- Simulation ---
dt = t(2)-t(1);
[aa,bb] = c2d(a,b,t(2)-t(1));
n = length(t);
[nx,nu] = size(b);
x = ltitr(aa,bb,zeros(n,nu),x0(:));
y=x*c.';
tout = t;

if nargout==0,		% If no output arguments, plot graph
  plot(t, y, [t(1),t(length(t))],[0;0],'w:')
  xlabel('Time (secs)'), ylabel('Amplitude')
  return % Suppress output
end
yout = y;