k-Wave

1. kcox
   

   Member
   Joined: May '17
   Posts: 10

   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)

   Posted 6 years ago #

2. kcox
   

   Member
   Joined: May '17
   Posts: 10

   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)

   Posted 6 years ago #

3. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   The 'SignalOffset' input for the toneBurst function is given in units of time samples, i.e., actual time divided by dt. There should be details in the help for toneBurst. 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).

   Posted 6 years ago #

4. robou
   

   Member
   Joined: Sep '20
   Posts: 1

   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

   Posted 3 years ago #

RSS feed for this topic

Reply

You must log in to post.