function [xM, xH, xV, xD] = haar2d(X)


N = size(X,1)/2;

xM = zeros(N,N);
xH = zeros(N,N);
xV = zeros(N,N);
xD = zeros(N,N);

for j=1:N
    for i=1:N
        
        I = 2*i-1;
        J = 2*j-1;
        
        xM(i,j) = 0.5*(X(I,J)+X(I,J+1)+X(I+1,J)+X(I+1,J+1));
        
        xH(i,j) = 0.5*(X(I,J)+X(I+1,J)-X(I+1,J)-X(I+1,J+1));
        
        xV(i,j) = 0.5*(X(I,J)+X(I,j+1)-X(I+1,J)-X(I+1,J+1));
        
        xD(i,j) = 0.5*(X(I,J)+X(I+1,J+1)-X(I+1,J)-X(I,J+1));
        
    end
end
