How to plot temporal frequency as a function of spatial frequency from a MATLAB FFT2 output of a time-space image?

Question

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:-

If you need more info you can ask so we can try to answer in a better way. However, if you appreciate the answer, you can accept the answer and to upvote if you agree. — Leos313, Sep 27 '16 at 14:38

score 1 · Accepted Answer · answered Sep 27 '16 at 12:42

1

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:

You can try to combine the different image and observe the different pattern in the 2Dfft:

I strongly recommend you to have a look on the related link at the beginnig of the answer.

answered Sep 27 '16 at 12:42

Leos313

5,152
6
40
69

1

Hey thanks for these. I'd just want to clarify, in the context of MATLAB, I'm assuming the 2D fft of these sinusoidal gratings have been fftshifted and log transformed? Or are they the raw fft2 output from the stimulus? Because I tried to recreate the fourier plots but the dots started off at the edges and not the centre. Also, suppose that the fftshifted image now produces the output at the centre, this is the DC term of the 2D plot? (Ie. as the image spans out further into the horizontal periphery, it reflects increasing spatial frequency and temporal frequency on the vertical periphery?) – Ang Jit Wei Aaron Sep 28 '16 at 06:10
fftshifed of course!!! And the little dot in the middle is the commonly-called DC term of the FT, you are right. I didn't understand the last question in the brackets. Try to explain better and I will try to answer :) and please accept or vote up if the answer is useful! :) – Leos313 Sep 28 '16 at 07:27
Yes it was indeed helpful thus far! Regarding the latter part of my question, refer to the edited main post where I added an image at the end of it. I am trying to recreate plots (C), which were generated from the FFT of (A) windowed by the gaussian window of (B). A brief explanation of why the range of -30 to 30 for plots (C) is that human vision cannot see anything higher than -30 to 30Hz and -30 to 30 degree visual angle, thus it has to be windowed within this range. What I am puzzled about, is how do you extract values from the FFT output to be able to plot values along such axis range? – Ang Jit Wei Aaron Sep 28 '16 at 10:15

How to plot temporal frequency as a function of spatial frequency from a MATLAB FFT2 output of a time-space image?

1 Answers1

Linked