2D FFT Frequency-wavenumber Dispersion Curves

Question

Lawrence Yule 2022-4-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1703815-2d-fft-frequency-wavenumber-dispersion-curves

回答： Brahmadev 2024-2-9

10responsesV2.zip

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;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Brahmadev 2024-2-9

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1703815-2d-fft-frequency-wavenumber-dispersion-curves#answer_1405596

在 MATLAB Online 中打开

Hi @Lawrence Yule,

To address the questions you have mentioned:

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.
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".
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.
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;