k-Wave User Forum » Topic: Velocity and pressure sources and the returned value

k-Wave User Forum » Topic: Velocity and pressure sources and the returned value http://www.k-wave.org/forum/topic/velocity-and-pressure-sources-and-the-returned-value Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:45:00 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Farnaz on "Velocity and pressure sources and the returned value" http://www.k-wave.org/forum/topic/velocity-and-pressure-sources-and-the-returned-value#post-7579 Mon, 08 Jun 2020 11:54:34 +0000 Farnaz 7579@http://www.k-wave.org/forum/ Dear Brad, Thank a lot for your response. For the artifact I was suspicious of PML size for a while, since in a post you mentioned bigger PML is not necessarily better... I realized there is a problem with speed of sound definition of the medium. I could remove the artifact with modifying the speed of sound, the variation of speed of sound in the medium affects the artifact. If I have an inclusion in which speed of sound drastically changes compared to background such artifact appears e.g. 1300 vs. 1700. About the inclusion shape, changing the distribution of speckles worked to get the whole shape. Best, Farnaz Bradley Treeby on "Velocity and pressure sources and the returned value" http://www.k-wave.org/forum/topic/velocity-and-pressure-sources-and-the-returned-value#post-7549 Sat, 06 Jun 2020 13:19:10 +0000 Bradley Treeby 7549@http://www.k-wave.org/forum/ Hi Farnaz, There is a more detailed description of pressure and velocity sources in the k-Wave user manual, as well as simple examples in the first two <a href="http://www.k-wave.org/documentation/k-wave_time_varying_source_problems.php">Time Varying Source Problems</a> examples (you can open these in MATLAB and run the code). I assume don't see the left and right edges of the inclusion, as the incident plane wave and the detector are placed above it, and the density of the scatterers is constant throughout your phantom (giving uniform speckle). This means there is no energy reflected from these edges back to the transducer. You could try changing the density (or magnitude) of the scatterers within the ellipse. Not sure about the artefact I'm afraid without studying your reconstruction more closely. Does it change if you run the simulation for longer? Brad Farnaz on "Velocity and pressure sources and the returned value" http://www.k-wave.org/forum/topic/velocity-and-pressure-sources-and-the-returned-value#post-7484 Wed, 06 May 2020 10:07:17 +0000 Farnaz 7484@http://www.k-wave.org/forum/ Hi, Since I am struggling with setting a correct setup for my simulation for very long time, I would like to rephrase my question and post the code. I added a b-mode reconstructed image from my simulation here: <img src="https://imagizer.imageshack.com/img922/916/UPZeHn.png" /> There is medium with specific speed of sound and an inclusion with a lower speed of sound on the center which after reconstruction I can only see a reflection from above and below edges in the b-mode image while horizontal reflections are missing and there is an artifact which can be removed if I set the speed of sound in the inclusion to a higher value compared to background. However, in real setup you can see a similar inclusion even with single planwave shot. What am I doing wrong? I have transducer with 384 elements (only 192 elements in middle are activated). I use a plane wave. The pitch of the transducer is 200 micrometer. My medium is a medium of size 7.6*3.8 Cm. Here is the setup of my code: <pre><code>clearvars; tic % ========================================================================= % Define Grid % ========================================================================= % set the size of the perfectly matched layer (PML) pml_y_size = 64; % [grid points] pml_x_size = 96; % [grid points] % set total number of grid points not including the PML Ny = 3072;% Nx = 1536;% % [grid points] dy = 2*3.84e-2/Ny; % grid point spacing in the y direction [m] dx = dy; % grid point spacing in the x direction [m] kgrid = kWaveGrid(Nx, dx, Ny, dy); % ========================================================================= % Define Medium % ========================================================================= % define the properties of the propagation medium medium.sound_speed = 1540; % [m/s] medium.alpha_coeff = 0.5; % [dB/(MHz^y cm)] medium.alpha_power = 1.05; medium.BonA = 10; rho0 = 900; c0 = 1300 ; % create the time array t_end = (Nx * dx) * 2.2 / c0; % [s] kgrid.makeTime(c0, [], t_end); % ========================================================================= % Define Medium % ========================================================================= % define a random distribution of scatterers for the medium background_map_mean = 1; SD = 1 ; num_run = 1 ; medium = repmat(medium,num_run,1) ; rng('shuffle') background_map = (0.9 - 0.03) + (0.06).* sprand(Nx, Ny,0.01); density_map = rho0 * ones(Nx, Ny).* background_map; sound_speed_map = 1540 * ones(Nx,Ny) ; theta = 0; R_x = 150 ; R_y = 150 ; C_x = Nx/2 ; C_y = Ny/2; scattering_region = makeEllipsoid2D([Nx,Ny],[C_x , C_y ], [R_x, R_y],theta) ; % assign region sound_speed_map(scattering_region == 1) = 1300 ; medium.sound_speed = sound_speed_map(:, :); medium.density = density_map(:,:); % ========================================================================= % Define Source middle % ========================================================================= kerf = 1; groupspacing = 7; element_num = 192*2; source_shape = reshape((1:groupspacing)' + (0:element_num-1)*(kerf+groupspacing), 1, []); % % define source mask for a linear transducer with an odd number of elements x_offset = 1; % [grid points] source_m.u_mask = zeros(Nx, Ny); source_m.u_mask(x_offset , source_shape ) = 1 ; %activate middle elements active_elements = element_num/2 ; sensor_m_offset_mask = groupspacing*element_num/4+kerf*(element_num/4)-1 ; %groupspacing*element_num/7+kerf*(element_num/7)-1 ; source_m.u_mask(x_offset , 1:sensor_m_offset_mask) = 0 ; source_m.u_mask(x_offset , Ny-sensor_m_offset_mask:Ny)=0 ; source_strength = 1e6 ; % define the properties of the tone burst used to drive the transducer sampling_freq = 1/kgrid.dt; % [Hz] steering_angle_m = 0; % [deg] tone_burst_freq = 5e6; % [Hz] [rw,t] = ricker(tone_burst_freq,120, kgrid.dt); input_signal = rw ; input_signal = (source_strength ./ (c0 * rho0)) .* input_signal ; source_m.ux = input_signal ; % ========================================================================= % Define Sensor % ========================================================================= sensor_m.mask = source_m.u_mask ; %sensor_m.directivity_angle = ones(Nx,Ny)*pi/2 ; sensor_m.directivity_size = kgrid.dx; input_args = {... 'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size], ... 'DataRecast', true, 'PlotSim', true, 'PlotPML', true , 'PlotScale', 'auto','DataCast', 'gpuArray-single'} ; %,'RecordMovie', true , 'MovieName', 'Final_simulation_1phantome', 'MovieProfile', 'MPEG-4'}; % run the simulation filename_m = sprintf('%s_%d.h5','test') ; kspaceFirstOrder2D(kgrid, medium, source_m,sensor_m, input_args{:},'SaveToDisk',filename_m ) ;</code></pre> I run the c++ code and save the output and again load it here do postprocessing and resample it and reconstruct the 'transducer_data_resample' with a standard delay and sum : <pre><code>%% sensor_data = h5read('test_out.h5','/p') ; size_data = size(sensor_data) ; transducer_data = zeros(active_elements , size_data(2)) ; element = 1 ; % % ========================================================================= % Extract transducer data % ========================================================================= sumnum = groupspacing-1 ; for i = 1:groupspacing:size_data(1)- sumnum transducer_data(element , :) = (sensor_data(i,:)+ sensor_data(i+1,:)+ sensor_data(i+2,:)+ sensor_data(i+3,:)+sensor_data(i+4,:)+ sensor_data(i+5,:)+ sensor_data(i+6,:)) ; element = element+1 ; end originalFs = 1/kgrid.dt ; desiredFs = 40e6; [p,q] = rat(desiredFs / originalFs); %originalFs * p / q ; for i = 1 : active_elements transducer_data_resample(i,:) = resample(transducer_data(i,:),p,q); end</code></pre> additional functions to run the code: ========================================================================= <pre><code>function [rw,t] = ricker(f,n,dt,t0,t1) %RICKER creates an causal ricker wavelet signal % % RICKER creates and plots a default causal ricker wavelet with: % % peak frequency = 20 Hz % sampling time = 0.001 seconds % number of points = 100; % peak location = 1/F = 1/20Hz % % RW = RICKER(...) returns the default wavelet in RW. % % [RW,T] = RICKER(...) returns the time vector in T. % % Specifying parameters of peak frequency (F, Hz), number of points (N), % and sampling time (DT) are specified by the syntax: % % [RW,T] = RICKER(F) % [RW,T] = RICKER(F,N) % [RW,T] = RICKER(F,N,DT) % % [RW,T] = RICKER(F,N,DT,T0) creates a ricker wavelet with peak centered % at T0. % % [RW,T] = RICKER(F,N,DT,T0,T1) creates a 2 dimensional symmetric % ricker wavelet with sift in 1st dimension of T0 and second dimension of % T1. % % Example 1: % ricker % plots a 20 Hz Ricker Wavelet over 0.1 seconds % % Example 2: % % create a ricker wavelet with 40 Hz, 200 points, and 0.02 s between % % samples % [rw,t] = ricker(40,200,0.002); % plot(t,rw), xlabel('Time'), ylabel('Amplitude') % Define inputs if needed switch nargin case 0 f = 20; n = 100; dt = 0.001; t0 = 1/f; is2d = false; case 1 n = 100; dt = 0.001; t0 = 1/f; is2d = false; case 2 dt = 0.001; t0 = 1/f; is2d = false; case 3 t0 = 1/f; is2d = false; case 4 % use all values is2d = false; case 5 % use all inputs is2d = true; otherwise warning('RICKER:tooManyInputs','Ignoring additional inputs') end % Create the wavelet and shift in time if needed T = dt*(n-1); t = 0:dt:T; tau = t-t0; if ~is2d s = (1-tau.*tau*f^2*pi^2).*exp(-tau.^2*pi^2*f^2); else [t1,t2] = meshgrid(tau,t-t1); s = (1-(t1.^2+t2.^2)*f^2*pi^2).*exp(-(t1.^2+t2.^2)*pi^2*f^2); end if nargout == 0 plot(t,s) xlabel('Time (s)') ylabel('Amplitude') title('Ricker Wavelet') else rw = s; end</code></pre> ============= <pre><code>function ellipsoid = makeEllipsoid2D(grid_size, position, radii,theta, plot_ellipsoid, binary) A = theta ; if nargin < 5 || isempty(plot_ellipsoid) plot_ellipsoid = false; end % check for binary input if nargin < 6 || isempty(binary) binary = false; end % force integer grid values Nx = round(grid_size(1)); Ny = round(grid_size(2)); % assign center and radius cx = position(1); cy = position(2); rx = radii(1); ry = radii(2); % check for zero values for center position, and set to middle of grid if cx == 0 cx = floor(Nx / 2) + 1; end if cy == 0 cy = floor(Ny / 2) + 1; end % create grid axes [x_index, y_index] = ndgrid(1:Nx, 1:Ny) ; % create ellipsoid based on code from <a href="https://stackoverflow.com/questions/36420023/generating-a-3d-binary-mask-of-geometric-shapes-in-matlab" rel="nofollow">https://stackoverflow.com/questions/36420023/generating-a-3d-binary-mask-of-geometric-shapes-in-matlab</a> %ellipsoid = ((x_index - cx).^2 / (rx.^2)) + ((y_index - cy).^2 / (ry.^2)) <= 1; ellipsoid = (((x_index-cx)*cos(A) + (y_index-cy)*sin(A)).^2/(rx.^2))+(((x_index-cx)*sin(A) - (y_index-cy)*cos(A)).^2 /ry.^2)<= 1; % convert to double precision if not binary if ~binary ellipsoid = double(ellipsoid); end % plot results if plot_ellipsoid voxelPlot(double(ellipsoid)); end</code></pre> Farnaz on "Velocity and pressure sources and the returned value" http://www.k-wave.org/forum/topic/velocity-and-pressure-sources-and-the-returned-value#post-7476 Mon, 04 May 2020 12:16:16 +0000 Farnaz 7476@http://www.k-wave.org/forum/ Hi, First of all thanks for responding to the questions in the forum, it's really helpful. secondly, I am a bit confused about pressure and velocity sources. How these two are related? and what is the difference between value 'p', 'p_final' that sensor records. I simulate a 2D transducer for plane wave using sensor and source functions. I modified the input shape to a ricker wavelet but the result somehow looks strange. If I use a velocity source and don't set sensor.record, then load '/p' from .h5 file what am I exactly looking at?