2D FFT Frequency-wavenumber Dispersion Curves

I would like to generate a frequency-wavenumber plot (2D FFT) from a series of time signals, in a similar manner to this post: https://uk.mathworks.com/matlabcentral/answers/342127-spatial-fft-of-a-series-of-wave-signals?s_tid=srchtitle

However, this code is not quite complete and I am struggling to understand what I need to do to finish it off.

I am exciting both A0 and S0 at 1 MHz, with a hamming windowed pulse, measuring y-axis displacement at 10 points, spaced 1mm apart. On a wavenumber-frequency plot I expect to see A0 at around 2.7 rad/mm, and S0 at around 1.2 rad/mm (both at 1 MHz).

• How can I get my wavenumber axis into these units?
• What should I do to increase wavenumber resolution? Increase the number of spatial measurement points while placing them closer together? Do I also need to increase sampling rate to improve frequency resolution?
• How do I plot the result in such a way that both frequency and wavenumber are positive values?

Code (data attached):

myfile = '10responsesV2.csv';
M1 = readtable(myfile,'ReadVariableNames', false);
A1 = table2array(M1);
A1(1,:) = [];

Fs = 1/((A1(2,1))-(A1(1,1)));
St = 1/Fs;     % Sampling period

Nx = 10; % Number of samples collected along first dimension
Nt = 28000; % Number of samples collected along second dimension
table = A1(:,2:11);
dx = 0.001;  % Distance increment (i.e., Spacing between each column)
dt = St; % Time increment (i.e., Spacing between each row)
x = 0 : dx : (Nx-1)*dx; % distance
t = 0 : dt : (Nt-1)*dt; % time
Nyq_k = 1/dx; % Nyquist of data in first dimension
Nyq_f = 1/dt; % Nyquist of data in second dimension
dk = 1/(Nx*dx);   % Wavenumber increment
df = 1/(Nt*dt);   % Frequency increment
k = -Nyq_k/2:dk:Nyq_k/2-dk; % wavenumber
f = -Nyq_f/2:df:Nyq_f/2-df; % frequency

figure(1);
fft2result = fftshift(fft2(table))*dx*dt;
imagesc(k,f,abs(fft2result));
colorbar;

Answer 1

[Brahmadev]

Hi @Lawrence Yule,

To address the questions you have mentioned:

1. Getting wavenumber axis into the correct units: The wavenumber k is typically in units of radians per unit length and it seems you have the correct units already, radians per millimeter in your case. However, when you create the wavenumber vector k, you should ensure that it is centered around zero to represent both positive and negative wavenumbers correctly. The "fftshift" function is used for this purpose.
2. Increasing wavenumber resolution: To increase the wavenumber resolution, you can either increase the number of spatial measurement points (Nx) or decrease the distance between them (dx). Both will result in a larger spatial aperture and thus a finer "dk".
3. Improving frequency resolution: Increasing the total measurement time (Nt*dt) will improve the frequency resolution, as df will become smaller. This could mean collecting data for a longer time period or increasing the number of temporal samples (Nt) while keeping the sampling period (dt) constant.
4. Plotting positive frequency and wavenumber values: To plot only positive values of frequency and wavenumber, you would typically use only the second half of the f and k vectors after performing fftshift.

Here is how you could modify your code to address these points:

% Read the data from the CSV file
myfile = '10responsesV2.csv';
M1 = readtable(myfile, 'ReadVariableNames', false);
A1 = table2array(M1);
A1(1,:) = []; % Assuming the first row is not part of the data

% Define the sampling frequency and other parameters
Fs = 1/((A1(2,1))-(A1(1,1)));
St = 1/Fs;     % Sampling period
Nx = 10;       % Number of samples collected along the x-axis (spatial)
Nt = size(A1, 1);    % Number of samples collected along the t-axis (temporal)
table = A1(:,2:11);
dx = 0.001;    % Distance increment in meters (1 mm)
dt = St;       % Time increment

% Define Nyquist limits
Nyq_k = pi/dx; % Nyquist wavenumber in rad/mm
Nyq_f = Fs/2;  % Nyquist frequency in Hz

% Calculate the wavenumber and frequency increments
dk = 2*pi/(Nx*dx);   % Wavenumber increment in rad/mm
df = Fs/Nt;          % Frequency increment in Hz

% Create wavenumber and frequency vectors
k = -Nyq_k:dk:Nyq_k-dk; % Wavenumber vector in rad/mm
f = -Nyq_f:df:Nyq_f-df; % Frequency vector in Hz

% Perform 2D FFT and shift zero frequency/wavenumber components to the center
fft2result = fftshift(fft2(table, Nt, Nx))*dx*dt;

% Plot the magnitude spectrum
figure(1);
imagesc(k, f, 20*log10(abs(fft2result)));
xlabel('Wavenumber (rad/mm)');
ylabel('Frequency (Hz)');
colorbar;
title('Magnitude Spectrum (dB)');

% Adjust color axis limits if necessary
caxis([min(20*log10(abs(fft2result(:)))), max(20*log10(abs(fft2result(:))))]);

% Ensure the axes are scaled correctly
axis xy;
axis tight;

Hope this helps in resolving your query!