How can I extract a single component from a Fourier decomposition of a 2D matrix and plot it as a 2D map?

I have a scalar field (potential vorticity PV) on a 2-D square grid (599x599 cells).

I need to do a Fourier decomposition of the field and extract the contribution of the most important components.

Then, I have to plot the map of each component alone.

I would like to know how to extract the single components of the field and plot them as 2D (599x599) maps.

So far:

    % Spatial sampling intervals 
    % (in km per pixel):
    dx = 15;
    dy = 15;  

    F1 = fft2(PV);  % 2D Fourier Transform of the matrix PV

    F2 = fftshift(F1);  % zero-frequency component in the middle of the spectrum

    imagesc(abs(F2)),colorbar

Not sure what to do after this.. any help would be appreciated!

[Rick Rosson]

It would be helpful if you could provide the spatial sampling intervals in the x- and y directions (in units of distance per pixel, e.g. meters, cm, feet, or inches per pixel). For example:

   % Spatial sampling intervals 
   % (in meters per pixel):
   dx = 0.02;
   dy = 0.03;

These parameters would allow you to calibrate the spatial frequency in the Fourier domain, which you will need to do to determine which spectral components are most prominent.

[Carlo]

Thanks Rick, I have added the spatial sampling intervals in the x- and y- directions

Answer 1

[Rick Rosson]

Here's a start:

[ M, N ] = size(PV);

Fs_x = 1/dx;
Fs_y = 1/dy;

dF_x = Fs_x/M;
dF_y = Fs_y/N;

Fx = -Fs_x/2:dF_x:Fs_x/2-dF_x + mod(M,2)*dF_x/2;
Fy = -Fs_y/2:dF_y:Fs_y/2-dF_y + mod(N,2)*dF_y/2;

F = meshgrid(Fx,Fy);

[Carlo]

Thanks! Is F the matrix of wavenumbers? How do I relate it to the matrix PV?

[Rick Rosson]

F is the matrix of spatial frequencies (in cycles per km). You can easily derive the matrix of wavenumbers as:

   K = 2*pi*F;

in radians per km.

[Rick Rosson]

Slight mistake: meshgrid returns two matrices, one for each spatial direction:

   [ Fxx, Fyy ] = meshgrid(Fx,Fy);

Then:

   Kxx = 2*pi*Fxx;
   Kyy = 2*pi*Fyy;

[Carlo]

Do I need these wavenumber matrices? To display the spectral image in the wavenumber domain I am using Kx and Ky to represent the X- and Y- axes, respectively:

Kx = 2*pi*Fx;
Ky = 2*pi*Fy;
figure;pcolor(Kx,Ky,abs(F2)),shading flat,colorbar

I am still not sure about how to have a 599x599 map in the x,y domain with only a specific wavenumber component of the original map. Should I do some filtering on the 2D FDT map?

[Rick Rosson]

Yes.