# fft2 function in matlab

244 views (last 30 days)
ddd on 30 Dec 2011
Commented: juntian chen on 16 Dec 2021
Does anyone have idea or good tutorials for manually programming a fft2 function (discrete fast fourier transform)?
The function already exists in MATLAB, I just want to be able to understand how it works

Dr. Seis on 30 Dec 2011
Here is a quick example. I assume a random 2D image where the horizontal axis (columns) represents data collected in space domain and the vertical axis (rows) represents data collected in time domain. The two domains could easily be collected in the same domain. I also chose the number of samples collected in either domain to be different (i.e., a rectangular image), but the number of samples could easily be the same for both dimensions.
Nx = 64; % Number of samples collected along first dimension
Nt = 32; % Number of samples collected along second dimension
dx = 1; % Distance increment (i.e., Spacing between each column)
dt = .1; % Time increment (i.e., Spacing between each row)
x = 0 : dx : (Nx-1)*dx; % distance
t = 0 : dt : (Nt-1)*dt; % time
data_spacetime = randn(Nt,Nx); % random 2D matrix
Nyq_k = 1/(2*dx); % Nyquist of data in first dimension
Nyq_f = 1/(2*dt); % Nyquist of data in second dimension
dk = 1/(Nx*dx); % Wavenumber increment
df = 1/(Nt*dt); % Frequency increment
k = -Nyq_k : dk : Nyq_k-dk; % wavenumber
f = -Nyq_f : df : Nyq_f-df; % frequency
data_wavenumberfrequency = zeros(size(data_spacetime)); % initialize data
% Compute 2D Discrete Fourier Transform
for i1 = 1 : Nx
for j1 = 1 : Nt
for i2 = 1 : Nx
for j2 = 1 : Nt
data_wavenumberfrequency(j1,i1) = data_wavenumberfrequency(j1,i1) + ...
data_spacetime(j2,i2)*exp(-1i*(2*pi)*(k(i1)*x(i2)+f(j1)*t(j2)))*dx*dt;
end
end
end
end
figure;
subplot(3,1,1);
imagesc(k,f,abs(data_wavenumberfrequency));
colorbar; v = caxis;
title('2D DFT');
fft2result = fftshift(fft2(data_spacetime))*dx*dt;
subplot(3,1,2);
imagesc(k,f,abs(fft2result));
colorbar; caxis(v);
title('FFT2');
subplot(3,1,3);
imagesc(k,f,abs(fft2result-data_wavenumberfrequency));
colorbar; caxis(v);
title('Difference');
juntian chen on 16 Dec 2021
Is there matlab tool like pwelch in wavenumber spectrum analysi?

Dr. Seis on 30 Dec 2011
The 2D case is associated with double integration... so the 1D case is single integration. The above answer will simplify down to the 1D case if you only have 1 sample location (e.g., only one "x" or only one "t"). So, for example, if you only had one spatial location (i.e., x = 0), then you would not need any of the variables associated with "Nx" or "dx". You would also not need the FOR loops associated with "Nx" either. Your FOR loop above would simply to:
data_time = randn(Nt,1);
data_frequency = zeros(size(data_time));
% Compute 1D Discrete Fourier Transform
for j1 = 1 : Nt
for j2 = 1 : Nt
data_frequency(j1) = data_frequency(j1) + ...
data_time(j2)*exp(-1i*(2*pi)*(f(j1)*t(j2)))*dt;
end
end
Kristian Torres on 2 Dec 2019
Hi. Does this mean that we also need to multiply df*dk if we are applying the Inverse fourier transform (e.g. ifft2(ifftshift(dummy))*df*dk )??
EDIT: I realized that the way I get the right amplitudes back through the inverse fft is with ifft2(ifftshift(dummy))/(dx*dt).
Thank you!

Bill Wood on 29 Sep 2021
It seems that for any square matrix A, fft2(A) = fft(fft(A).').' although I do get small differences when I compute this both ways. Differences in implementation, I guess.