Interpretting a 2D fourier transform

Hello all

Sorry for basic question, my understanding of fourier transforms is apparently quite limited.

I've got an image with some spatial oscillations hidden in it, and I'm trying to use a 2D fourier transform to find components of their wavevectors (in both the x and y directions).

The image is 80 pixels high and 100 wide, and for now I'm just assuming that 1 pixel = 1 unit length, so the spatial sampling frequency is 1.

I've done a 2D fourier transform of the image, but I can't figure out how to work out the spatial frequencies of the oscillations from the resulting plot.

The code below is a minimal working example, which produces the image and the 2D FT.

vidHeight = 80; 
vidWidth = 100; 

% These are the wavevectors I want to find
k1 = 0.1;
k2 = 0.6; 
k3 = -0.8;
k4 = 0.7;

generatedOscillation = zeros(vidHeight, vidWidth);
for x = 1:vidWidth
    for y = 1:vidHeight
        % here i add a zero frequency component, which I get rid of by
        % subtracting the average value of the image
        generatedOscillation(y,x) = generatedOscillation(y,x) + 1.2;

          % add one oscillation
          generatedOscillation(y,x) = generatedOscillation(y,x) + sin(k1*x + k2*y);

          % add another oscillation
          generatedOscillation(y,x) = generatedOscillation(y,x) + sin(k3*x + k4*y);
      end
  end

% show resulting image
imagesc(generatedOscillation)

% dft 2d
NFFTY = 2^nextpow2(vidHeight);
NFFTX = 2^nextpow2(vidWidth);

% 'detrend' data to eliminate zero frequency component
av = sum(generatedOscillation(:)) / length(generatedOscillation(:));
generatedOscillation = generatedOscillation - av;

% this section I'm not too sure about, I pretty much copied the example
% from the matlab 1D fft and adapted it slightly for 2D
samplingFreq = 1;
spatialFreqsX = samplingFreq/2*linspace(0,1,NFFTX/2+1);
spatialFreqsY = samplingFreq/2*linspace(0,1,NFFTY/2+1);

spectrum2D = fft2(generatedOscillation, NFFTY,NFFTX);

% shift the 2D spectrum. I'm not entirely sure why, but all the examples
% I've found seem to do it.
spectrum2D = fftshift(spectrum2D);

imagesc(abs(spectrum2D))

I can see four peaks on the 2D FT, but I can't figure out how to translate these peaks into the wavevectors I'm looking for (ie find k1, k2, k3 and k4). I'm also confused about what the x and y axis of the 2D FT mean, and their units. Comparing this with the 1D FT, I would guess that they have units of spatial cycles per unit length (or in this case, spatial cycles per pixel, since I've assumed 1 pixel = 1 unit length).

Any help would be most appreciated.

Many thanks

Dave

[David]

I've found in a textbook that for image size (Nx, Ny), the 2D fourier transform of sin(2*pi*(x*f_x+ y*f_y)) is NxNy/2 when (kx, ky) = -(f_x,f_y), and -iNxNy/2 when (kx,ky) = (f_x,f_y).

Based on this I expect to see spikes at locations (2*pi*k1, 2*pi*k2) and (2*pi*k3, 2*pi*k4). But there are four spikes not two, and also they are not where I expect them to be.

Is there something wrong with the code? Or is it simply my interpretation which is wrong?

Many thanks for your help

Dave

Answer 1

[David]

Hello all

I've since solved the problem, it was a simple coding issue, which caused the scales on the FFT to be wrong. The peaks were exactly where I expected them to be after the fix.

The code is posted below

vidHeight = 80; 
vidWidth = 100; 

% These are the wavevectors I want to find
k1 = 2*pi*0.1;
k2 = 2*pi*0.1; 
k3 = 2*pi*0.2;
k4 = 2*pi*0.2;

generatedOscillation = zeros(vidHeight, vidWidth);
for x = 1:vidWidth
    for y = 1:vidHeight
        % here i add a zero frequency component, which I later get rid of by
        % subtracting the average value of the image
        generatedOscillation(y,x) = generatedOscillation(y,x) + 1.2;

          % add one oscillation
          generatedOscillation(y,x) = generatedOscillation(y,x) + sin(k1*x + k2*y);

          % add another oscillation
          generatedOscillation(y,x) = generatedOscillation(y,x) + sin(k3*x + k4*y);
      end
  end

% show resulting image
% imagesc(generatedOscillation)

% dft 2d
NFFTY = 2^nextpow2(vidHeight);
NFFTX = 2^nextpow2(vidWidth);

% 'detrend' data to eliminate zero frequency component
av = sum(generatedOscillation(:)) / length(generatedOscillation(:));
generatedOscillation = generatedOscillation - av;

% Find X and Y frequency spaces, assuming sampling rate of 1
samplingFreq = 1;
spatialFreqsX = samplingFreq/2*linspace(0,1,NFFTX/2+1);
spatialFreqsY = samplingFreq/2*linspace(0,1,NFFTY/2+1);

spectrum2D = fft2(generatedOscillation, NFFTY,NFFTX);

contourf(spatialFreqsX, spatialFreqsY, abs(spectrum2D(1:NFFTY/2+1, 1:NFFTX/2+1)))

Maybe this will help someone one day.

(ps I'm not sure about the rules regarding accepting my own answer, feel free to take this post down if it's not ok)

[Rollo]

Thanks for this. It helped in interpretation of a 2D FFT

[Alon Shapira]

Thanks a lot!