如何将 2D FFT 图归一化到正确的频率(Matlab)? [英] How to normalise a 2D FFT plot to the right frequency (Matlab)?

查看:64
本文介绍了如何将 2D FFT 图归一化到正确的频率(Matlab)?的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

首先参考

以下内容

但是,我正在努力找出一种方法来规范化 2D FFT 图的 x 和 y 轴值(这些图的图像可在本文第一句的上面给出的链接中找到.)

有人知道我应该怎么做吗?

我的代码工作部分的片段是:-

clear;deg_speed = 15.35;%度视角/秒max_speed = deg_speed/5.15;% 以帧为单位转换所需的 deg_speednr_of_dots = 10;%点数sin_cycle_dur = 80;完成正弦波所需的帧数(沿 Nt).正弦总= 0;Nx = 160;% 沿 x 轴的帧数.1 帧 = 0.1 dvaNt = 200;% 沿 y 轴的帧数.1 帧 = 10msstart_dot_pos = round(rand(1,nr_of_dots) .* Nx);%产生点的随机起始位置dot_pos = zeros(Nt, nr_of_dots);%初始化二维刺激数组dot_pos(1,:) = start_dot_pos;%用点的起始位置填充二维数组的第一行dot_pos_sim = zeros(Nt, nr_of_dots);% 设置模拟数组,以便最终的 dot_pos 可以缩放为外部条件的平均速度dot_pos_sim(1,:) = start_dot_pos;%用点的起始位置填充二维数组的第一行对于 a = 2:Ntsine_speed = max_speed .* sin((a-1)/sin_cycle_dur *2*pi);%正弦公式sineTOTAL = sineTOTAL + abs(sine_speed);%从正弦公式中添加所有正弦生成的值以获得平均值计算的总和dot_pos_sim(a,:) = dot_pos_sim(a-1,:) + max_speed .* sin((a-1)/sin_cycle_dur *2*pi);%Sine 模拟矩阵(缩放前)结尾%暂时忽略这个 for 循环.稍后需要对模拟进行归一化%array 到其他条件下的平均速度.对于 b = 1:Ntdot_pos(b,:) = dot_pos_sim(b,:);结尾点位置 = 圆(点位置);%Frames 是整数,因此所有浮点值都需要四舍五入.dot_pos = mod(dot_pos,Nx)+1;%将超出边缘的点包裹到图的另一侧%对于所有用点填充的插槽,将值设置为 0 到 1.对于 c = 1:Nt刺激(c,dot_pos(c,:))= 1;结尾图1)x=linspace(0,16,5);y=linspace(0,2,10);图像c(x,y,stim);xlabel('度数');ylabel('秒');颜色图('灰色');X = abs(fft2(stim));X = fftshift(X);%标准化数据X = 日志(1+X);图(2)图像c(X);颜色图('灰色');

我一直在尝试在线查找指南和帮助,但到目前为止都无济于事.任何帮助将不胜感激!

解决方案

每当我不确定轴和缩放比例时,我都会回到基础:复指数(复正弦曲线,实数=cos 和虚数=sin).

我知道 exp(j * 2 * pi * f * t) 的一维 FFT,对于时间样本向量 t(以秒为单位)和一个频率 f(以赫兹为单位)将在 f 处有一个峰值,只要 fmin <;f<fmax,其中 fmax = 1/diff(t(1:2))/2fmin = -1.0 * fmax,并且峰值将具有值 1.0.

完全相同的事情适用于 2D 情况.频率为 (fx, fy) 的二维复指数将在各自轴的 fxfy 处出现峰值,峰值将为1.0.

这是一个完整的 Matlab 示例,它通过详细信息和约定来获得此已知结果.它在矩形 2D 网格上模拟 2D 复指数,x 频率为 2 Hz,y 频率为 -3 Hz.然后在补零后进行 FFT.

clearvarsx = linspace(-2, 2, 100);% 秒y = linspace(-3, 3, 200);% 秒xFreq = 2;%赫兹yFreq = -3;%赫兹im = exp(2j * pi * y(:) * yFreq) * exp(2j * pi * x(:)' * xFreq);图;imagesc(x, y, real(im))xlabel('x (秒)');ylabel('y (秒)');title('时域数据(真实)')颜色条;颜色图(翻转(灰色))Nfft = 4 * 2 .^ nextpow2(size(im));imF = fftshift(fft2(im, Nfft(1), Nfft(2)))/numel(im);fx = ([0: Nfft(2) - 1]/Nfft(2) - 0.5)/diff(x(1:2));fy = ([0 : Nfft(1) - 1]/Nfft(1) - 0.5)/diff(y(1:2));数字;图像c(fx,fy,abs(imF));颜色条;颜色图(翻转(灰色))xlabel('f_x (Hz)');ylabel('f_y (Hz)')title('频域数据(abs)')网格;xy轴

这是输入的时域数据:

确认您在 x 维度上看到 两个 峰峰值循环,在 y 维度上看到 三个 循环——如果您研究图的下边缘和左边缘.

这是 2D FFT,适当移动(使用 fftshift),轴正确缩放(参见 fxfy),峰值缩放正确(看看我如何将 fft2 的输出与 numel(im) 分开).

确认[fx, fy]对应的峰值在(2, -3)Hz处,峰值接近1.0(小一点,因为量化网格).

所以,我做了三件事来完成这一切:

  • fftshift fft2 的输出,
  • 正确生成fxfy 以匹配fftshift
  • fft2 的输出按在补零之前操作的元素数量进行缩放.

希望您可以将此完整示例扩展到您自己的案例中.

First, refer to How to plot temporal frequency as a function of spatial frequency from a MATLAB FFT2 output of a time-space image? for a bit more of a background to this question.

Assuming in the case of this sample signal:-

n = [0:1024];
signal = sin(2*pi*n/10) + sin(2*pi*n/20) + sin(2*pi*n/30);

N = 2048; %At least twice of the n value
X = abs(fft(signal,N));
X = fftshift(X); %normalise data
F = [-N/2:N/2-1]/N; %normalise data - shift it to the correct frequency
plot(F,X);

The variable F range here is what sorts out the normalisation of the x-axis from

to the following

However, I'm struggling to figure out a way to normalise the x and y-axis values for a 2D FFT plot (The image for the plots are available on the above given link at the first sentence of this post.)

Does anyone have a clue as to how I should go about doing this?

A snippet of a working portion of my codes are:-

clear;

deg_speed = 15.35; %degrees visual angle/sec
max_speed = deg_speed/5.15; %converting the required deg_speed in terms of frames
nr_of_dots = 10; %number of dots
sin_cycle_dur = 80; %number of frames (along Nt) required to complete a sin wave.
sineTOTAL = 0;

Nx = 160; % Frames along x-axis. 1 frame = 0.1 dva
Nt = 200; % Frames along y-asis. 1 frame = 10ms

start_dot_pos = round(rand(1,nr_of_dots) .* Nx); %spawn random starting positions of dots
dot_pos = zeros(Nt, nr_of_dots); %Initialise 2D stimulus array
dot_pos(1,:) = start_dot_pos; %Fill up first line of 2D array with the starting position of dots

dot_pos_sim = zeros(Nt, nr_of_dots); %Setup simulated array so the final dot_pos can be scaled to mean speed of outher condition
dot_pos_sim(1,:) = start_dot_pos; %Fill up first line of 2D array with the starting position of dots

for a = 2:Nt
    sine_speed = max_speed .* sin((a-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(a,:) = dot_pos_sim(a-1,:) + max_speed .* sin((a-1) / sin_cycle_dur *2*pi); %Sine simulated matrix (before scaling)
end

%Ignore this for loop for now. This is later required for normalising simulated
%array to the mean speed across other conditions.
for b = 1:Nt
    dot_pos(b,:) = dot_pos_sim(b,:);
end

dot_pos = round(dot_pos); %Frames are in integers, therefore all float values needed to be rounded up.
dot_pos = mod(dot_pos,Nx)+1; %Wrap the dots the go beyond the edges to the other side of the plot

%For all of the slots filled with dots, set the value from 0 to 1.
for c = 1:Nt
    stim(c,dot_pos(c,:)) = 1;
end

figure (1)
x=linspace(0,16,5);
y=linspace(0,2,10);
imagesc(x,y,stim); 
xlabel('degrees');
ylabel('seconds');
colormap('gray');

X = abs(fft2(stim));
X = fftshift(X); %normalise data
X = log(1+X);
figure (2)
imagesc(X);
colormap('gray');

I have been trying to find guides and help online but to no avail so far. Any help would be greatly appreciated!

解决方案

Whenever I'm not sure about axes and scalings, I go back to basics: the complex exponential (a complex sinusoid, with real=cos and imaginary=sin).

I know that the 1D FFT of exp(j * 2 * pi * f * t), for a vector of time samples t (in seconds) and a frequency f (in Hz) will have a peak at f, as long as fmin < f < fmax, where fmax = 1 / diff(t(1:2)) / 2 and fmin = -1.0 * fmax, and that the peak will have value 1.0.

The exact same thing applies in the 2D case. A 2D complex exponential with frequency (fx, fy) will have peak at fx and fy in the respective axes, and the peak value will be 1.0.

Here's a complete Matlab example that works through the details and conventions to get this known result. It simulates a 2D complex exponential with x-frequency at 2 Hz and y-frequency at -3 Hz over a rectangular 2D grid. Then it takes the FFT after zero-padding.

clearvars

x = linspace(-2, 2, 100); % seconds
y = linspace(-3, 3, 200); % seconds

xFreq = 2; % Hz
yFreq = -3; % Hz

im = exp(2j * pi * y(:) * yFreq) * exp(2j * pi * x(:)' * xFreq);

figure;imagesc(x, y, real(im))
xlabel('x (seconds)'); ylabel('y (seconds)');
title('time-domain data (real)')
colorbar; colormap(flipud(gray))

Nfft = 4 * 2 .^ nextpow2(size(im));
imF = fftshift(fft2(im, Nfft(1), Nfft(2))) / numel(im);

fx = ([0 : Nfft(2) - 1] / Nfft(2) - 0.5) / diff(x(1:2));
fy = ([0 : Nfft(1) - 1] / Nfft(1) - 0.5) / diff(y(1:2));

figure; imagesc(fx, fy, abs(imF));
colorbar; colormap(flipud(gray))
xlabel('f_x (Hz)'); ylabel('f_y (Hz)')
title('Frequency-domain data (abs)')
grid; axis xy

Here's the input time-domain data:

Confirm that you see two peak-to-peak cycles in the x dimension and three cycles in the y-dimension—it's easy to see these if you study the bottom and left edges of the figure.

Here's the 2D FFT, appropriately shifted (with fftshift), with axes scaled correctly (see fx and fy), and peak scaled correctly (see how I divide the output of fft2 with numel(im)).

Confirm that the peak is at (2, -3) Hz corresponding to [fx, fy], and that the value of the peak is almost 1.0 (it's a little smaller because of of the quantized grid).

So, there’s three things I did to make all this work:

  • fftshift the output of fft2,
  • generate fx and fy correctly to match the fftshift, and
  • scale the output of fft2 by the number of elements it operates on before zero-padding.

Hopefully you can extend this complete example to your own case.

这篇关于如何将 2D FFT 图归一化到正确的频率(Matlab)?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆