k-Wave User Forum » Recent Posts

k-Wave User Forum » Recent Posts http://www.k-wave.org/forum/ Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 19:15:44 +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! atandri 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-9239 Tue, 30 Sep 2025 16:06:13 +0000 atandri 9239@http://www.k-wave.org/forum/ It seems as though you should be able to do this similar to the "Simulating Ultrasound Beam Patterns" example- essentially by setting the whole grid as sensors. 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. 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-9234 Tue, 16 Sep 2025 10:32:02 +0000 antonio.marzoa 9234@http://www.k-wave.org/forum/ Hi again, I have sorted the issue of the zeroing of the pressure (it was just an issue regarding reaching the steady state - thus, not having put enough computation time! -). I am still looking for a way of obtaining the axial pressure plots from the example_diff_focused_ultrasound_heating.m, any ideas? 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 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 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. atandri on "Concave ring transducer" http://www.k-wave.org/forum/topic/concave-ring-transducer#post-9230 Thu, 04 Sep 2025 20:25:39 +0000 atandri 9230@http://www.k-wave.org/forum/ the kWaveArray addAnnularArray seems to exactly fit what you're looking for nmin 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-9229 Mon, 01 Sep 2025 09:41:34 +0000 nmin 9229@http://www.k-wave.org/forum/ Hi ujjal, I am experiencing the same issue. Could you kindly share how you resolved it? Thanks a lot! 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? mhopeng on "2D linear array transducer with karray.addRectElement" http://www.k-wave.org/forum/topic/2d-linear-array-transducer-with-karrayaddrectelement#post-9227 Sat, 02 Aug 2025 01:08:13 +0000 mhopeng 9227@http://www.k-wave.org/forum/ You need to define the elements of the Array using one of the methods for adding elements, such as <code>addRectElement()</code>. I don't see that in the code you have posted here. 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-9226 Wed, 30 Jul 2025 07:40:21 +0000 ayitime 9226@http://www.k-wave.org/forum/ This question has been resolved. frezaei 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-9225 Tue, 29 Jul 2025 02:59:04 +0000 frezaei 9225@http://www.k-wave.org/forum/ I have a follow-up question: when I change the frequency, I observe that the sensor signal is higher for a higher frequency than lower frequency. Why does this happen? The material has an attenuation coefficient that depends on frequency, so for higher frequencies, we should see a lower amplitude at the sensor. Why do I see the reverse behavior? I also tried to use the same number of cycles for each case, but again seeing the same behavior, I am interested to learn how the following papers simulated with reasonable results? paper1.K‑wave modelling of ultrasound wave propagation in aerogels and the effect of physical parameters on attenuation and loss paper2.Simultaneous use of pulse-echo and through-transmission methods in determining a combined reflection coefficient frezaei 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-9224 Mon, 28 Jul 2025 19:05:11 +0000 frezaei 9224@http://www.k-wave.org/forum/ Hello, I am new to k-wave simulation, and I am trying to simulate the attenuation through media. First of all, for medium.alpha_coeff, what should the unit be? dB/MHz.cm or Np/MHz.cm? Secondly, how do I make sure the attenuation through the medium is correct? I placed two sensors, one near the source and the other 1.1mm away, and now by looking at sensor signals, I see there is an attenuation, but if I consider the peak-to-peak value of the sensor signal, it does not follow this relationship: P=P0*exp(-alpha*f*x), why? Here is my code: % Define medium properties medium.sound_speed = zeros(Nx, Ny); medium.density = zeros(Nx, Ny); medium.sound_speed(:,:) = 5555; medium.density(:,:) = 2440; % Define attenuation coefficients (Np/(MHz^y cm)) medium.alpha_coeff = zeros(Nx, Ny); % in dB/(MHz^y cm) but K-Wave expects Np/(MHz^y m) medium.alpha_coeff(:,:) = 0.1; % low loss (converted to Np/(MHz^y m)) medium.alpha_power = 1.0; % frequency power medium.alpha_mode = 'no_dispersion'; % Simulation time cfl=0.3; t_end = 16e-6; % [s] sound_speed=5555; dt = cfl * dx / sound_speed; kgrid.makeTime(sound_speed, cfl,t_end) % Tone burst source settings sample_freq = 20e6; % [Hz] signal_freq = 1e6; % [Hz] num_cycles = 3; signal=toneBurst(sample_freq,signal_freq,num_cycles,'Plot',true); %Defining a single source point source.p_mask = zeros(Nx, Ny); source.p_mask(round(Nx/2), 1) = 1; % define a time varying sinusoidal source source_mag = 1e2; % Amplitude in Pascals source.p=source_mag *signal; % Initialize sensor mask with unique labels sensor.mask = zeros(Nx, Ny); % Sensor 1: label '1' in second row (just below source) sensor.mask(round(Nx/2), 2) = 1; % Sensor 2: label '2' in the row right after the glass layer sensor.mask(round(Nx/2), 24) = 1; sensor.record = {'p', 'p_final','I'}; % Run the simulation with memory-safe options input_args = {'PMLInside', false, 'PMLSize', 20, 'PMLAlpha',3, 'PlotPML', true, 'DataCast', 'single', 'PlotLayout', true,'RecordMovie', true}; % input_args = {'DisplayMask', display_mask,'PMLInside', false, 'PMLSize', 20, 'PMLAlpha',2, 'PlotPML', true, 'DataCast', 'single', 'PlotLayout', true}; sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,input_args{:}); % Extract pressure signals p_signals = sensor_data.p; % [2 x Nt] Nt = length(kgrid.t_array); t_micro = (0:Nt-1) * kgrid.dt * 1e6; % Plot recorded signals figure; plot(t_micro, p_signals(1,:), 'b-', 'LineWidth', 1.5); xlabel('Time (µs)'); ylabel('Pressure (Pa)'); legend('Sensor 1 (near source)'); title('Pressure1 Time‑Series Recorded by Sensors'); figure; plot(t_micro, p_signals(2,:), 'r-', 'LineWidth', 1.5); xlabel('Time (µs)'); ylabel('Pressure (Pa)'); legend('Sensor 2 (after 1.1mm)'); title('Pressure2 Time‑Series Recorded by Sensors'); grid on; [row_sens, col_sens] = find(sensor.mask); figure; imagesc(x_mm, y_mm, medium.sound_speed'); axis image tight; colormap(jet); colorbar; xlabel('x (mm)'); ylabel('y (mm)'); title('Medium with Sensor and Source Locations'); hold on; % Plot source in red [row_src, col_src] = find(source.p_mask); plot(x_mm(row_src), y_mm(col_src), 'rs', 'MarkerSize', 8, 'LineWidth', 1.5); % Plot sensors in green plot(x_mm(row_sens), y_mm(col_sens), 'go', 'MarkerSize', 8, 'LineWidth', 1.5); legend('Source', 'Sensor'); 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); 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 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> mcdannold on "Issue with XXX_output.h5 files" http://www.k-wave.org/forum/topic/issue-with-xxx_outputh5-files#post-9219 Mon, 07 Jul 2025 15:35:39 +0000 mcdannold 9219@http://www.k-wave.org/forum/ Hello, The computing cluster I run kWave code on was recently migrated to RedHat from centOS. The compiled kspaceFirstOrder3DC code seems to be running fine, but there seems to be an issue with the output temporary files. The XXX_output.h5 files are only 1 kB, and there are additional files starting with '.nfs' that are generated (such as '.nfs7eed9ff995d7a06200000304'). These files don't seem to be deleted. I am using matlab2024b. I used to use an older version, but that version was not installed after the OS migration. Any idea what is going on? Thanks, Nathan ejoa1234561 on "Linux vs windows GPU memory allocation" http://www.k-wave.org/forum/topic/linux-vs-windows-gpu-memory-allocation#post-9218 Wed, 28 May 2025 04:03:52 +0000 ejoa1234561 9218@http://www.k-wave.org/forum/ Hello, I am running this code in 2 workstations. The first workstation is a Windows 10 with a GPU ADA A4000 of 12gb VRAM. The second is a cluster linux with 2 GPUs A100 of 48gb of VRAM each. In the W10 the code can allocate until t_end = 0.12; (11.48gb of 12gb) without problems. However, when I try to do the same in the cluster the VRAM only accepts (t_end=0.1) using 39.8gb of 40gb. What is happening? Linux manages different the VRAM? Is there a way to fix this? I installed the binaries in both cases and the problem is on pstdElastic2D(kgrid, medium, source, sensor, input_args{:}) function. clc; clear; close all; gpuDevice(1) format compact; set(0,'defaultAxesFontSize',18); set(groot, 'defaultAxesTickLabelInterpreter','tex'); set(groot, 'defaultLegendInterpreter','tex'); set(0,'defaulttextInterpreter','tex'); DATA_CAST = 'gpuArray-single'; %% Create the computational grid Nx = 60; % number of grid points in the x direction Ny = 100; % number of grid points in the y direction dx = 2e-3; % grid spacing in x [m] dy = 2e-3; % grid spacing in y [m] kgrid = kWaveGrid(Nx, dx, Ny, dy); %% Define medium properties medium.sound_speed_compression = 1500 * ones(Nx, Ny); % [m/s] realistic compressional speed medium.sound_speed_shear = 2.5 * ones(Nx, Ny); % [m/s] for gelatin-like shear medium.density = 1000 * ones(Nx, Ny); % [kg/m^3] % Add boundary rigid_thickness = 5; %AIR medium.density(:,1:rigid_thickness) = 12; % medium.density(:,end-rigid_thickness:end) = 1.2; medium.density(1:rigid_thickness,:) = 12; medium.density(end-rigid_thickness:end,:) = 12; medium.sound_speed_shear(:, 1:rigid_thickness) = 1; % medium.sound_speed_shear(:,end-rigid_thickness:end) = 0.1; medium.sound_speed_shear(1:rigid_thickness,:) = 1; medium.sound_speed_shear(end-rigid_thickness:end,:) = 1; medium.sound_speed_compression(:, 1:rigid_thickness) = 1500; % medium.sound_speed_compression(:,end-rigid_thickness:end) = 340; medium.sound_speed_compression(1:rigid_thickness,:) = 1500; medium.sound_speed_compression(end-rigid_thickness:end,:) = 1500; %TX medium.density(1:rigid_thickness,40:60) = 12; medium.sound_speed_shear(1:rigid_thickness,40:60) = 1000; medium.sound_speed_compression(1:rigid_thickness,40:60) = 1500; % Attenuation medium.alpha_coeff_compression = 0.05; % [dB/(MHz^2 cm)] medium.alpha_coeff_shear = 1000; % [dB/(MHz^2 cm)] % medium.alpha_power_shear = 1.3; %% Time array (safe dt based on c_max) t_end = 0.1; %0.12; % [s] to allow steady-state observation c_max = max([max(medium.sound_speed_compression(:)), max(medium.sound_speed_shear(:))]); cfl = 0.3; kgrid.makeTime(c_max, cfl, t_end); %% Define source cx1 = Nx/2; cy1 = Ny-14; radius = 10; plot_disc = 'true'; disc1 = makeDisc(Nx, Ny, cx1, cy1, radius, plot_disc); source.u_mask = disc1; source_freq = 200; % [Hz] source_mag = 1e-6; source_signal = source_mag * sin(2 * pi * source_freq * kgrid.t_array); source.uy= source_signal; % source_points = find(source.u_mask); % num_sources = length(source_points); % source.ux = zeros(length(kgrid.t_array), num_sources); % for i = 1:num_sources % source.ux(:, i) = source_signal; % end %% Define sensor mask = zeros(Nx, Ny); mask(:, 10:end) = 1; sensor.mask = mask; sensor.record = {'u', 'u_split_field'}; %% Run simulation input_args = {'PMLAlpha', 2, 'PlotPML', false, 'PMLInside', false, ... 'DisplayMask', 'off', 'DataCast', DATA_CAST}; sensor_data = pstdElastic2D(kgrid, medium, source, sensor, input_args{:}); ux_s = gather(reshape(sensor_data.ux_split_s, Nx, Ny-9, [])); ux_p = gather(reshape(sensor_data.ux_split_p, Nx, Ny-9, [])); ux_total = gather(reshape(sensor_data.ux, Nx, Ny-9, [])); dinf.dx = dx; dinf.dy = dy; dinf.dt = kgrid.dt; %% Visualization figure xx = dinf.dx * (0:size(ux_total,1)-1); yy = dinf.dy * (0:size(ux_total,2)-1); for ii = 1:500:20000 im = real(squeeze(ux_total(:,:,ii))); % amplitude of shear wave imagesc(yy(10:end),xx, im); xlabel('x (m)'); ylabel('y (m)'); axis image; clim([-1 1]); title(['Frame ' num2str(ii)]); grid on; colormap('jet'); colorbar; drawnow pause(0.1) end % Optional: save filename = ['ShearReflection_' num2str(source_freq) 'Hz']; save(filename, "ux_total", "ux_p", "ux_s", "dinf", "source_freq", '-v7.3'); energy_t = zeros(1, size(ux_s,3)); for ii = 1:size(ux_s,3) frame = ux_s(:,:,ii); energy_t(ii) = sum(frame(:).^2); % proportional to elastic energy end % 2. Difference between consecutive frames diff_energy = zeros(1, size(ux_s,3)-1); for ii = 1:(size(ux_s,3)-1) diff = ux_s(:,:,ii+1) - ux_s(:,:,ii); diff_energy(ii) = sum(diff(:).^2); end % 3. Plot energy over time figure; subplot(2,1,1); plot(kgrid.t_array, energy_t, 'b-', 'LineWidth', 2); xlabel('Time [s]'); ylabel('Total Shear Energy'); title('Shear Field Energy vs Time'); grid on; % 4. Plot frame-to-frame difference subplot(2,1,2); plot(kgrid.t_array(2:end), diff_energy, 'r-', 'LineWidth', 2); xlabel('Time [s]'); ylabel('Frame-to-Frame Energy Change'); title('Difference Between Consecutive Shear Frames'); grid on; Error output: Running k-Wave elastic simulation... start time: 27-May-2025 22:34:42 reference sound speed: 1500m/s [Warning: Support for ver('distcomp') will be removed in a future release. Use ver('parallel') instead.] [> In ver>locGetSingleToolboxInfo (line 283) In ver (line 56) In verLessThan (line 39) In kspaceFirstOrder_inputChecking (line 1306) In pstdElastic2D (line 385) In sphere (line 104)] dt: 400ns, t_end: 150ms, time steps: 375001 input grid size: 60 by 100 grid points (120 by 200mm) maximum supported compressional frequency: 375kHz maximum supported shear frequency: 250Hz expanding computational grid... computational grid size: 100 by 140 grid points {Error using gpuArray.zeros Out of memory on device. To view more detail about available memory on the GPU, use 'gpuDevice()'. If the problem persists, reset the GPU by calling 'gpuDevice(1)'. Error in kspaceFirstOrder_inputChecking>@(sz)gpuArray.zeros(sz,'single') (line 1314) castZeros = @(sz) gpuArray.zeros(sz, 'single'); Error in kspaceFirstOrder_createStorageVariables (line 336) sensor_data.uy_split_s = castZeros([num_sensor_points, num_recorded_time_points]); Error in kspaceFirstOrder_inputChecking (line 1663) kspaceFirstOrder_createStorageVariables; Error in pstdElastic2D (line 385) kspaceFirstOrder_inputChecking; Error in sphere (line 104) sensor_data = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});} Mekdes Bezabh on "Regarding Homogenous B-mode image simulation" http://www.k-wave.org/forum/topic/regarding-homogenous-b-mode-image-simulation#post-9217 Mon, 26 May 2025 21:37:15 +0000 Mekdes Bezabh 9217@http://www.k-wave.org/forum/ Please, I wonder if anyone can suggest any solution. Thank you in advance. yanyan on "超声传感器阵列显示The defined transducer is too large or positioned outside the grid in" http://www.k-wave.org/forum/topic/%e8%b6%85%e5%a3%b0%e4%bc%a0%e6%84%9f%e5%99%a8%e9%98%b5%e5%88%97%e6%98%be%e7%a4%bathe-defined-transducer-is-too-large-or-positioned-outside-the-grid-in#post-9216 Mon, 26 May 2025 11:25:23 +0000 yanyan 9216@http://www.k-wave.org/forum/ clear dx =5e-6; % x方向分辨率 dy =5e-6; % y方向分辨率 PML_size=20;%设置完美匹配层 % 阵列参数 element_size = [12, 12]; % x/y方向各6网格点，z方向1点 element_spacing = [1, 1]; % x/y间距1点，z方向无间隔 array_dims = [80, 1]; % x/y各80阵元，单层 % 总网格数 Nx=array_dims(1)*(element_size(1)+element_spacing(1)) - element_spacing(1)+2*PML_size;%x方向总网格数量 Ny=Nx;%y方向总网格数量 kgrid=kWaveGrid(Nx, dx, Ny, dy); %设置介质材料属性 sound_speed_glass = 5600; % 玻璃声速 [m/s] density_glass = 2400; % 玻璃密度 [kg/m^3] sound_speed_water = 1500; % 水声速 density_water = 1000; % 水密度 sound_speed_air =343; % 空气声速 density_air = 1; % 空气密度 medium.sound_speed = sound_speed_glass*ones(Nx,Ny);%声速初始化 medium.density = density_glass*ones(Nx,Ny); %密度初始化 reflector_y_pos = 1+(PML_size + round(1e-3 / kgrid.dy));%反射面设置为在21个，避免反射波对计算区域造成影响 % 定义交替层的位置 x_start = PML_size + 1; %定义开始位置 x_end = Nx-PML_size; %定义结束位置 interval_num = 20; %分成20个区域 interval_grid_points = floor((Nx-2*PML_size+1) / interval_num); % 计算每个区域的网格点数，并使用floor函数取整 remaining_points =Nx-2*PML_size - interval_num * interval_grid_points; % 剩余的网格点数 target_column =x_start; interval_num=20; i = zeros(interval_num,2); %网格单元的长度 for interval_idx = 0:(interval_num-1) current_x_start = x_start + interval_idx * interval_grid_points; current_x_end = current_x_start + interval_grid_points - 1; i(interval_idx+1,1)=current_x_start; i(interval_idx+1,2)=current_x_end; % 最后一层分配剩余点数 if interval_idx == interval_num - 1 current_x_end = current_x_end + remaining_points; end % 边界保护 if current_x_end > x_end current_x_end = x_end; end % 交替材料设置 if mod(interval_idx, 2) == 0 medium.sound_speed(current_x_start:current_x_end, :, reflector_y_pos:end) = sound_speed_water; medium.density(current_x_start:current_x_end, :, reflector_y_pos:end) = density_water; else medium.sound_speed(current_x_start:current_x_end, :, reflector_y_pos:end) = sound_speed_air; medium.density(current_x_start:current_x_end, :, reflector_y_pos:end) = density_air; end end imagesc(medium.sound_speed), %设置换能器的位置和检测传感器的位置 tr.element_width = 12; % 换能器宽度所占网格数 tr.element_length= 1;% 换能器长度所占网格数 tr.element_spacing =1; % 换能器间隔所占网格数 tr.number_elements= 80; % 换能器数量 tr.radius=Inf; transducer_width = tr.number_elements*tr.element_width+(tr.number_elements - 1) * tr.element_spacing; tr.position = round([Nx/2 - transducer_width/2,Ny-PML_size]); tr.active_elements=zeros(tr.number_elements,1); tr.active_elements(1:80)=1; % 定义激励超声源信号 %transducer.mask = transducer_mask; source_strength = 1e6; % [Pa] tone_burst_freq = 10e6; % [Hz] tone_burst_cycles = 5; t_end = 1e-6; % 设置仿真结束时间为10个周期 kgrid.makeTime(medium.sound_speed, [], t_end); % 使用kgrid的makeTime方法生成时间数组 % create the input signal using toneBurst input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles); %input_signal=(source_strength./(medium.sound_speed* medium.density)).*input_signal; tr.input_signal=input_signal; transducer=kWaveTransducer(kgrid, tr); sensor_data= kspaceFirstOrder2D(kgrid, medium, transducer,transducer); Mekdes Bezabh on "Regarding Homogenous B-mode image simulation" http://www.k-wave.org/forum/topic/regarding-homogenous-b-mode-image-simulation#post-9215 Wed, 21 May 2025 18:15:24 +0000 Mekdes Bezabh 9215@http://www.k-wave.org/forum/ Dear All, I have some questions regarding the code "example_us_bmode_linear_transducer." I am trying to create a homogeneous B-mode image without any artifacts. I have tested different values for the tone burst frequency—3.47 MHz, 7.8 MHz, and 15.62 MHz—along with varying values for the medium's alpha coefficient and alpha power, while keeping the sound speed and density uniform. I am using a plane wave setup with the focus set to infinity. However, I encountered some unexpected results: 1. At a frequency of 3.47 MHz, the resulting B-mode image was not as expected. 2. At 7.8 MHz, I observed a bright line in the middle of the final B-mode image. 3. The same issue occurred at 15.62 MHz. I would appreciate your assistance in determining whether my code is correct, as this simulation software is crucial for my project. I would also like to share my code for your review if you have some time to look at it and provide your valuable suggestions. Thank you for your attention. The code: 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 = 1; % [grid points] % set total number of grid points not including the PML Nx = 256*2 - 2 * pml_x_size; % [grid points] Ny = 128*2 - 2 * pml_y_size; % [grid points] Nz = 2; % [grid points] %128 - 2 * pml_z_size; x= 40e-3; dx = x / Nx; dy = dx; dz = dx; % [m] % create the k-space grid kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); % ========================================================================= % DEFINE THE MEDIUM PARAMETERS % ========================================================================= % define the properties of the propagation medium c0 = 1540; % [m/s] rho0 = 1000; % [kg/m^3] % medium.alpha_mode = 'no_dispersion'; medium.alpha_coeff = 0.5; % [dB/(MHz^y cm)] medium.alpha_power = 1.0; % create the time array t_end = (Nx * dx) * 2.2 / c0; % [s] kgrid.makeTime(c0, [], t_end); % ========================================================================= % DEFINE THE INPUT SIGNAL % ========================================================================= % define properties of the input signal source_strength = 1e6; % [Pa] tone_burst_freq = 15.62e6; % [Hz] 7.8 MHz, 3.47 MHz tone_burst_cycles = 4; % 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 = 2; % width of each element [grid points] transducer.element_length = 1; %24; % 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 = Inf; %40e-3; % focus distance [m] % transducer.elevation_focus_distance = Inf; %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 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 = kWaveTransducer(kgrid, transducer); % print out transducer properties transducer.properties; % ========================================================================= % DEFINE THE MEDIUM PROPERTIES % ========================================================================= % define a large image size to move across number_scan_lines = 96; 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; background_map = background_map_mean + background_map_std * randn([Nx_tot, Ny_tot, Nz_tot]); sound_speed_map = c0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map; density_map = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map; % RUN THE SIMULATION % ========================================================================= % preallocate the storage scan_lines = zeros(number_scan_lines, kgrid.Nt); % set the input settings input_args = {... 'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size], ... 'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false}; % run the simulation if set to true, otherwise, load previous results from % disk h = waitbar(0,'Please wait...'); if RUN_SIMULATION % set medium position medium_position = 1; % loop through the scan lines for scan_line_index = 1:number_scan_lines scan_line_index waitbar(scan_line_index/number_scan_lines,h) % update the command line status disp(''); disp(['Computing scan line ' num2str(scan_line_index) ' of ' num2str(number_scan_lines)]); % load the current section of the medium medium.sound_speed = sound_speed_map(:, medium_position:medium_position + Ny - 1, :); medium.density = density_map(:, medium_position:medium_position + Ny - 1, :); % run the simulation sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:}); % extract the scan line from the sensor data scan_lines(scan_line_index, :) = transducer.scan_line(sensor_data); % update medium position medium_position = medium_position + transducer.element_width; end % save the scan lines to disk save HomogeneousTriaNew_scan_lines.mat; else % load the scan lines from disk load ('HomogeneousReferenceNew_scan_lines.mat'); end scan_lines_saved = scan_lines; scan_lines = scan_lines_saved; jgreen19 on "kspaceFirstOrder2D - function not recognised as valid" http://www.k-wave.org/forum/topic/kspacefirstorder2d-function-not-recognised-as-valid#post-9214 Sat, 10 May 2025 14:41:13 +0000 jgreen19 9214@http://www.k-wave.org/forum/ have managed to solve somehow, no need to respond! jgreen19 on "kspaceFirstOrder2D - function not recognised as valid" http://www.k-wave.org/forum/topic/kspacefirstorder2d-function-not-recognised-as-valid#post-9213 Fri, 09 May 2025 11:35:02 +0000 jgreen19 9213@http://www.k-wave.org/forum/ Hi, I've been able to run simple functions using k-Wave, then somehow without touching it, MATLAB is now failing to recognize: sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor); function. Even when using the code directly copy and pasted from chapter 3 of the user manual. I wonder if there is a solution to this? Perhaps I have installed the add-on incorrectly, but I feel I have tried all methods, but as a simple fact, no code (no matter how simple) is able to be run at the moment. Sun_hh on "Cuda code kspaceFirstOrder3DG not running on Windows 11" http://www.k-wave.org/forum/topic/cuda-code-kspacefirstorder3dg-not-running-on-windows-11#post-9212 Wed, 07 May 2025 17:39:53 +0000 Sun_hh 9212@http://www.k-wave.org/forum/ the same problem，need help ClaireH on "pstdElastic3D in layered fluid and elastic media" http://www.k-wave.org/forum/topic/pstdelastic3d-in-layered-fluid-and-elastic-media#post-9211 Mon, 28 Apr 2025 14:23:22 +0000 ClaireH 9211@http://www.k-wave.org/forum/ Hi, I'm a new k-Wave user interested in simulating ultrasound propagation through a layered medium that includes both fluid and elastic materials. The layers, in order of wave propagation, are: water → solid → soft tissue → solid → thin water layer → air I’m considering using the pstdElastic3D model for this simulation, and I’d like to know: Can pstdElastic3D handle a mixture of elastic and fluid media? Especially, can it accurately model the interface between water and air, given the large impedance mismatch and strong reflections at that boundary? Thanks so much for your help! 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.