Hi,
I am trying to model the reflection and refraction of ultrasound at a interface between soft tissue and air. The sound speed and density of air are 340 m/s and 1.29 Kg/m3 respectively. the size of the interface is much larger than the wavelength of ultrasound source. ultrasound frequency is on the order of 10 MHz.
However, the ultrasound echo received by the linear transducer is a matrix of NaN (not a number). The Matlab code is shown below. I don't know why it didn't work out.
clc
clear all
close all
% simulation settings
DATA_CAST = 'double'; % set to 'single' or 'gpuArray-single' to speed up computations
%%
========================================================================
=
% DEFINE THE K-WAVE GRID
%
========================================================================
=
% compute grid spacing based on a desired x_size and f_max
x_size = 50e-3; % [m] desired domain size on x direction (depth...)
c0_min = 1500; % [m/s]
f_max = 15e6; % [Hz] maximum frequency of ultrasound interested
points_per_wavelength = 2;
dx = c0_min/(points_per_wavelength*f_max); % [m]
dy = dx; % [m]
dz = dx; % [m]
pml_x_size = 20; % [grid points]
pml_y_size = 20; % [grid points]
pml_z_size = 20; % [grid points]
Nx_at_smallest_prime_factor = 256;
Ny_at_smallest_prime_factor = 128;
Nz_at_smallest_prime_factor = 128;
% set number of grid points not including the PML
Nx = Nx_at_smallest_prime_factor - 2*pml_x_size; % [grid points]
Ny = Ny_at_smallest_prime_factor - 2*pml_y_size; % [grid points]
Nz = Nz_at_smallest_prime_factor - 2*pml_z_size; % [grid points]
% create the k-space grid
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
% create the time array
t_end = (Nx*dx)*1.5/c0_min; % [s]
kgrid.t_array = makeTime(kgrid, c0_min, [], t_end);
%%
========================================================================
=
% DEFINE THE MEDIUM PARAMETERS
%
========================================================================
=
% define the properties of the air
c0_air = 340; % [m/s]
rho0_air = 1.29; % [kg/m^3]
% medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
% medium.alpha_power = 1.5;
% medium.BonA = 6;
medium.sound_speed = c0_air*ones(Nx, Ny, Nz);
medium.density = rho0_air*ones(Nx, Ny, Nz);
% define the properties of the soft tissue
c0_tissue = 1540; % [m/s]
rho0_tissue = 1000; % [kg/m^3]
std = 0.005;
speed_map_tissue = c0_tissue*(ones(Nx, Ny, Nz) + std*randn(Nx, Ny, Nz));
density_map_tissue = rho0_tissue*(ones(Nx, Ny, Nz) + std*randn(Nx, Ny, Nz));
% define the region of tissue
region_tissue = zeros(Nx, Ny, Nz);
radius = Nx/2 -10;
for i = 1:Nz
region_tissue(:,:,i) = makeDisc(Nx, Ny, 1, Ny/2, radius);
end
% define the properties of tissue
medium.sound_speed(region_tissue == 1) = speed_map_tissue(region_tissue == 1);
medium.density(region_tissue == 1) = density_map_tissue(region_tissue == 1);
% show the sound speed map of compound medium, the air and tissue
figure
imagesc(medium.sound_speed(:,:,44))
%%
========================================================================
=
% DEFINE THE INPUT SIGNAL
%
========================================================================
=
% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 12e6; % [Hz]
tone_burst_cycles = 2;
% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles,...
'Envelope', 'Gaussian', 'Plot',true);
% 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_tissue*rho0_tissue)).*input_signal;
%%
========================================================================
=
% DEFINE THE TX TRANSDUCER
%
========================================================================
=
% the linear array transducer for TX
transducer_TX.number_elements = 32; % total number of transducer elements
transducer_TX.element_width = 2 ; % width of each element [grid points]
transducer_TX.element_length = 20 ; % length of each element [grid points]
transducer_TX.element_spacing = 0; % spacing (kerf width) between the elements [grid points]
transducer_TX.radius = inf; % radius of curvature of the transducer [m]
transducer_width_TX = transducer_TX.number_elements*transducer_TX.element_width ...
+ (transducer_TX.number_elements - 1)*transducer_TX.element_spacing;
% transducer_TX.position: X = 1
transducer_TX.position = round([1, Ny/2 - transducer_width_TX/2, Nz/2 -
transducer_TX.element_length/2]);
% properties used to derive the beamforming delays
transducer_TX.sound_speed = 1540; % sound speed [m/s]
transducer_TX.focus_distance = inf; % focus distance [m]
transducer_TX.elevation_focus_distance = 20e-3;% focus distance in the elevation plane [m]
transducer_TX.steering_angle = 0; % steering angle [degrees]
% apodization
transducer_TX.transmit_apodization = 'Rectangular';
transducer_TX.receive_apodization = 'Rectangular';
% use all the element to TX
transducer_TX.active_elements = ones(transducer_TX.number_elements, 1);
% append input signal used to drive the transducer
transducer_TX.input_signal = input_signal;
% create the transducer using the defined settings
transducer_TX_class = makeTransducer(kgrid, transducer_TX);
transducer_TX_class.properties
transducer_TX_class.plot
%%
========================================================================
=
% RUN THE SIMULATION
%
========================================================================
=
% set the input settings
input_args = {'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size],...
'DataCast', DATA_CAST};
[sensor_data] = kspaceFirstOrder3D(kgrid, medium, transducer_TX_class, transducer_TX_class,
input_args{:});
element_data = transducer_TX_class.combine_sensor_data(sensor_data);
save element_data.mat element_data