k-Wave User Forum » Topic: Simulating US Beam by shifting the transducer

k-Wave User Forum » Topic: Simulating US Beam by shifting the transducer http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-transducer Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 02:32:02 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php maktjik on "Simulating US Beam by shifting the transducer" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-transducer#post-1747 Tue, 26 Mar 2013 14:42:06 +0000 maktjik 1747@http://www.k-wave.org/forum/ Here I post again the code with better arrangement: Initial transducer position: <pre><code>transducer.position = round([1, (((Ny-128)/2) + (128-transducer_width)), Nz/2 - transducer.element_length/2]); Complete 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</code></pre> maktjik on "Simulating US Beam by shifting the transducer" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-transducer#post-1740 Tue, 26 Mar 2013 10:47:43 +0000 maktjik 1740@http://www.k-wave.org/forum/ 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