k-Wave

1. POORNIMA
   

   Member
   Joined: Jul '19
   Posts: 6

   Hi,

   Is it possible to calculate the spatial average time average intensity(ISATA) from the simulation values. I have tried to simulate nonlinear ultrasound beam pattern using the example "simulating ultrasound beam pattern". As the input signal for the transducer is defined in units(pa), the matrix obtained from that is also pressure value. In my real time experiments, the input signal for the transducer is defined in unit (voltage).

   Hence I need to calculate the ISATA values of the simulated beam pattern. Using these simulated ISATA values, I would compare the values with the actual ISATA values obtained from measuring the intensity of the transducer using a hydrophone.

   Thank you very much for your help.

   This is my code

   % Simulating Ultrasound Beam Patterns Example
   %
   % This example shows how the nonlinear beam pattern from an ultrasound
   % transducer can be modelled. It builds on the Defining An Ultrasound
   % Transducer and Simulating Transducer Field Patterns examples.
   %
   % author: Bradley Treeby
   % date: 27th July 2011
   % last update: 25th September 2012
   %
   % This function is part of the k-Wave Toolbox (http://www.k-wave.org)
   % Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

   % This file is part of k-Wave. k-Wave is free software: you can
   % redistribute it and/or modify it under the terms of the GNU Lesser
   % General Public License as published by the Free Software Foundation,
   % either version 3 of the License, or (at your option) any later version.
   %
   % k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
   % WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
   % FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
   % more details.
   %
   % You should have received a copy of the GNU Lesser General Public License
   % along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

   depth = 50e-3; %50mm [m]

   % simulation settings
   DATA_CAST = 'single'; % set to 'single' or 'gpuArray-single' to speed up computations
   MASK_PLANE = 'yz'; % set to 'xy' or 'xz' to generate the beam pattern in different planes
   USE_STATISTICS = false; % set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns

   % =========================================================================
   % 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 = 10; % [grid points]

   % set total number of grid points not including the PML
   Nx = 128 - 2*PML_X_SIZE; % [grid points]
   Ny = 64 - 2*PML_Y_SIZE; % [grid points]
   Nz = 64 - 2*PML_Z_SIZE; % [grid points]

   % % set desired grid size in the x-direction not including the PML
   % x = 80e-3; % [m]

   % calculate the spacing between the grid points
   dx = 1e-3; % [m]
   dy = dx; % [m]
   dz = dx; % [m]

   % create the k-space grid
   kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

   % =========================================================================
   % DEFINE THE MEDIUM PARAMETERS
   % =========================================================================

   % define the properties of the propagation medium (skeletal muscle)
   medium.sound_speed = 1580; % [m/s]
   medium.density = 1050; % [kg/m^3]
   medium.alpha_coeff = 0.74; % [dB/(MHz^y cm)]
   medium.alpha_power = 1.5;
   medium.BonA = 6.6;

   % create the time array
   t_end = 100e-6; % [s]
   kgrid.t_array = makeTime(kgrid, medium.sound_speed, [], t_end);

   % =========================================================================
   % DEFINE THE INPUT SIGNAL
   % =========================================================================

   % define properties of the input signal
   source_strength = 1e6; % [Pa]
   tone_burst_freq = 1.25e6; % [Hz]

   % create the input signal using toneBurst
   input_signal = sin(2.*pi.*tone_burst_freq.*kgrid.t_array/2);

   % scale the source magnitude by the source_strength divided by the
   % impedance (the source is assigned to the particle velocity)
   input_signal = (source_strength./(medium.sound_speed*medium.density)).*input_signal;
   %%
   % =========================================================================
   % DEFINE THE ULTRASOUND TRANSDUCER
   % =========================================================================

   % physical properties of the transducer
   transducer.number_elements = 1; % total number of transducer elements
   transducer.element_width = 12; % width of each element [grid points]
   transducer.element_length = 12; % 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 = 1540; % sound speed [m/s]
   transducer.focus_distance = 0; % focus distance [m]
   transducer.elevation_focus_distance = 0;% focus distance in the elevation plane [m]
   transducer.steering_angle = 0; % steering angle [degrees]

   % apodization
   transducer.transmit_apodization = 'Rectangular';
   transducer.receive_apodization = 'Rectangular';

   % define the transducer elements that are currently active
   transducer.active_elements = 1;

   % append input signal used to drive the transducer
   transducer.input_signal = input_signal;

   % create the transducer using the defined settings
   transducer = makeTransducer(kgrid, transducer);

   % print out transducer properties
   transducer.properties;

   RADIUS = inf;
   DIAMETER = 22E-3 / dx;

   grid_size = [Nx; Ny; Nz];
   centre_pos = [2;round(Ny/2);round(Nz/2)];
   focus_pos = [round(Nx/2);round(Ny/2);round(Nz/2)];
   plot_bowl = true;

   source.p_mask = makeBowl(grid_size, centre_pos, RADIUS, DIAMETER, ...
   focus_pos, plot_bowl);

   source.p = input_signal;

   % =========================================================================
   % DEFINE SENSOR MASK
   % =========================================================================

   % define a sensor mask through the central plane
   sensor.mask = zeros(Nx, Ny, Nz);
   switch MASK_PLANE
   case 'xy'
   % define mask
   sensor.mask(:, :, Nz/2) = 1;

   % store y axis properties
   Nj = Ny;
   j_vec = kgrid.y_vec;
   j_label = 'y';

   case 'xz'
   % define mask
   sensor.mask(:, Ny/2, :) = 1;

   % store z axis properties
   Nj = Nz;
   j_vec = kgrid.z_vec;
   j_label = 'z';

   case 'yz'
   % define mask
   sensor.mask(round(depth/dx),:, :) = 1;

   % store z axis properties
   Nj = Nz;
   j_vec = kgrid.x_vec;
   j_label = 'x';

   end

   % set the record mode such that only the rms and peak values are stored
   if USE_STATISTICS
   sensor.record = {'p_rms', 'p_max'};
   end

   % =========================================================================
   % RUN THE SIMULATION
   % =========================================================================

   % set the input settings
   input_args = {'DisplayMask', source.p_mask, ...
   'PMLInside', false, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
   'DataCast', DATA_CAST, 'DataRecast', true, 'PlotScale', [-source_strength/2, source_strength/2]};

   % stream the data to disk in blocks of 100 if storing the complete time
   % history
   if ~USE_STATISTICS
   input_args = [input_args {'StreamToDisk', 100}];
   end

   % run the simulation
   sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

   % =========================================================================
   % COMPUTE THE BEAM PATTERN USING SIMULATION STATISTICS
   % =========================================================================

   if USE_STATISTICS

   % reshape the returned rms and max fields to their original position
   sensor_data.p_rms = reshape(sensor_data.p_rms, [Nx, Nj]);
   sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Nj]);

   % plot the beam pattern using the pressure maximum
   figure;
   imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, sensor_data.p_max/1e6);
   xlabel([j_label '-position [mm]']);
   ylabel('x-position [mm]');
   title('Total Beam Pattern Using Maximum Of Recorded Pressure');
   colormap(jet(256));
   c = colorbar;
   ylabel(c, 'Pressure [MPa]');
   axis image;

   % plot the beam pattern using the pressure rms
   figure;
   imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, sensor_data.p_rms/1e6);
   xlabel([j_label '-position [mm]']);
   ylabel('x-position [mm]');
   title('Total Beam Pattern Using RMS Of Recorded Pressure');
   colormap(jet(256));
   c = colorbar;
   ylabel(c, 'Pressure [MPa]');
   axis image;

   % end the example
   return

   end

   voxelPlot(source.p_mask+sensor.mask)
   % =========================================================================
   % COMPUTE THE BEAM PATTERN FROM THE AMPLITUDE SPECTRUM
   % =========================================================================
   %%
   % reshape the sensor data to its original position so that it can be
   % indexed as sensor_data(x, j, t)
   sensor_data = reshape(sensor_data, [Ny, Nz, kgrid.Nt]);

   % compute the amplitude spectrum
   [freq, amp_spect] = spect(sensor_data, 1/kgrid.dt, 'Dim', 3);

   % compute the index at which the source frequency and its harmonics occur
   [f1_value, f1_index] = findClosest(freq, tone_burst_freq);
   [f2_value, f2_index] = findClosest(freq, tone_burst_freq*2);

   % extract the amplitude at the source frequency and store
   beam_pattern_f1 = amp_spect(:, :, f1_index);

   % extract the amplitude at the second harmonic and store
   beam_pattern_f2 = amp_spect(:, :, f2_index);

   % extract the integral of the total amplitude spectrum
   beam_pattern_total = sum(amp_spect, 3);

   % plot the beam patterns
   figure;
   imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, beam_pattern_f1/1e6);
   xlabel([j_label '-position [mm]']);
   ylabel('x-position [mm]');
   title('Beam Pattern At Source Fundamental');
   colormap(jet(256));
   c = colorbar;
   ylabel(c, 'Pressure [MPa]');
   axis image;

   figure;
   imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, beam_pattern_f2/1e3);
   xlabel([j_label '-position [mm]']);
   ylabel('x-position [mm]');
   title('Beam Pattern At Second Harmonic');
   colormap(jet(256));
   c = colorbar;
   ylabel(c, 'Pressure [kPa]');
   axis image;

   figure;
   imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, beam_pattern_total/1e6);
   xlabel([j_label '-position [mm]']);
   ylabel('x-position [mm]');
   title('Total Beam Pattern Using Integral Of Recorded Pressure');
   colormap(jet(256));
   c = colorbar;
   ylabel(c, 'Pressure [MPa]');
   axis image;
   hold on

   % =========================================================================
   % PLOT DIRECTIVITY PATTERN AT FOCUS
   % =========================================================================

   % compute the directivity at each of the harmonics
   directivity_f1 = squeeze(beam_pattern_f1(round(transducer.focus_distance/dx), :));
   directivity_f2 = squeeze(beam_pattern_f2(round(transducer.focus_distance/dx), :));

   % normalise
   directivity_f1 = directivity_f1./max(directivity_f1(:));
   directivity_f2 = directivity_f2./max(directivity_f2(:));

   % compute relative angles from transducer
   if strcmp(MASK_PLANE, 'xy')
   horz_axis = ((1:Ny) - Ny/2)*dy;
   else
   horz_axis = ((1:Nz) - Nz/2)*dz;
   end
   angles = 180*atan2(horz_axis, transducer.focus_distance)/pi;

   % plot the directivity
   figure;
   plot(angles, directivity_f1, 'k-', angles, directivity_f2, 'k--');
   axis tight;
   set(gca, 'FontSize', 12);
   xlabel('Angle [deg]');
   ylabel('Normalised Amplitude');
   legend('Fundamental', 'Second Harmonic', 'Location', 'NorthWest');

   Posted 4 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi POORNIMA,

   In general, you can average the time-average intensity over the beam area (e.g., as defined in BS EN 62127-1:2007 +A1:2013).

   Why not use the same approach you have used to calculate Isata in the experiment?

   Brad.

   Posted 4 years ago #

3. POORNIMA
   

   Member
   Joined: Jul '19
   Posts: 6

   Hi Brad,

   Thank you for the reply. But in the simulation the values are expressed in pressure. In actual experiments, the Isata values are calculated in terms of voltage.

   I would like to compare the simulated intensity with the experimental intensity. How can I convert the simulated pressure values to voltage.

   Thanks once again

   Posted 4 years ago #

4. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi POORNIMA,

   You can't calculate acoustic intensity directly from voltage, you will need to convert your measured voltages to pressure. To convert between voltage and pressure for your experiment, you will need to use the hydrophone calibration.

   Hope that helps,

   Brad.

   Posted 4 years ago #

RSS feed for this topic

Reply

You must log in to post.