Matlab: 2D Discrete Fourier Transform and Inverse

I'm trying to run a program in matlab to obtain the direct and inverse DFT for a grey scale image, but I'm not able to recover the original image after applying the inverse. I'm getting complex numbers as my inverse output. Is like i'm losing information. Any ideas on this? Here is my code:

%2D discrete Fourier transform
%Image Dimension

M=3;
N=3;
f=zeros(M,N);
f(2,1:3)=1;
f(3,1:3)=0.5;
f(1,2)=0.5;
f(3,2)=1;
f(2,2)=0;

figure;imshow(f,[0 1],'InitialMagnification','fit')


%Direct transform


for u=0:1:M-1
   for v=0:1:N-1
       for x=1:1:M
           for y=1:1:N

             F(u+1,v+1)=f(x,y)*exp(-2*pi*(1i)*((u*(x-1)/M)+(v*(y-1)/N)));


            end
        end
    end
end


Fab=abs(F);

figure;imshow(Fab,[0 1],'InitialMagnification','fit')



%Inverse Transform

for x=0:1:M-1
    for y=0:1:N-1
       for u=1:1:M
            for v=1:1:N

                z(x+1,y+1)=(1/M*N)*F(u,v)*exp(2*pi*(1i)*(((u-1)*x/M)+((v-1)*y/N)));


            end
        end
    end
end

figure;imshow(real(z),[0 1],'InitialMagnification','fit')

There are a couple of issues with your code:

1. You are not applying the definition of the DFT (or IDFT) correctly: you need to sum over the original variable(s) to obtain the transform. See the formula here; notice the sum.

2. In the IDFT the normalization constant should be 1/(M*N) (not 1/M*N).

Note also that the code could be made mucho more compact by vectorization, avoiding the loops; or just using the fft2 and ifft2 functions. I assume you want to compute it manually and "low-level" to verify the results.

The code, with the two corrections, is as follows. The modifications are marked with comments.

M=3;
N=3;
f=zeros(M,N);
f(2,1:3)=1;
f(3,1:3)=0.5;
f(1,2)=0.5;
f(3,2)=1;
f(2,2)=0;

figure;imshow(f,[0 1],'InitialMagnification','fit')

%Direct transform

F = zeros(M,N); % initiallize to 0
for u=0:1:M-1
   for v=0:1:N-1
       for x=1:1:M
           for y=1:1:N
               F(u+1,v+1) = F(u+1,v+1) + ...
                   f(x,y)*exp(-2*pi*(1i)*((u*(x-1)/M)+(v*(y-1)/N))); % add term
            end
        end
    end
end

Fab=abs(F);
figure;imshow(Fab,[0 1],'InitialMagnification','fit')

%Inverse Transform

z = zeros(M,N);
for x=0:1:M-1
    for y=0:1:N-1
       for u=1:1:M
            for v=1:1:N
                z(x+1,y+1) = z(x+1,y+1) + (1/(M*N)) * ... % corrected scale factor
                    F(u,v)*exp(2*pi*(1i)*(((u-1)*x/M)+((v-1)*y/N))); % add term
            end
        end
    end
end

figure;imshow(real(z),[0 1],'InitialMagnification','fit')

Now the original and recovered image differ only by very small values, of the order of eps, due to the usual floating-point inaccuacies:

>> f-z
ans =
   1.0e-15 *
  Columns 1 through 2
  0.180411241501588 + 0.666133814775094i -0.111022302462516 - 0.027755575615629i
  0.000000000000000 + 0.027755575615629i  0.277555756156289 + 0.212603775716506i
  0.000000000000000 - 0.194289029309402i  0.000000000000000 + 0.027755575615629i
  Column 3
 -0.194289029309402 - 0.027755575615629i
 -0.222044604925031 - 0.055511151231258i
  0.111022302462516 - 0.111022302462516i