问如何从时空图像的MATLAB FFT2输出中将时间频率绘制为空间频率的函数？
EN

Stack Overflow用户

提问于 2016-09-27 18:39:55

回答 1查看 615关注 0票数 1

我是一个业余的FFT爱好者(没有受过物理训练！)所以我希望周围的人有专业知识来给我一个提示，告诉我下一步应该如何做。

因此，我正在尝试通过MATLAB从视觉刺激中生成时空模式的功率谱，如下所示。这基本上是在2秒的时间范围内10个点(正弦波)的运动轨迹的曲线图，距离以度标记。(200x160矩阵-y轴上每帧10ms，x轴上每帧0.1度)。

我对这个刺激做了fft2，fftshift和对数变换，结果输出是这样的。

首先，我有点困惑，这张变形的图像到底代表了什么？中心显示的是刺激的高频数据还是低频数据？在这个转换后的图中，x轴和y轴现在代表什么？

我实际上希望转换转换后的图像，使y轴反映-30到30 to之间的时间频率和x轴，空间频率在-30度/周期到30度/周期之间。也许有人可以给我一个想法，我应该怎么做？(即有没有一个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周/度和-30到30 the范围内的功率谱图：-