1.10 二维离散拉普拉斯逆变换

二维离散拉普拉斯逆变换在解决描述线性动态系统动态行为的部分微分方程上有着很重要的应用。本节将给出两种实现二维离散拉普拉斯逆变换的算法，其中第一种方法是根据文献提出的二维离散拉普拉斯逆变换算法改进的基于二维FFT的；第二种方法是我们自行提出的，将一种我们改进的基于FFT的一维离散拉普拉斯逆变换拓展到二维实现二维离散拉普拉斯逆变换法。下面将分两部分分别讨论。

通过二维系统函数分析、描述线性系统特性，包括分析描述线性动态系统的瞬间线性行为特性正逐渐成为一个新的研究热点。而将拉普拉斯变换算法，特别是二维离散拉普拉斯逆变换实现就显得尤为重要。国外许多学者已经提出越来越多的各种实现二维离散拉普拉斯逆变换的算法，相应国内这方面工作还很少有人涉及。在描述非线性动态系统的瞬间动态行为上，特别是描述线性系统行为上，离散拉普拉斯变换，特别是二维离散拉普拉斯变换具有极其重要的作用。提出较为理想的算法，以进一步提高计算精度，对二维系统函数分析将具有一定实际意义。

拉普拉斯变换本身可以从傅里叶变换导出，它们的表示式相似，性质也有很多相似之处，拉普拉斯变换只是将后者拓展到复频域。因此我们自然联想到，如果想实现离散拉普拉斯逆变换，可以从FFT算法中获得一些有益的启发。由此，文献[30]提出了两种算法：一种是基于二维FFT的二维离散拉普拉斯逆算法；另一种是通过基于一维FFT的一维拉普拉斯逆变换拓展到二维的离散拉普拉斯逆算法。

二维离散拉普拉斯逆变换的实现可以将基于FFT的一维离散拉普拉斯逆变换拓展到二维。该方法对二维离散拉普拉斯逆变换的定义是：对f（t1，t2）不同变量分别进行部分一维离散拉普拉斯逆变换。

1.10.1 基于一维FFT的一维离散拉普拉斯逆变换

首先设定关于t的实函数f（t），当t<0，f（t）=0。

定义1.14 f（t）的拉普拉斯变换定义为

式中，s是复频率变量。

复函数F（s）的拉普拉斯逆变换定义为

式中，a>0是一个正实数；j表示虚数单位。

依照文献[24-29]已经提出的算法，通过下列公式计算离散拉普拉斯逆变换：

fN（t）：0<t<2T

式中，t=hn，n=0，1，…，N-1；。

上述方法实际应用时，我们发现误差较大，逆变换效果不是很理想。为了获得更好的计算结果，参考文献[24，25]的改进办法，文献[30]对上述方法系数进行了调整，得到了如下公式：

下面分别给出运用该方法对正弦函数和阶跃函数的拉普拉斯变换进行拉普拉斯逆变换的处理效果（注：因为选取的点数越大，误差越小，选取无穷个点，误差趋向于零；为了减小误差，达到更好的效果，我们用选取点数的10倍进行离散拉普拉斯逆变换计算）。

例1.1 对正弦函数的拉普拉斯变换进行拉普拉斯逆变换实验结果如图1.1所示，选取了128个点画图，选用的正弦函数的拉普拉斯变换函数模型为

F（s）=1/（s2+1） （1.114）

可以看到，逆变换的正弦曲线光滑对称，计算结果与理论结果基本相符。

例1.2 对阶跃函数的拉普拉斯变换进行拉普拉斯逆变换实验结果如图1.2所示，选取了128个点画图。选用的阶跃函数的拉普拉斯变换函数模型为

F（s）=1/（s2+1.414s+1） （1.115）

可以看到，除了最开始有少量突变，后面基本上是平滑的阶跃曲线，结果较为理想。有了上述的工作铺垫，就可以将其拓展到下面的对二维信号的拉普拉斯逆变换处理。

1.10.2 基于一维FFT的二维离散拉普拉斯逆变换

通过前面提出的对于一维信号基于FFT的拉普拉斯逆变换，文献[24]把它拓展到二维信号处理。对二维信号的二维拉普拉斯变换的结果进行二维离散拉普拉斯逆变换恢复原信号。该方法的思想实际上是对系统输出的对二维信号f（t1，t2）的二维拉普拉斯变换F（n1，n2）沿两维分别进行一维离散拉普拉斯逆变换，以实现二维离散拉普拉斯逆变换，实现对二维信号F（n1，n2）的二维离散拉普拉斯逆变换，恢复成原信号f（t1，t2），并以此来实现对非线性系统特性的分析。

二维离散拉普拉斯逆变换的定义为：对f（t1，t2）不同变量分别进行部分一维离散拉普拉斯逆变换。

令二维信号f（t1，t2）沿两维分别进行部分拉普拉斯变换的结果为F（n1，n2）。依照前面提出的基于FFT的一维拉普拉斯逆变换算法，本节给出该变换算法的定义。

定义1.15 二维信号F（n1，n2）偏n1的部分一维离散拉普拉斯逆变换为

记为

定义1.16 二维信号F（n1，n2）偏n2的部分一维离散拉普拉斯逆变换为

记为

根据定义，通过分别沿两维进行一维离散拉普拉斯逆变换的方法实现二维离散拉普拉斯逆算法。由定义1.15和定义1.16，可以推出二维离散拉普拉斯逆变换定义如下。

定义1.17 二维离散拉普拉斯逆变换定义为

记为

式（1.120）中

式中，Re[*]表示取实部；N、M分别代表两维取的点数。

1.10.3 基于一维FFT的二维离散拉普拉斯逆变换实验效果

下面分别给出运用该方法对二维正弦函数、二维阶跃函数以及一维正弦函数一维阶跃函数的二维函数的拉普拉斯变换进行二维离散拉普拉斯逆变换的处理结果。因为选取的计算点数越多，误差越小，若选取无穷个计算点，虽误差趋向于零，但计算量也趋于无穷。为了减小计算误差，同时达到更好的逼近理论值的效果，每一维用选取点数的10倍进行离散拉普拉斯逆变换计算。

例1.3 对二维阶跃函数拉普拉斯变换进行二维离散拉普拉斯逆变换实验结果如图1.3所示，选取了个128×128点画图。我们选用的阶跃函数的拉普拉斯变换函数模型为

F（s1，s2）=1/（（s1+1.1）×（s2+0.1）） （1.122）

例1.4 对二维正弦函数拉普拉斯变换进行二维离散拉普拉斯逆变换实验结果如图1.4所示，选取了128×128个点画图。

选用的阶跃函数的拉普拉斯变换函数模型为

F（s1，s2）=1/（（s12+1.1）×（s22+1.1）） （1.123）

例1.5 对横向一维沿正弦曲线纵向一维沿阶跃曲线的二维函数的拉普拉斯变换进行二维离散拉普拉斯逆变换实验结果如图1.5所示，选取了个128×128点画图。

选用的阶跃函数的拉普拉斯变换函数模型为

F（s1，s2）=1/（s1×（s22+1.1）） （1.124）

通过上面的例子，可以看到文献[30]提出的算法实验效果理想。从而可以得出结论：通过上述改进的基于FFT的一维离散拉普拉斯逆变换可以拓展到二维实现二维离散拉普拉斯逆变换。

1.10.4 二维拉普拉斯逆变换程序实现

二维拉普拉斯逆变换算法有两种：算法1是直接使用二维FFT逆变换进行计算，算法2是使用两次一维FFT逆变换进行计算。附录1.2给出了算法2的二维拉普拉斯逆变换程序。

算法1：

（1）利用二维拉普拉斯变换公式，赋值10N×10M的二维函数F（s1，s2），其中s1=a+jπ（k1-1）/N，s2=a+jπ（k2-1）/M；

（2）对上面赋值的二维函数F（s1，s2）进行二维FFT逆变换；

（3）将上一步的结果乘以调节系数：；

（4）取f（t1，t2）中的前N×M个点画图。

算法2：

（1）先将沿第一维对10N×10M的二维函数F（s1，s2）的第一变量进行一维离散拉普拉斯逆变换，其中s1=a+jπ（k1-1）/N，得到F（t1，s2）；

（2）再沿第二维对二维函数F（t1，s2）的第二变量进行一维离散拉普拉斯逆变换，其中s2=a+jπ（k2-1）/M，得到f（t1，t2）；

（3）取f（t1，t2）中的前N×M个点画图。

附录1.1 二维多项式的Schur稳定性检验程序

% Schur Stability Test of 2-D Polynomials
% 2-D Polynomial： B（z1，z2）=[1 z1^（-1） z1^（-2）]*B*[1 z2^（-1） z2^（-2）]'
clear；
i=sqrt（-1）；
M=512；             % The number of discrete frequency points on unit circle
w=0；
for n=1：M
w=pi*n/M； z=exp（j*2*pi*（n-1）/M）； %
%---------------------------
% stable
%  a=[12 6 0；
%     10 5 0；
%      2 0 1]
%-----------------------------
%stable
%  B=[1 -1.2 .5；
%      -1.5 1.8 -.75；
%      .6 -.72 .2719]；
%---------------------
%unstable
B=[1 -1.2 .5；           % The coefficients of a 2-D polynomial
-1.5 1.8 -.75；
.6 -.72 .269]；
%---------------------
bb=B*[1 z^（-1） z^（-2）]'；
%aa=a*[1 z z^2 z^3]'
B0=bb（1）； B1=bb（2）； B2=bb（3）；
%A=abs（b0）-abs（b1）；
BZ=[B0 B1 B2]；         % The complex varible polynomial coefficients B0（z1），
B1（z1） and B2（z1） of the 2-D  polynomial
x（n）=max（abs（roots（BZ）））；
if x（n）>1
e=x（n）；
end
end
xmax=max（x）
plot（x）
fprintf（'Display the maximum absolute value of the roots of the given polynomial\n'）
if xmax>1
fprintf（'The 2-D polynomial is not Schur stable'）
else
   fprintf（'The 2-D polynomial is Schur stable'）
end
plot（x）
ylabel（'The maximum absolute value of the roots of the polynomial'）；
xlabel（'The discrete frequency points on unit circle'）；

附录1.2 二维拉普拉斯逆变换程序

clear；
a=1e-6；
T=1e-3；
M=64；
N=64；
h=2*T/N；
j=sqrt（-1）；
I=40；
for i=1：I
   t1=i*h；
   i
for n=1：N
   % inverse Laplace transform_s1
   xL=0；
      for k=1：N
      s1（k）=a+8*j*pi*（k-1）/N；
     % s（k）=a+2*j*pi*（k-1）*t；
   % inverse z transform
      x2=0；
      xa=0；
         for m=1：M
         s2（m）=a+8*j*2*pi*（m-1）/M；
      %  Example 1
        x_s2=s2（m）*s1（k）/（s1（k）^2+.25）/（s2（m）^2+1.1）；
        x_a2=s2（m）*a/（a^2+.25）/（s2（m）^2+1.1）；
      %  x_sz=z（m）/（z（m）+1.1）；
        %  x_az=0；
         %  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+6*s（k））*z（m）^2+（10+5*s（k））*z（m）+2+s（k））；
     %    x_az=1/（（12+6*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）；
           x2=x2+x_s2*exp（j*2*pi*（m-1）*（n-1）/M）；
           xa=xa+x_a2*exp（j*2*pi*（m-1）*（n-1）/M）；
           end
       %  end of inverse Laplace transform_s2
         %  xL=xL+xz*exp（j*pi*（k-1）*i/N）；
          xL=xL+x2*exp（j*pi*（k-1）*t1/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*t1）/M/N/T；
    end
end
mesh（x）
ylabel（'time：t1*1e-6'）；
xlabel（'digital value： n'）；
zlabel（'h（t，n）'）；
title（'Inverse 2-D L Transform of given 2-D function H（s1，s2）'）；