Inverse_2D_Laplace_z_transform.m

% The program can get spatial-time response of 2-D Continuous-Discrete systems by taking inverse 2-D Laplace-z transform [1]. 
% The detailed algorithm is provided in Ref. [1].
% Copyright (C) Yang XIAO, BJTU, July 28, 2007, E-Mail: [e-mail]. 
% Based on recent results for 2-D continuous-discrete systems, Ref. [1] develops 2-D Laplace-z transform, which can be used 
% to analyze 2-D continuous-discrete signals and system in Laplace-z hybrid domain. Current 1-D Laplace transformation and
% z transform can be combined into the new 2-D s-z transform. However, 2-D s-z transformation is not a simple extension of
% 1-D transform, in 2-D case, we need consider the 2-D boundary conditions which dont occur in 1-D case. The hybrid 2-D 
% definitions and theorems are given in Ref. [1]. This program is derived from the numerical inverse 2-D Laplace-z transform,
% it shows the 2-D pulse response of a stable  2-D  continuous-discrete system.
% Ref: 
% [1] Y. Xiao and M. CH. Lee, 2-D Laplace-Z Transformation, IEICE TRANS. FUNDAMENTALS, VOL. E89-A, No. 5, May, 2006, pp.1500-1504.
% The paper [1] can be downloaded from the Web Site of IEICE.

clear;
a=1e-6;
T=1e-3;
M=50;
N=50;
h=2*T/N;
j=sqrt(-1);
I=N;

for i=1:I
    t=i*h;
    i
    for n=1:N
       % inverse Laplace transform 
       xL=0;
          for k=1:N
          s(k)=a+j*pi*(k-1)/N;
         % s(k)=a+2*j*pi*(k-1)*t;
       
       % inverse z transform 
          xz=0;
          xa=0;
              for m=1:M
              z(m)=j*2*pi*(m-1)/M;
          %  Example 1
              x_sz=z(m)*s(k)/(s(k)^2+.25)/(z(m)-.5);
              x_az=z(m)*a/(a^2+.25)/(z(m)-.5);
         %  Example 3
           %   x_sz=z(m)/(s(k)+1)/(z(m)-.5);
           %   x_az=z(m)/(a+1)/(z(m)-.5);
              %  x_sz=1/(1+.5*s(k)+.4*z(m));
         %  x_az=1/(1+.5*a+.4*z(m));
          %  Example 2
          % x_sz=z(m)^2/((12+s(k))*z(m)^2+(10+5*s(k))*z(m)+2+s(k));
          % x_az=1/((12+a)*z(m)^2+(10+5*a)*z(m)+2+a);
          %------------
          % x_sz=1/(s(k)+1)/(z(m)-.9);
            % x_az=1/(a+1)/(z(m)-.9);
              xz=xz+x_sz*exp(j*2*pi*(m-1)*(n-1)/M);
              xa=xa+x_az*exp(j*2*pi*(m-1)*(n-1)/M);
              end
       %  end of inverse z transform 
          %  xL=xL+xz*exp(j*pi*(k-1)*i/N);
           xL=xL+xz*exp(j*pi*(k-1)*t/T);
       end 
       %  end of Laplace transform  
       %x(i,n)=real(xL-xa/2)*exp(-a*t)/T;
       x(i,n)=real(xL-xa/2)*exp(a*t)/M/N;
   end
end
mesh(x)
ylabel('time:t*1e-6');
xlabel('digital value: n');
zlabel('h(t,n)');
title('Inverse L-z Transform of given 2-D function H(s,z)');