How to plot temporal frequency as a function of sp

I'm a little bit of a FFT amateur (not trained in physics!) so I'm hoping someone around here has the expertise to give me a hint as to how I should go about doing this next step.

So I'm trying to generate the power spectra of time-space pattern via MATLAB from a visual stimulus as shown below. This is basically a plot of the movement trajectory of 10 dots (sine wave) within a time frame of 2 seconds with the distance labelled in degrees. (200x160 matrix - 10ms per frame on the y-axis and 0.1 degrees per frame on the x-axis).

I have done fft2, fftshift and a log transform on this stimulus and the resulting output is this.

First off, I am a little confused as to what this transformed image exactly represent? Is the centre displaying the high or low frequency data of the stimulus? And what do the x and y-axis now represents in this transformed plot?

I am actually hoping to convert the transformed image such that the y axis reflects temporal frequency between -30 to 30Hz and the x axis, spatial frequency between -30deg/cycle to 30deg/cycle. Perhaps someone could give me an idea of how I should go about doing this? (ie. is there a MATLAB function that is able to handle this sort of conversion?)

A sample of the codes to reproduce the plots are:-

function STotal = playINTOdotty (varargin)

deg_speed = 15.35; %dva/s
nr_of_dots = 10;
motion_type = 'const';

%Number of iterations
runs = 1;

stim_x = 160; %1 frame = 0.1d
stim_t = 200; %1 frame = 10ms
sin_cycle_dur = 80; %80; 

max_speed = deg_speed/5.15; %This is very, very abstract. Basically plot out stim image and you'll see 5.15 is the best value.

sd = (sin_cycle_dur/2)/6;
mu = (sin_cycle_dur/2)/2;

sineTOTAL = 0;
counter = 1;

if nargin > 0
    nr_of_dots = varargin{1};
end
if nargin > 1
    deg_speed = varargin{2};
end
if nargin > 2
    motion_type = varargin{3};
end

thisFTTOTAL = zeros(stim_t,stim_x);
stimTOTAL = zeros(stim_t,stim_x);

% initialize stim
stim = zeros(stim_t, stim_x) + .5;

%% define random dots for simulation/generation of position (before scaling to mean speed)

start_dot_pos = round(rand(1,nr_of_dots) .* stim_x);
dot_pos = zeros(stim_t, nr_of_dots);
dot_pos(1,:) = start_dot_pos;
%dot_pos(1,:) = 0;

dot_pos_sim = zeros(stim_t, nr_of_dots);
dot_pos_sim(1,:) = start_dot_pos;
%dot_pos_sim(1,:) = 0;

%% define random dots for neutral condition. dot_pos1 is for Sine and dot_pos2 for Constant 

start_dot_pos1 = round(rand(1,nr_of_dots/2) .* stim_x);
dot_pos1 = zeros(stim_t, nr_of_dots/2);
dot_pos1(1,:) = start_dot_pos1;

dot_pos_sim1 = zeros(stim_t, nr_of_dots/2);
dot_pos_sim1(1,:) = start_dot_pos1;

start_dot_pos2 = round(rand(1,nr_of_dots/2) .* stim_x);
dot_pos2 = zeros(stim_t, nr_of_dots/2);
dot_pos2(1,:) = start_dot_pos2;

dot_pos_sim2 = zeros(stim_t, nr_of_dots/2);
dot_pos_sim2(1,:) = start_dot_pos2;

%% Mean of Constant speed

CTotal = max_speed*sin_cycle_dur;
Cmean = max_speed/2;

for q = 1:runs
    %% Calculate position list to allow calculation of Gmean and Smean for scaling
    for t = 2:stim_t
        switch motion_type
            case 'sine'
                sine_speed = max_speed .* sin((t-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(t,:) = dot_pos_sim(t-1,:) + max_speed .* sin((t-1) / sin_cycle_dur *2*pi); %Sine simulated matrix (before scaling)
            case 'gaussian'
                x = linspace((mu-4*sd),(mu+4*sd),sin_cycle_dur/2); %Gaussian formula part 1
                y = 1/(2*pi*sd)*exp(-(x-mu).^2/(2*sd^2)); %Gaussian formula part 2
                scalefactor = max_speed / (1/(2*pi*sd)); 
                y = y*scalefactor;
                y1 = y;
                y2 = -y;
                yTOTAL = [y,y2,y,y2,y,y2,y,y2,y,y2]; %y and y2 forms a full gaussian cycle. Two cycles here (80+80 frames) + 1 (Because stim_t is 161)
                dot_pos_sim(t,:) = dot_pos_sim(t-1,:) + yTOTAL(:,t); %Gaussian simulated matrix (before scaling)
            case 'const'
                if t > 10 && t <= 30 %This is hard coding at its best. Need to change this some time. Basically definding dot positions based on the specified stim_t range. 
                    con_speed = max_speed;
                    dot_pos_sim(t,:) = dot_pos_sim(t-1,:) + con_speed; 
                elseif t > 50 && t <= 70
                    con_speed = -max_speed;
                    dot_pos_sim(t,:) = dot_pos_sim(t-1,:) + con_speed; 
                elseif t > 90 && t <= 110
                    con_speed = max_speed;
                    dot_pos_sim(t,:) = dot_pos_sim(t-1,:) + con_speed;
                elseif t > 130 && t <= 150
                    con_speed = -max_speed;
                    dot_pos_sim(t,:) = dot_pos_sim(t-1,:) + con_speed; 
                elseif t > 170 && t <= 190
                    con_speed = max_speed;
                    dot_pos_sim(t,:) = dot_pos_sim(t-1,:) + con_speed;                
                else
                    con_speed = 0;
                    dot_pos_sim(t,:) = dot_pos_sim(t-1,:) + con_speed;           
                end
            case 'neutral' %Fusion of Sine + Const codes (similar to above) to generate neutral.
                sine_speed = max_speed .* sin((t-1) / sin_cycle_dur *2*pi);
                sineTOTAL = sineTOTAL + abs(sine_speed); 
                dot_pos_sim1(t,:) = dot_pos_sim1(t-1,:) + max_speed .* sin((t-1) / sin_cycle_dur *2*pi); 

                if t > 10 && t <= 30 
                    con_speed = max_speed;
                    dot_pos_sim2(t,:) = dot_pos_sim2(t-1,:) + con_speed; 
                elseif t > 50 && t <= 70
                    con_speed = -max_speed;
                    dot_pos_sim2(t,:) = dot_pos_sim2(t-1,:) + con_speed; 
                elseif t > 90 && t <= 110
                    con_speed = max_speed;
                    dot_pos_sim2(t,:) = dot_pos_sim2(t-1,:) + con_speed;
                elseif t > 130 && t <= 150
                    con_speed = -max_speed;
                    dot_pos_sim2(t,:) = dot_pos_sim2(t-1,:) + con_speed; 
                elseif t > 170 && t <= 190
                    con_speed = max_speed;
                    dot_pos_sim2(t,:) = dot_pos_sim2(t-1,:) + con_speed;   
                else
                    con_speed = 0;
                    dot_pos_sim2(t,:) = dot_pos_sim2(t-1,:) + con_speed;           
                end
        end     
    end 

    yT = 0; %counter to sum up all of gaussian's speed to form a total from all frames

    %% Calculate means
    for y = 1:stim_t
        switch motion_type
            case 'sine'
                Smean = sineTOTAL/stim_t;
            case 'gaussian'
                yT = sum(y1) + sum(abs(y2)) * 5; %5 cycles of y,y2
                Gmean = yT/stim_t;
            case 'neutral'
                Smean = sineTOTAL/stim_t;
        end
    end

    %% Scale positions to Cmean
    for t = 1:stim_t
        switch motion_type
            case 'sine'
                dot_pos(t,:) = dot_pos_sim(t,:) .* (Cmean/Smean); 
            case 'gaussian'
                dot_pos(t,:) = dot_pos_sim(t,:) .* (Cmean/Gmean);
            case 'const'
                dot_pos(t,:) = dot_pos_sim(t,:);
            case 'neutral'
                dot_pos1(t,:) = dot_pos_sim1(t,:) .* (Cmean/Smean); %For Sine
                dot_pos2(t,:) = dot_pos_sim2(t,:); %For Constant
        end     
    end 

    %rounding
    dot_pos = round(dot_pos);
    dot_pos1 = round(dot_pos1);
    dot_pos2 = round(dot_pos2);
    %wrapping
    dot_pos = mod(dot_pos,stim_x)+1;
    dot_pos1 = mod(dot_pos1,stim_x)+1;
    dot_pos2 = mod(dot_pos2,stim_x)+1;

    %Dots given a value of 1 to the 0.5 stim matrix
    for t = 1:stim_t
        switch motion_type
            case 'sine'
                stim(t,dot_pos(t,:)) = 1;
            case 'gaussian'
                stim(t,dot_pos(t,:)) = 1;
            case 'const'
                stim(t,dot_pos(t,:)) = 1;
            case 'neutral'
                stim(t,dot_pos1(t,:)) = 1;
                stim(t,dot_pos2(t,:)) = 1;
        end
    end

    F = fft2(stim);
    S = abs(F);

    Fc = (fftshift(F));
    S2 = abs(Fc); %If without log transform within iteration
    %S2 = log(1+abs(Fc)); %Log transform within iteration

    thisFTTOTAL = thisFTTOTAL + S2;
end

thisFTTOTAL = thisFTTOTAL/runs;
S2 = log(1+abs(thisFTTOTAL)); %If without log transform within iteration
%S2 = thisFTTOTAL; %If log transform within iteration

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

figure (2)
colormap('gray');
imagesc(S2);

**EDIT : Trying to recreate something along the lines of the following, where I only want the power-spectra plots within the range of -30 to 30 cycle/degree and -30 to 30Hz:-

Just to have an idea on how the fft works on a 2D space,you can have a look here and,more useful, here.

In other words, if you do an 2D fft of an image like this (please note that a row it is just a sin function, really easy to implement in matlab):

corresponds to:

Now, if you build a similar image but with a different period you will obtain a similar result but the point in the 2D fft will be closer. For example:

where the fft will be:

The orientation of the sinusoid correlates with the orientation of the peaks in the Fourier image relative to the central DC point. In this case a tilted sinusoidal pattern creates a tilted pair of peaks in the Fourier image: