Как нормализовать 2D-график FFT на правильную частоту (Matlab)?

Question

Прежде всего, обратитесь к How to plot temporal frequency as a function of spatial frequency from a MATLAB FFT2 output of a time-space image?, чтобы узнать больше об этом вопросе.

Предполагая, что в случае этого сигнала образца: -

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); 

Переменная диапазон F вот что перебирает нормализацию оси х от

к следующему

Тем не менее, я изо всех сил пытаюсь найти способ нормализовать значения оси x и y для графика 2D FFT (изображение для сюжетов доступно по указанной выше ссылке в первом предложении этого сообщения.)

Кто-нибудь знает, как я должен это делать?

Фрагмент рабочей части моих кодов: -

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'); 

Я пытался найти руководства и помочь в Интернете, но безрезультатно до сих пор. Любая помощь будет принята с благодарностью!

Answer 1

Всякий раз, когда я не уверен в осях и масштабировании, я возвращаюсь к основам: сложная экспоненциальная (сложная синусоидальная, с вещественным = cos и мнимой = sin).

Я знаю, что 1D БПФ exp(j * 2 * pi * f * t), для вектора временных выборок t (в секундах), и частоту f (в Гц), будет иметь пик при f, до тех пор, как fmin < f < fmax, где fmax = 1/diff(t(1:2))/2 и fmin = -1.0 * fmax, и что пик будет иметь значение 1.0.

То же самое применяется в 2D-корпусе. 2D-комплексная экспоненциальная с частотой (fx, fy) будет иметь пик при fx и fy в соответствующих осях, а пиковое значение будет 1.0.

Вот полный пример Matlab, который работает через детали и соглашения, чтобы получить этот известный результат. Он моделирует двумерную сложную экспоненту с частотой x на частоте 2 Гц и y-частотой при -3 Гц по прямоугольной 2D-сетке. Затем он принимает БПФ после нулевого заполнения.

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 

Вот входные данные во временной области:

Убедитесь, что вы видите два пик-пик цикла в х размерности и три циклов в y- измерение - их легко увидеть, если вы изучите нижний и левый края фигуры.

Вот 2D FFT, соответственно смещаются (с fftshift), с осями масштабируется правильно (см fx и fy), и пик масштабируется правильно (посмотреть, как я делю выход fft2 с numel(im)).

Убедитесь, что пик находится в точке (2, -3) Гц, соответствующих [fx, fy], и что значение пики почти 1,0 (это немного меньше из-за квантованную сетку).

Итак, есть три вещи, которые я сделал, чтобы сделать всю эту работу:

• fftshift выход fft2,
• генерировать fx и fy правильно, чтобы соответствовать fftshift и
• шкала выходной fft2 по количеству элементов, которые он использует на до нулевого заполнения.

Надеюсь, вы можете расширить этот полный пример до своего собственного случая.

Answer 2

Я думаю, что это очень простой ответ - то, что вы ищете, - это масштаб по осям X и Y, который представляет нормированную частоту. Поскольку вы использовали fftshift, термин DC будет находиться в середине графика, поэтому ваши шкалы должны начинаться от -0,5 до 0,5. Это легко сделать с помощью команды imagesc(x,y,C), а не только imagesc(C), как вы сейчас делаете. (Вы используете X, но в помощи Matlab используется C, поскольку он представляет собой цветовой код).

Измените свою вторую последнюю линию:

imagesc([-0.5,0.5],[-0.5,0.5],X); 

, и это дает вам то, что я думаю, что вы просите.