function xdaub2 = dwtdaub2(x2, J, Hord)
% Daubechies wavelet transform in two dimensions
%  xdaub2 = dwtdaub2(x2, J, Hord)
%
% Inputs:
%  x2      2-dimensional signal (array)
%  J       number of stages (resolution levels)
%  Hord    positive integer. Order of Daubechies filter is 2*Hord-1. 
% Output:
%  xdaub2  Daubechies wavelet transform of x2
%

sz = size(x2);
nrows = sz(1); % number of rows in x2
ncols = sz(2); % number of columns in x2

% check whether nrows and ncols are divisable with 2^J
if ( rem(nrows, 2^J) ~= 0 ) || ( rem(ncols, 2^J) ~= 0 )
    error('The dimensions of x2 must be multiples of 2^J');
end

% compute Daubechies filter coefficients
H = hdaub(Hord);

xdaub1 = zeros(nrows,ncols);
xdaub2 = zeros(nrows,ncols);

% Transform first row-wise, then column-wise
for i=1:nrows
    xdaub1(i,:) = dwtdaub( x2(i,:), J, H );
end

for j=1:ncols
    xdaub2(:,j) = dwtdaub( xdaub1(:,j), J, H );
end