Inverze_2D_z.m

% The program can get the 2-D impulse response of a 2-D recursive discrete system in 2-D digital domain. 
% For a 2-D recursive discrete system, the inverse 2-D z transform can be implemented by a
% 2-D IIR filter[1-3]. 
% Copyright (C) Yang XIAO, Beijing Jiaotong University, Aug.1, 2007, E-Mail: [e-mail].
%
% Resort to the theorems in Ref. [1-3], we can get a stable 2-D recursive discrete system, with stable
% transfer function in 2-D z-domain. Since the transfer function % H(z1,z2)=A(z1,z2) is stable, 
% the 2-D impulse response  h(n1,n2) of the system equals the solution of a 2-D IIR filter to process a unit impulse
% in 2-D digital domain [3]. Then this program can get the 2-D impulse response of a given 2-D recursive
% discrete system in 2-D digital domain if the system is stable. 
% The stability test for the two examples is provided in Ref [1-3]. In the two examples, the program
% verifies the convergence of the 2-D impulse response, here, one is of stable denominator polynomial B(z1,z2), 
% and the other is of an unstable one. According to the theorems in [1-3], the stability of the systems
% is depended on that of the denominator polynomial, the program let us to observe the stability from
% the 2-D impulse responses of the two systems.
% Ref:
% [1] XIAO Yang, Unbehauen R. New stability test algorithm for two-dimensional digital filters,
% IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, Volume 45,  Issue 7,
% July 1998 Page(s):739 - 741
% [2] Yang XIAO, Unbehauen R. Schur stability of polytopes of bivariate polynomials, IEEE Transactions 
% on Circuits and Systems I: Fundamental Theory and Applications, Volume 49, Issue 7,  July 2002 Page(s):
% 1020 C 1023.
% [3] XIAO Yang, Stability Analysis of Multidimensional Systems, Shanghai Science and Technology Press,
% Shanghai, 2003.
%
% The papers [1] and [2] can be downloaded from Web site of IEEE Explore.
%
clear;
NP1=64
NP2=64
N1=3;
N2=3;
X=zeros(NP1,NP1);
Y=X;
X(N1,N2)=1;
%Y(1,1)=1;
%Example 1: Stable system
% B(z1,z2)=[1 z1^(-1) z1^(-2)]*[b00 b01 b02;b10 b11 b12;b20 b21 b22]*[1 z2^(-1) z2^(-2)]';
% A(z1,z2)=[1 z1^(-1) z1^(-2)]*[a00 a01 a02;a10 a11 a12;a20 a21 a22]*[1 z2^(-1) z2^(-2)]';

%  Contruct transfer function: H(z1,z2)=A(z1,z2)/B(z1,z2)
%---stable
% B=[1 -1.2 .5;
%  -1.5 1.8 -.72;
%  .6 -.75 .29]
% A=[1 -1.2 .5;
%  -1.5 1.8 -.72;
%  .6 -.75 .25]

%Example 2:unstable system
B=[1 -1.2 .5;
  -1.5 1.8 -.72;
  .6 -.75 .25]
A=[1 -1.2 .5;
 -1.5 1.8 -.72;
 .6 -.75 .29]

% Calculate h(n1,n2)=inver_2D_z(H(z1,z2))

for I=N1:NP1
for J=N2:NP2
%     Y0=0;
%     for m=1:N1-1
%         for n=1:N2-1
%         Y0=Y0+B(m+1,n+1)*Y(I-m,J-n);    
%          end
%     end
% Y(I,J)=X(I,J)-Y0;

Y(I,J)=A(1,1)*X(I,J)+A(1,2)*X(I,J-1)+A(1,3)*X(I,J-2)...
       +A(2,1)*X(I-1,J)+A(2,2)*X(I-1,J-1)+A(2,3)*X(I-1,J-2)...
       +A(3,1)*X(I-2,J)+A(3,2)*X(I-2,J-1)+A(3,3)*X(I-2,J-2)...    
                    -B(1,2)*Y(I,J-1)-B(1,3)*Y(I,J-2)...
                    -B(2,1)*Y(I-1,J)-B(2,2)*Y(I-1,J-1)-B(2,3)*Y(I-1,J-2)...
                    -B(3,1)*Y(I-2,J)-B(3,2)*Y(I-2,J-1)-B(3,3)*Y(I-2,J-2);    
end
end
% Display the 2-D impulse response of the given system
mesh(Y)
zlabel('2-D impulse response of the given system');
xlabel('x-direction:n1');
ylabel('y-direction:n1');