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

What I try to implement is Discrete Fourier Transform: 2 Dimensional DFT

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 my own function

The transform result of MATLAB fft2() function: enter image description here

mehmet
  • 1,631
  • 16
  • 21
  • 1
    For the "black image" and the darker results, try to add the square brackets as in `imshow(image,[])`. This will scale the result accordingly. That might do it. – kkuilla Dec 05 '14 at 11:55
  • @kkuilla I have tried it but unfortunately it hasn't worked :/ – mehmet Dec 05 '14 at 12:02
  • 1
    you have a scaling problem indeed. Try doing one divided to the other, it may give you some insight of what is there missing. It looks like the fft itself is rigth, but the scale is worng – Ander Biguri Dec 05 '14 at 12:56
  • @AnderBiguri Yess, you are right! The scale difference is 255. This is the max value of color (for uint8). What is that supposed to mean? Should I always multiply the value 255 with the matrix while I am using im2double() function? – mehmet Dec 05 '14 at 13:26
  • 1
    Usually an image is described from 0-255 in uint8 or from 0-1 in double. If you want to conver to double, try doing `double()` instead of `im2double()` – Ander Biguri Dec 05 '14 at 13:46

2 Answers2

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().

Ander Biguri
  • 35,140
  • 11
  • 74
  • 120
0

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.

arbiterr
  • 1