How can I write a Matlab function haar_dec.m for computing the J-level, J = 1, 2, . . . , log2N, Haar wavelet transform of an N x N image, where N is a power of 2? The function can repeatedly call the function haar.m (see below) which computes 1-level of the 2-D Haar wavelet transform. The function should take two arguments: an input image im and the number of levels J, and output an array of N x N wavelet coefficients C in the arrangement as illustrated in the figure attached.
function [w1,w2,w3,w4] = haar(x)
n = length(x);
for i=1:n/2
for j=1:n/2
i0 = 2*(i-1);
j0 = 2*(j-1);
w1(i,j) = 1/2*(x(i0+1,j0+1)+x(i0+1,j0+2)+x(i0+2,j0+1)+x(i0+2,j0+2));
w2(i,j) = 1/2*(x(i0+1,j0+1)+x(i0+1,j0+2)-x(i0+2,j0+1)-x(i0+2,j0+2));
w3(i,j) = 1/2*(x(i0+1,j0+1)-x(i0+1,j0+2)+x(i0+2,j0+1)-x(i0+2,j0+2));
w4(i,j) = 1/2*(x(i0+1,j0+1)-x(i0+1,j0+2)-x(i0+2,j0+1)+x(i0+2,j0+2));
end
end
function [w1,w2,w3,w4] = haar(x)
n = length(x);
for i=1:n/2
for j=1:n/2
i0 = 2*(i-1);
j0 = 2*(j-1);
w1(i,j) = 1/2*(x(i0+1,j0+1)+x(i0+1,j0+2)+x(i0+2,j0+1)+x(i0+2,j0+2));
w2(i,j) = 1/2*(x(i0+1,j0+1)+x(i0+1,j0+2)-x(i0+2,j0+1)-x(i0+2,j0+2));
w3(i,j) = 1/2*(x(i0+1,j0+1)-x(i0+1,j0+2)+x(i0+2,j0+1)-x(i0+2,j0+2));
w4(i,j) = 1/2*(x(i0+1,j0+1)-x(i0+1,j0+2)-x(i0+2,j0+1)+x(i0+2,j0+2));
end
end