Hi all,
I would like to simulate US beam by shifting the transducer. I have 32 elements transducer (with 32 active elements) and I shift it 1 element_width in each iteration. My looping iterates from 0:(97-1), so with this iteration I intended to simulate 'pseudo' transducer with 128 elements. Is it valid?
I generate the final output in 3D matrix with dimensions 128x1078x97 = (desired_num_channels, total_time_index, number_scan_lines). So I need to do zero padding in every iteration.
Initial position of my transducer is in the edge of my 'pseudo' 128 elements transducer, if that 'pseudo' transducer placed in the middle of computational grid.
transducer.position = round([1, (((Ny-128)/2) + (128-transducer_width)), Nz/2 - transducer.element_length/2]);
My medium is 4 layers phantom (with different soundspeed and density in each layer) above water.
This is the copy of the script:
% =========================================================================
% ULTRASOUND SIMULATION (TRANSDUCER SHIFTING)
% Linear array transducer (32 element)
% Transmit: Soundbeam
% Medium : Layer phantom (4 layers) in water
% =========================================================================
clear all;
%simulation setting
data_cast = 'single';
% run the simulation
run_simulation = true;
% =========================================================================
% DEFINE THE K-WAVE GRID
% =========================================================================
% set the size of the perfectly matched layer (PML)
PML_X_SIZE = 20; % [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 = 256 - 2*PML_X_SIZE; % [grid points]
Ny = 256 - 2*PML_Y_SIZE; % [grid points]
Nz = 128 - 2*PML_Z_SIZE; % [grid points]
% set desired grid size in the x-direction not including the PML
x = 40e-3; % [m]
y = 50e-3; % [m]
z = 60e-3; % [m]
% calculate the spacing between the grid points
dx = x/Nx; % [m]
dy = y/Ny; % [m]
dz = z/Nz; % [m]
% create the k-space grid
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================
c0 = 1497; % c in water [m/s]
rho0 = 1000; % density of water [kg/m^3]
% create the time array
t_end = 40e-6; % [s]
kgrid.t_array = makeTime(kgrid, c0, [], t_end);
% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 5e6; % [Hz]
tone_burst_cycles = 10;
% 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./(c0*rho0)).*input_signal;
% =========================================================================
% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================
% physical properties of the transducer
transducer.number_elements = 32; % total number of transducer elements
transducer.element_width = 1; % width of each element [grid points]
transducer.element_length = 8; % length of each element [grid points]
transducer.element_spacing = 0; % 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 edge of ('pseudo' transducer 128 elements if it is located in the middle of computational grid)
transducer.position = round([1, (((Ny-128)/2) + (128-transducer_width)), Nz/2 - transducer.element_length/2]);
% properties used to derive the beamforming delays
transducer.sound_speed = 1530; % sound speed [m/s]
transducer.focus_distance = 10.5e-3; % focus distance [m]
transducer.elevation_focus_distance = 19e-3; % focus distance in the elevation plane [m]
transducer.steering_angle = 0; % steering angle [degrees]
% apodization
transducer.transmit_apodization = 'Hanning';
transducer.receive_apodization = 'Rectangular';
% define the transducer elements that are currently active
number_active_elements = 32;
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 = makeTransducer(kgrid, transducer);
% print out transducer properties
transducer.properties;
% break
% =========================================================================
% DEFINE THE MEDIUM PROPERTIES
% =========================================================================
% define a random distribution of scatterers for the medium
background_map_mean = 1;
background_map_std = 0.008;
% background_map = background_map_mean + background_map_std*randn([Nx, Ny, Nz]);
% define the thicknesses of layer phantom
d_layer_1 = 5.4e-3; % [m]
d_layer_2 = 4.9e-3; % [m]
d_layer_3 = 4.4e-3; % [m]
d_layer_4 = 6.3e-3; % [m]
% define the properties of the propagation medium
% water
medium.sound_speed = c0*ones(Nx, Ny, Nz); % [m/s]
medium.density = rho0*ones(Nx, Ny, Nz); % [kg/m^3]
% layer_1
background_map_layer_1 = background_map_mean + background_map_std*randn([round(Nx*(d_layer_1/x)), Ny, Nz]);
medium.sound_speed(1:round(Nx*(d_layer_1/x)), :, :) = 1520.*background_map_layer_1; % [m/s]
medium.density(1:round(Nx*(d_layer_1/x)), :, :) = 1281.*background_map_layer_1; % [kg/m^3]
% layer_2
background_map_layer_2 = background_map_mean + background_map_std*randn([round(Nx*((d_layer_1+d_layer_2)/x))-(round(Nx*(d_layer_1/x))), Ny, Nz]);
medium.sound_speed(round(Nx*(d_layer_1/x))+1:round(Nx*((d_layer_1+d_layer_2)/x)), :, :) = 1218.*background_map_layer_2; % [m/s]
medium.density(round(Nx*(d_layer_1/x))+1:round(Nx*((d_layer_1+d_layer_2)/x)), :, :) = 1510.*background_map_layer_2; % [kg/m^3]
% layer_3
background_map_layer_3 = background_map_mean + background_map_std*randn([round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))-(round(Nx*((d_layer_1+d_layer_2)/x))), Ny, Nz]);
medium.sound_speed(round(Nx*((d_layer_1+d_layer_2)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x)), :, :) = 1325.*background_map_layer_3; % [m/s]
medium.density(round(Nx*((d_layer_1+d_layer_2)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x)), :, :) = 1643.*background_map_layer_3; % [kg/m^3]
% layer_4
background_map_layer_4 = background_map_mean + background_map_std*randn([round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x))-(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))), Ny, Nz]);
medium.sound_speed(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x)), :, :) = 1429.*background_map_layer_4; % [m/s]
medium.density(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x)), :, :) = 1684.*background_map_layer_4; % [kg/m^3]
% =========================================================================
% RUN THE SIMULATION
% =========================================================================
% set the input settings
% plot medium map
input_args = {...
'PlotLayout', true, 'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', data_cast};
% run the simulation if set to true, otherwise, load previous results from disk
if run_simulation
desired_num_channels = 128;
total_time_index = 1078;
total_scan_lines = 97;
sensor_data_output = zeros(desired_num_channels, total_time_index, total_scan_lines);
for scan_line = 0 : total_scan_lines-1
fprintf('Calculating scan line %d of %d\n', scan_line+1, total_scan_lines);
transducer_position = transducer.position(:,:,:);
sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
% Zero padding to fill 3D matrix output
padding_before = zeros(scan_line, total_time_index);
padding_after = zeros(desired_num_channels - number_active_elements - scan_line, total_time_index);
% Populate 3D matrix output with the output of each iteration
sensor_data_output(:, :, scan_line + 1) = [padding_before; sensor_data; padding_after];
transducer_position = transducer_position + transducer.element_width;
end
% save the scan lines to disk
save sensor_data_output sensor_data_output;
else
% load the scan lines from disk
load sensor_data_output;
end
% =========================================================================
Thanks,
Makcik