k-Wave User Forum » Recent Topics

k-Wave User Forum » Recent Topics http://www.k-wave.org/forum/ Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 19:08:03 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php nakterm on "Getting RF data in k-wave" http://www.k-wave.org/forum/topic/getting-rf-data-in-k-wave#post-9249 Tue, 03 Mar 2026 15:38:06 +0000 nakterm 9249@http://www.k-wave.org/forum/ Hello What is the way to convert the pressure field from k-wave simulation to RF data? Can anyone please let me know. Thanks in advance. xf20228 on "Complete Beginner Question" http://www.k-wave.org/forum/topic/complete-beginner-question#post-9247 Tue, 17 Feb 2026 12:48:38 +0000 xf20228 9247@http://www.k-wave.org/forum/ Hi, I am working through the Homogeneous Propagation Medium Example but get this error when I run the example .m file: Arrays have incompatible sizes for this operation. Error in kspaceFirstOrder2D (line 920) ux_sgx = dt .* rho0_sgx_inv .* real(ifft2( bsxfun(@times, ddx_k_shift_pos, kappa .* fft2(p)) )) / 2; ^ Error in example_ivp_homogeneous_medium (line 67) sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor); ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ So far, I have installed MATLAB R2025a, installed k-wave from add-ons, opened the example_ivp_homogeneous_medium.m file in the MATLAB editor, and run the code. If someone could tell me what I am doing wrong, I would appreciate it. GioF on "Axysimmetric simulations with kspaceFirstOrder-CUDA" http://www.k-wave.org/forum/topic/axysimmetric-simulations-with-kspacefirstorder-cuda#post-9246 Mon, 09 Feb 2026 13:31:57 +0000 GioF 9246@http://www.k-wave.org/forum/ I'm trying to move my simulations from MATLAB to C++ exploiting the kspaceFirstOrder executable. I managed to obtain comparable results with my axysimmetric simulation in MATLAB and the one from kspaceFirstOrder-OMP. I would like now to exploit the CUDA executable to speed up the computations, but I noticed that the axysimmetric simulation is not available for this executable. As moving to 3D will increase the computational burden, I was wondering if you are planning to maka the axysimmetric simulation available soon. Thank you for your feedback! guojia on "The different values of the "cycle" in function " toneBurst" in HIFU simulation" http://www.k-wave.org/forum/topic/the-different-values-of-the-cycle-in-function-toneburst-in-hifu-simulation#post-9245 Thu, 29 Jan 2026 09:16:11 +0000 guojia 9245@http://www.k-wave.org/forum/ Hi everyone, recently, I have been using kWaveArray to simulate the phased array ultrasound focusing sound field. However, when using the function " toneBurst" , I found that the value of "cycle" must be large enough (≥500) to obtain reasonable results; otherwise, the gain of the focused sound field would be too high, which is not reasonable. Could you please tell me what the reason is? The code is as follows: input_args={'DisplayMask',source.p_mask + sensor.mask,... 'DataCast','double',... 'CartInterp','linear',...% 'CartInterp','nearest',... 'PMLSize',[20 20 20],... 'PMLAlpha',5,... 'PMLInside',false}; sensor_data=kspaceFirstOrder3DG(kgrid,medium,source,sensor,input_args{:}); % 求取各阵元延时 % [max_val,max_pos]=max(sensor_data,[],2); % travel_time=max_pos; % time_delays=max(travel_time)-travel_time; sig_rev=fliplr(sensor_data); [max_val,max_pos]=max(sig_rev,[],2); time_delays=max_pos-min(max_pos); focus_mask=source.p_mask; source.p_mask=sensor.mask; p0=0.1e6; % source_cycle2=490; % 能够聚焦，但是结果会发散，声压增益500倍 source_cycle2=500; source.p=1e6*toneBurst(1/kgrid.dt,f0,source_cycle2,'SignalOffset',time_delays,'Envelope','Rectangular'); sensor.mask=zeros(nx,ny,nz); sensor.mask(Index_x-30:Index_x+30,Index_y-30:Index_y+30,:)=1; sensor.record={'p','p_max_all'}; nT=100; sensor.record_start_index=kgrid.Nt-T0/kgrid.dt*nT+1; It is currently unclear what specific reason caused this. I would greatly appreciate it if friends could help figure out exactly what the cause is. Thank you very much. toehead on "Thanks" http://www.k-wave.org/forum/topic/thanks#post-9244 Wed, 28 Jan 2026 17:12:35 +0000 toehead 9244@http://www.k-wave.org/forum/ I just wanted to say thank you for such an amazing toolbox. It was just what I needed and is very user-friendly. Oskar Sadowski on "Time Reversal image reconstruction with kWaveArray elements" http://www.k-wave.org/forum/topic/time-reversal-image-reconstruction-with-kwavearray-elements#post-9243 Tue, 27 Jan 2026 11:00:38 +0000 Oskar Sadowski 9243@http://www.k-wave.org/forum/ Hello, I'm still somewhat new to kWave and am trying to work out Time Reversal image reconstruction for a segmented ring array, created using kWaveArray (i.e. <code>addArcElement</code> and <code>addCustomElement</code>). During testing around, I noticed a unexpected behaviour, and would be grateful for some insights on the matter. I've put together a simple minimal example issued below to explain my problem. The geometry is simply a couple arc shaped elements equally distributed in a circular pattern, all with the same parameters and focus on (0,0). Identical signals (here just <code>toneBurst</code>, the exact signal doesn't matter) are assigned to all elements using <code>sensor.time_reversal_boundary_data</code> I would expect each elements to send out the same signal, as they are identical and the geometry is symmetrical. However, the simulation process shows that the signal intensity differs quite heavily between different elements, with (as far as I can tell) the ones closer to the cartesian axes having a higher intensity assigned than the ones off-axis. The number of elements does not affect this behavior While I'm sure there is a proper explanation for this, that I'm just missing, I can not make sense of this behavior currently, and would be grateful for any advise on how to get the proper results. I should note that a identical result can be achieved when constructing the sensor as a source and setting <code>source.p_mode = 'dirichlet'</code> Minimal example: clear; clearvars; % ========================================================================= % PARAMETERS % ========================================================================= numElements = 7; radius = 1E-2; sizeArc = 1E-3; focus = [0,0]; gridSize = [2.5E-2, 2.5E-2]; gridN = [512, 512]; gridStep = [gridSize(1)/gridN(1), gridSize(2)/gridN(2)]; % ========================================================================= % ELEMENT DEFINITION % ========================================================================= karray = kWaveArray; for i=1:numElements pos = [radius*cosd(i*360/numElements), ... radius*sind(i*360/numElements)]; karray.addArcElement(pos, radius, sizeArc, focus); end % ========================================================================= % GRID & MEDIUM PROPERTIES % ========================================================================= medium.sound_speed = 1520; medium.density = 1040; kgrid = kWaveGrid(gridN(1), gridStep(1), gridN(2), gridStep(2)); kgrid.makeTime(medium.sound_speed) % ========================================================================= % ASSIGN SOURCE & SENSOR % ========================================================================= signal = repmat(toneBurst(1/kgrid.dt, 100E4, 24).*15, numElements, 1); source.p0 = 0; sensor.mask = karray.getArrayBinaryMask(kgrid); sensor.time_reversal_boundary_data = karray.getDistributedSourceSignal(kgrid, signal); % ========================================================================= % RUN SIMULATION % ========================================================================= simulation = kspaceFirstOrder2D(kgrid, medium, source, sensor, ... 'PMLInside', false, 'DataCast', 'gpuArray-single'); figure; imagesc(kgrid.y_vec, kgrid.x_vec, simulation) axis image; colormap("gray") hold on; karray.plotArray(false); scmmw on "HIFU Wave through a heterogenous medium" http://www.k-wave.org/forum/topic/hifu-wave-through-a-heterogenous-medium#post-9242 Tue, 13 Jan 2026 16:35:52 +0000 scmmw 9242@http://www.k-wave.org/forum/ Hi, I'm trying to simulate a HIFU wave propagating through a heterogeneous medium (water, soft tissue, muscle, thyroid, and nodule). My simulation works as expected at 2 MHz, producing a lesion at the focus point. However, at 3 MHz, side lesions appear instead of a lesion at the focus. I'm unsure how to resolve this issue. Any help would be greatly appreciated! Here is my script (I have a lot of visualisations, this is just to help me understand the process) clearvars; % ========================================================================= % ACOUSTIC SIMULATION % ========================================================================= %% define the PML size pml_size = 20; % [grid points] % define the grid parameters Nx = 524 - 2 * pml_size; % [grid points] Ny = 427 - 2 * pml_size; % [grid points] dx = 206e-6; % [m] dy = 206e-6; % [m] % create the computational grid kgrid = kWaveGrid(Nx, dx, Ny, dy); % define the properties of the propagation medium medium.sound_speed = zeros(Nx,Ny); % [m/s] medium.density = zeros(Nx,Ny); % [kg/m^3] medium.alpha_coeff = zeros(Nx,Ny); % [dB/(MHz^y cm)] water.sound_speed = 1480; water.density = 1000; % [kg/m^3] water.alpha_coeff = 0.0022; % [dB/(MHz^y cm)] water.alpha_power = 1; tissue.sound_speed = 1540; tissue.density = 1020; % [kg/m^3] tissue.alpha_coeff = 0.75; % [dB/(MHz^y cm)] tissue.alpha_power = 1.5; muscle.sound_speed = 1588; muscle.density = 1090; % [kg/m^3] muscle.alpha_coeff = 0.63; % [dB/(MHz^y cm)] muscle.alpha_power = 1.3; thyroid.sound_speed = 1585; thyroid.density = 1050; % [kg/m^3] thyroid.alpha_coeff = 1.34; % [dB/(MHz^y cm)] thyroid.alpha_power = 1.5; nodule.sound_speed = 1583; nodule.density = 1080; % [kg/m^3] nodule.alpha_coeff = 0.31; % [dB/(MHz^y cm)] nodule.alpha_power = 1.1; % now break down the grid % --- Water region ------------------------------------------------------- medium.sound_speed(1:156, :) = water.sound_speed; medium.density(1:156, :) = water.density; medium.alpha_coeff(1:156, :) = water.alpha_coeff; % --- Tissue region ------------------------------------------------------- medium.sound_speed(157:341, :) = tissue.sound_speed; medium.density(157:341, :) = tissue.density; medium.alpha_coeff(157:341, :) = tissue.alpha_coeff; %medium.alpha_power(1:185, :) = tissue.alpha_power; % --- Muscle region ------------------------------------------------------- medium.sound_speed(342:416, :) = muscle.sound_speed; medium.density(342:416, :) = muscle.density; medium.alpha_coeff(342:416, :) = muscle.alpha_coeff; %medium.alpha_power(186:260, :) = muscle.alpha_power; % --- Thyroid region ------------------------------------------------------ medium.sound_speed(417:427, :) = thyroid.sound_speed; medium.density(417:427, :) = thyroid.density; medium.alpha_coeff(417:427, :) = thyroid.alpha_coeff; %medium.alpha_power(261:270, :) = thyroid.alpha_power; % --- Nodule region ------------------------------------------------------- medium.sound_speed(427:end, :) = nodule.sound_speed; medium.density(427:end, :) = nodule.density; medium.alpha_coeff(427:end, :) = nodule.alpha_coeff; %medium.alpha_power(271:end, :) = nodule.alpha_power; medium.sound_speed = smooth(medium.sound_speed, true); medium.density = smooth(medium.density, true); medium.alpha_coeff = smooth(medium.alpha_coeff, true); % define the source parameters diameter = 64e-3; % [m] radius = 63.2e-3; % [m] aperture = 20e-3; % [m] freq = 2e6; % [Hz] amp = 0.5e6; % [Pa] %% This may need running a few times to optimise alphafit to get best fit power value i.e. lines overlap as best as possible. alpha_coeff_vec = [water.alpha_coeff, tissue.alpha_coeff, muscle.alpha_coeff, thyroid.alpha_coeff, nodule.alpha_coeff]; alpha_power_vec = [water.alpha_power,tissue.alpha_power, muscle.alpha_power, thyroid.alpha_power, nodule.alpha_power]; c_vec = [water.sound_speed, tissue.sound_speed, muscle.sound_speed, thyroid.sound_speed, nodule.sound_speed]; alphaFit = 1.5; alpha_coeff_fit_vec = fitPowerLawParamsMulti(alpha_coeff_vec, alpha_power_vec, c_vec, freq, alphaFit, true); medium.alpha_power = alphaFit; %% % --- Tissue --- medium.alpha_coeff(1:156, :) = alpha_coeff_fit_vec(1); % --- Tissue --- medium.alpha_coeff(157:341, :) = alpha_coeff_fit_vec(2); % --- Muscle --- medium.alpha_coeff(342:416, :) = alpha_coeff_fit_vec(3); % --- Thyroid --- medium.alpha_coeff(417:427, :) = alpha_coeff_fit_vec(4); % --- Nodule --- medium.alpha_coeff(427:end, :) = alpha_coeff_fit_vec(5); %% close all source.p_mask = zeros(Nx,Ny); % define a focused ultrasound transducer null_mask = makeArc([Nx, Ny], [1, Ny/2], round(radius / dx), round(aperture / dx), [Nx/2, Ny/2]); % define aperture in transducer source.p_mask = source.p_mask + makeArc([Nx, Ny], [1, Ny/2], round(radius / dx), round(diameter / dx), [Nx/2, Ny/2]) - null_mask; % calculate the time step using an integer number of points per period ppw = mean(mean(medium.sound_speed)) / (freq * dx); % points per wavelength cfl = 0.3; % cfl number ppp = ceil(ppw / cfl); % points per period T = 1 / freq; % period [s] dt = T / ppp; % time step [s] % calculate the number of time steps to reach steady state t_end = sqrt( kgrid.x_size.^2 + kgrid.y_size.^2 ) / mean(mean(medium.sound_speed)); Nt = round(t_end / dt); % create the time array kgrid.setTime(Nt, dt); % define the input signal source.p = createCWSignals(kgrid.t_array, freq, amp, 0); % set the sensor mask to cover the entire grid sensor.mask = ones(Nx, Ny); sensor.record = {'p', 'p_max_all'}; % record the last 3 cycles in steady state num_periods = 3; T_points = round(num_periods * T / kgrid.dt); sensor.record_start_index = Nt - T_points + 1; % set the input arguements input_args = {'PMLInside', false, 'PlotPML', false, 'DisplayMask', ... 'off', 'PlotScale', [-1, 1] * amp}; % run the acoustic simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % ========================================================================= % CALCULATE HEATING % ========================================================================= % convert the absorption coefficient to nepers/m alpha_np = db2neper(medium.alpha_coeff, medium.alpha_power) * ... (2 * pi * freq).^medium.alpha_power; % extract the pressure amplitude at each position p = extractAmpPhase(sensor_data.p, 1/kgrid.dt, freq); % reshape the data, and calculate the volume rate of heat deposition p = reshape(p, Nx, Ny); Q = alpha_np .* p.^2 ./ (medium.density .* medium.sound_speed); % Find the location of maximum heat deposition [maxQ, idxMaxQ] = max(Q(:)); [rowMax, colMax] = ind2sub(size(Q), idxMaxQ); % Place the thermal sensor at that location sensor.mask = zeros(Nx, Ny); sensor.mask(rowMax, colMax) = 1; % Report coordinates fprintf('Thermal sensor placed at row=%d, col=%d (x=%.3f m, y=%.3f m)\n', ... rowMax, colMax, kgrid.x_vec(rowMax), kgrid.y_vec(colMax)); % ========================================================================= % THERMAL SIMULATION % ========================================================================= % clear the input structures %clear source sensor; %medium % set the background temperature and heating term source.Q = Q; source.T0 = 37; medium.thermal_conductivity = zeros(Nx,Ny); % [W/(m.K)] medium.specific_heat = zeros(Nx,Ny); % [J/(kg.K)] water.thermal_conductivity = 0.598; % [W/(m.K)] water.specific_heat = 4200; % [J/(kg.K)] tissue.thermal_conductivity = 0.5; % [W/(m.K)] tissue.specific_heat = 3600; % [J/(kg.K)] muscle.thermal_conductivity = 0.49; muscle.specific_heat = 3421; thyroid.thermal_conductivity = 0.5; thyroid.specific_heat = 3826; nodule.thermal_conductivity = 0.5; nodule.specific_heat = 3826; % --- Water region ------------------------------------------------------- medium.thermal_conductivity(1:156, :) = water.thermal_conductivity; medium.specific_heat(1:156, :) = water.specific_heat; % --- Tissue region ------------------------------------------------------- medium.thermal_conductivity(157:341, :) = tissue.thermal_conductivity; medium.specific_heat(157:341, :) = tissue.specific_heat; % --- Muscle region ------------------------------------------------------- medium.thermal_conductivity(342:416, :) = muscle.thermal_conductivity; medium.specific_heat(342:416, :) = muscle.specific_heat; % --- Thyroid region ------------------------------------------------------ medium.thermal_conductivity(417:427, :) = thyroid.thermal_conductivity; medium.specific_heat(417:427, :) = thyroid.specific_heat; % --- Nodule region ------------------------------------------------------- medium.thermal_conductivity(427:end, :) = nodule.thermal_conductivity; medium.specific_heat(427:end, :) = nodule.specific_heat; % create kWaveDiffusion object kdiff = kWaveDiffusion(kgrid, medium, source, sensor); % set source on time and off time on_time = 10; % [s] off_time = 200; % [s] % set time step size dt = 0.1; % Number of timesteps for the entire heating period num_steps = round(on_time / dt); % Initialize a container for CEM43 values at each timestep cem43_values = zeros(1, num_steps); % Initialize variable to store the time when CEM43 exceeds 240 time_cem43_exceeded = NaN; % Default to NaN (not exceeded yet) % Iterate over each timestep for step = 1:num_steps % Take one timestep kdiff.takeTimeStep(1, dt); % Record CEM43 value at the sensor location cem43_values(step) = max(kdiff.cem43(sensor.mask == 1)); % Assuming you want the max CEM43 at the sensor position % Display the CEM43 value and time %fprintf('Time: %.2f s, CEM43: %.2f\n', step * dt, cem43_values(step)); % Check if CEM43 exceeds 240 if cem43_values(step) >= 240 && isnan(time_cem43_exceeded) time_cem43_exceeded = step * dt; % Record the time when it exceeds 240 fprintf('CEM43 exceeded 240 at time %.2f s.\n', time_cem43_exceeded); end end % After the loop, you can use time_cem43_exceeded to check when the threshold was met if ~isnan(time_cem43_exceeded) fprintf('CEM43 exceeded 240 at %.2f seconds during the simulation.\n', time_cem43_exceeded); else disp('CEM43 did not exceed 240 during the simulation.'); end % store the current temperature field T1 = kdiff.T; if max(max(T1)) >= 100 disp('Boiling occured'); end % turn off heat source and take time steps kdiff.Q = 0; kdiff.takeTimeStep(round(off_time / dt), dt); % store the current temperature field T2 = kdiff.T; [max_dose, idx_dose] = max(kdiff.cem43(:)); % Find max thermal dose and index [row_dose, col_dose] = ind2sub(size(kdiff.cem43), idx_dose); % Get coordinates % Get the x and y coordinates of the maximum thermal dose center_x_dose = kgrid.x_vec(row_dose); center_y_dose = kgrid.y_vec(col_dose); % Display the result disp(['Center of Thermal Dose: x = ', num2str(center_x_dose), ' m, y = ', num2str(center_y_dose), ' m']); disp(['Center of Thermal Dose: x = ', num2str(row_dose), ' , y = ', num2str(col_dose), ' ']); %% break % ========================================================================= % VISUALISATION % ========================================================================= % Define the intensity threshold for ablated tissue threshold = 0.99; % Values close to 1 indicate ablated tissue % Create a binary mask of the ablated region binary_mask = kdiff.lesion_map >= threshold; % Find the rows and columns where the ablated tissue exists [row_idx, col_idx] = find(binary_mask); % Add half the grid spacing to account for pixel edges min_x = min(kgrid.x_vec(row_idx)) - abs(kgrid.x_vec(2) - kgrid.x_vec(1)) / 2; max_x = max(kgrid.x_vec(row_idx)) + abs(kgrid.x_vec(2) - kgrid.x_vec(1)) / 2; min_y = min(kgrid.y_vec(col_idx)) - abs(kgrid.y_vec(2) - kgrid.y_vec(1)) / 2; max_y = max(kgrid.y_vec(col_idx)) + abs(kgrid.y_vec(2) - kgrid.y_vec(1)) / 2; % Compute the adjusted dimensions length_ablated = (max_x - min_x) * 1e3; % in mm width_ablated = (max_y - min_y) * 1e3; % in mm % Display the results fprintf('Ablated Tissue Dimensions:\n'); fprintf('Length: %.2f mm\n', length_ablated); fprintf('Width: %.2f mm\n', width_ablated); % Optional: Visualize the lesion map with bounding box figure; imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, kdiff.lesion_map, [0, 1]); colorbar; xlabel('y-position [mm]'); ylabel('x-position [mm]'); axis image; title('Ablated Tissue'); % Use adjusted bounds in the rectangle function rectangle('Position', [min_y * 1e3, min_x * 1e3, ... width_ablated, length_ablated], ... 'EdgeColor', 'r', 'LineWidth', 1.5); % Create time axis for the recorded data t_axis = (0:kdiff.time_steps_taken - 1) * dt; % Plot temperature at every grid point (optional: you can visualize data at a specific point) figure; imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, kdiff.sensor_data(:,:,end)); % Plot temperature at last time step colorbar; xlabel('y-position [mm]'); ylabel('x-position [mm]'); title('Temperature at Last Time Step'); % Plot temperature time profile figure; plot(t_axis, kdiff.sensor_data, 'LineWidth', 1.5); % Enhanced line width for better visibility xlabel('Time [s]'); ylabel('Temperature [^\circC]'); % Add grid lines grid on; % Turns on the grid % Customize y-axis ticks ylim([floor(min(kdiff.sensor_data(:))), ceil(max(kdiff.sensor_data(:)))]); % Ensure y-axis limits encompass the data yticks(floor(min(kdiff.sensor_data(:))):1:ceil(max(kdiff.sensor_data(:)))); % Set ticks with a step of 1 % Optionally add a title title('Temperature Time Profile'); % plot the thermal dose and lesion map figure; % plot the acoustic pressure subplot(2, 3, 1); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, p * 1e-6); h = colorbar; xlabel(h, '[MPa]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Acoustic Pressure Amplitude'); % plot the volume rate of heat deposition subplot(2, 3, 2); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, Q * 1e-7); h = colorbar; xlabel(h, '[kW/cm^2]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Volume Rate Of Heat Deposition'); % plot the temperature after heating subplot(2, 3, 3); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, T1); h = colorbar; xlabel(h, '[degC]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Temperature After Heating'); % plot the temperature after cooling subplot(2, 3, 4); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, T2); h = colorbar; xlabel(h, '[degC]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Temperature After Cooling'); % plot thermal dose subplot(2, 3, 5); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, kdiff.cem43, [0, 1000]); h = colorbar; xlabel(h, '[CEM43]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Thermal Dose'); % plot lesion map subplot(2, 3, 6); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, kdiff.lesion_map, [0, 1]); colorbar; ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Ablated Tissue'); % set colormap and enlarge figure window colormap(jet(256)); scaleFig(1.5, 1); amadou on "Artifact pattern (chevron-like) in 2D (and 3D) k-Wave simulation with linear arr" http://www.k-wave.org/forum/topic/artifact-pattern-chevron-like-in-2d-and-3d-k-wave-simulation-with-linear-arr#post-9241 Thu, 16 Oct 2025 16:39:01 +0000 amadou 9241@http://www.k-wave.org/forum/ Dear k-Wave team, I’m currently running a 2D simulation using kspaceFirstOrder2DG to model ultrasound propagation of a plane wave at 0° from a linear array source in water (L11-5), and I’m observing a strong “chevron-like” interference pattern in the reference medium (homogeneous water). This pattern looks like crossed diagonal bands, even though the medium is fully homogeneous.I define the sources using karraysources class. My grid size is 1/15 of the wavelengt The same simulation were also conducted in 3D considering probe elevation of 5 mm and using kwave_transducer. If the simulatios are run with a kerf (both 2D and 3D) of some grid points, the issue is amplificated and the artifacts become stronger. At lower frequencies, these artifacts are still present but at a smaller scale. Here are the main parameters of my setup: fc = 7.6e6; % Central frequency c_water = 1540; % Sound speed in water rho_water = 1000; % Density ppw = 15; % Points per wavelength dx = floor((c_water/fc/ppw)*1e5)*1e-5; dy = dx; N_elements = 64; pitch_m = 300e-6; elem_width_m = pitch_m; ix_probe = round(1e-3/dx)+1; tone_burst_cycles=3; win = hanning(N_elements); source_signal = toneBurst(1/kgrid.dt, fc, tone_burst_cycles); PML_size = 30; PML_alpha =4; % higher values were also used but no change cfl = 0.15; ` The transducer and sensor are both linear arrays located at ix_probe (which is 1mm away from the boundary). I’m running two simulations (one with scatterers, one reference), then subtracting their pressure fields. However, even in the reference case (homogeneous water), I see this strong interference pattern (see attached figure). My questions are: Is this interference pattern expected given my pitch (300 µm) and frequency (7.6 MHz)? Could it be a spatial aliasing or grating-lobe effect related to the array discretization? Are there recommended values for pitch, element_width, or PML parameters to avoid these artifacts in k-Wave 2D simulations? Would reducing the pitch below λ/2 or adjusting the apodization help mitigate this? Any insight or guidance on proper transducer discretization or grid setup to avoid this pattern would be greatly appreciated. Best regards, Amadou antonio.marzoa on "Issues with beam propagation when comparing two different examples" http://www.k-wave.org/forum/topic/issues-with-beam-propagation-when-comparing-two-different-examples#post-9240 Wed, 08 Oct 2025 01:06:40 +0000 antonio.marzoa 9240@http://www.k-wave.org/forum/ Hi all, I am simulating the propagation in water of a beam generated by a single focused transducer (focal length = 100 mm, aperture diameter = 64 mm, operating frequency = 1.1 MHz, initial pressure = 100 kPa) with different examples from the toolbox. I've done with the "Heating By A Focused Ultrasound Transducer" script and by using the "example_at_focused_bowl_3D" script. I am using the same parameters for the transducer and for the medium in both examples, however, the peak pressure at focus is of 0.5 MPa in the former and of 2.5 MPa when running the later. Does anybody had encounter with this issue before? What could be happening? Thanks! antonio.marzoa on "Obtaining pressure axial plots and issues using example_at_focused_bowl_3D.m" http://www.k-wave.org/forum/topic/obtaining-pressure-axial-plots-and-issues-using-example_at_focused_bowl_3dm#post-9232 Mon, 15 Sep 2025 18:30:04 +0000 antonio.marzoa 9232@http://www.k-wave.org/forum/ Hi everyone, I am trying to simulate the behaviour of a focused spherical ultrasound transducer (64 mm diameter, 100 mm focal length, operating at 1.1 MHz) in water. I am particularly interested in obtaining an axial plot of the pressure and I was using the example_at_focused_bowl_3D.m file to do so, however, it seems that the example does not complete simulating my transducer, as the pressure field map obtained is constant (zero) before reaching the focal position (the simulation is completed, taking almost 2h to complete). Then the axial plot is obtained in which the exact O’Neil solution is perfectly plotted, but the k-wave numerical computation goes to be zero before reaching the focal position. I tried to compare the results with other software’s (like HIFU Beam) but since the k-wave computation is null for most of the region, I couldn’t. The issue didn’t occur with the original example parameters and I am using these for my simulation: % medium parameters c0 = 1500; % sound speed [m/s] rho0 = 1000; % density [kg/m^3] % source parameters source_f0 = 1.1e6; % source frequency [Hz] source_roc = 100e-3; % bowl radius of curvature [m] source_diameter = 64e-3; % bowl aperture diameter [m] source_amp = 0.879e6; % source pressure [Pa] % grid parameters axial_size = 120e-3; % total grid size in the axial dimension [m] lateral_size = 120e-3; % total grid size in the lateral dimension [m] Does anybody have come across something like this before? As an alternative I am running the example_diff_focused_ultrasound_heating.m and I was trying to extract the 2D plot of the axial pressure. Any ideas to do so? Thank you in advance. All the best, Tony sgb2155 on "Instability Simulating a Constant Volume Force" http://www.k-wave.org/forum/topic/instability-simulating-a-constant-volume-force#post-9236 Thu, 18 Sep 2025 22:12:39 +0000 sgb2155 9236@http://www.k-wave.org/forum/ I am running a simulation based on the elastic wave propagation example for pstd3dElastic(). I'm first running my simulation as described in the example to reach a steady state. I am then simulating a longer period, where instead of the velocity input as described in the example, I am simulating a volumetric velocity input based on the steady state (a constant velocity field spanning my whole volume, derived from the acoustic radiation force). I am able to simulate for exactly one time step with this volumetric input before reaching an instability and getting only NaN outputs. I've experimented with lowering my frequency (0.5e6) and reducing my CFL (0.08), but the former didn't get me any further than the first time step with my velocity field, and the latter led to issues with RAM. I have confirmed that my velocity field input contains real values across my simulated time range, so it does appear to me to be a convergence issue to me, but I want to make sure I am headed towards a solution and not completely in the wrong direction. Artyom on "Hydroacoustic, boundary conditions and sound propagation" http://www.k-wave.org/forum/topic/hydroacoustic-boundary-conditions-and-sound-propagation#post-9233 Mon, 15 Sep 2025 20:25:47 +0000 Artyom 9233@http://www.k-wave.org/forum/ Hello everyone, First, let me introduce myself: my name is Artem, and I am a student at the Faculty of Physics. In my research work at the Department of Acoustics, I have been using k-Wave for about six months now to numerically model sound propagation in a coastal wedge. There are a few nuances I would like to clarify, as they are preventing me from fully converging my results with the analytical solution. Here’s the issue: I am currently trying to numerically solve the problem of sound propagation in an ideal fluid wedge with acoustically soft boundaries using k-Wave. The analytical solution for this scenario is presented in the paper by Buckingham and Tolstoy (1990). In k-Wave, I’ve encountered a problem where I cannot set boundary conditions on the wedge edges—specifically, the requirement that the acoustic pressure must be zero at the boundaries. Most likely, this is the reason I am unable to achieve full convergence with the analytical solution. Therefore, I would like to know if it is possible to set such boundary conditions in k-Wave, or whether k-Wave automatically accounts for these boundary conditions based on the defined medium model. Thank you to anyone who responds and helps me resolve this issue! Best regards, Artyom atandri on "Adding elements in kArray" http://www.k-wave.org/forum/topic/adding-elements-in-karray#post-9231 Thu, 04 Sep 2025 20:52:23 +0000 atandri 9231@http://www.k-wave.org/forum/ I am having some trouble creating sensical shapes when rotating kArray 2D elements in a 3D simulation. For example, the addRectElement makes a flat rectangle structure when the angle of rotation is [0,0,0], but with any other angle of rotation, it forms a 3D rectangular prism with + shaped structures on the end when the binary mask is visualized in volumeViewer or similar. I have found this to be the case with other element shapes as well when rotated. Is there a fix for this to obtain the expected rotated structures? Sample code: % ========================================================================= % 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] PML_size = [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE]; % 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] grid_size = [Nx, Ny, Nz]; % calculate the spacing between the grid points dx = 2e-4; % [m] dy = dx; % [m] dz = dx; % [m] % create the k-space grid kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); %% ========================================================================= % DEFINE KARRAY % ========================================================================= karray = kWaveArray; karray.addRectElement([-6e-3,1.5e-3,0],5e-3, 2e-4, [90, 0, 0]); source_mask = karray.getArrayBinaryMask(kgrid); slice(single(source_mask), 127,88,118); % ========================================================================= Edit: My apologies, it seems I have submitted this under the wrong forum subject but there is no way to delete. ukilic on "Concave ring transducer" http://www.k-wave.org/forum/topic/concave-ring-transducer#post-9222 Tue, 15 Jul 2025 22:19:15 +0000 ukilic 9222@http://www.k-wave.org/forum/ Hi everyone, I try to simulate a concave ring transducer. addDisc is good for flat transducers. My transducer is better fit for addBowl but then I cannot put more input, eg. inner diameter. Any ideas? Thanks a lot! BTW- it is single element sonic concept H233 ujjal on "Issue with Accessing k-Wave GPU Code on Ubuntu 24" http://www.k-wave.org/forum/topic/issue-with-accessing-k-wave-gpu-code-on-ubuntu-24#post-9204 Sat, 05 Apr 2025 06:40:53 +0000 ujjal 9204@http://www.k-wave.org/forum/ Hi everyone, We are trying to access MATLAB (k-Wave) GPU code on Ubuntu 24. Our system configuration is as follows: CUDA Device: NVIDIA RTX A4000 CUDA Toolkit Version: 12.2.0 NVIDIA-SMI: 535.230.02 Driver Version: 535.230.02 CUDA Version: 12.2 MATLAB Version: R2024b k-Wave Toolbox Version: 1.4 k-Wave Toolbox Version (CUDA): 1.3 (CPP Linux Binaries) When we run the following command in MATLAB: sensor_data = kspaceFirstOrder2DG(kgrid, medium, source, sensor, input_args{:}, 'RecordMovie', false, 'DataCast', 'single'); We encounter the following error message: ┌───────────────────────────────────────────────────────────────┐ │ kspaceFirstOrder-CUDA v1.3 │ ├───────────────────────────────────────────────────────────────┤ │ Reading simulation configuration: Done │ │ Selected GPU device id: Failed │ └───────────────────────────────────────────────────────────────┘ ┌───────────────────────────────────────────────────────────────┐ │ !!! K-Wave experienced a fatal error !!! │ ├───────────────────────────────────────────────────────────────┤ │ Error: All CUDA-capable devices are busy or unavailable. │ ├───────────────────────────────────────────────────────────────┤ │ Execution terminated │ └───────────────────────────────────────────────────────────────┘ Error using h5readc Unable to open '/tmp/kwave_output_data05-Apr-2025-16-29-44.h5'. File or folder not found. Error in h5read (line 95) [data,var_class] = h5readc(Filename,Dataset,start,count,stride); ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Error in kspaceFirstOrder3DC (line 569) Nx = h5read(output_filename, '/Nx'); ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Error in kspaceFirstOrder2DG (line 76) sensor_data = kspaceFirstOrder3DC(varargin{:}); ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Error in untitled (line 51) sensor_data = kspaceFirstOrder2DG(kgrid, medium, source, sensor, input_args{:}, 'RecordMovie', false, 'DataCast', 'single'); ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Is it not possible to run the k-Wave binaries (related to the Linux version) with CUDA support? If so, how can we resolve this issue? Thank you for your help! abbasK on "Transducer object for elastic simulations" http://www.k-wave.org/forum/topic/transducer-object-for-elastic-simulations#post-9228 Fri, 15 Aug 2025 15:27:02 +0000 abbasK 9228@http://www.k-wave.org/forum/ Hello, can the kWaveTransducer be used with the pstdElastic3D function? masud407 on "2D linear array transducer with karray.addRectElement" http://www.k-wave.org/forum/topic/2d-linear-array-transducer-with-karrayaddrectelement#post-9210 Thu, 24 Apr 2025 17:25:58 +0000 masud407 9210@http://www.k-wave.org/forum/ I am trying to simulate a linear array which I can do for 3D case using the following example: <a href="http://www.k-wave.org/documentation/example_at_linear_array_transducer.php" rel="nofollow">http://www.k-wave.org/documentation/example_at_linear_array_transducer.php</a> However, I was trying to modify it for 2D case as following: clearvars; % ========================================================================= % DEFINE LITERALS % ========================================================================= % select which k-Wave code to run % 1: MATLAB CPU code % 2: MATLAB GPU code % 3: C++ code % 4: CUDA code model = 2; % medium parameters c0 =1540; % sound speed [m/s] cs1 = 5; % shear wave speed alpha0_p1 = 0.002; % absorption coefficient for water % [dB/(MHz^2 cm)] alpha0_s1 = 0.002; % absorption coefficient for water (shear) rho0 = 1000; % density [kg/m^3] %define acrylic layer_tissue cp2 = 2750; % Longitudinal wave speed cs2 = 1150; % shear wave speed rho2 = 1180; %density %alpha0_p2 = 0.5; % absorption coefficient %alpha0_s2 = 0.5; % absorption coefficient alpha0_p2 = 1.58; % absorption coefficient % [dB/(MHz^2 cm)] alpha0_s2 = 0.016; % absorption coefficient % source parameters source_f0 =1e6; % source frequency [Hz] source_amp = 1e6; % source pressure [Pa] source_cycles = 12; % number of toneburst cycles %12 source_focus = 35e-3; % focal length [m] element_num = 90; % number of elements element_width = 0.25e-3; % width [m] element_length = 7e-3; % elevation height [m] element_pitch = 0.3e-3; % pitch [m]0.03mm % transducer position translation = [0, 0]; rotation = [0, 0]; % grid parameters grid_size_x = 40e-3; % [m] grid_size_y = 20e-3; % [m] % grid_size_z = 40e-3; % [m] % computational parameters ppw = 4; % number of points per wavelength t_end = 35e-6; % total compute time [s] cfl = 0.1; % CFL number % ========================================================================= % RUN SIMULATION % ========================================================================= % -------------------- % GRID % -------------------- % calculate the grid spacing based on the PPW and F0 % dx = 2*c0 / (ppw * source_f0); % [m] dx= 1.9250e-04; % compute the size of the grid Nx = roundEven(grid_size_x / dx); Ny = roundEven(grid_size_y / dx); % Nz = roundEven(grid_size_z / dx); % create the computational grid % kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx); kgrid = kWaveGrid(Nx, dx, Ny, dx); % create the time array kgrid.makeTime(c0, cfl, t_end); % -------------------- % SOURCE % -------------------- % set indices for each element if rem(element_num, 2) ids = (1:element_num) - ceil(element_num/2); else ids = (1:element_num) - (element_num + 1)/2; end % set time delays for each element to focus at source_focus time_delays = -(sqrt((ids .* element_pitch).^2 + source_focus.^2) - source_focus) ./ c0; time_delays = time_delays - min(time_delays); time_delays = 1*time_delays; % create time varying source signals (one for each physical element) source_sig = source_amp .* toneBurst(1/kgrid.dt, source_f0, source_cycles, 'SignalOffset', round(time_delays / kgrid.dt)); % create empty kWaveArray karray = kWaveArray('BLITolerance', 0.05, 'UpsamplingRate', 10); Unfortunately, I am getting the following error: Expected input number 4, theta, to be an array with number of elements equal to 1. Error in kWaveArray/addRectElement (line 715) validateattributes(theta, {'numeric'}, {'finite', 'real', 'numel', theta_length}, 'addRectElement', 'theta', 4); Any suggestion on resolving it is highly appreciated. ayitime on "Delay Between Symmetric Sensors Despite Symmetric Setup in k-Wave" http://www.k-wave.org/forum/topic/delay-between-symmetric-sensors-despite-symmetric-setup-in-k-wave#post-9223 Wed, 16 Jul 2025 10:25:29 +0000 ayitime 9223@http://www.k-wave.org/forum/ Hi, I'm running a 2D simulation in k-Wave using a line source placed symmetrically along the x-axis and a symmetrical arc-shaped sensor array (e.g., 6 elements). Both the source and sensor are symmetric in geometry. However, although symmetric sensor elements visibly receive the wavefront at the same time in the simulation movie, their recorded signals in combined_sensor_data show a delay of about 200 time steps (out of 1200), which is unexpected. I’m wondering: What might cause this delay between symmetric receivers? Could it be due to the use of kWaveArray, combineSensorData(), or the element size? Any ideas would be appreciated. Thanks! =================== code ================== clear; clc; Nx = 256; Ny = 256; dx = 0.5e-3; dy = dx; %% grid kgrid = kWaveGrid(Nx, dx, Ny, dy); %% medium medium.sound_speed = 1500; kgrid.makeTime(medium.sound_speed); %% source [array type] karrayy = kWaveArray; x_bias = -20e-3; y_width = 50e-3; start_pos = [x_bias, -y_width / 2]; end_pos = [x_bias, +y_width / 2]; karrayy.addLineElement(start_pos, end_pos); source.p_mask = karrayy.getArrayBinaryMask(kgrid); imagesc(kgrid.x_vec*1000, kgrid.y_vec*1000, source.p_mask); xlabel("y/mm"); ylabel("x/mm"); axis equal; axis tight; colorbar; f1 = 100e3; % 0.1MHz sig1 = toneBurst(1/kgrid.dt, f1, 3); source.p = sig1; % plot(sig1); %% sensor [array type] karray = kWaveArray; radius = 50e-3; num_elements = 6; elem_pos = makeCartCircle(radius, num_elements, [0, 0]); radius_of_curv = 10e-3; diameter = 1e-3; focus_pos = [0, 0]; for ind = 1:num_elements karray.addArcElement(elem_pos(:, ind), radius_of_curv, diameter, focus_pos); end sensor.mask = karray.getArrayBinaryMask(kgrid); imagesc(sensor.mask); xlabel("y/mm"); ylabel("x/mm"); axis equal; axis tight; colorbar; input_args = {'DisplayMask',(source.p_mask | sensor.mask),... 'RecordMovie', true,... 'MovieName','video2',... 'DataCast', 'gpuArray-single',... 'PlotSim',true}; %% run sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); combined_sensor_data = karray.combineSensorData(kgrid, sensor_data); stackedPlot(combined_sensor_data); Crusher on "I want to do a 2D elastic wave simulation on step geometry" http://www.k-wave.org/forum/topic/i-want-to-do-a-2d-elastic-wave-simulation-on-step-geometry#post-8920 Tue, 03 Oct 2023 11:53:34 +0000 Crusher 8920@http://www.k-wave.org/forum/ I want to do a 2D elastic wave simulation on step geometry, like a small rectangle of different property placed on a big rectangle of different property. How to do that? NicholasHorton on "Seems attenuation settings (following power law) affect signals' time-of-arrival" http://www.k-wave.org/forum/topic/seems-attenuation-settings-following-power-law-affect-signals-time-of-arrival#post-9221 Mon, 14 Jul 2025 07:59:55 +0000 NicholasHorton 9221@http://www.k-wave.org/forum/ Hi everyone, I was working on the simulation of ultrasound propagation. I designed three homogeneous mediums. Their speed-of-sounds were 1500 m/s. Their densities were 1000 kg/m^3. Their alpha_power were set as 1.1. Their alpha_coeff were 0.002, 0.7, and 3.0, respectively. A sensor was set to record transmitted signal. I just found that with different alpha_coeff, the time-of-arrivals of the signals were different. The larger alpha_coeff was, the earlier the signal arrived. The change was so significant that it was counter-intuitive. I wonder if dispersion can explain this phenomenon. I would like to ask you, is this simulation result consistent with the physical phenomenon? If not, how should I correct my code? The code has been attached below. Thanks a lot. <pre><code>% To compare TOFs of transmitted signals, with different [constant] att.s or wvfroms. % * The 2-D numerical simulations are based on k-Wave toolbox % * Emitted pulse: a 3-cycle tone-burst whose centre-frequency is 2 MHz % * Emitter and receiver: two points about 25 cm apart, recording transmitted signal % * Medium: homogeneous, whose attenuation power is 1.1, and attenuation coefficient % is 0.002, 0.7 or 3.0 dB/cm/MHz^1.1 clc;clear close all %% Define Params % default medium params sos_ref = 1500; % [m/s] density_ref = 1000; % [kg/m^3] att_power = 1.1; % y % define k-Wave Grid Nx = 3072; dx = 1e-4; kgrid = kWaveGrid(Nx, dx, Nx, dx); t_end = Nx * dx / sos_ref * 2; kgrid.makeTime(sos_ref, 0.5, t_end); % CFL = 0.5 % different wvforms and att.s small_att = 0.002; % dB/cm/MHz^y big_att = 0.7; bigger_att = 3.0; wv_tb = toneBurst(1./kgrid.dt, 2e6, 3, 'SignalLength', kgrid.Nt); %% Define numerical phantom %------------define Medium One, no attenuation at all medium_noAtt.sound_speed = sos_ref * ones(Nx); medium_noAtt.density = density_ref * ones(Nx); %------------define Medium Two, with a small atten. coeff., medium_smallAtt = medium_noAtt; medium_smallAtt.alpha_power = att_power; medium_smallAtt.alpha_coeff = small_att * ones(Nx); %------------define Medium Three, with a big atten. coeff., medium_bigAtt = medium_smallAtt; medium_bigAtt.alpha_coeff = big_att * ones(Nx); %------------define Medium Three, with a big atten. coeff., medium_biggerAtt = medium_smallAtt; medium_biggerAtt.alpha_coeff = bigger_att * ones(Nx); %% Define source & sensor source.p_mask = zeros(Nx, Nx); source.p_mask(Nx / 2, 250) = 1; sensor.mask = zeros(Nx, Nx) sensor.mask(Nx / 2, Nx - 250) = 1; % source.p not defined %% Execute the simulations input_args = {'PMLSize', 'auto',... 'PMLInside', false,... 'PlotPML', false,... 'DisplayMask', 'off',... 'DataPath', pwd,... 'DeleteData', false}; for medium_name = {'smallAtt', 'bigAtt', 'biggerAtt'} cmd = ['medium = medium_', medium_name{1}, ';']; eval(cmd) source.p = wv_tb; rf_data = kspaceFirstOrder2DC(kgrid, medium, source, sensor, input_args{:}); cmd = ['sig_', medium_name{1}, ' = rf_data;']; eval(cmd) end %% Visulaization f1 = figure('Units','centimeters', 'Position', [2, 2, 25, 40]); plt_cnt = 0; plt_c = 3; plt_r = 1; range_ = [160, 175] * 1e-6; % unit: s for medium_name = { 'smallAtt', 'bigAtt', 'biggerAtt'} cmd = ['sig = sig_', medium_name{1}, ';']; eval(cmd) env = abs(hilbert(sig)); % ---------- plt_cnt = plt_cnt + 1; subplot(plt_c, plt_r, plt_cnt) plot(kgrid.t_array * 1e6, env, 'color','r', 'lineWidth', 1.5) grid on grid minor xlabel('Time [\mus]') xlim(range_ * 1e6) title(medium_name{1}, 'Interpreter', 'none') end sgtitle('Signal Envelope')</code></pre>