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