2D Fourier transform- problems with the shift
显示 更早的评论
I have a 2D matrix representing a fluid's density field, I am trying to see if I can present it in an intelegent way with less parameters while maintaining key features. One of my thoughts was to represent it as a 2D fourier series and maintain only the leading spatial frequencies. I have written a code (attached) that performs an fft, sorts the terms based on their absolute values, and then adds them up from high to low to reconstruct the original field.
However, my reconstructed field is flipped. the corners and center are flipped. I am not sure why but it might be related to the fact that I used fftshit on the shifted data. I tried flipping it back by using fftshift again on the reconstructed field (or ifftshift), and it almost works, but not quite. you can see this by uncommenting line
would be thankful for any assistance!
Nathan
采纳的回答
David Goodmanson
2024-11-1
编辑：David Goodmanson
2024-11-2
Hi nathan,
It is a bit more inconvenent but you can reconstruct using a shifted grid as you did. I kept that method in in the code below.
I am not sure why using fftshift on rho_reconstructed did not work for you since it does seem to work fine.
The only real change made below is that for true accuracy the xy grid should not be
% Define the sample 2D density field 'rho'
Nx = 100;
Ny = 100;
xx = linspace(-5, 5, Nx);
yy = linspace(-5, 5, Nx);
[x, y] = meshgrid(xx,yy); (1)
but rather
xx = linspace(-5, 5, Nx+1);
xx(end) = [];
yy = linspace(-5, 5, Ny+1);
yy(end) = [];
[x, y] = meshgrid(xx,yy); (2)
i.e. the grid comes up one point short of going from -5 to 5, exactly in parallel with what you are doing to create the frequency grid with
fx = (-Nx/2:Nx/2-1) / Lx; % etc
Using (1) and reconstructing with all 10000 points results in
offby = max(max(abs((rho_adj)-rho))) % rho_adj is the final reconstructed result
offby = 0.0302
whereas using (2) results in
offby = 1.0568e-14
The plot below uses 100 terms for which
offby = 0.0451
% Define the sample 2D density field 'rho'
Nx = 100;
Ny = 100;
% xx = linspace(-5, 5, Nx);
% yy = linspace(-5, 5, Nx);
xx = linspace(-5, 5, Nx+1);
xx(end) = [];
yy = linspace(-5, 5, Ny+1);
yy(end) = [];
[x, y] = meshgrid(xx,yy);
rho = exp(-x.^2 - y.^2) .* sin(2 * pi * 0.2 * x) .* sin(2 * pi * 0.2 * y);
% Compute the 2D Fourier Transform and shift zero frequency to the center
F_rho = fft2(rho);
F_rho_shifted = fftshift(F_rho);
% Flatten F_rho_shifted and calculate magnitudes
F_flat = F_rho_shifted(:);
magnitudes = abs(F_flat);
% Sort the indices of Fourier coefficients by magnitude in descending order
[~, sorted_indices] = sort(magnitudes, 'descend');
% Initialize the reconstructed density field with zeros
rho_recon = zeros(Nx, Ny);
% Set up frequency ranges based on sampling intervals
dx = x(1,2) - x(1,1);
dy = y(2,1) - y(1,1);
Lx = Nx * dx;
Ly = Ny * dy;
fx = (-Nx/2:Nx/2-1) / Lx;
fy = (-Ny/2:Ny/2-1) / Ly;
% Generate the frequency grids for Fourier reconstruction
[FX, FY] = meshgrid(fx, fy);
% Loop through each Fourier coefficient in sorted order, adding one term at a time
%for idx = 1:(Nx*Ny)
for idx = 1:100
% Get the original 2D indices of the current sorted Fourier coefficient
[i, j] = ind2sub([Nx, Ny], sorted_indices(idx));
% Get the current Fourier coefficient (complex value)
F_coeff = F_rho_shifted(i, j);
% Add the Fourier term to the reconstruction
rho_recon = rho_recon + F_coeff/Nx/Ny*exp(2i*pi*(FX(i,j)*x + FY(i,j)*y));
end
% fftshift for final result
rho_adj = fftshift(real(rho_recon));
offby = max(max(abs((rho_adj)-rho)))
figure(1)
subplot(1,2,1)
contourf(x,y,rho,20,'LineStyle','none')
colorbar
title('original density field')
subplot(1,2,2)
contourf(x,y,rho_adj,20,'LineStyle','none')
title('100 terms')
colorbar
更多回答（1 个）
William Rose
2024-11-1
0 个投票
Reconstructing the 2d DFT is a little tricky. The attached script uses the Nk largest components of the 2d DFT to reconstuct ther image. When Nk=10, it makes this:
When Nk=25, it makes this:
The inverse 2d DFT has very small non-zero imaginary parts, due to round-off error. Therefore I take the real part.
The code includes lines that reconstruct the full 2d DFT of the original image from the left half of its DFT. The image reconstructed this way is called rho1. You can edit the script to display rho1 instead of rho2. rho1 should exactly match the original, rho, and it appears to do so. The lines to make rho1 may help illustrate the symmetry of the 2d DFT.
I expect the script will fail if Nx or Ny is odd. You could add some if(...) lines to deal with this possibility, if necessary.
2 个评论
William Rose
2024-11-1
@nathan blanc, notice that my script does not use shift(). I think shift() is an unnecessary step that makes it harder to get the indexing right. By the way, due to symmetry of the 2d DFT of a real image, you can divide the 2d DFT top/bottom or left/right, and reconstruct one half from the other half. I chose to split left/right. I find the Nk largest components in one half. Then I build the other half, using the known symmetry of the 2d DFT of a real sequence.
Nathan Blanc
2024-11-3
类别
在 帮助中心 和 File Exchange 中查找有关 Bartlett 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
