http://www.k-wave.org/forum/topic/simulation-of-ultrasound-propagating-from-soft-tissue-to-air Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 01:52:59 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/simulation-of-ultrasound-propagating-from-soft-tissue-to-air#post-5450 Tue, 05 Apr 2016 08:55:34 +0000 Bradley Treeby 5450@http://www.k-wave.org/forum/ <p>Hi Jingba,</p> <p>First, it's worth pointing out that the situation you are simulating has an impedance mismatch on the order of 1e6, which makes it challenging to deal with computationally. The NaNs are because the code is unstable with the grid parameters you have used. </p> <p>There are a few problems with your code:</p> <p>(1) The definition of dx should use the minimum sound speed in the medium, which in your case will be air (not 1500 m/s). You can still use the sound speed in tissue to define your value of <code>t_end</code>.</p> <p>(2) The CFL needs to be reduced for the code to be stable. Something like the following should work</p> <pre><code>% create the time array t_end = (Nx*dx)*1.5/c0_tissue; % [s] CFL = 0.05; kgrid.t_array = makeTime(kgrid, c0_tissue, CFL, t_end);</code></pre> <p>(3) You should be aware of the frequency content of your source, and the maximum frequency you have chosen for the grid. If you plot the frequency spectrum of <code>transducer_TX.input_signal</code>, you will see that it has energy above your chosen maximum frequency. You could either increase the maximum frequency, reduce the source frequency, or remove the unsupported frequencies using <code>filterTimeSeries</code>.</p> <p>Hope that helps,</p> <p>Brad </p> http://www.k-wave.org/forum/topic/simulation-of-ultrasound-propagating-from-soft-tissue-to-air#post-5431 Fri, 18 Mar 2016 08:36:15 +0000 Bowen Jing 5431@http://www.k-wave.org/forum/ <p>Hi,</p> <p>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.<br /> 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.</p> <pre><code>clc clear all close all % simulation settings DATA_CAST = &#39;double&#39;; % set to &#39;single&#39; or &#39;gpuArray-single&#39; 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,... &#39;Envelope&#39;, &#39;Gaussian&#39;, &#39;Plot&#39;,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 = &#39;Rectangular&#39;; transducer_TX.receive_apodization = &#39;Rectangular&#39;; % 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 = {&#39;PMLInside&#39;, false, &#39;PMLSize&#39;, [pml_x_size, pml_y_size, pml_z_size],... &#39;DataCast&#39;, 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</code></pre>