Hi,
I’m trying to determine the medium property (speed of sound) based on the measured Time-of-Flight. When I started simulating some cases in k-wave I noticed, that the results I get differ from the theory especially close to the source. I’m using a time varying source and dirichlet boundary conditions. The magnitude of the error seems to depend on the grid resolution (Nx) and is independent from the cfl-number. Can anyone tell me what might be the cause?
Thanks
Michael
Here is a simple case that demonstrates the problem:
clearvars;
% =========================================================================
% SIMULATION
% =========================================================================
% create the computational grid
Nx = 512; % number of grid points in the x (row) direction
Ny = 64; % number of grid points in the y (column) direction
dx = 15/Nx; % grid point spacing in the x direction [m]
dy = dx; % grid point spacing in the y direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy);
% define the properties of the propagation medium
medium.sound_speed = 1500; % [m/s]
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
cfl = 0.1;
endtime = 0.013; % [s]
% create the time array
kgrid.makeTime(medium.sound_speed, cfl, endtime);
% define a row of source points
source.p_mask = zeros(Nx, Ny);
source.p_mask(end - Nx/4, Ny/2) = 1;
% define the properties of the tone burst used to drive the transducer
sampling_freq = 1/kgrid.dt; % [Hz]
tone_burst_freq = 800; % [Hz]
tone_burst_cycles = 6;
source_mag = 2; % [Pa]
% create the tone burst signals
source.p = source_mag * toneBurst(sampling_freq, tone_burst_freq, tone_burst_cycles);
padsz = length(kgrid.t_array)-length(source.p);
source.p = padarray(source.p, [0 padsz], 0, 'post');
% filter the source to remove high frequencies not supported by the grid
source.p = filterTimeSeries(kgrid, medium, source.p);
source.p_mode = 'dirichlet';
% define a single sensor point
sensor.mask = zeros(Nx, Ny);
sensor.mask(Nx/4:end - Nx/4, Ny/2) = 1;
% define the acoustic parameters to record
sensor.record = {'p', 'p_final'};
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);
%% calculate Time of Flight between two sensor points
[NS, ~] = size(sensor_data.p);
tof = zeros(1, NS);
x2 = source.p; % source signal
for j=1:NS
x1 = squeeze(sensor_data.p(j, :)); % receiver signal
[y1, lags] = xcorr(x1,x2);
[~,idx] = max(y1);
tof(j) = lags(idx) * kgrid.dt;
end
dist = (kgrid.x(logical(source.p_mask))-kgrid.x(logical(sensor.mask)))';
c = dist./tof; % mean speed of sound between the source an the sensor
c(end) = nan;
tofn = tof/kgrid.dt; % Normalized ToF
deltatn = -diff(tofn); % Number of Timesteps it takes the signal to get from one grid cell to the next
% theoretical values
c_theo = ones(1, length(c))* medium.sound_speed;
deltatn_theo = ones(1, length(deltatn))*(kgrid.dx / (medium.sound_speed * kgrid.dt));
% =========================================================================
%% VISUALISATION
% =========================================================================
% plot calulated speed of Sound
figure;
[~, scale, prefix] = scaleSI(max(dist(:)));
subplot(2, 1, 1);
plot(dist * scale, c, 'k-', 'DisplayName', 'simulation');
hold on;
plot(dist * scale, c_theo, 'r-', 'DisplayName', 'theoretical');
hold off;
xlabel(['distance [' prefix 'm]']);
ylabel('speed of sound [m/s]');
axis tight;
title('effective SoS between source and sensor');
legend('simulation','theoretical');
subplot(2, 1, 2);
plot(dist * scale, [deltatn 0], 'k-', 'DisplayName', 'simulation');
hold on;
plot(dist * scale, [deltatn_theo 0], 'r-', 'DisplayName', 'theoretical');
hold off;
xlabel(['distance [' prefix 'm]']);
ylabel('number of time steps');
axis tight;
title('number of time steps for each grid step');
legend('simulation','theoretical');