How to normalise a 2D FFT plot to the right frequency (Matlab)?

Question

First, refer to How to plot temporal frequency as a function of spatial frequency from a MATLAB FFT2 output of a time-space image? for a bit more of a background to this question.

Assuming in the case of this sample signal:-

n = [0:1024];
signal = sin(2*pi*n/10) + sin(2*pi*n/20) + sin(2*pi*n/30);

N = 2048; %At least twice of the n value
X = abs(fft(signal,N));
X = fftshift(X); %normalise data
F = [-N/2:N/2-1]/N; %normalise data - shift it to the correct frequency
plot(F,X);

The variable F range here is what sorts out the normalisation of the x-axis from

to the following

However, I'm struggling to figure out a way to normalise the x and y-axis values for a 2D FFT plot (The image for the plots are available on the above given link at the first sentence of this post.)

Does anyone have a clue as to how I should go about doing this?

A snippet of a working portion of my codes are:-

clear;

deg_speed = 15.35; %degrees visual angle/sec
max_speed = deg_speed/5.15; %converting the required deg_speed in terms of frames
nr_of_dots = 10; %number of dots
sin_cycle_dur = 80; %number of frames (along Nt) required to complete a sin wave.
sineTOTAL = 0;

Nx = 160; % Frames along x-axis. 1 frame = 0.1 dva
Nt = 200; % Frames along y-asis. 1 frame = 10ms

start_dot_pos = round(rand(1,nr_of_dots) .* Nx); %spawn random starting positions of dots
dot_pos = zeros(Nt, nr_of_dots); %Initialise 2D stimulus array
dot_pos(1,:) = start_dot_pos; %Fill up first line of 2D array with the starting position of dots

dot_pos_sim = zeros(Nt, nr_of_dots); %Setup simulated array so the final dot_pos can be scaled to mean speed of outher condition
dot_pos_sim(1,:) = start_dot_pos; %Fill up first line of 2D array with the starting position of dots

for a = 2:Nt
    sine_speed = max_speed .* sin((a-1) / sin_cycle_dur *2*pi); %Sine formula
    sineTOTAL = sineTOTAL + abs(sine_speed); %Add all sine generated values from Sine formula to get an overall total for mean calculation
    dot_pos_sim(a,:) = dot_pos_sim(a-1,:) + max_speed .* sin((a-1) / sin_cycle_dur *2*pi); %Sine simulated matrix (before scaling)
end

%Ignore this for loop for now. This is later required for normalising simulated
%array to the mean speed across other conditions.
for b = 1:Nt
    dot_pos(b,:) = dot_pos_sim(b,:);
end

dot_pos = round(dot_pos); %Frames are in integers, therefore all float values needed to be rounded up.
dot_pos = mod(dot_pos,Nx)+1; %Wrap the dots the go beyond the edges to the other side of the plot

%For all of the slots filled with dots, set the value from 0 to 1.
for c = 1:Nt
    stim(c,dot_pos(c,:)) = 1;
end

figure (1)
x=linspace(0,16,5);
y=linspace(0,2,10);
imagesc(x,y,stim); 
xlabel('degrees');
ylabel('seconds');
colormap('gray');

X = abs(fft2(stim));
X = fftshift(X); %normalise data
X = log(1+X);
figure (2)
imagesc(X);
colormap('gray');

I have been trying to find guides and help online but to no avail so far. Any help would be greatly appreciated!

Answer 1

Whenever I'm not sure about axes and scalings, I go back to basics: the complex exponential (a complex sinusoid, with real=cos and imaginary=sin).

I know that the 1D FFT of exp(j * 2 * pi * f * t), for a vector of time samples t (in seconds) and a frequency f (in Hz) will have a peak at f, as long as fmin < f < fmax, where fmax = 1 / diff(t(1:2)) / 2 and fmin = -1.0 * fmax, and that the peak will have value 1.0.

The exact same thing applies in the 2D case. A 2D complex exponential with frequency (fx, fy) will have peak at fx and fy in the respective axes, and the peak value will be 1.0.

Here's a complete Matlab example that works through the details and conventions to get this known result. It simulates a 2D complex exponential with x-frequency at 2 Hz and y-frequency at -3 Hz over a rectangular 2D grid. Then it takes the FFT after zero-padding.

clearvars

x = linspace(-2, 2, 100); % seconds
y = linspace(-3, 3, 200); % seconds

xFreq = 2; % Hz
yFreq = -3; % Hz

im = exp(2j * pi * y(:) * yFreq) * exp(2j * pi * x(:)' * xFreq);

figure;imagesc(x, y, real(im))
xlabel('x (seconds)'); ylabel('y (seconds)');
title('time-domain data (real)')
colorbar; colormap(flipud(gray))

Nfft = 4 * 2 .^ nextpow2(size(im));
imF = fftshift(fft2(im, Nfft(1), Nfft(2))) / numel(im);

fx = ([0 : Nfft(2) - 1] / Nfft(2) - 0.5) / diff(x(1:2));
fy = ([0 : Nfft(1) - 1] / Nfft(1) - 0.5) / diff(y(1:2));

figure; imagesc(fx, fy, abs(imF));
colorbar; colormap(flipud(gray))
xlabel('f_x (Hz)'); ylabel('f_y (Hz)')
title('Frequency-domain data (abs)')
grid; axis xy

Here's the input time-domain data:

Confirm that you see two peak-to-peak cycles in the x dimension and three cycles in the y-dimension—it's easy to see these if you study the bottom and left edges of the figure.

Here's the 2D FFT, appropriately shifted (with fftshift), with axes scaled correctly (see fx and fy), and peak scaled correctly (see how I divide the output of fft2 with numel(im)).

Confirm that the peak is at (2, -3) Hz corresponding to [fx, fy], and that the value of the peak is almost 1.0 (it's a little smaller because of of the quantized grid).

So, there’s three things I did to make all this work:

• fftshift the output of fft2,
• generate fx and fy correctly to match the fftshift, and
• scale the output of fft2 by the number of elements it operates on before zero-padding.

Hopefully you can extend this complete example to your own case.

Answer 2

I think this is a very simple answer - what you appear to be looking for is a scale along the X and Y axes that represents the normalised frequency. As you've used fftshift, the DC term will be in the middle of the plot, so your scales should run from -0.5 to 0.5. This is easily done with the imagesc(x,y,C) command, instead of just imagesc(C) as you are currently doing. (You use X, but the Matlab help uses C as it represents a colormap).

Change your second last line to:

imagesc([-0.5,0.5],[-0.5,0.5],X);

and it gives you what I think you are asking for.