Dear Mr Treeby,
I am trying to simulate a "simple" set up, and I don't really understand the signals I receive.
My phantom is homogenous in density but for a 2x2x2 pixels highly reflective scatterer in the middle (high density). In terms of speed of sound, it is homogenous, with a bulk at 1540m/s and a layer, placed before the scatterer, at higher speed of sound. The scatterer gets an extra reflective punch from having also a higher speed of sound. There is no attenuation set.
I have a linear transducer, with only 10 elements, focused on my scatterer in azimuth, and not focused in elevation.
My first question is regarding the signal I get when I do not inforce the layer.
I get this:
https://www.dropbox.com/s/cwymc0crk3tn9cc/Signal_no_layer_after_Tukey_before_TGC.PNG
This, as the name says, is without any other post-processing step than the Tukey windowing to remove the trace of the input signal, but cropped to a the depth of my phantom.
I do not understand the first, very big, echo, as my scatterer shows up where it is expected, at 2cm*. Do you know where this comes from?
* about the scatterer position: I actually found out that if I wanted it to be at the expected depth, I had to build up the r vector not with t0 at half the input signal, but at a quarter of it. Any idea why? I had the feeling that the half-signal correction was more like a compromise than a definite rule.
Then, when I put up the layer, I see echos within it that I'm not expecting.
https://www.dropbox.com/s/mb8o57iayypaaty/Signal_layer_after_Tukey_before_TGC.PNG
The fact that my scatterer is drowned in the tail is not so much of a problem if I can explain where this "tail" comes from. But the presence of some signal between the two interfaces of the layer is strange to me.
My hypothesis at the moment is that the grid spacing is too coarse. Do you have a better idea?
Thank you very much,
Margot
Code:
% define heterogenous medium
% phantom size
X = 0.04; %[m]
Y = 0.02; %[m]
Z = 0.02; %[m]
% SOS bulk
c = 1540; % [m/s]
delta_c = 18; % [m/s]
% density bulk
rho = 1000; % [kg/m^3]
% highly reflective scatterer properties
scatterer.c = 1590;
scatterer.rho = 2000; % [kg/m^3]
% grid spacing
dx = c/(3*2.85e6); % [m]
dy = dx; % better for stability
dz = dx; % better for stability
% define the transducer
% physical properties
transducer.number_elements = 10;
transducer.element_width = round(9.905e-4/dx);
transducer.element_length = transducer.element_width;
transducer.element_spacing = 0;
transducer.radius = Inf;
% effective grid size
Nx = round(X/dx);
Nz = round(Z/dz);
n_scan_lines = 1;
Ny_tot = round(Y/dy);
Ny = Ny_tot - (n_scan_lines-1)*transducer.element_width;
% layer depth
L_top = 0.01; %[m]
LTH = 0.005; %[m] layer thickness
% highly reflective scatterer position
scatterer.x = X/2; %[m]
scatterer.y = Y/2; %[m]
scatterer.z = Z/2; %[m]
% PML size for border conditions
PML_X = 30;
PML_Y = 15;
PML_Z = 15;
% make grid
% space
kgrid = makeGrid(Nx,dx,Ny,dy,Nz,dz);
% time
t_end = (Nx*dx)*2.2/c;
[kgrid.t_array,dt] = makeTime(kgrid, c, [], t_end);
% medium physical properties
SOS_map = repmat(c,[Nx,Ny_tot,Nz]);
rho_map = repmat(rho,[Nx,Ny_tot,Nz]);
% layer mask
layer_mask = zeros(size(SOS_map));
layer_mask(round(L_top/dx):round(L_top/dx)+round(LTH/dx),:,:) = 1;
% neutral background map
SOS_map(layer_mask == 1) = SOS_map(layer_mask == 1)+ delta_c;
% highly reflective scatterer
SOS_map(round(scatterer.x/dx):round(scatterer.x/dx)+1,...
round(scatterer.y/dy):round(scatterer.y/dy)+1,...
round(scatterer.z/dz):round(scatterer.z/dz)+1) = scatterer.c;
rho_map(round(scatterer.x/dx):round(scatterer.x/dx)+1,...
round(scatterer.y/dy):round(scatterer.y/dy)+1,...
round(scatterer.z/dz):round(scatterer.z/dz)+1) = scatterer.rho;
% Transducer geometry
% position and size
element_size = [transducer.element_width transducer.element_length]*dx;
transducer_aperture = element_size(1)*transducer.number_elements + (transducer.number_elements-1)*transducer.element_spacing*dx;
transducer_width = transducer.number_elements*transducer.element_width ...
+ (transducer.number_elements - 1)*transducer.element_spacing;
transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]);
% beamforming properties
transducer.sound_speed = c;
transducer.focus_distance = X/2 + dx; % in the middle of the phantom
transducer.steering_angle = 0;
transducer.active_elements = ones(transducer.number_elements,1);
transducer.transmit_apodization = 'Hanning';
transducer.receive_apodization = 'Hanning';
% input signal
% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 2.85e6; % [Hz]
tone_burst_cycles = 4;
% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
% scale the source magnitude by the source_strength divided by the
% impedance (the source is assigned to the particle velocity)
input_signal = (source_strength./(c*rho)).*input_signal;
% assign
transducer.input_signal = input_signal;
% create
transducer = makeTransducer(kgrid,transducer);
% Simulation scan
scan_line = zeros(n_scan_lines,kgrid.Nt);
medium_position = 1;
for l = 1:n_scan_lines
medium.sound_speed = SOS_map(:,medium_position:medium_position + Ny - 1 ,:);
medium.density = rho_map(:,medium_position:medium_position + Ny - 1,:);
input_args = {...
'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X, PML_Y, PML_Z], ...
'DataCast', 'single'};
sensor_data = kspaceFirstOrder3D(kgrid,medium,transducer,transducer,input_args{:});
[rf,ch_data] = transducer.scan_line(sensor_data);
