http://www.k-wave.org/forum/topic/understanding-speed-of-sound-of-the-medium Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 02:34:06 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/understanding-speed-of-sound-of-the-medium#post-6211 Fri, 10 Nov 2017 15:14:04 +0000 Bradley Treeby 6211@http://www.k-wave.org/forum/ <p>Hi maktijk,</p> <p>1. When you use <code>makeTime</code>, you can directly pass the heterogeneous sound speed. This will use the highest sound speed to calculate the time step, and the lowest sound speed to calculate the number of time steps (although in your example you set this manually).</p> <p>Something to note, by default, k-Wave will use the highest sound speed in the k-space correction term to make the simulation unconditionally stable (at least in the case of no absorption). However, this means dispersion will accumulate in other regions of the domain where the sound speed is lower. In your case, I would suggest setting <code>medium.sound_speed_ref</code> to the background sound speed (water), so the propagation is dispersion free in these regions. In this case, the simulation will no longer be unconditionally stable. You can query the largest stable time step using the <code>checkStability</code> function.</p> <p>2. For the input signal, use the sound speed where the transducer is (water in your case).</p> <p>3. The transducer sound speed is only used to set the focus delays, so I would use the sound speed in water.</p> <p>Hope that helps,</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/understanding-speed-of-sound-of-the-medium#post-6190 Tue, 07 Nov 2017 04:41:34 +0000 maktjik 6190@http://www.k-wave.org/forum/ <p>Dear Dr. Treeby and Dr. Cox,</p> <p>I have some questions regarding my simulation. In my case, I want to simulate plane wave ultrasound (128 elements array) propagation through heterogeneous medium. My medium is steel cylinder in water (3D) that I generate by integrating 2D images (see the 2nd section of my code and the link for the image source). </p> <p>From several example I got confuse if I should use c in steel, c in water or c average of the medium to define:<br /> 1. time array<br /> 2. input_signal<br /> 3. transducer sound speed</p> <p>Thank you.</p> <p>Here is the code:<br /> <pre><code>% ========================================================================= % ULTRASOUND SIMULATION % Linear array transducer (128 element) % Transmit: Planewave, 128 active elements % Receive : 128 active elements % Medium : steel cylinder in water (256 x 512 x 256) % ========================================================================= clear all; %simulation setting data_cast = &#39;single&#39;; % 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 = 1498; % c in water [m/s] rho0 = 1000; % density of water [kg/m^3] c_needle = 5940; % c in steel [m/s] rho_needle = 7800; % density of steel [kg/m^3] c_avg = 3719; % average c of the medium (needle and water) [m/s] rho_avg = 4400; % average density of the medium (needle and water) [kg/m^3] % create the time array t_end = (Nx*dx)*2.2/c0; % [s] kgrid.t_array = makeTime(kgrid, c_needle, [], 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 [SOURCE AND SENSOR] % ========================================================================= % 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]); transducer.sound_speed = c0; % sound speed [m/s] transducer.elevation_focus_distance = 16e-3; % focus distance in the elevation plane [m] % apodization transducer.transmit_apodization = &#39;Hanning&#39;; transducer.receive_apodization = &#39;Rectangular&#39;; % 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 MEDIUM PROPERTIES % ========================================================================= % define a random distribution of scatterers for the medium background_map_mean = 1; background_map_std = 0.008; % % define the properties of the propagation medium load c_needle_water_171103; medium.sound_speed = c_needle_water_171103; % [m/s] load rho_needle_water_171103; medium.density = rho_needle_water_171103; % [kg/m^3] % ========================================================================= % RUN THE SIMULATION % ========================================================================= % set the input settings % plot medium map input_args = {... &#39;PlotLayout&#39;, true, &#39;PMLInside&#39;, false, &#39;PlotSim&#39;, false, &#39;PMLSize&#39;, [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ... &#39;DataCast&#39;, data_cast}; % run the simulation if set to true, otherwise, load previous results from disk if run_simulation RF_data_PW_needle_171104 = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:}); % save the scan lines to disk save RF_data_PW_needle_171104 RF_data_PW_needle_171104; else % % load the scan lines from disk load RF_data_PW_needle_171104 end</code></pre> <p>=======================================================================<br /> and this is the code for creating the medium:</p> <pre><code>clear all close all [data]=imread(&#39;needle_2D.bmp&#39;); data_bw=double(data(:,:,1)); [xx,yy]=size(data_bw); clear data for i=1:xx for j=1:yy if data_bw(i,j)==1 c(i,j)=5940; % average c of the medium (phantom and water) [m/s] rho(i,j)=7800; % average density of the medium (phantom and water) [kg/m^3] elseif data_bw(i,j)==0 c(i,j)=1498; % c in water [m/s] rho(i,j)=1000; % density of water [kg/m^3] end; end; end; cit_3d=zeros(xx,yy,492); for i=1:492 c_3d(:,:,i)=c; rho_3d(:,:,i)=rho; end; n=zeros(236,492,236); m=zeros(236,492,236); [x,y,z]=size(n); for i=1:x for j=1:y for k=1:z c_needle_water_171103(i,j,k)=c_3d(i,k,j); rho_needle_water_171103(i,j,k)=rho_3d(i,k,j); end end end save c_needle_water_171103 c_needle_water_171103 save rho_needle_water_171103 rho_needle_water_171103</code></pre> <p>=============================================================<br /> and here is the image:</p> <p><a href="https://drive.google.com/file/d/1nLgtWFzVAlHxV_gxr_dJjuELNnMxi38h/view?usp=sharing" rel="nofollow">https://drive.google.com/file/d/1nLgtWFzVAlHxV_gxr_dJjuELNnMxi38h/view?usp=sharing</a> </p>