k-Wave User Forum » Topic: Creating a waveform with custom time offsets

k-Wave User Forum » Topic: Creating a waveform with custom time offsets http://www.k-wave.org/forum/topic/creating-a-waveform-with-custom-time-offsets Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:15:42 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php robou on "Creating a waveform with custom time offsets" http://www.k-wave.org/forum/topic/creating-a-waveform-with-custom-time-offsets#post-8121 Fri, 16 Apr 2021 00:46:27 +0000 robou 8121@http://www.k-wave.org/forum/ Hello, I am... blocked... is the term. I am sure that it´s a simple error that i do not catch. I trying to create a multiring trandsucer, It´s seems that geomotry is correct.Then I create signals and calculate needed delays for each rings. I understand that if there are 4 rings, we must have 4 differents toneburst vectors, isn't it? When I do not use offset everything is working, but that's not the case for multirings. Rings value index go from 1 to 4. Code is more explicit, or not:). Thank you Romain %% Creation of transducer geometry source.p_mask = zeros(Nz,Ny,Nx); % source mask % if plane transducer if(type == 0) radius = round(aperture*1e-3/2/ds); % radius of aperture [grid points] cy = round(Ny/2); % Y-axis center of the disk cx = round(Nx/2); % X-axis center of the disk geom = makeDisc(Nx,Ny,cx,cy,radius); % disk geometry % add the geometry to the source domain source.p_mask(Nz_trans_pos+PML_size, 1:Ny, 1:Nx) = geom; % if curved transducer elseif(type == 1) Nfd = round(fd/ds); Nap = round(aperture*1e-3/ds); % Redefine aperture to make it an odd number for makeBowl to work if rem(Nap,2) == 0 Nap = Nap + 1; end trans_pos = [PML_size+Nz_trans_pos, round(Ny/2), round(Nx/2)]; focus_pos = [PML_size+Nz_trans_pos+Nfd, round(Ny/2), round(Nx/2)] ; % makeBowl(grid_size, origin bowl_pos, radius curvature, aperture diameter, location of focus_pos, 'Plot', true,'RemoveOverlap',true); geom = makeBowl(size(source.p_mask), trans_pos, Nfd+1, Nap, focus_pos, 'Plot', false,'RemoveOverlap',true); %----------------------- Version Romain ------------------------% % Multirings configuration with identification % 2D rings creation NradioInnerDisc = round((1e-3 * ri)/ds - ds); % in grid number NradioOuterDisc = round((1e-3 * ro)/ds); innerDisc = zeros (Ny,Nx); % in grid number outerDisc = zeros (Ny,Nx); multiRings2D = zeros (Ny,Nx); for i=1:length(ri) innerDisc = makeDisc(Ny,Nx,round(Ny/2),round(Nx/2),NradioInnerDisc(i)); outerDisc = makeDisc(Ny,Nx,round(Ny/2),round(Nx/2),NradioOuterDisc(i)); ring2D = outerDisc - innerDisc; multiRings2D = multiRings2D + i.* ring2D; % index identification grow with radius figure; imagesc(multiRings2D); axis image; end % 3D Projection multiRings3D = repmat(multiRings2D,1,1,Nz); multiRings3D = permute(multiRings3D,[3,2,1]); % Apply rings to 3D bowl geomMultiRings3D = geom.*multiRings3D; % Add the geometry to the source domain source.p_mask = geomMultiRings3D; % Verification transductor = squeeze(sum(geomMultiRings3D,1)); figure; imagesc(squeeze(transductor)); axis image; %flyThrough(geomMultiRings3D(1:50,:,:),1,50,1); % figure % for i=1:floor(Ny/2) % subplot(3,1,1) % imagesc(squeeze(geomMultiRings3D(i,:,:)));axis image % subplot(3,1,2) % imagesc(squeeze(geomMultiRings3D(:,round(Ny/2)+i,:)));axis image % subplot(3,1,3) % imagesc(squeeze(geomMultiRings3D(:,:,round(Nx/2)+i)));axis image % drawnow; pause; % end %% Create time array and emitted signal signal % Function 'makeTime' automatically specifies Nt, dt, and t_array based on the % Courant-Friedrichs-Lewy (CFL) number and the grid size. % Default CFL=0.3, where dt = cfl * dx / sound_speed. c_max = max(medium.sound_speed(:)); % maximum sound speed needed for "makeTime" function [m/s] c_min = min(medium.sound_speed(:)); % minimum sound speed needed for calculation of "t_end" [m/s] cfl = 0.3; % Courant-Friedrich-Lewy stability number tone_burst_cycles = 10; % number of emitted cycles working at "source_freq" t_end = 3*aperture*1e-3/c_min+tone_burst_cycles/source_freq; % total time of wave propagation [s] kgrid.t_array = makeTime(kgrid, c_max, cfl, t_end); % calculate time steps input_signal = p0 * toneBurst(1/kgrid.dt, source_freq, tone_burst_cycles,'Envelope',[1 1],'SignalLength',150); % obtain the emitted burst signal waveform [Pa] %% Beamforming: Calculation of delays %-------------------- Version Romain -------------------------- % Set the configuration if no beanforming used if strcmp(beamforming, 'no') %Set same identification to all rings source.p_mask( source.p_mask > 0 ) = 1; %Set same signal for all transducer elements source.p = input_signal; else % Create CenterRings for delays calculation % rc is the radius that cuts a ring into two subrings of equal surface area % rc, distance from the center of the rings to Z axis knowed from datasheet CenterLines = zeros(Nz,Ny); rc = [18.6; 22.98; 26.59; 29.67]; radioCenterCircle = (1e-3 *rc)/ds - ds; % in grid number for i=1:length(apertureIn) CenterLines = CenterLines + ... makeLine( Nz,Ny,... [1 , round(Ny/2 + radioCenterCircle(i))],... [Nz, round(Ny/2 + radioCenterCircle(i))]... ); end figure;imagesc(CenterLines); axis image; % Delays calculation % Create a ZY axial plane of geom geom_ZYplane = geom(:,:,round(Nx/2)); % Find index where CenterLines intercept geometry source ZY plane [intercept_z,intercept_y] = find((geom_ZYplane==CenterLines).*(geom_ZYplane>0)); intercept_x = intercept_y NcenterRings = cat(2,intercept_z,intercept_y,intercept_x) % Calculate propagation times and delays for each rings Nzf = round(zf/ds) % Desired Z axis focus distance Nfocal = [PML_size+Nz_trans_pos+Nzf, round(Ny/2), round(Nx/2)] ; Nfocal = repmat(Nfocal,length(ri),1) NcenterRings2focus = Nfocal - NcenterRings centerRings2focus = abs(NcenterRings2focus) * ds % convert to distancee centerRings2focus_distance = sqrt(centerRings2focus(:,1).^2 + centerRings2focus(:,2).^2) centerRings2focus_propaTime = centerRings2focus_distance ./ c0(1) centerRings2focus_delays = max(centerRings2focus_propaTime) - centerRings2focus_propaTime % Convert delays in (s) to number of samples tone_burst_offset = round(centerRings2focus_delays ./ kgrid.dt); % transpose to obtain a row vector % Create signals input_signals = p0*toneBurst(1/kgrid.dt, source_freq, tone_burst_cycles,... 'Envelope',[1 , 1],'SignalLength',150, 'SignalOffset', tone_burst_offset); % Plot signals figure; t_axis = (1:size(input_signals, 2)) * kgrid.dt; stackedPlot(1e6 * t_axis, input_signals); xlabel('Time [\mus]'); ylabel('Signal Number'); % Set signals for all transducer elements source.p = input_signals; unning k-Wave simulation... start time: 16-Apr-2021 01:13:06 reference sound speed: 1482.4595m/s dt: 86.3636ns, t_end: 134.6409us, time steps: 1560 input grid size: 224 by 225 by 225 grid points (95.5962 by 96.0229 by 96.0229mm) maximum supported frequency: 1.7368MHz by 1.7291MHz by 1.7291MHz precomputation completed in 0.62563s saving input files to disk... completed in 0.62381s +---------------------------------------------------------------+ | kspaceFirstOrder-OMP v1.3 | +---------------------------------------------------------------+ | Reading simulation configuration: Done | | Number of CPU threads: 40 | | Processor name: Intel(R) Xeon(R) CPU E5-2680 v2 @ 2.80GHz | +---------------------------------------------------------------+ | Simulation details | +---------------------------------------------------------------+ | Domain dimensions: 224 x 225 x 225 | | Medium type: 3D | | Simulation time steps: 1560 | +---------------------------------------------------------------+ | Initialization | +---------------------------------------------------------------+ | Memory allocation: Done | | Data loading: Failed | +---------------------------------------------------------------+ +---------------------------------------------------------------+ | !!! K-Wave experienced a fatal error !!! | +---------------------------------------------------------------+ | Error: Dataset p_source_input has wrong dimension sizes.: | | iostream stream error | +---------------------------------------------------------------+ | Execution terminated | +---------------------------------------------------------------+ Error using h5readc Unable to open file. Filename may be corrupt or have unsupported characters. Error in h5read (line 58) [data,var_class] = h5readc(Filename,Dataset,start,count,stride); Error in kspaceFirstOrder3DC (line 569) Nx = h5read(output_filename, '/Nx'); Error in eye_acoustic_sim_Ta_multiRings (line 552) sensor_data = kspaceFirstOrder3DC(kgrid, medium, source, sensor, input_args{:}, 'DataPath',tempPath); % simulation Bradley Treeby on "Creating a waveform with custom time offsets" http://www.k-wave.org/forum/topic/creating-a-waveform-with-custom-time-offsets#post-5965 Mon, 29 May 2017 08:58:55 +0000 Bradley Treeby 5965@http://www.k-wave.org/forum/ The <code>'SignalOffset'</code> input for the <code>toneBurst</code> function is given in units of time samples, i.e., actual time divided by dt. There should be details in the help for <code>toneBurst</code>. The 40 is used to add a constant offset to prevent negative values. This number just delays the source (try changing it to 20 and running the example). kcox on "Creating a waveform with custom time offsets" http://www.k-wave.org/forum/topic/creating-a-waveform-with-custom-time-offsets#post-5962 Fri, 26 May 2017 18:24:05 +0000 kcox 5962@http://www.k-wave.org/forum/ Also, where does the '40' in the following equation for the tone_burst_offset come from? tone_burst_offset = 40 + element_spacing*element_index*sin(steering_angle*pi/180)/(medium.sound_speed*dt) kcox on "Creating a waveform with custom time offsets" http://www.k-wave.org/forum/topic/creating-a-waveform-with-custom-time-offsets#post-5961 Fri, 26 May 2017 18:14:14 +0000 kcox 5961@http://www.k-wave.org/forum/ Hey there, I modified the steering the linear array example to have all elements fire at the appropriate times in order to meet at some desired focus. However, to get it to work, I had to apply a mock signal for a set amount of time, as I was not using the toneBurst function. I was wondering what units the tone_burst_offset is, so that I can't apply a similar concept to my function and get waveforms representative of actual acoustic sources. My code is shown below: % ========================================================================= % SIMULATION % ========================================================================= clear all close all % create the computational grid Nx = 216; % number of grid points in the x (row) direction Ny = Nx; % number of grid points in the y (column) direction dx = 50e-3/Nx; % grid point spacing in the x direction [m] dy = dx; % grid point spacing in the y direction [m] kgrid = makeGrid(Nx, dx, Ny, dy); % define the properties of the propagation medium medium.sound_speed = 1500; % [m/s] medium.alpha_power = 1.5; % [dB/(MHz^y cm)] medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)] medium.density = 1; % create the time array [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed); kgrid.t_array = 0:dt:1000*dt; % define sensor for entire grid sensor.mask = zeros(Nx, Ny); sensor.mask = ones(size(sensor.mask)); % define source mask for a linear transducer num_elements = 216; % [grid points] x_offset = 10; % [grid points] source.p_mask = zeros(Nx, Ny); start_index = Ny/2 - round(num_elements/2) + 1; source.p_mask(x_offset, start_index:start_index + num_elements - 1) = 1; % define the properties of the tone burst used to drive the transducer sampling_freq = 1/dt; % [Hz] steering_angle = 0; % [deg] element_spacing = dx; % [m] tone_burst_freq = 12e6; % [Hz] tone_burst_cycles = 10; % finds x and y positions of the elements [x_pos,y_pos] = ind2sub([Nx,Ny],find(source.p_mask == 1)); % define x and y position of the desired focus location x_focus = 120; y_focus = Nx/2; % calculate distances between elements and focus location dists = sqrt((x_focus - x_pos).^2 + (y_focus - y_pos).^2); % calculate number of actual steps steps = dists.*kgrid.dx; % calculate the times times = steps./medium.sound_speed; % calculate the actual time steps timesteps = round(times./kgrid.dt); % calculate time delays istarts = round(max(timesteps)+1-timesteps); % signal = toneBurstCopy(sampling_freq, tone_burst_freq, tone_burst_cycles); for ii = 1:num_elements source.p(ii,istarts(ii):istarts(ii)+3) = 5*sin(2*pi*tone_burst_freq*dt); end input_args = {'DisplayMask', source.p_mask}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % ========================================================================= % VISUALISATION % ========================================================================= % plot the input time series figure; num_source_time_points = length(source.p(1,:)); [t_sc, scale, prefix] = scaleSI(kgrid.t_array(num_source_time_points)); stackedPlot(kgrid.t_array(1:num_source_time_points)*scale, source.p); xlabel(['Time [' prefix 's]']); ylabel('Input Signals'); % reshaping of sensor_data to show acoustic field and focus point size_sensor_data = size(sensor_data); dimx = sqrt(size_sensor_data(1)); dimy = dimx; dimz = size_sensor_data(2); tmp = reshape(sensor_data,dimx,dimy,dimz); tmp2 = max(tmp,[],3); figure imagesc(tmp2)