How to calculate wavenumber, for fft2 of a 2D array?

I have perfomed fft2 on an image/2D-array (361x211). When I plot the result, using imagesc, I would like to label the axes with wavenumbers. However, for the fft2 output, I do not know how to calculate the wavenumbers (I commented out my this part of my code). Is it related to the image resolution (dpi)?

My code and output plots (amplitude and phase) are attached below. Thanks!

data = load ('86tr.txt'); % load xyz file

x = data(:,1); % longitude

y = data(:,2); % latitude

z = data(:,3); % anomaly (nano Tesla)

% Find the size (row,col) of the xyz file

row = length(unique(x));

col = length(unique(y));

% Reshape and transpose z values to get correct data orientaion

zm = reshape(z,row,col);

zm = fliplr(zm);

zm = zm';

% Perform fft2 on z-values (nT)

fft2_zm = fftshift(log(abs(fft2(zm)))); % amplitude

% phase_zm = angle(fftshift(fft2(zm))); % phase

% % Calculate wavenumbers

% Nc = col; Nr = row;

% dp = 1/300; % dpi = 300?

% kx = (-Nr/2:Nr/2-1)/(Nr*dp);

% ky = (-Nc/2:Nc/2-1)/(Nc*dp);

figure(1);

imagesc(fft2_zm); colorbar;

[Star Strider]

The wavenumber values would be along the z axis and would be in terms of the (probably spatial) frequency (although there are no units listed for z in the file). You would need to calculate and show them as different colours, likely using the colorbar function. The other option would be to use surf, and then display them on the z-axis, likely also using colorbar as well, for clarity.

% T1 = readtable('86tr.txt')

.

[Jay Ghosh]

@Star Strider, The unit for the z-values is nano Tesla (nT). Wavenumbers should be along the z-axis? I was under the impression that the colorbar here (now added to the figure) represents the z-axis, i.e the amplitudes (energy), which is the fft2 values of z (zm). Should the wavenumbers not be labelled along x, and y axis, since this is 2D array?

Also, how do I calculate wavenumbers? Is it 2π/N, where N = length(row) or length(col)? Thank you.

[Star Strider]

I looked at the original data by ‘slicing’ it horizontally at various levels. There are patterns that show up in the east-west direction, however nothing in the north-south direction. I would simply go with spatial frequency rather than wavenumbers.

[Jay Ghosh]

@Star Strider, What would be the spatial frequency range then, (-pi to pi)? I think I am still lost. Moreover, as I asked before, for the fft2 of lena512.bmp image, was my axes labeling correct? I do not see my post (and the codes) related to that anymore. Thanks.

Also, did you delete some of your previous responses (including the 3D (surf) plot? Is it possible to reinstate those? I wanted to take a look at those codes and your respoinses again. Thank you. I appreciate it.

[Star Strider]

The spatial frequency would be in terms of . I am not certain that it would be possible to expres anything about your data in terms of wavenumber.

I deleted my Answer because it became obvious to me that it was not going to provide the result you want, essentially because I am not certain that is possible.

[Jay Ghosh]

@Star Strider, Thank you for your comments. I appreciate it.

Answer 1

[Hiro Yoshino]

We use "normalized frequency" with FFT and this is a linear transformation between the time space and the frequency space. So the number of points along the x-axis and the y-axis corresponds to the wavenumber.

<< This is what I called "normalized frequency".

If you ask me, I would rewrite your code as follows:

data = load ('86tr.txt');  % load xyz file
x = data(:,1);   % longitude 
y = data(:,2);   % latitude
z = data(:,3);   % anomaly (nano Tesla)
% Find the size (row,col) of the xyz file
row = length(unique(x));  
col = length(unique(y));
% Reshape and transpose z values to get correct data orientaion
zm = reshape(z,row,col);
zm = fliplr(zm);
zm = zm';
% Perform fft2 on z-values (nT)
fft2_zm = fftshift(log(abs(fft2(zm))));  % amplitude
% phase_zm = angle(fftshift(fft2(zm)));  % phase
% % Calculate wavenumbers
% Nc = col; Nr = row;
% dp = 1/300; % dpi = 300?
% kx = (-Nr/2:Nr/2-1)/(Nr*dp);
% ky = (-Nc/2:Nc/2-1)/(Nc*dp);
figure(1); 
%imagesc(fft2_zm); colorbar;
x_axis = linspace(-pi,pi,size(fft2_zm,1));
y_axis = linspace(-pi,pi,size(fft2_zm,2));
imagesc('XData',x_axis,'YData',y_axis,'CData',fft2_zm);
xlim([x_axis(1),x_axis(end)]);
ylim([y_axis(1),y_axis(end)]);
xlabel('Rad');
ylabel('Rad');
ax = gca;
ax.XTick = [-pi, -pi/2, 0, pi/2, pi];
ax.YTick = [-pi, -pi/2, 0, pi/2, pi];
ax.XTickLabel = {"-\pi", "-\pi/2", "0","\pi/2", "/pi"};
ax.YTickLabel = {"-\pi", "-\pi/2", "0","\pi/2", "/pi"};

[Jay Ghosh]

@Hiro Yoshino, Thank you very mcuh for the explanation. I am going to go with your suggestion. However, I still have an unanswered question that is related. I will be posting the script and the question below. (as the previous post was removed).

The following is a script to read an image file (lena512.bmp), and apply fft2 on the grayscale version of it. For the fft2 of the image (subplot 2), how would you label the axes, in terms of wavenumber or spatial frequency? In another word, for fft of an image how would you labe the axes. Thank you so much, in advance.

%% fft of Lena512
lena512rgb = imread('Lena512.bmp');
lena512gray = rgb2gray(lena512rgb);
fft2_lena512 = fftshift(log(abs(fft2(lena512gray))));

figure(2)
subplot (1,2,1)
imagesc(lena512gray)
subplot (1,2,2)
imagesc(fft2_lena512)

[Hiro Yoshino]

@Jay Ghosh

You should use "tiledlayout" rather than subplot and the same principle is applied to this too:

lena512rgb = imread('Lena512.bmp');
lena512gray = rgb2gray(lena512rgb);
fft2_lena512 = fftshift(log(abs(fft2(lena512gray))));

t = tiledlayout(1,2);
nexttile
imagesc(lena512gray)
xlabel("X");
ylabel("Y");

nexttile
imagesc(fft2_lena512)
xlabel("Normalized Frequency X");
ylabel("Normalized Frequency Y");
title(t,"Lena FFT");

Also, when you want to convert the normalized frequency into spacial frequency you need to use the sampling frequency at which you sampled the data points from the picture.

For example. if your pixel is sampled at every 2mm, your sampling period T = 2 [mm/Sample].

Your normalized frequency has the unit of w [rad/Sample]. So if you work out w/T [rad/mm] then this is so-called the spacial frequency.

[Jay Ghosh]

@Hiro Yoshino, For this image (lena512.bmp), how would I calculate the sampling period? The image dimension is 512x512. So, does this say anything about the pixel sampling frequency? If not, how do I calculate it? May be it is simple, but I am not able to understand it. Thanks.

[Hiro Yoshino]

@Jay Ghosh

Yes. Sampling frequency = 1/(Sampling Period) is hold.

You would never know what the sampling period is unless you know the actual size of the picture. So 512 x 512 is 512 x 512, that's it

[Jay Ghosh]

@Hiro Yoshino, Thank you!

In your first response/comment, you said, 'the number of points along the x-axis and the y-axis corresponds to the wavenumber'. So, if I want to express the axes as wavenumbers, could I do it as following (where wavenumbers range from -250 to 250)? Thanks.

%% fft of Lena512
lena512rgb = imread('Lena512.bmp');
lena512gray = rgb2gray(lena512rgb);
fft2_lena512 = fftshift(log(abs(fft2(lena512gray))));

% Calculate wavenumbers
Nr = length(lena512gray);
Nc = Nr;
dp = 1/512;
kx = (-Nr/2:Nr/2-1)/(Nr*dp);
ky = (-Nc/2:Nc/2-1)/(Nc*dp);

t = tiledlayout(1,2);
nexttile
imagesc(lena512gray)

nexttile
imagesc(kx,ky, fft2_lena512)
xlabel("Wavenumber (kx)");
ylabel("Wavenumber (ky)");
title(t,"Lena FFT");

[Hiro Yoshino]

If you consider 2pi/N as the base frequency, the way you show wavenumber in the figure is right.

[Jay Ghosh]

Thank you @Hiro Yoshino. I have a follow up question/comment. Why is it that here we see that the highest energy corrsponds to the lowest spatial frequency? I thought the higher the energy in kx-ky space, the higher the frequency. In another word, high energy is associated with high frequency. But, for both the transforations above (tr86 and lena512), we see that the highest energy (in the centre) is associated with lowest frequency. Thanks.

[Hiro Yoshino]

@Jay Ghosh

The power spectrum we got is typical. We often see the highest power at lower frequencies.

Unless we have specific high frequencies in the data, we normally observe similar powe spectrums.

The energy at 0 frequency corresponds to the bias. So you can remove it by subtracting the mean value from the data.