izt.m

很多matlab的源代码

function ilt = izt(ty1,u,v,k,ty,tol)
% IZT Inverse z-transform in SYMBOLIC FORM.
%
%	I = IZT('tf',NN,DD,TY,TOL) Inverse z-transform of H(z) = N(z)/D(z).
%	If TY=0, NN,DD are the numerator and denominator of H(z).
%	If TY=1, DEN contain roots of D(z) and leading coeff of D(z) = 1.
%	TOL = tolerance for repeated roots [Default: TOL = 0.002]
%
%	I = IZT('pf',R,P,K,TY,TOL) IZT from PF expansion of H(z)/z.
%	R, P, and K are provided by the PF expansion M-file TF2PF.M
%	TY = 0 for this form
%
%	I is the SYMBOLIC inverse as a string function of n.
%	WARNING: IZT may be unstable, especially for repeated poles.
%
%	USAGE: 	To find the IZT of H(z)=(z-1)/(z*z-z+.5), use 
%		y = izt('tf',[1 -1],[1 -1 .5])
%
%     	IZT (with no input arguments) invokes the following example:
%
%	% Find the IZT of H(z)=(z*z-1)/[z(z-.5)(z*z+.25)]
%	% Since the den roots are known, we use TY=1 and the denominator roots
%	  >>hinv = izt('tf',[1 0 -1],[0;0.5;-j*0.5;j*0.5],1)


% 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 izt,disp('Strike a key to see results of the example')
pause,j=sqrt(-1);yin=izt('tf',[1 0 -1],[0;0.5;-j*0.5;j*0.5],1),return,end


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,lu=length(u);lv=length(v);
if u(lu)==0,u(lu)=[];lu=lu-1;else,v=[v 0];lv=lv+1;end  %Form H(z)/z
[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


[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=[]; %%% 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; %Separate into complex and real
lc=length(pc);lr=length(pr);ilt=[];

i=find(pr==0);if ~isempty(i),li=length(i);lk=lk+li;
k=[k rr(i).'];ik=[ik 1-ir(i)];pr(i)=[];rr(i)=[];ir(i)=[];lr=lr-li;end

for j=1:lr  %Now for the real terms
if abs(rr(j))>50*eps
p=num2str(pr(j));if eval(p)<0,p=['(' p ')'];else,p=[p ' '];end
m=ir(j);r=num2str(rr(j));
ex=['*(' p '.^n)'];
if eval(p)==1,ex=[];end
if rr(j)<0,a1='+';add=[];else,a1=[];add='+';end
if isempty(ilt),a1=[];add=[];end
if m==1,lt=[add r ex];ilt=[ilt lt];end
if m>2,fr=1;for j=1:m-1,fr=fr*j;end
   rr1=['(' r '/' num2str(fr) ')'];
   rre=eval(rr1);r=[a1 '(' num2str(rre) ')'];
   %if rre>0,r=[a1 r];end   
end
if m>1,
pow=['*n'];for i=1:m-2,pow=[pow '.*(n-' num2str(i) ')'];end
if eval(p)~=1,pow=[pow '.*(' p '.^(n-' num2str(m-1) '))'];end
ilt=[ilt add r pow];end
end,end
for j=1:lc  %Now for the complex terms
if abs(rc(j))>50*eps
ex=[];pa=abs(pc(j));p=num2str(pa);m=ic(j);thp=angle(pc(j))/pi;
ra=2*abs(rc(j));thr=angle(rc(j))/pi;
rm=num2str(ra);
a1='+';if isempty(ilt),a1=[];end
if thr>0,th=['+' num2str(thr) '*pi)'];
elseif thr<0,th=[num2str(thr) '*pi)'];else,th=[')'];end
aa=num2str(thp);cs=['.*cos(' aa '*n*pi'];
if eval(p) ~= 1,ex=['*(' p ' .^n)'];
else,if m==1,cs=['*cos(' aa '*n*pi'];end,end
if m==1,lt=[a1 rm ex cs th];ilt=[ilt lt];end
if m>2,fr=1;for j=1:m-1,fr=fr*j;end
   rm=['(' num2str(ra) '/' num2str(fr) ')'];
rm=eval(rm);rm=num2str(rm);rm=['(' rm ')'];   
end
if m>1,
mm=int2str(m-1);
cs=['.*cos(' aa '*(n-' mm ')*pi']; %%NEWNEWNEW
pow=['*n'];for i=1:m-2,pow=[pow '.*(n-' int2str(i) ')'];end
if eval(p) ~= 1,pow=[pow '.*(' p ' .^(n-' mm '))'];end
lt=[a1 rm pow cs th];ilt=[ilt lt];end
end,end
%if lk>0,if ~isempty(ilt),ilt=['(' ilt ').*ustep(n)'];end,end %%%NEW
if ~isempty(ilt),ilt=['(' ilt ').*ustep(n)'];end
for j=1:lk  %Now for the direct terms
if abs(k(j))>50*eps
if k(j)>=0,add='+';else,add='';end,if isempty(ilt),add='';end
nk=num2str(k(j));nn=[];if ik(j)>0,nn=['+' num2str(ik(j))];end
if ik(j)<0,nn=num2str(ik(j));end
del=['*udelta(n' nn ')'];ilt=[ilt add nk del];
end,end