如何根据 time-space 图像的 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 从视觉刺激生成时间-space 模式的功率谱，如下所示。这基本上是 2 秒时间范围内 10 个点（正弦波）的运动轨迹图，距离以度为单位。 （200x160 矩阵 - y 轴上每帧 10 毫秒，x 轴上每帧 0.1 度）。
我已经对这个刺激进行了 fft2、fftshift 和对数变换，结果输出是这样的。
首先，我对这个转换后的图像到底代表什么感到有点困惑？中心显示的是刺激的高频数据还是低频数据？在此转换后的图中，x 轴和 y 轴现在代表什么？
我实际上希望转换转换后的图像，使 y 轴反映 -30 到 30Hz 之间的时间频率，x 轴反映 -30deg/cycle 到 30deg/cycle 之间的空间频率。也许有人可以告诉我应该如何去做这件事？ （即是否有能够处理此类转换的 MATLAB 函数？）
重现情节的代码示例是：-
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);
**编辑：尝试按照以下方式重新创建一些东西，我只想要 -30 到 30 cycle/degree 和 -30 到 30Hz 范围内的功率谱图：-
只是想了解 fft 如何在 2D 上工作 space，您可以看看 here and,more useful, here。
换句话说，如果你像这样对图像进行 2D fft（请注意，一行只是一个 sin 函数，在 matlab 中非常容易实现）：
对应于：
现在，如果您构建相似的图像但周期不同，您将获得相似的结果，但 2D fft 中的点会更近。例如：
fft 的位置：
正弦曲线的方向与傅立叶图像中峰值相对于中心 DC 点的方向相关。在这种情况下，倾斜的正弦曲线图案会在傅里叶图像中产生一对倾斜的峰：
您可以尝试将不同的图像组合起来，观察 2Dfft 中的不同图案：
强烈建议你看一下答案开头的相关link
我是 FFT 的业余爱好者（没有受过物理学方面的训练！）所以我希望这里的人有专业知识可以给我一些提示，告诉我应该如何进行下一步。
所以我正在尝试通过 MATLAB 从视觉刺激生成时间-space 模式的功率谱，如下所示。这基本上是 2 秒时间范围内 10 个点（正弦波）的运动轨迹图，距离以度为单位。 （200x160 矩阵 - y 轴上每帧 10 毫秒，x 轴上每帧 0.1 度）。
我已经对这个刺激进行了 fft2、fftshift 和对数变换，结果输出是这样的。
首先，我对这个转换后的图像到底代表什么感到有点困惑？中心显示的是刺激的高频数据还是低频数据？在此转换后的图中，x 轴和 y 轴现在代表什么？
我实际上希望转换转换后的图像，使 y 轴反映 -30 到 30Hz 之间的时间频率，x 轴反映 -30deg/cycle 到 30deg/cycle 之间的空间频率。也许有人可以告诉我应该如何去做这件事？ （即是否有能够处理此类转换的 MATLAB 函数？）
重现情节的代码示例是：-
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);
**编辑：尝试按照以下方式重新创建一些东西，我只想要 -30 到 30 cycle/degree 和 -30 到 30Hz 范围内的功率谱图：-
只是想了解 fft 如何在 2D 上工作 space，您可以看看 here and,more useful, here。
换句话说，如果你像这样对图像进行 2D fft（请注意，一行只是一个 sin 函数，在 matlab 中非常容易实现）：
对应于：
现在，如果您构建相似的图像但周期不同，您将获得相似的结果，但 2D fft 中的点会更近。例如：
fft 的位置：
正弦曲线的方向与傅立叶图像中峰值相对于中心 DC 点的方向相关。在这种情况下，倾斜的正弦曲线图案会在傅里叶图像中产生一对倾斜的峰：
您可以尝试将不同的图像组合起来，观察 2Dfft 中的不同图案：
强烈建议你看一下答案开头的相关link