如何根据时空图像的MATLAB FFT2输出绘制时间频率作为空间频率的函数？ [英] How to plot temporal frequency as a function of spatial frequency from a MATLAB FFT2 output of a time-space image?
问题描述
我是一名FFT业余爱好者（没有受过物理学训练！）所以我希望周围的人有专业知识给我一个关于如何进行下一步的提示。
所以我试图通过MATLAB从视觉刺激生成时空模式的功率谱，如下所示。这基本上是在2秒的时间范围内以10度（正弦波）的运动轨迹的曲线图，其中距离以度为单位。 （200x160矩阵 - y轴每帧10ms，x轴每帧0.1度）。
首先，我对这个转换后的图像究竟代表什么感到困惑？中心是否显示刺激的高频或低频数据？ x和y轴现在在这个转换后的图中表示什么？
我实际上希望转换变换后的图像，使y轴反映之间的时间频率 - 30至30Hz和x轴，空间频率在-30度/周期至30度/周期之间。也许有人可以让我知道我应该怎么做呢？ （即。是否有能够处理这种转换的MATLAB函数？）
重现这些图的代码示例如下： -
function STotal = playINTOdotty（varargin）
deg_speed = 15.35; ％dva / s
nr_of_dots = 10;
motion_type ='const';
％迭代次数
runs = 1;
stim_x = 160; ％1帧= 0.1d
stim_t = 200; ％1帧= 10ms
sin_cycle_dur = 80; ％80;
max_speed = deg_speed / 5.15; ％这非常非常抽象。基本上绘制出刺激图像，你会看到5.15是最好的价值。
sd =（sin_cycle_dur / 2）/ 6;
mu =（sin_cycle_dur / 2）/ 2;
sineTOTAL = 0;
counter = 1;
if nargin> 0
nr_of_dots = varargin {1};
结束
如果是nargin> 1
deg_speed = varargin {2};
结束
如果是nargin> 2
motion_type = varargin {3};
end
thisFTTOTAL = zeros（stim_t，stim_x）;
stimTOTAL =零（stim_t，stim_x）;
％initialize stim
stim = zeros（stim_t，stim_x）+ .5;
%%定义模拟/生成位置的随机点（在缩放到平均速度之前）
start_dot_pos = round（rand（1，nr_of_dots）。* stim_x）;
dot_pos =零（stim_t，nr_of_dots）;
dot_pos（1，:) = start_dot_pos;
％dot_pos（1，:) = 0;
dot_pos_sim =零（stim_t，nr_of_dots）;
dot_pos_sim（1，:) = start_dot_pos;
％dot_pos_sim（1，:) = 0;
%%为中性条件定义随机点。 dot_pos1用于Sine，dot_pos2用于Constant
start_dot_pos1 = round（rand（1，nr_of_dots / 2）。* stim_x）;
dot_pos1 =零（stim_t，nr_of_dots / 2）;
dot_pos1（1，:) = start_dot_pos1;
dot_pos_sim1 =零（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 =零（stim_t，nr_of_dots / 2）;
dot_pos2（1，:) = start_dot_pos2;
dot_pos_sim2 =零（stim_t，nr_of_dots / 2）;
dot_pos_sim2（1，:) = start_dot_pos2;
%%恒速平均值
CTotal = max_speed * sin_cycle_dur;
Cmean = max_speed / 2;
for q = 1：运行
%%计算位置列表以允许计算Gmean和Smean以缩放
t = 2：stim_t
switch motion_type
case'sine'
sine_speed = max_speed。* sin（（t-1）/ sin_cycle_dur * 2 * pi）; ％正弦公式
sineTOTAL = sineTOTAL + abs（sine_speed）; ％从Sine公式中添加所有正弦生成的值以获得平均值计算的总计
dot_pos_sim（t，:) = dot_pos_sim（t-1，:) + max_speed。* sin（（t-1）/ sin_cycle_dur * 2 * PI）; ％正弦模拟矩阵（缩放前）
case'高斯'
x = linspace（（mu-4 * sd），（mu + 4 * sd），sin_cycle_dur / 2）; ％高斯公式第1部分
y = 1 /（2 * pi * sd）* exp（ - （x-mu）。^ 2 /（2 * sd ^ 2））; ％高斯公式第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和y2形成全高斯循环。这里有两个周期（80 + 80帧）+ 1（因为stim_t是161）
dot_pos_sim（t，:) = dot_pos_sim（t-1，:) + yTOTAL（：，t）; ％高斯模拟矩阵（缩放前）
case'const'
如果t> 10&& t <= 30％这是最好的硬编码。需要改变一段时间。基本上根据指定的stim_t范围来定位点位置。
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'inutral'％Sine + Const代码融合（类似于上面）生成中性。
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）;如果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;
结束
结束
结束
yT = 0; ％计数器总结所有高斯的速度以形成所有帧的总数
%%计算意味着
为y = 1：stim_t
switch motion_type
case'正弦'
Smean = sineTOTAL / stim_t;
case'aussian'
yT = sum（y1）+ sum（abs（y2））* 5; ％y周期为y，y2
Gmean = yT / stim_t;
case'inutral'
Smean = sineTOTAL / stim_t;
结束
结束
%%比例位置为Cmean
为t = 1：stim_t
开关motion_type
case'sine'
dot_pos（t，:) = dot_pos_sim（t，:)。*（Cmean / Smean）;
case'augssian'
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，:); ％对于常数
结束
结束
％舍入
dot_pos = round（dot_pos）;
dot_pos1 = round（dot_pos1）;
dot_pos2 = round（dot_pos2）;
％包裹
dot_pos = mod（dot_pos，stim_x）+1;
dot_pos1 = mod（dot_pos1，stim_x）+1;
dot_pos2 = mod（dot_pos2，stim_x）+1;
％点数值为1到0.5刺激矩阵
表示t = 1：stim_t
开关motion_type
case'sine'
stim（ t，dot_pos（t，:)）= 1;
case'aussian'
stim（t，dot_pos（t，:)）= 1;
case'const'
stim（t，dot_pos（t，:)）= 1;
case'inutral'
stim（t，dot_pos1（t，:)）= 1;
stim（t，dot_pos2（t，:)）= 1;
结束
结束
F = fft2（刺激）;
S = abs（F）;
Fc =（fftshift（F））;
S2 = abs（Fc）; ％如果迭代内没有对数变换，则
％S2 = log（1 + abs（Fc））;迭代中的％Log变换
thisFTTOTAL = thisFTTOTAL + S2;
结束
thisFTTOTAL = thisFTTOTAL / runs;
S2 = log（1 + abs（thisFTTOTAL））; ％如果迭代内没有log变换
％S2 = thisFTTOTAL; ％如果迭代内的log变换
figure（1）
colormap（'gray'）;
x = linspace（0,16,5）;
y = linspace（0,2,10）;
imagesc（x，y，stim）;
xlabel（'degrees'）;
ylabel（'秒'）;
xlim（[0 16]）
figure（2）
colormap（'grey'）;
imagesc（S2）;
**编辑：尝试重新创建以下内容，我只想要力量-spectra绘图在-30到30周期/度和-30到30Hz的范围内： -
只是想知道fft如何在2D空间上工作，你可以看看
对应于：
现在，如果你建立一个类似的图像但是在不同的时期你将获得类似的结果，但2D fft中的点将更接近。例如：
其中fft将是：
正弦曲线的方向与傅立叶图像中的峰值相对于中心DC点的方向。在这种情况下，倾斜的正弦曲线模式在傅立叶图像中创建一对倾斜的峰值：
您可以尝试组合不同的图像并观察不同的图案在2Dfft中：
我强烈建议您查看一下相关链接nnig的答案。
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:
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.
这篇关于如何根据时空图像的MATLAB FFT2输出绘制时间频率作为空间频率的函数？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
