Hello,
I am trying to generate simulated pulse-echos created/detected by a single element concave transducer in a heterogeneous medium for purposes of attempting different methods of estimating attenuation coefficient. As the simulation seems a bit generic, could someone point me towards an example program if it already exists?
I've created a script to run this simulation, but the signals are rather damped despite making many adjustments. So I was curious if someone could perhaps indicate something I'm not doing quite right. Any tips would be most appreciated!
% Focussed Detector In 3D Example
%
% This example shows how k-Wave can be used to model the output of a
% focussed bowl detector in 3D where the directionality arises from
% spatially averaging across the detector surface. It builds on the
% Focussed Detector In 2D Example.
%
clearvars;
% =========================================================================
% SIMULATION
% =========================================================================
% create the computational grid
% dx, dy, dz dictate the size of the grid, as their units is meters. i.e.
% the size of the grid is [1:Nx]*dx in the depth dimension
Nx = 256; %10*64; % number of grid points in the x direction
Ny = 64; %64; % number of grid points in the y direction
Nz = 64; % number of grid points in the z direction
res = 1e-4;
dx = res; %10e-3/Nx; % grid point spacing in the x direction [m]
dy = res; % grid point spacing in the y direction [m]
dz = res; % grid point spacing in the z direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
%%
% define the properties of the propagation medium
c0 = 1540; % [m/s]
rho0 = 1000; % [kg/m^3]
medium.alpha_coeff = 0.5; % [dB/(MHz-cm)]
medium.alpha_power = 1.01;
medium.BonA = 6; % parameter of nonlinearity
medium.sound_speed = c0; % [m/s]
% define a random distribution of scatterers for the medium
background_map_mean = 1;
background_map_std = 0.008;
background_map = background_map_mean + background_map_std * randn([Nx, Ny, Nz]);
% define properties
sound_speed_map = c0 * ones(Nx, Ny, Nz) .* background_map;
density_map = rho0 * ones(Nx, Ny, Nz) .* background_map;
%%
% define the properties of the propagation medium
% medium.sound_speed = 1500; % [m/s]
% create the time array
kgrid.makeTime(medium.sound_speed);
% create a concave sensor
sphere_offset = 0; % [grid points]
diameter = Ny/4 + 1; % [grid points]
radius = 100*Nx; %Nx/2; % radius of curvature [grid points]
bowl_pos = [1 + sphere_offset, Ny/2, Nz/2]; % [grid points]
% x is depth dimension
focus_pos = [100*Nx, Ny/2, Nz/2]; %[Nx/2, Ny/2, Nz/2]; % [grid points]
sensor.mask = makeBowl([Nx, Ny, Nz], bowl_pos, radius, diameter, focus_pos);
% % define a time varying sinusoidal source
% source_freq = 0.25e6; % [Hz]
% source_mag = 1; % [Pa]
% source.p = source_mag * sin(2 * pi * source_freq * kgrid.t_array);
% source.p = filterTimeSeries(kgrid, medium, source.p);
%
%%
% define properties of the input signal
source_strength = 100e6; % [Pa]
tone_burst_freq = 5e6;%0.25e6; %1.5e6; % [Hz]
tone_burst_cycles = 4;
% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
% scale the source magnitude by the source_strength divided by the
% impedance (the source is assigned to the particle velocity)
input_signal = (source_strength ./ (c0 * rho0)) .* input_signal;
% filter the source to remove any high frequencies not supported by the grid
source.p = filterTimeSeries(kgrid, medium, input_signal);
% give sensor finite bandwidth
bandwidth_pct = 50;
sensor.frequency_response = [tone_burst_freq, bandwidth_pct];
source.p_mask = sensor.mask; %source1;
input_args = {'PMLSize', 30, 'DataCast', 'off', 'PlotSim', false, ...
'DataRecast', true};
nIterations = 2;
for nnn = 1:nIterations
% shuffle scatterers to simulate another instance of media with the
% same properties.
id1 = randperm(size(sound_speed_map,1));
id2 = randperm(size(sound_speed_map,2));
id3 = randperm(size(sound_speed_map,3));
sound_speed_map_shuffled = sound_speed_map(id1, id2, id3);
id1 = randperm(size(density_map,1));
id2 = randperm(size(density_map,2));
id3 = randperm(size(density_map,3));
density_map_shuffled = density_map(id1, id2, id3);
medium.sound_speed = sound_speed_map_shuffled;
medium.density = density_map_shuffled;
%% run the first simulation
sensor_data1 = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
% average the data recorded at each grid point to simulate the measured
% signal from a single element focussed detector
sensor_data1_sum(:, nnn) = sum(sensor_data1, 1);
end
% =========================================================================
% VISUALISATION
% =========================================================================
% plot the detector and on-axis and off-axis point sources
% voxelPlot(sensor.mask + source1);
voxelPlot(sensor.mask)
view([130, 40]);
% get a suitable scaling factor for the time axis
[~, t_scale, t_prefix] = scaleSI(kgrid.t_array(end));
% plot the time series corresponding to the on-axis and off-axis sources
figure
plot(kgrid.t_array * t_scale, sensor_data1_sum, '-');
% hold on
% plot(kgrid.t_array * t_scale, sensor_data2, 'r-');
xlabel(['Time [' t_prefix 's]']);
ylabel('Average Pressure Measured By Focussed Detector [Pa]');
% legend('Source on axis', 'Source off axis');
%save