Matlab: 2D Discrete Fourier Transform and Inverse

Question

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')

Answer 1

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  

Answer 2

Firstly, the biggest error is that you are computing the Fourier transform incorrectly. When computing F, you need to be summing over x and y, which you are not doing. Here's how to rectify that:

F = zeros(M, N);

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)));
            end
        end
    end
end

---

Secondly, in the inverse transform, your bracketing is incorrect. It should be 1/(M*N) not (1/M*N).

---

As an aside, at the cost of a bit more memory, you can speed up the computation by not nesting so many loops. Namely, when computing the FFT, do the following instead

x = (1:1:M)'; % x is a column vector
y = (1:1:N) ; % y is a row vector

for u = 0:1:M-1
    for v = 0:1:N-1
        F2(u+1,v+1) = sum(f .* exp(-2i * pi * (u*(x-1)/M + v*(y-1)/N)), 'all');
    end
end

To take this method to the extreme, i.e. not using any loops at all, you would do the following (though this is not recommended, since you would lose code readability and the memory cost would increase exponentially)

x = (1:1:M)'; % x is in dimension 1
y = (1:1:N) ; % y is in dimension 2
u = permute(0:1:M-1, [1, 3, 2]); % x-freqs in dimension 3
v = permute(0:1:N-1, [1, 4, 3, 2]); % y-freqs in dimension 4
% sum the exponential terms in x and y, which are in dimensions 1 and 2. 
% If you are using r2018a or older, the below summation should be 
%    sum(sum(..., 1), 2) 
% instead of 
%    sum(..., [1,2])
F3 = sum(f .* exp(-2i * pi * (u.*(x-1)/M + v.*(y-1)/N)), [1, 2]); 
% The resulting array F3 is 1 x 1 x M x N, to make it M x N, simply shiftdim or squeeze
F3 = squeeze(F3);