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

Answer 1

[David Goodmanson]

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    

[Nathan Bblanc]

Thank you david, this seems to work!

Answer 2

[William Rose]

@nathan blanc,

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.

[William Rose]

@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 Bblanc]

Thank you Will, this seems to work!

Nathan