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
