sysresp2.m

很多matlab的源代码

function [yt,yzs,yzi] = sysresp2(ty,tfn,tfd,inn,ind,ic)
% SYSRESP2 System time response.
%
%	[YT,YZS,YZI] = SYSRESP2(TY,NH,DH,NX,DX,IC) System time response
%	TY = 's' for CT systems, TY = 'z' for DT systems. 
%	NH, DH are num and den of H(s) or H(z) in descending powers.
%	NX, DX are the num and den of the LT of the input x(t) or ZT of x[n]   
%	For N inputs, NX and DX must have as many rows as inputs
%	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 = SYSRESP2(TY,NUM,DEN) computes the system IMPULSE RESPONSE.
%
%	NOTE: ALL OUTPUTS are in SYMBOLIC form (string functions of t or n)
%
%	SYSRESP2 (with no input arguments) invokes the following example:
%	% Find and plot the response of H(s)=(s+2)/(s*s+4s+3) to 2texp(-2t)u(t) 
%	% with y(0) = 3, y'(0) = -1. Note that X(s) =2/(s*s+4s+4)
%	  >>[ytot,yzsr,yzir]=sysresp2('s',[1 2],[1 4 3],2,[1 4 4],[3 -1])
%	  >>t=0: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 sysresp2,disp('Strike a key to see results of the example')
pause,[ytot,yzsr,yzir]=sysresp2('s',[1 2],[1 4 3],2,[1 4 4],[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=[];nzs=[];rd=[];

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<6,ic=0;end,
if ~isempty(inn)
[m,n]=size(inn);
emp=[];for i=1:m,if ~any(inn(i,:)),emp=[emp i];end,end
inn(emp,:)=[];ind(emp,:)=[];
end

if ~isempty(ind)
[m,n]=size(ind);
emp=[];for i=1:m,if ~any(ind(i,:)),emp=[emp i];end,end
inn(emp,:)=[];ind(emp,:)=[];
end

if any(inn)

tfin=inn(1,:);tfid=ind(1,:);
[m,n]=size(inn);
if m>1,
for j=2:m
t1=inn(j,:);t2=ind(j,:);
den=conv(tfid,t2);n1=conv(tfin,t2);n2=conv(tfid,t1);
l1=length(n1);l2=length(n2);ld=length(den);
n1=[zeros(1,ld-l1) n1];n2=[zeros(1,ld-l2) n2];
tfin=n1+n2;tfid=den;
end,
end
inn=tfin;ind=tfid;


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<6 | ~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 any(inn)
nic=conv(nzi,ind);
else
nic=0*nzi;
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