Hi all,
I am trying to simulate a real transducer, I came across an artifact that I am not sure where it is coming from and whether it could be solved or it is just a physical limitation. I could not find my answer in the forum or the manual. The artifact seems to be some sort of cross talk. It would be a great help if you could share your insights.
I am trying to simulate a 128 element transducer (64 elements of the transducer are activated during TX and RX), transducer width is 37.5 mm therefore only the 64 middle elements are active(18.75mm). I use the transducer class and CUDA code. My intention is to transmit a plane wave trough a medium and receive the echoes with same transducer. My grid size is 1024*1024*64 with pml size 64*64*16(outside).
After extracting the scan-lines using transducer.combine_sensor_data(sensor_data) and ploting it, the area near transducer receives what appears to me as a horizontal wave propagated from the side elements. When I remove the kerf the artifact to some extend disappears (not completely though and I need kerf)
-What can possibly cause such artifact? (e.g. can this be a result of small number of grid points in z directions?(I have restriction in computational resources right now))
-Can I avoid it in any way (increasing pml size, pml alpha, placing the transducer inside a highly absorbing material, increase or decrease dt, etc.)
-What are the acoustic properties of the transducer and the kerf?
-As the grid is too large, I would like to increase time step but the artifact seems to be even worse when I do so (I use t_min_stable = checkStability(kgrid,medium) and set t_end to this value, changing cfl value doesn't help too). How can I increase time step (or decreasing frequency or in general any solution that results in smaller data size)
my code is as follows:
% simulation settings
DATA_CAST = 'gpuArray-single';
RUN_SIMULATION = true;
f_max = 5e6 ;
% set the size of the perfectly matched layer (PML)
pml_x_size = 64; % [grid points]
pml_y_size = 64; % [grid points]
pml_z_size = 16 ;
pml_x_alpha = 3;
pml_y_alpha = 3 ;
pml_z_alpha = 3 ;
%total size
x_size = 4.24e-2 ; %[m]
y_size = 4.24e-2 ;
dy = 0.0375/(128*4+127*4 ); %[m]
dx = dy ; %
dz = dy ;
% Ny = y_size/dy grid points
Ny = 1024 ; %ceil(y_size/dy)-2*pml_x_size ;
Nx = 1024 ; %ceil(x_size/dx)-2*pml_y_size ;
Nz = 32;
% create the k-space grid
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================
% define the properties of the propagation medium
medium.sound_speed = 1540 ;
rho0 = 900; % [kg/m^3]
medium.alpha_coeff = 0.5 ;
medium.alpha_power = 1.05 ;
% create the time array
c0_min = 1300 ; %m/s
t_end = (Nx * dx) * 2.2 / c0_min; % [s]
kgrid.makeTime(c0_min,0.3, t_end);
% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% define properties of the input signal
source_strength = 1e6; % [Pa]
sampling_freq = 20e6 ; % [Hz]
tone_burst_freq = 5e6; % [Hz]
tone_burst_cycles = 2;
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
input_signal = (source_strength ./ (c0_min * rho0)) .* input_signal;
% =========================================================================
% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================
% physical properties of the transducer
transducer.number_elements = 128; % total number of transducer elements
transducer.element_width = 4; % width of each element [grid points]
transducer.element_length =28; % length of each element [grid points]
transducer.element_spacing = 4; % spacing (kerf width) between the elements [grid points]
transducer.radius = inf; % radius of curvature of the transducer [m]
% calculate the width of the transducer in grid points
transducer_width = transducer.number_elements * transducer.element_width ...
+ (transducer.number_elements - 1) * transducer.element_spacing;
% use this to position the transducer in the middle of the computational grid
transducer.position = round([2, 2, 1]);
% properties used to derive the beamforming delays
transducer.sound_speed = c0_min; % sound speed [m/s]
transducer.focus_distance = inf; % focus distance [m]
transducer.elevation_focus_distance = inf; % focus distance in the elevation plane [m]
transducer.steering_angle = 0; % steering angle [degrees]
% apodization
transducer.transmit_apodization = 'Hanning';
transducer.receive_apodization = 'Hanning';
% define the transducer elements that are currently active
transducer.active_elements = zeros(transducer.number_elements, 1);
transducer.active_elements(33:96) = 1;
% append input signal used to drive the transducer
transducer.input_signal = input_signal;
% create the transducer using the defined settings
transducer = kWaveTransducer(kgrid, transducer);
% print out transducer properties
transducer.properties;
% =========================================================================
% DEFINE THE MEDIUM PROPERTIES
% =========================================================================
% define a large image size to move across
Nx_tot = Nx;
Ny_tot = Ny ;
Nz_tot = Nz;
% define a random distribution of scatterers for the medium
background_map_mean = 1;
background_map = background_map_mean ;
scattering_c0 = c0_min ;
scattering_rho0 = scattering_c0 / 0.9 ;
% define properties
sound_speed_map = c0_min * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;
density_map = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;
medium.sound_speed = sound_speed_map(:,:, :);
medium.density = density_map(:, :, :);
% % =========================================================================
% % RUN THE SIMULATION
% % =========================================================================
% set the input settings
input_args = {...
'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size], ...
'PMLAlpha', [pml_x_alpha, pml_y_alpha, pml_z_alpha], 'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false};
% save input files to disk
% run the simulation
[sensor_data] = kspaceFirstOrder3DG(kgrid, medium, transducer, transducer, input_args{:});
scan_lines = transducer.combine_sensor_data(sensor_data) ;
%
save sensor_data
%
It would be really helpful if you can share your insights.
Thanks in advance,
Farnaz
