I just tried to play with k-Wave using different parameters and see what happened. Then when I was simulating a linear medium with a small scatter, which has a radius of 1 (i know this small radius does not make much sense and is not practical), I got some strange sensor data.
Following is my code:
'%% define medium, scatter and grid points
clc;
clearvars;
% set the size of the perfectly matched layer (PML)
pml_x_size = 10; % [grid points]
pml_y_size = 10; % [grid points]
pml_z_size = 10; % [grid points]
% set total number of grid points not including the PML
Nx = 190 - 2*pml_x_size; % [grid points]
Ny = 100 - 2*pml_y_size; % [grid points]
Nz = 600 - 2*pml_z_size; % [grid points]
% define the properties of the propagation medium
c0=1500;
rho0=1040;
tone_burst_freq = 3e6;
tone_burst_cycles = 4;
source_strength = 1e6; % [Pa]
% at least 3 points per wavelength is recommended
dx = c0/tone_burst_freq/4; % grid point spacing in the x direction [m]
dy = dx;
dz = dx;
kgrid = kWaveGrid(Nz,dz,Nx, dx, Ny, dy);
kgrid.makeTime(c0, 0.3, (Nz * dz) * 2 / c0);
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
input_signal = (source_strength ./ (c0 * rho0)) .* input_signal;
%% Create The Transducer
% physical properties of the transducer
transducer.number_elements = 80; % total number of transducer elements
transducer.element_width = round(0.23e-3/dx); % width of each element [grid points]
transducer.element_length = round(6.5e-3/dx); % length of each element [grid points]
transducer.element_spacing = round(0.04e-3/dx); % 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([1, Nx/2-transducer_width/2, Ny/2-transducer.element_length/2]);
%transducer_central.position = [z x y];
% properties used to derive the beamforming delays
transducer.sound_speed = c0; % sound speed [m/s]
transducer.focus_distance = inf; % focus distance [m]
transducer.elevation_focus_distance = 77e-3; % 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
number_active_elements = transducer.number_elements;
transducer.active_elements = ones(transducer.number_elements, 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 properties
sound_speed_map = c0 * ones(Nz, Nx, Ny) ;
density_map = rho0 * ones(Nz, Nx, Ny) ;
scatterRadius = 1;
scatterCoef = 1.06;
% define a sphere for a highly scattering region
scattering_region1 = makeBall(Nz, Nx, Ny, 4*Nz/5, ...
Nx/2, Ny/2, scatterRadius(1));
% assign region
sound_speed_map(scattering_region1 == 1) = c0*scatterCoef(1);
density_map(scattering_region1 == 1) = rho0*scatterCoef(1);
medium.sound_speed = sound_speed_map;
medium.density = density_map;
voxelPlot(single(transducer.mask | scattering_region1));
view(142,25);
%% simulation
input_args = {...
'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size], ...
'DataCast', 'gpuArray-single', 'DataRecast', true};
% run the simulation
tic
sensor_data = kspaceFirstOrder3DG(kgrid, medium, transducer, transducer, input_args{:});
toc
%% sensor_data
sensor_data_com = transducer.combine_sensor_data(sensor_data);
data_win = getWin(kgrid.Nt * 2, 'Tukey', 'Param', 0.05).';
data_win = [zeros(1, ceil(length(transducer.input_signal) * 2)), data_win(1:end/2 - ceil(length(transducer.input_signal) * 2))];
figure;
subplot(211);
plot(sensor_data_com(40,:));title('before windowing');
subplot(212);
sensor_data_com = bsxfun(@times, data_win, sensor_data_com);
plot(sensor_data_com(40,:));title('after windowing');'
If you have tried the code above, you may have already seen the sensor_data after the windowing operation, some strong signals appeared in the near field and around the half of the time series, but there were no scatter points, can you please tell me is there anything wrong with my code or is it hard to simulate such a small scatter in k-Wave? If it is the latter, why is it?
And I noticed that if the scan_line() method is used, I will get one smooth sensor_data plot, does anyone know what operation was done beneath scan_line() function?
Thank you for your responses in advance.
Kind regards,
Bing
