Writing Own fft2() Function on MATLAB

Question

• I want to write my own 2 Dimensional DFT function with reduced loops.

What I try to implement is Discrete Fourier Transform:

Using the separability property of transform (actually exponential function), we can write this as multiplication of two 1 dimensional DFT. Then, we can calculate the exponential terms for rows (the matrix wM below) and columns (the matrix wN below) of transform. Then, for summation process we can multiply them as "F = wM * original_matrix * wN"

Here is the code I wrote:

f = imread('cameraman.tif');

[M, N, ~] = size(f);
wM        = zeros(M, M);
wN        = zeros(N, N);

for u = 0 : (M - 1)
    for x = 0 : (M - 1)
        wM(u+1, x+1) = exp(-2 * pi * 1i / M * x * u);
    end    
end

for v = 0 : (N - 1)
    for y = 0 : (N - 1)
        wN(y+1, v+1) = exp(-2 * pi * 1i / N * y * v);
    end    
end

F = wM * im2double(f) * wN;

The first thing is I dont want to use 2 loops which are MxM and NxN times running. If I used a huge matrix (or image), that would be a problem. Is there any chance to make this code faster (for example eliminating the loops)?

The second thing is displaying the Fourier Transform result. I use the codes below to display the transform:

% // "log" method
fl = log(1 + abs(F));
fm = max(fl(:));
imshow(im2uint8(fl / fm))

and

% // "abs" method
fa = abs(F);
fm = max(fa(:));
imshow(fa / fm)

When I use the "abs" method, I see only black figure, nothing else. What is wrong with "abs" method you think?

And the last thing is when I compare the transform result of my own function with MATLAB' s fft2() function', mine displays darker figure than MATLAB' s result. What am I missing here? Implementation misktake?

The transform result of my own function:

The transform result of MATLAB fft2() function:

Answer 1

I am happy you solved your problem but unfortunately you answer is not completely right. Indeed it does the job, but as I commented, im2double will normalize everything to 1, therefore showing the scaled result you have. What you want (if you are looking for performance) is not doing im2doubleand then multiply by 255, but directly casting to double().

Answer 2

You can eliminate loops by using meshgrid. For example:

M = 1024;
tic
[ mX, mY ] = meshgrid( 0 : M - 1, 0 : M - 1 );
wM1 = exp( -2 * pi * 1i / M .* mX .* mY );
toc
tic
for u = 0 : (M - 1)
    for x = 0 : (M - 1)
        wM2( u + 1, x + 1 ) = exp( -2 * pi * 1i / M * x * u );
    end    
end
toc
all( wM1( : ) == wM2( : ) )

The timing on my system was: Elapsed time is 0.130923 seconds. Elapsed time is 0.493163 seconds.