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
采纳的回答
Dr. Seis
2011-12-30
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');
9 个评论
ddd
2011-12-30
yebin tong
2017-7-5
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
2017-7-5
"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
2018-9-6
Thanks Dr. seis, very instructive!
Addy
2021-3-20
Well written code. But this method is very slow for large data. (I tried for 3200X200 matrix)
juntian chen
2021-12-14
If dt=0.5，then comes the error.
How can I solve this?
@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
2021-12-16
Is there matlab tool like pwelch in wavenumber spectrum analysi?
Justin Kin Jun Hew
2022-9-12
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
更多回答（2 个）
Dr. Seis
2011-12-30
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
5 个评论
Dr. Seis
2011-12-30
Incidentally, do you understand why I am multiplying things by "dt" (in the 1D example) and "dx*dt" (in the 2D example)?
ddd
2011-12-30
Dr. Seis
2011-12-30
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
2018-4-4
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
2019-12-2
编辑：Kristian Torres
2019-12-2
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
2021-9-29
0 个投票
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Resizing and Reshaping Matrices 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
