ilt.m

很多matlab的源代码

function [iltr,ilt0] = ilt(ty1,u,v,k,ty,tol)
% ILT Inverse Laplace transform in SYMBOLIC FORM.
%
%	I = ILT('tf',NN,DD,TY,TOL) Finds the inverse LT from H(s) = N(s)/D(s). 
%	If TY=0, NN, DD are the num and den of H(s) in descending powers of s
%	If TY = 1, DD contains roots of D(s) and leading coeff of D(s) =1.
%	TOL  =  tolerance for repeated roots [Default: 0.002].
%
%     	I = ILT('pf',R,P,K,TY,TOL) Finds ILT from PF expansion of H(s) 
%	R, P are residues and poles and K are the direct terms in the form 
%	K = [...s^2, s^1, s^0] (DEFAULT K = []). (See TF2PF)
%
%	I is returned in SYMBOLIC form as a string function of t. Impulses are 
%       included as udelta(t,0) or udelta(t,n) (nth derivative) (see UDELTA)
%     	WARNING: ILT may be unstable, especially for repeated poles.
%
%	USAGE: To find ILT of (s+2)/(s*s+4s+3), use yt=ilt('tf',[1 2],[1 4 3]) 
%
%       ILT (with no input arguments) invokes the following example:
%
%     	% Find and plot the ILT of H(s)=(s+2)/(s*s+4)(s+1)(s). 
%     	  >>yit = ilt('tf',[1 2],[-j*2;j*2;-1;0],1) %Uses den roots and TY=1 
%	  >>t=0:0.1:10;yi=eval(yit);plot(t,yi),grid


% 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 ilt,disp('Strike a key to see results of last example')
pause,j=sqrt(-1);yit=ilt('tf',[1 2],[-j*2;j*2;-1;0],1),
t=0:0.1:10;yi=eval(yit);plot(t,yi),grid,pause,return,end

rty=0;if length(ty1)>2,rty=1;ty1(1)=[];end
%EXAMPLE: if ty1='rtf', result is in rational fractions.
lv=[]; % Initialize
if ty1=='tf',
if nargin==5,tol=ty;ty=k;end
if nargin==4,ty=k;tol=0.002;end,
if nargin==3,k=[];ty=0;tol=0.002;end
if ty==1,
if length(v)==1,v=[1 -v];ty=0;elseif length(v)==0,v=1;ty=0;end;
end

%%%%%%%%%%%%NEW SEP 17,2003 Make sure it works
addik=0;
if ty==0,
   lu=length(u);lv=length(v);
   while u(lu)==0 & v(lv)==0,
      u(lu)=[];v(lv)=[];lu=lu-1;lv=lv-1;
   end
   if lv==1,
      while u(lu)==0
         u(lu)=[];lu=lu-1;addik=addik+1;
      end
   end
end

v=v(:).';
%if length(v)<2,u=[u 0];v=[v 0];end %%OLD
if length(v)<2,v=[v 0];end %%NEW

[u,v,k]=tf2pf(u,v,ty,tol);
else
if nargin<6,tol=0.002;end,if nargin<5,ty=0;end,if nargin<4,k=[];end
end


%lk=length(k);ik=[];
%if lk>0,ik=[lk:-1:1];end  %%OLD
%if addik>0 & ~isempty(u),
%   k=[k u];u=0;%lk=lk+1;
%   ik=addik+[lk:-1:0];
%end  %%NEW

if lv<2 & ~isempty(u),
   k=[k u];u=0;
end


[p,i]=sort(-imag(v));p=v(i);r=u(i);pp=p;rr=r;
i=find(imag(p)<0);if ~isempty(i),pp(i)=[];rr(i)=[];end
n=length(pp);qq=ones(1,n);
pc=[];pr=[];%%% New for v5

for j=2:n,if abs(pp(j)-pp(j-1))<tol,qq(j)=qq(j-1)+1;end,end
i=find(imag(pp)>0);if ~isempty(i),pc=pp(i);ic=qq(i);rc=rr(i);
pp(i)=[];qq(i)=[];rr(i)=[];end,pr=pp;ir=qq; 
lc=length(pc);lr=length(pr);iltr=[];

for j=1:lr  %Now for the real terms
if abs(rr(j))>60*eps
ex=[];if rty==0,p=num2str(pr(j));else,
[ee,ff]=rat(pr(j));ee=int2str(ee);ff=int2str(ff);
if ff=='1',p=ee;else,p=['(' ee '/' ff ')'];end,end

m=ir(j);if rty==0,r=num2str(abs(rr(j)));else,
[ee,ff]=rat(abs(rr(j)));ee=int2str(ee);ff=int2str(ff);
if ff=='1',r=ee;else,r=['(' ee '/' ff ')'];end,end
if eval(p)~=0,ex=['*exp(' p '*t)'];end
if eval(p)==1,ex=['*exp(t)'];end,if eval(p)==-1,ex=['*exp(-t)'];end 
if rr(j)<0,add='-';else,add='+';end,
r1=r;
if isempty(iltr),if rr(j)>=0,add=[];end,end
if eval(p)==0,ex=[];end 
if m==1,lt=[add r1 ex];if eval(r)==0,lt=[];end,iltr=[iltr lt];end
if m>2,fr=1;for i=1:m-1,fr=fr*i;end
   r1=['(' r '/' num2str(fr) ')'];
if rty==0,rr1=eval(r1);r1=num2str(rr1);r1=['(' r1 ')'];end %NEW Sep 18,2003   
end
if m>1,pow=['*(t.^' num2str(m-1) ')'];
if m==2,pow='*t';end 
if ~isempty(ex),pow=[pow '.'];end
lt=[add r1 pow ex];if eval(r)==0,lt=[];end
iltr=[iltr lt];end
end,end

for j=1:lc  %Now for the complex terms
if abs(rc(j))>50*eps
pa=real(pc(j));if rty==0,p=num2str(pa);else,
[ee,ff]=rat(pa);ee=int2str(ee);ff=int2str(ff);
if ff=='1',p=ee;else,p=['(' ee '/' ff ')'];end,end
m=ic(j);ra=2*abs(rc(j));thc=angle(rc(j))/pi;if abs(thc)<50*eps,thc=0;end
if rty==0,rm=[num2str(ra)];else,
[ee,ff]=rat(ra);ee=int2str(ee);ff=int2str(ff);
if ff=='1',rm=ee;else,rm=['(' ee '/' ff ')'];end,end
if isempty(iltr),if j==1,
if rty==0,rm=num2str(ra);else,
[ee,ff]=rat(ra);ee=int2str(ee);ff=int2str(ff);
if ff=='1',rm=ee;else,rm=['(' ee '/' ff ')'];end,end
end,end

if thc>0,su='+';end,if thc<0,su='-';end
if thc==0,th=')';else
if rty==0,th=[su num2str(abs(thc)) '*pi)'];
else,[ee,ff]=rat(thc);ee=int2str(abs(ee));ff=int2str(ff);
if ff=='1',th=[su ee '*pi)'];else,th=[su '(' ee '/' ff ')*pi)'];end
end,end

pca=imag(pc(j));if rty==0,aa=num2str(pca);else,
[ee,ff]=rat(pca);ee=int2str(ee);ff=int2str(ff);
if ff=='1',aa=ee;else,aa=['(' ee '/' ff ')'];end,end
cs1=['cos('];cs2=[aa '*t'];
if eval(p)~=0,ex=['exp(' p '*t).*'];else, ex=[];end
if eval(p)==1,ex=['exp(t).*'];end,if eval(p)==-1,ex=['exp(-t).*'];end 
if eval(aa)==1,cs2=['t'];end,if eval(aa)==-1,cs2=['-t'];end 
if isempty(iltr),a1=[];else,a1='+';end,tk=eval(num2str(thc));
if tk==1,th=[')'];a1='-';end,if tk==-1,th=[')'];a1='-';end
if tk==.5,a1='-';th=[')'];cs1=['sin('];end,if tk==-.5,th=[')'];cs1=['sin('];end
if tk==0,th=[')'];end,if eval(rm)==1,rm1=[];else,rm1=[rm '*'];end 
if m==1,lt=[a1 rm1 ex cs1 cs2 th];iltr=[iltr lt];end
if m>2,fr=1;for i=1:m-1,fr=fr*i;end,rm1=['(' rm '/' num2str(fr) ')'];
   if rty==0,rmm1=eval(rm1);rm1=['(' num2str(rmm1) ')'];end
   rm1=[rm1 '*'];%NEW Sep 18,2003   
end
if m>1,pow=['(t.^' num2str(m-1) ').*'];if m==2,pow='t.*';end
lt=[a1 rm1 pow ex cs1 cs2 th];iltr=[iltr lt];end
end,end
if ~isempty(iltr),iltr=['(' iltr ').*ustep(t)'];end

ilt0=iltr;if ~isempty(k)
l=length(k);k=k(l:-1:1);m=0:l-1;i=find(k==0);
if lv<2 & addik>0,m=addik+m;end

if ~isempty(i),k(i)=[];m(i)=[];l=length(k);end
for j=1:l,if abs(k(j))>50*eps,if k(j)>=0,add='+';else,add='-';end
if rty==0,nk=num2str(abs(k(j)));else,
[ee,ff]=rat(abs(k(j)));ee=int2str(ee);ff=int2str(ff);
if ff=='1',nk=ee;else,nk=['(' ee '/' ff ')'];end,end
if eval(nk)==1,nk=[];else,nk=[nk '*'];end,del=['udelta(t,' num2str(m(j)) ')'];
if isempty(iltr),if add=='+';add=[];end,end
iltr=[iltr add nk del];
end,end,end