fft2 function in matlab

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

Answer 1

[Dr. Seis]

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

[ddd]

can you tell me in case 1 D same thing can be applied ?

[yebin tong]

Hello Why the timespace data is random, I think it should be a function relate to x and t. I tried this program on my own data, but I found there is difference between two images

[Walter Roberson]

"Why the timespace data is random"

In order to have some data to work with for example purposes.

"I think it should be a function relate to x and t"

No, it should be whatever it is.

[Aitor Garcia]

Thanks Dr. seis, very instructive!

[Addy]

Well written code. But this method is very slow for large data. (I tried for 3200X200 matrix)

[juntian chen]

If dt=0.5，then comes the error.

How can I solve this?

[Walter Roberson]

@juntian chen I am not seeing an error with dt=0.5 ?

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 = .5; % 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]

Is there matlab tool like pwelch in wavenumber spectrum analysi?

[Justin Kin Jun Hew]

There is no 2D welch method implemented in any package, you will have to do the power spectrum directly, i.e. summing the squares of the fourier coefficients through FFT algorithm

Answer 2

[Dr. Seis]

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

[Dr. Seis]

Incidentally, do you understand why I am multiplying things by "dt" (in the 1D example) and "dx*dt" (in the 2D example)?

[ddd]

Thanks for ur help . I tried to understand what you wrote,but I want to ask is dx*dt for know the size of each sample or for what ? and in case 2D did you mean to show that your DFT fuction and fft2 have the same result?

[Dr. Seis]

1. "dx*dt" or "dt" is included in the calculation for determining the discrete integral (i.e., summation). Above, I am computing the area (length x width) under each data point and summing all those areas together. If you do not include these values, it is implicitly assumed that width between samples (i.e., dx or dt) is 1 unit. However, it is likely that your data were not sampled at 1 unit per sample (e.g., most quality seismic data is rarely collected at 1 second per sample). If you simply want to look at the shape of your data in the frequency domain, then it is not necessary to include this as you will be multiplying your data by the same constant. However, if you plan to do any sort of calculations from the amplitude information in the frequency domain, then it is important that you include this information because your amplitudes will be off by 1/(dx*dt) or 1/dt (in the examples above).

2. Yes, the plots are simply showing that you get the same result (to double precision) in either case.

[Umar Farooq Ghumman]

Can you elaborate how the space changes to frequency domain? Initially you had distances on x-axis and time on y-axis. What will be the result of FFt2 on the space of this data? it will most definitely not be in time and distance space.

[Kristian Torres]

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!

Answer 3

[Bill Wood]

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.