sysresp1.m

很多matlab的源代码

function [yt,yzs,yzi] = sysresp1(ty,tfn,tfd,x,ic)
% SYSRESP1 System time response.
%
%	[YT,YZS,YZI] = SYSRESP1(TY,N,D,X,IC) Computes system time response
%	TY = 's' for CT systems, TY = 'z' for DT systems. 
%	N and D are num and den of the operational transfer function 
%	OR numerator and dennominator of H=N/D in descending powers.
%	X is the TIME-DOMAIN INPUT as an array with the following form
%	CT: X = [K, a, p, w, r] for X(t) = Kexp(-at)(t^p)cos(wt+r)
%	DT: X=[K, a, p, w, r,n0] for X(n)=K[a^(n-n0)][(n-n0)^p]cos[w(n-n0)n+r] 
%	For N inputs, X is an Nx5 matrix for CT and Nx6 matrix for DT
%	IC = [y(0), y'(0), y"(0),... ] for CT systems (DEFAULT: IC = 0)
%	IC = {y[-1], y[-2],....,y[-N]} for DT systems (DEFAULT: IC = 0)
%	YT = total, YZS = zero-state response, and YZI = zero-input response.
%
%       YI = SYSRESP1(TY,NUM,DEN) computes the system IMPULSE RESPONSE.
%	NOTE: ALL OUTPUTS are in SYMBOLIC form (string functions of t or n)
%
%	SYSRESP1 (with no input arguments) invokes the following example:
%	% Find and plot the response of [s*s+4s+3]y=[s+2]x to 2*t*exp(-2t)u(t) 
%	% with y(0) = 3, y'(0) = -1.
%	  >>[ytot,yzsr,yzir]=sysresp1('s',[1 2],[1 4 3],[2 2 1 0 0],[3 -1])
%	  >>t=0:.1:3;>>ye=eval(ytot);plot(t,ye)


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


if nargin==0,help sysresp1,disp('Strike a key to see results of the example')
pause,[ytot,yzsr,yzir]=sysresp1('s',[1 2],[1 4 3],[2 2 1 0 0],[3 -1])
t=0:.1:3;ye=eval(ytot);
vx=matverch;
if vx < 4, eval('clg');else,eval('clf');end
plot(t,ye),grid,return,end

yt=[];yzs=[];yzi=[];ind=[];inn=[];nzs=[];rd=[]; %Initialize variables
while tfd(1)==0,tfd(1)=[];end,
while tfn(1)==0,tfn(1)=[];end
tfn=tfn/tfd(1);tfd=tfd/tfd(1);
ln=length(tfn);ld=length(tfd);

if nargin==3
if ty=='s',yt=ilt('tf',tfn,tfd);else,yt=izt('tf',tfn,tfd);end
yzs=yt;return,end

if nargin<5,ic=0;end

if ~isempty(x)
xxa=x(:,1);i=find(xxa==0);
if ~isempty(i),x(i,:)=[];end
end

if ~isempty(x)
%%% NEW ERROR CHECK
if ty~='s'
if any(x(:,2)==0)
error('Exponential base must be non-zero');
return
end
end


if ty=='s',[inn,ind]=ltr(x);else,[inn,ind]=ztr(x);end

%if ~isempty(inn)
while inn(1)==0,inn(1)=[];end,
while ind(1)==0,ind(1)=[];end
inn=inn/ind(1);ind=ind/ind(1);
if length(inn)>1,nzs=conv(tfn,inn);else,nzs=tfn*inn;end
rd=roots(tfd);
if length(ind)>1,rd=[rd;roots(ind)];end
if ty=='s',yzs=ilt('tf',nzs,rd,1);else,
yzs=izt('tf',nzs,rd,1);
%if ln<=ld,yzs=['(' yzs ').*ustep(n)'];end
end
else
yzs=[];
end

if nargin<5 | ~any(ic),yt=yzs;return,end
l=length(tfd);ic=[ic zeros(1,l-1-length(ic))];

nzi=zir(tfd,ic,ty);l1=length(nzs);
if ~isempty(ind)
nic=conv(nzi,ind);
else
nic=0*nzi;
%nic=conv(nzi,0);
end
l2=length(nic);
m=max(l1,l2);nt=[zeros(1,m-l1) nzs]+[zeros(1,m-l2) nic];
%if ty=='s',yt=ilt('tf',nt,rd,1);else,yt=izt('tf',nt,rd,1);end  %see below
%if nargout<3,return,end   %%see below
if ty=='s',yzi=ilt('tf',nzi,tfd);
else,
yzi=izt('tf',nzi,tfd);
%yzi=['(' yzi ').*ustep(n+' int2str(length(ic)) ')'];
end
if isempty(rd),yt=yzi;return,end

if ty=='s',yt=ilt('tf',nt,rd,1);else
if ln>ld
yzi=['(' yzi ').*ustep(n+' int2str(length(ic)) ')'];
yt=[yzi '+ (' yzs ')'];
else
yt=izt('tf',nt,rd,1);
%yt=['(' yt ').*ustep(n+' int2str(length(ic)) ')'];
end
end