Hello Forum,
I am trying to perform 2D photoacoustic simulation of point target sources in a homogeneous medium. To record sensor data, a 128-element linear array with center frequency of 21 MHz and 78% fractional bandwidth is being used. Delay-and-sum beamforming is used for image reconstruction. I have noticed that K-wave allows us to set two sampling frequencies - spatial and temporal and changing either of them changes the reconstruction results. Higher spatial and temporal frequcy seems to provide better results (image enclosed here: https://drive.google.com/file/d/1hD-wVWQgKJN7id1qKPJQ0_My4JnVbg23/view?usp=sharing). Can anyone help me to understand why this is happening? I am also providing my code below. Thank you very much for the help.
%% Linear Array Photoacoustic Imaging for Point Initial Source
% Define Probe Properties
pitch = 90e-6; % pitch
TransducerElem = 128; % Number of transducer elements
f_cen = 21e6; % Transducer center frequency
f_BW = 78; % Transducer Bandwidth
imaging_aperture = pitch*TransducerElem; % aperture length for single quadrant imaging
c0 = 1540;
lambda=c0/f_cen; % wavelengsth [m]
% Define K-space spatial grids
Ny = 2500;
% Ny = 256;
Nx = TransducerElem;
imagingdepth = 25e-3; % [m]
imagingwidth = imaging_aperture; % [m]
dy = imagingdepth/Ny; % grid point spacing in the y direction [m] --> Simulation Depth
dx = imagingwidth/Nx; % grid point spacing in the x direction [m] --> Aperture is placed in this direction
PML_size=[20,200]; % size of the PML in grid points
kgrid = kWaveGrid(Nx, dx, Ny, dy);
% Define Medium properties
medium.density = 1000;
medium.sound_speed = 1540; % [m/s]
% medium.alpha_coeff = 0.3; % [dB/(MHz^y cm)]
% medium.alpha_power = 1.5;
% medium.alpha_mode = 'no_dispersion';
% Define K-wave Temporal Grid
f_s = 42e6*2;
dt = 1/f_s; % 1/sampling frequency
dt_stability_limit = checkStability(kgrid, medium);
while dt>dt_stability_limit
f_s = f_s*2;
dt = 1/f_s;
end
Nt = ceil((imagingdepth*f_s)/c0); % time steps to simulate
kgrid.t_array = dt*(0:Nt);
% Define K-wave sensor
sensor.mask = zeros(Nx, Ny);
sensor.mask(:,1) = 1;
sensor.frequency_response = [21e6,78]; % [center frequency = 21Mz with 78% fractional BW]
% Define initial pressure distrbution
rad = 0.0001/dx; % radius = 0.1mm
mag = 3; % magnitude in PA
x_index = [round(Nx/2-Nx/6) round(Nx/2) round(Nx/2+Nx/6)];
xPos = kgrid.x_vec(x_index) + kgrid.x_vec(end); % x-location of point targets [m]
yPos = [8e-3 12e-3 16e-3 20e-3]; % y-location of point targets [m]
p0sum = zeros(Nx,Ny);
for j = 1:length(xPos)
x = round(xPos(j)/dx);
for k=1:length(yPos)
y = round(yPos(k)/dy);
disc = mag*makeDisc(Nx, Ny, x, y, rad);
p0 = smooth(kgrid, disc, true);
p0sum = p0sum + p0;
end
end
%% Run the Simualtion on GPU
%
input_args = {'PMLSize',PML_size,'PMLAlpha',4,'PlotPML',false,'PMLInside', false,'Smooth', false,'DataCast','gpuArray-single', 'PlotSim', true,'PlotLayout',false};
sensorData = kspaceFirstOrder2D(kgrid, medium, ...
struct('p0',p0sum), sensor, input_args{:});
sensorData = gather(sensorData);