k-Wave User Forum » Topic: Simulating US Beam by 'Shifting' the medium

k-Wave User Forum » Topic: Simulating US Beam by 'Shifting' the medium http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 01:49:55 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2575 Mon, 22 Apr 2013 17:38:14 +0000 Bradley Treeby 2575@http://www.k-wave.org/forum/ Hi Makcik, Your approach sounds right for simulating plane wave imaging using a conventional diagnostic ultrasound transducer. In this case, you could use the same transducer object for both source and sensor as the active elements will be 128 in both cases. Brad. maktjik on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2343 Wed, 17 Apr 2013 15:17:31 +0000 maktjik 2343@http://www.k-wave.org/forum/ Hi Brad, Thanks for your suggestion, now I can simply save the final output in 3D matrix. I also reduce tone_burst center frequency, because I just check the spectrum exceed the maximum supported frequency. And what about the plane wave case I asked before? - in transmit: each element (total: 128 elements), propagate the signal in the same time, focus_distance is infinite, and elevation_focus_distance is finite. - in receive: total 128 elements also receive the echo in the same time. Can I simulate it by simply define 128 total elements transducer (with also 128 active elements), define focus_distance to be inf, and set them as KSpaceFirstOrder3D input arguments in the places of source and sensor? Another parameters, such as computational grid, medium parameters, medium properties, input signal and time array are still the same as before. Thanks, Makcik Bradley Treeby on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2293 Tue, 16 Apr 2013 16:05:46 +0000 Bradley Treeby 2293@http://www.k-wave.org/forum/ Hi Makcik, For each individual simulation you run, you will obtain 128 time signals. These are equivalent to the raw RF data from an ultrasound machine. If you run 97 simulations changing the transmit elements each time, in the end you will finally obtain 97*128 time series. You could store these using something like: <code>RF_data(sim_index, :, :) = kspaceFirstOrder3D(...)</code>. Glancing at your code, it looks like you set 128 elements on both transmit and receive. Also, see my previous content about the frequency of the source and the frequency supported by the grid. Brad. maktjik on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2285 Tue, 16 Apr 2013 11:24:44 +0000 maktjik 2285@http://www.k-wave.org/forum/ PS: I just get the result from my simulation (transducer shifting, 32 active elements transmit and 128 active elements receive), but I just try it for 1 scanline and 2 scanlines. Why the first 2 columns and the last 2 columns of my data are always filled by zero? Here is my latest code: <pre><code>% ========================================================================= % ULTRASOUND SIMULATION % Linear array transducer (128 element) % Transmit: Sound beam (Focusing), 32 active elements % Receive : 128 active elements % Medium : Layer phantom (4 layers) in water % ========================================================================= clear all; %simulation setting data_cast = 'single'; % run the simulation % true : run the simulation % false : load previous results from disk run_simulation = true; % ========================================================================= % DEFINE THE K-WAVE GRID % ========================================================================= % set the size of the perfectly matched layer (PML) PML_X_SIZE = 10; % [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 = 512 - 2*PML_Y_SIZE; % [grid points] Nz = 256 - 2*PML_Z_SIZE; % [grid points] % set desired grid size in the x-direction not including the PML x = 38e-3; % [m] y = 50e-3; % [m] z = 38e-3; % [m] % calculate the spacing between the grid points dx = z/Nz; % [m] dy = dx; % [m] dz = dx; % [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] medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)] medium.alpha_power = 1.5; medium.BonA = 6; % represent nonlinearity properties of the medium % 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 = 3.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 [SOURCE] % ========================================================================= % physical properties of the transducer transducer.number_elements = 128; % total number of transducer elements transducer.element_width = 2; % width of each element [grid points] transducer.element_length = 25; % 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 middle of the computational grid transducer.position = round([1, Ny/2 - transducer_width/2, 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 = 16e-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 = 128; % 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; % ========================================================================= % DEFINE THE ULTRASOUND TRANSDUCER [SENSOR] % ========================================================================= % physical properties of the transducer transducer2.number_elements = 128; % total number of transducer elements transducer2.element_width = 2; % width of each element [grid points] transducer2.element_length = 25; % length of each element [grid points] transducer2.element_spacing = 0; % spacing (kerf width) between the elements [grid points] transducer2.radius = inf; % radius of curvature of the transducer [m] % calculate the width of the transducer in grid points transducer_width2 = transducer2.number_elements*transducer2.element_width ... + (transducer2.number_elements - 1)*transducer2.element_spacing; % use this to position the transducer in the middle of the computational grid transducer2.position = round([1, Ny/2 - transducer_width2/2, Nz/2 - transducer2.element_length/2]); % properties used to derive the beamforming delays transducer2.sound_speed = 1530; % sound speed [m/s] transducer2.focus_distance = 10.5e-3; % focus distance [m] transducer2.elevation_focus_distance = 16e-3; % focus distance in the elevation plane [m] transducer2.steering_angle = 0; % steering angle [degrees] % apodization transducer2.transmit_apodization = 'Hanning'; transducer2.receive_apodization = 'Rectangular'; % define the transducer elements that are currently active number_active_elements2 = 128; transducer2.active_elements = ones(transducer2.number_elements, 1); % create the transducer using the defined settings transducer2 = makeTransducer(kgrid, transducer2); % print out transducer properties transducer2.properties; % ========================================================================= % DEFINE THE MEDIUM PROPERTIES % ========================================================================= % define a random distribution of scatterers for the medium background_map_mean = 1; background_map_std = 0.008; % define the thicknesses of layer phantom in x direction 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 thicknesses of layer phantom in x direction position_x1 = 1:round(Nx*(d_layer_1/x)); position_x2 = round(Nx*(d_layer_1/x))+1:round(Nx*((d_layer_1+d_layer_2)/x)); position_x3 = round(Nx*((d_layer_1+d_layer_2)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x)); position_x4 = 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)); % define the width of phantom in y and z direction width_y = 29e-3; % width_phantom_in y direction [grid points]*dy = [m] width_z = 35e-3; % width_phantom_in z direction [grid points]*dz = [m] % define position of each layer in y and z direction position_y = round(Ny*(((y-width_y)/2)/y))+1:round(Ny*(y-(y-width_y)/2)/y); position_z = round(Nz*(((z-width_z)/2)/z))+1:round(Nz*(z-(z-width_z)/2)/z); % define the dimension of each layer (equal to its background map) width_backgroundmap_layer_1 = [round(Nx*(d_layer_1/x)), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))]; width_backgroundmap_layer_2 = [round(Nx*((d_layer_1+d_layer_2)/x))-(round(Nx*(d_layer_1/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))]; width_backgroundmap_layer_3 = [round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))-(round(Nx*((d_layer_1+d_layer_2)/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))]; width_backgroundmap_layer_4 = [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))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))]; % % 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(width_backgroundmap_layer_1); medium.sound_speed(position_x1, position_y, position_z) = 1520.*background_map_layer_1; % [m/s] medium.density(position_x1, position_y, position_z) = 1281.*background_map_layer_1; % [kg/m^3] % layer_2 background_map_layer_2 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_2); medium.sound_speed(position_x2, position_y, position_z) = 1218.*background_map_layer_2; % [m/s] medium.density(position_x2, position_y, position_z) = 1510.*background_map_layer_2; % [kg/m^3] % layer_3 background_map_layer_3 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_3); medium.sound_speed(position_x3, position_y, position_z) = 1325.*background_map_layer_3; % [m/s] medium.density(position_x3, position_y, position_z) = 1643.*background_map_layer_3; % [kg/m^3] % layer_4 background_map_layer_4 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_4); medium.sound_speed(position_x4, position_y, position_z) = 1429.*background_map_layer_4; % [m/s] medium.density(position_x4, position_y, position_z) = 1684.*background_map_layer_4; % [kg/m^3] % ========================================================================= % RUN THE SIMULATION % ========================================================================= % set the input settings % plot medium map input_args = {... '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 total_scan_lines = 97; total_scan_lines_test = 2; for scan_line = 1 : total_scan_lines_test fprintf('Calculating scan line %d of %d\n', scan_line, total_scan_lines_test); % define the transducer elements that are currently active transducer.active_elements = zeros(transducer.number_elements, 1); transducer.active_elements(scan_line:scan_line+31)=1; RF_data_SB_transd_shift = kspaceFirstOrder3D(kgrid, medium, transducer, transducer2, input_args{:}); end end</code></pre> Thanks Makcik maktjik on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2278 Mon, 15 Apr 2013 10:54:45 +0000 maktjik 2278@http://www.k-wave.org/forum/ Hi Brad, Thanks for your suggestion, and currently my simulation is still running. For question number 2: I still don't really understand how it works. I read in the manual: "The way in which the signals across each element are calculated depends on the setting for transducer.elevation_focus_distance. If this is set to Inf, the signals across the grid points within each sensor element are averaged at each time step,and only the average is stored. If a finite elevation focus distance is defined, a buffer the length of the longest beamforming delay is filled using a FIFO queue (first-in, first-out). The elevation beamforming is then computed on the fly once the buffer is filled." In my case: I define the elevation focus distance to be finite and each element has 2 grid points. In the first iteration I transmit the signal from the first 32 active elements and receive in 128 active elements, and so on until the last 32 active elements transmit the 97th beam. So from the first iteration, I have first 128 data filled in each transducer element. My question: for the final <code>sensor_data</code> which I name <code>RF_data_SB_transd_shift</code>, is it the average or summing of the 97 set data? For question number 3: My MATLAB was simply crash and didn't work at all every time I call transducer.plot, I think it's the issue of memory because I have relatively large 3D computational grid. Another question: I have another plan to simulate another type of US propagation (plane wave), which: - in transmit: each element (total: 128 elements), propagate the signal in the same time, focus_distance is infinite, and elevation_focus_distance is finite. - in receive: total 128 elements also receive the echo in the same time. Can I simulate it by simply define 128 total elements transducer (with also 128 active elements), define focus_distance to be inf, and set them as KSpaceFirstOrder3D input arguments in the places of source and sensor? Another parameters, such as computational grid, medium parameters, medium properties, input signal and time array are still the same as before. Thank you very much, Makcik Bradley Treeby on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2245 Fri, 12 Apr 2013 10:58:46 +0000 Bradley Treeby 2245@http://www.k-wave.org/forum/ Hi Makcik, To answer your questions: 1. If you use an object of the <code>kWaveTransducer</code> class (created using <code>makeTransducer</code>) to replace the <code>sensor</code> input, the number of time signals returned will be equal to the number of active elements. With the current code, if you wish to transmit using 32 elements, and receive using a different number of elements, you will need to create a second copy of your transducer, set the desired number of active elements, and use this in place of the <code>sensor</code> input. 2. If your physical transducer elements are made up from multiple grid points, the signal received by the physical element is calculated at each time point by averaging over the grid points. There is more information about how this is done within the k-Wave manual. 3. What is the error your receive, is it a memory error? Unfortunately the function <code>voxelPlot</code> does not deal very well with very large 3D matrices. 4. The number displayed should be correct. Perhaps double check your grid spacing (<code>kgrid.dx</code>, <code>kgrid.dt</code>, <code>kgrid.dz</code>) and your sound speeds (<code>min(medium.sound_speed(:)</code>). Brad. maktjik on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2243 Fri, 12 Apr 2013 09:41:00 +0000 maktjik 2243@http://www.k-wave.org/forum/ Hi Brad, Regarding my 1st question above, I also have tried to run the example from the toolbox (example_us_bmode_linear_transducer.m). * I changed the total number of elements to 128, and number of active elements was still the same (32). * I run the loop, only until scanline number 2. * And I get the result <code>sensor_data</code> with dimension 128 x total time step. All of 128 element is not empty, eventhough the scanline was only 2. * But why in my case (in my latest code), I get the data with dimension only 32 x total time index, eventhough I have 128 total number of elements and any number of scanlines? Thank you very much, Makcik maktjik on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2213 Thu, 11 Apr 2013 10:23:36 +0000 maktjik 2213@http://www.k-wave.org/forum/ Hi Brad, I implemented your suggestion in my loop: <pre><code>% ========================================================================= % RUN THE SIMULATION % ========================================================================= % set the input settings % plot medium map input_args = {... '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 num_channels = transducer.number_elements; total_time_index = 1078; total_scan_lines = 97; RF_data_SB_transd_shift = zeros(num_channels, total_time_index); for scan_line = 1 : total_scan_lines fprintf('Calculating scan line %d of %d\n', scan_line, total_scan_lines); % define the transducer elements that are currently active transducer.active_elements = zeros(transducer.number_elements, 1); transducer.active_elements(scan_line:scan_line+31)=1; RF_data_SB_transd_shift = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:}); end save RF_data_SB_transd_shift RF_data_SB_transd_shift;</code></pre> My questions: 1. In this case, I have only 2D data, which previously I expected to have the information about total number of elements and time index, but, I get the result which only contains 32 number of elements instead of 128. Why? And also where I can find the information about scanline in my data (in this case 97 beams or scanlines)? Actually, in this simulation, I try to mimic my laboratory experiment, which produce 3D DAQ raw RF data, which contains the information about total number of elements, total number of sample, and number of scanlines (beams). 2. In KSpaceFirstOrder3D function, Is it in each iteration, the data which are filled in each transducer element being overwrited or being summed? 3. Why does it crash every time I call transducer.plot? I run this toolbox in MATLAB 2012a. 4. Why is the maximum supported frequency is not match with the formula: <code>f_max = min(c_0)/(2*dx)</code>. In my case: It's supposed to be <code>1218/(2*254.2373um) = 2.3954 MHz</code>, but on the command line when the simulation runs it was 2.3086MHz. Thank you very much. Makcik Bradley Treeby on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2039 Fri, 05 Apr 2013 17:27:00 +0000 Bradley Treeby 2039@http://www.k-wave.org/forum/ Hi Makcik, I haven't run your code, but it looks like you haven't updated the active transducer elements with the for loop. You could use something like: <pre><code>% define the transducer elements that are currently active transducer.active_elements = zeros(transducer.number_elements, 1); transducer.active_elements(scan_line:scan_line+31)=1;</code></pre> Here <code>transducer</code> is an object of the <code>kWaveTransducer</code> class which contains the properties of the transducer. There are more details in the k-Wave manual if you get stuck. Brad. maktjik on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2035 Fri, 05 Apr 2013 10:42:33 +0000 maktjik 2035@http://www.k-wave.org/forum/ Hi Brad, I try to implement your idea to simulate the whole medium by update which elements are active within each iteration. My total number of elements is 128. I also define new computational grid and its properties, and the position of my medium. Could you take a look at this? whether my assumptions is already valid or not. Here is the copy of my new script. Thank you very much Makcik <pre><code>% ========================================================================= % ULTRASOUND SIMULATION % Linear array transducer (128 element) % Transmit: Sound beam (Focusing) % Medium : Layer phantom (4 layers) in water % ========================================================================= clear all; %simulation setting data_cast = 'single'; % run the simulation % true : run the simulation % false : load previous results from disk run_simulation = true; % ========================================================================= % DEFINE THE K-WAVE GRID % ========================================================================= % set the size of the perfectly matched layer (PML) PML_X_SIZE = 10; % [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 = 256 - 2*PML_Z_SIZE; % [grid points] % set desired grid size in the x-direction not including the PML x = 60e-3; % [m] y = 60e-3; % [m] z = 60e-3; % [m] % calculate the spacing between the grid points dx = x/Nx; % [m] dy = dx; % [m] dz = dx; % [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] medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)] medium.alpha_power = 1.5; medium.BonA = 6; % represent nonlinearity properties of the medium % 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 = 2e6; % [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 = 128; % 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 middle of the computational grid transducer.position = round([1, Ny/2 - transducer_width/2, 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 transducer.active_elements = zeros(transducer.number_elements, 1); transducer.active_elements(1:32)=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; % ========================================================================= % DEFINE THE MEDIUM PROPERTIES % ========================================================================= % define a random distribution of scatterers for the medium background_map_mean = 1; background_map_std = 0.008; % 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 position of each layer in x direction (4 layers: 1 until 4) position_x1 = 1:round(Nx*(d_layer_1/x)); position_x2 = round(Nx*(d_layer_1/x))+1:round(Nx*((d_layer_1+d_layer_2)/x)); position_x3 = round(Nx*((d_layer_1+d_layer_2)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x)); position_x4 = 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)); % define the width of phantom in y and z direction (in this case it's equal to transducer width) width_y = 32.5424e-3; % [m] width_z = 32.5424e-3; % [m] % define position of each layer in y and z direction position_y = round(Ny*(((y-width_y)/2)/y))+1:round(Ny*(y-(y-width_y)/2)/y); position_z = round(Nz*(((z-width_z)/2)/z))+1:round(Nz*(z-(z-width_z)/2)/z); % define the dimension of each layer (and its background map) width_backgroundmap_layer_1 = [round(Nx*(d_layer_1/x)), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))]; width_backgroundmap_layer_2 = [round(Nx*((d_layer_1+d_layer_2)/x))-(round(Nx*(d_layer_1/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))]; width_backgroundmap_layer_3 = [round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))-(round(Nx*((d_layer_1+d_layer_2)/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))]; width_backgroundmap_layer_4 = [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))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))]; % % 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(width_backgroundmap_layer_1); medium.sound_speed(position_x1, position_y, position_z) = 1520.*background_map_layer_1; % [m/s] medium.density(position_x1, position_y, position_z) = 1281.*background_map_layer_1; % [kg/m^3] % layer_2 background_map_layer_2 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_2); medium.sound_speed(position_x2, position_y, position_z) = 1218.*background_map_layer_2; % [m/s] medium.density(position_x2, position_y, position_z) = 1510.*background_map_layer_2; % [kg/m^3] % layer_3 background_map_layer_3 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_3); medium.sound_speed(position_x3, position_y, position_z) = 1325.*background_map_layer_3; % [m/s] medium.density(position_x3, position_y, position_z) = 1643.*background_map_layer_3; % [kg/m^3] % layer_4 background_map_layer_4 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_4); medium.sound_speed(position_x4, position_y, position_z) = 1429.*background_map_layer_4; % [m/s] medium.density(position_x4, position_y, position_z) = 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 num_channels = transducer.number_elements; total_time_index = 1078; total_scan_lines = 97; RF_data_SB_transd_shift = zeros(num_channels, total_time_index, total_scan_lines); for scan_line = 1 : total_scan_lines fprintf('Calculating scan line %d of %d\n', scan_line, total_scan_lines); transducer_active_elements_position = transducer.active_elements(1:32); RF_data_SB_transd_shift = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:}); % update the position of transducer active elements transducer_active_elements_position = transducer_active_elements_position + transducer.element_width; end save RF_data_SB_transd_shift RF_data_SB_transd_shift; else % load the scan lines from disk load RF_data_SB_transd_shift; end</code></pre> Bradley Treeby on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2029 Thu, 04 Apr 2013 09:06:57 +0000 Bradley Treeby 2029@http://www.k-wave.org/forum/ Hi Makcik, The maximum supported temporal frequency in each direction can be calculated based on the grid spacing and the sound speed via the relation <code>f_max = min(c_0)/(2*dx)</code>. These frequencies are also reported on the command line when the simulation runs. There is some background information that you may find useful in the <a href="http://www.k-wave.org/documentation.php">k-Wave Manual</a>. The entire frequency content of your source (including nonlinear harmonics) should be less than the maximum frequency to avoid aliasing errors. In this particular example, the frequency is set by <code>tone_burst_freq</code> as you say, but keep in mind that a tone-burst will contain a range of frequencies (you can see this by looking at the frequency spectrum of the input signal). If you are moving the medium, the number of mimicked transducer elements will be N_elements + N_simulations - 1, so in the example you mention, this will be 127. If you want to simulate the entire medium, the alternative way would be to define the transducer with 128 elements, then within each iteration, update which elements are active. Hope that helps, Brad. maktjik on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-2017 Wed, 03 Apr 2013 08:54:43 +0000 maktjik 2017@http://www.k-wave.org/forum/ Hi Brad, Sorry for my late response. I am quite new to this area and this toolbox, so I didn't realize that the spacing grid affects the supported frequencies. How does it work (how to calculate those 3.2MHz by 2.8MHz by 1.1MHz frequencies)? And for the source frequency, is it : <code>tone_burst_freq = 5e6; % [Hz]</code> which you mean? Is it suppose to be less or more than supported frequencies? Sorry, if I repeat these questions: In the example I reference: 96 separate simulations are performed, mimicking the collection of scan lines from a transducer with a large number of elements. What is exactly the number of elements of the transducer that you mimic? Is it 128? In my case: if I want to simulate 'pseudo' transducer with 128 elements, is it valid to make iteration from 0:(97-1)? And do this step to rearrange the final output: <pre><code>% Zero padding to fill 3D matrix output padding_before = zeros(scan_line_index, total_time_index); padding_after = zeros(desired_num_channels - number_active_elements - scan_line_index, total_time_index); % Populate 3D matrix output with the output of each iteration sensor_data_output(:, :, scan_line_index + 1) = [padding_before; sensor_data; padding_after];</code></pre> so my final otput will be <code>3D matrix with dimensions 128x1078x97 = (desired_num_channels x total_time_index x number_scan_lines).</code> Actually my goal is to simulate ultrasound beam with 32 elements transmit aperture and it shifts (overlap 1 element in each iteration), so If I shift it 96 times I get 97 beam until the last element of my transducer (real or 'pseudo' 128 elements transducer). Therefore, in another post: <a href="http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-transducer" rel="nofollow">http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-transducer</a> I try to do it in different way, Is it also possible? Or do you have suggestion how to do it properly? I will do beamforming later on outside this toolbox, so in this case I just use this toolbox to generate (simulate) raw RF ultrasound data. Thank you very much. Makcik Bradley Treeby on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-1815 Wed, 27 Mar 2013 16:01:53 +0000 Bradley Treeby 1815@http://www.k-wave.org/forum/ Hi Makcik, For each simulation, it is only the total number of active elements that matters, as this defines where the acoustic waves will be added to the medium. Similarly for the beamforming step - it is only the number of active transducer elements that matters. In the example you reference, 96 separate simulations are performed, mimicking the collection of scan lines from a transducer with a large number of elements. These are only combined in the final scan conversion step. I noticed in your simulation you have a different grid spacing in X, Y, and Z, giving rise to different maximum supported frequencies (3.2MHz by 2.8MHz by 1.1MHz). Is there any reason for this? Also, your source frequency always be less than the maximum supported frequency. Hope that helps, Brad. maktjik on "Simulating US Beam by 'Shifting' the medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-1748 Tue, 26 Mar 2013 14:48:18 +0000 maktjik 1748@http://www.k-wave.org/forum/ Here I post again the code with better arrangement: <pre><code>% ========================================================================= % ULTRASOUND SIMULATION (MEDIUM 'SHIFTING') % Linear array transducer (32 element) % Transmit: Sound beam % 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 middle of the computational grid transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]); % properties used to derive the beamforming delays transducer.sound_speed = c0; % sound speed [m/s] transducer.focus_distance = 11e-3; % focus distance [m] transducer.elevation_focus_distance = 11e-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; % ========================================================================= % DEFINE THE MEDIUM PROPERTIES % ========================================================================= % defining a large image size to move across number_scan_lines = 97; Nx_tot = Nx; Ny_tot = Ny + number_scan_lines*transducer.element_width; Nz_tot = Nz; % define a random distribution of scatterers for the medium background_map_mean = 1; background_map_std = 0.008; % 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 global_medium.sound_speed = c0*ones(Nx, Ny_tot, Nz); % [m/s] global_medium.density = rho0*ones(Nx, Ny_tot, Nz); % [kg/m^3] % layer_1 background_map_layer_1 = background_map_mean + background_map_std*randn([round(Nx*(d_layer_1/x)), Ny_tot, Nz]); global_medium.sound_speed(1:round(Nx*(d_layer_1/x)), :, :) = 1520.*background_map_layer_1; % [m/s] global_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_tot, Nz]); global_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] global_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_tot, Nz]); global_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] global_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_tot, Nz]); global_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] global_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 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 if run_simulation desired_num_channels = 128; total_time_index = 1078; % set medium position medium_position = 1; sensor_data_output = zeros(desired_num_channels, total_time_index, number_scan_lines); for scan_line_index = 0:number_scan_lines-1 % update the command line status disp(''); disp(['Computing scan line ' num2str(scan_line_index+1) ' of ' num2str(number_scan_lines)]); % load the current section of the medium local_medium.sound_speed = global_medium.sound_speed(:, medium_position:medium_position + Ny - 1, :); local_medium.density = global_medium.density(:, medium_position:medium_position + Ny - 1, :); % run the simulation sensor_data = kspaceFirstOrder3D(kgrid, local_medium, transducer, transducer, input_args{:}); % Zero padding to fill 3D matrix output padding_before = zeros(scan_line_index, total_time_index); padding_after = zeros(desired_num_channels - number_active_elements - scan_line_index, total_time_index); % Populate 3D matrix output with the output of each iteration sensor_data_output(:, :, scan_line_index + 1) = [padding_before; sensor_data; padding_after]; % update medium position medium_position = medium_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 medium" http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-medium#post-1739 Tue, 26 Mar 2013 10:23:12 +0000 maktjik 1739@http://www.k-wave.org/forum/ Hi all, I would like to simulate US beam by 'shifting' the medium. I have 32 elements transducer (with 32 active elements). For each new scanline I updated the value of medium position by define local_medium.soundspeed and local_medium.density. My medium is 4 layers phantom (with different soundspeed and density in each layer) above water. I try to copy the idea from the 'Simulating B-mode Ultrasound Images Example'. In this example, scan_line_index iterates from 1:96. My question: it aims to model 'pseudo' transducer with 128 elements? Because it shifts 1 element_width in each iteration. Does the number of scanlines affect the output 'sensor_data' from the function 'kspaceFirstOrder3D'? or It only affects this beamforming step?: % extract the scan line from the sensor data scan_lines(scan_line_index, :) = transducer.scan_line(sensor_data); In my case, if I want to simulate 'pseudo' transducer with 128 elements, is it valid to make iteration from 0:(97-1)? Because it makes me easier to adjust the zero padding which I did later on. I just need the output from the function 'kspaceFirstOrder3D' (not the output after beamformed by using 'transducer.scan_line(sensor_data)'). 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. This is the copy of the script: % ========================================================================= % ULTRASOUND SIMULATION (MEDIUM 'SHIFTING') % Linear array transducer (32 element) % Transmit: Sound beam % 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 middle of the computational grid transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]); % properties used to derive the beamforming delays transducer.sound_speed = c0; % sound speed [m/s] transducer.focus_distance = 11e-3; % focus distance [m] transducer.elevation_focus_distance = 11e-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; % ========================================================================= % DEFINE THE MEDIUM PROPERTIES % ========================================================================= % defining a large image size to move across number_scan_lines = 97; Nx_tot = Nx; Ny_tot = Ny + number_scan_lines*transducer.element_width; Nz_tot = Nz; % define a random distribution of scatterers for the medium background_map_mean = 1; background_map_std = 0.008; % 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 global_medium.sound_speed = c0*ones(Nx, Ny_tot, Nz); % [m/s] global_medium.density = rho0*ones(Nx, Ny_tot, Nz); % [kg/m^3] % layer_1 background_map_layer_1 = background_map_mean + background_map_std*randn([round(Nx*(d_layer_1/x)), Ny_tot, Nz]); global_medium.sound_speed(1:round(Nx*(d_layer_1/x)), :, :) = 1520.*background_map_layer_1; % [m/s] global_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_tot, Nz]); global_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] global_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_tot, Nz]); global_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] global_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_tot, Nz]); global_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] global_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 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 if run_simulation desired_num_channels = 128; total_time_index = 1078; % set medium position medium_position = 1; sensor_data_output = zeros(desired_num_channels, total_time_index, number_scan_lines); for scan_line_index = 0:number_scan_lines-1 % update the command line status disp(''); disp(['Computing scan line ' num2str(scan_line_index+1) ' of ' num2str(number_scan_lines)]); % load the current section of the medium local_medium.sound_speed = global_medium.sound_speed(:, medium_position:medium_position + Ny - 1, :); local_medium.density = global_medium.density(:, medium_position:medium_position + Ny - 1, :); % run the simulation sensor_data = kspaceFirstOrder3D(kgrid, local_medium, transducer, transducer, input_args{:}); % Zero padding to fill 3D matrix output padding_before = zeros(scan_line_index, total_time_index); padding_after = zeros(desired_num_channels - number_active_elements - scan_line_index, total_time_index); % Populate 3D matrix output with the output of each iteration sensor_data_output(:, :, scan_line_index + 1) = [padding_before; sensor_data; padding_after]; % update medium position medium_position = medium_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