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;
Hope this helps in resolving your query!
