Hello. I am trying to run an ultrasound ring array simulation with a heterogeneous speed of sound medium and a time-varying pressure source. I can see the original transmitted signal in all sensor channels, but I am not seeing the reflection signals off of the object I have embedded in my medium. I cannot see where the mistake is. Any help would be greatly appreciated.
___________________________________________________________
load TransducerPositions_1024.txt
TransPos=TransducerPositions_1024(1:4:end,:); %Convert 1024 element array to 256
%===================================================
% User Input
%===================================================
% Select the simulation type
%
% 1. C++ Simulation from Matlab
% 2. CUDA Simulation from Matlab
% 3. MATLAB simulation
%
SIMULATION_TYPE = 1;
SIGNAL_AMPLITUDE = 40; % [Pa]
SIGNAL_BANDWIDTH = 0.80; % [ratio], e.g. 50% --> 0.50
SIGNAL_CENTER_FREQUENCY = 1.5; % [MHz]
SIGNAL_SAMPLING_FREQUENCY = 12; % [MHz]
%=================================
% Computational Grid
%=================================
Nx = 256; % Number of grid points in the x direction
Ny = 256; % Number of grid points in the y direction
Nz = 30; % Number of grid points in the z direction
dx = 1e-3; % Grid point spacing in the x direction [m]
dy = 1e-3; % Grid point spacing in the y direction [m]
dz = 1e-3; % Grid point spacing in the z direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
%=================================
% Medium
%=================================
%create cylinder embedded in medium
cx=128; cy=128; cz=3;
radius=50;
height=25;
medium.sound_speed=makeCylinder(Nx, Ny, Nz, cx, cy, cz, radius, height);
%create sphere inside cylinder inside medium
cx=128; cy=128; cz=15;
radius=5;
medium.sound_speed=medium.sound_speed+makeSphere(Nx, Ny, Nz, cx, cy, cz, radius);
medium.sound_speed(medium.sound_speed==0)=1500; %water bath sound speed
medium.sound_speed(medium.sound_speed==1)=1450; %cylinder object sound speed
medium.sound_speed(medium.sound_speed==2)=1600; %sphere object sound speed
medium.density=1000*ones(Nx,Ny,Nz);
%%
% Init time array
dt = 1/(SIGNAL_SAMPLING_FREQUENCY*10^6);
c0_peak = max(max(max(medium.sound_speed)));
CFL = (dt/dx)*c0_peak;
if CFL>0.3
disp("WARNING: CLF number is larger than 0.3, which could result in numerical errors!");
end
kgrid.makeTime(medium.sound_speed, CFL);
%=================================
% Source
%=================================
tc = gauspuls('cutoff',SIGNAL_CENTER_FREQUENCY*10^6 ...
,SIGNAL_BANDWIDTH,[],-40);
t = -tc : dt : tc;
signal = SIGNAL_AMPLITUDE*gauspuls(t,SIGNAL_CENTER_FREQUENCY*10^6,SIGNAL_BANDWIDTH);
p0 = zeros(size(kgrid.t_array));
p0(1:length(signal)) = signal;
source.p = p0;
%=================================
% Sensors
%=================================
ring.position = ceil((TransPos*1e-3./dx - 0.5));
ring.position(:,1) = ring.position(:,1) + Nx/2;
ring.position(:,2) = ring.position(:,2) + Ny/2;
% define a sensor mask through the central z-plane
sensor.mask = zeros(Nx, Ny, Nz);
for idx = 1:length(ring.position)
sensor.mask(ring.position(idx,1),ring.position(idx,2), ceil(Nz/2)) = 1;
end
%===================================================
% Run Simulation
%===================================================
input_args = {};
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(ring.position(eleNum,1),ring.position(eleNum,2),Nz/2) = 1;
switch SIMULATION_TYPE
case 1
% Run the C++ simulation
sensor_data = kspaceFirstOrder3DC(kgrid, medium, source, sensor, input_args{:});
case 2
% Run the CUDA simulation
sensor_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:});
case 3
% Run the MATLAB simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
end