k-Wave

1. ghost82
   

   Member
   Joined: Jan '12
   Posts: 4

   Hi,
   my name is Daniele. I found and downloaded k-wave just yestarday and the first thing I can say is that this is a great toolbox :)
   Now, if you can provide me some more information to solve my doubts I will be very greatful.
   Premise: I'm quite new to matlab..

   I want to simulate an ultrasound transducer, 1 element, in water medium; I studied the US example (example_us_beam_patterns) and I have a few questions:

   1) My input signal for the transducer, as shown in scientific literature is: P0-Pa*sin(2*pi*f*t), where P0 is the atmospheric pressure, Pa is the Pressure amplitude, f the frequence, t the time. In the example I think that I have understand that the signal is generated by toneburst, and the signal amplitude always range from -1 to 1; then this pressure signal is converted in velocity units;
   What is the source_strength? Is it my "Pa" "Pressure amplitude"?
   I think there is also an error in the comments as this "source_strength" maybe expressed in Pascal, not MPa.
   Is it for me, the following signal source input ok?

   input_signal = 0-sin(2*pi*freq*kgrid.t_array);
   input_signal = (source_strength./(medium.sound_speed*medium.density)).*input_signal;

   where freq is the ultrasonic frequency (25 kHz) and source_strength is my Pressure amplitude of 240000 (Pa)

   2) I don't understand why there are these lines in the script:
   % 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]

   So I commented them, modified the total number of grid points without subtracting PML grid points, and modified the input_args.

   Can you explain me the meaning of defining a PML?I thought to create only a domain with PML boundaries applied to the edges of my domain.

   3) When the total beam pattern charts are displayed at the end of the simulation I read on the coloured bar (MPa); however I think that the correct unit is MPa gauge: Am I right?So when I read 0 absolute pressure is atmospheric pressure in my case.

   Thank you very much for reading,
   and I still write that this is a very good toolbox.
   Thank you,

   Daniele

   Posted 12 years ago #

2. ghost82
   

   Member
   Joined: Jan '12
   Posts: 4

   Good morning all!

   ok, I think I understand a bit more...

   1)My input signal should be:
   source_strength = 113000;
   freq = 25000;
   input_signal = source_strength*sin(2*pi*freq*kgrid.t_array);
   input_signal = (1./(medium.sound_speed*medium.density)).*input_signal;

   where 113000 is my pressure amplitude calculated from (2*intensity*sounda velocity*rho water)^0.5 with the intensity equal to power/area of irradiation of the transducer.
   P0, my ambient pressure, hasn't to be included in the transducer signal, as all is referenced to 0 Pascal gauge.

   2)I think I understand also the PML lines: they generate a border around my domain as perfectly matched layer.

   3)Units in the final charts should be Mpag.

   Am I right?

   Now I have another doubt:
   I run a simulation with this code:

   clear all;

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

   % =========================================================================
   % 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]

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

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

   % calculate the spacing between the grid points
   dx = 1e-2; % [m]
   dy = 1e-2; % [m]
   dz = 1e-2; % [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
   medium.sound_speed = 1450; % [m/s]
   medium.density = 998.2; % [kg/m^3]
   medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
   medium.alpha_power = 1.5;
   medium.BonA = 6;

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

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

   % define properties of the input signal
   source_strength = 113000; % [Pa]
   freq = 25000; % [Hz]

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

   % scale the source magnitude by the source_strength divided by the
   % impedance (the source is assigned to the particle velocity)
   input_signal = (1./(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 = 4; % width of each element [grid points]
   transducer.element_length = 4; % 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 = 6000; % sound speed [m/s]
   transducer.focus_distance = 1e-3; % focus distance [m]
   transducer.elevation_focus_distance = 1e-3;% 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 = 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 = makeTransducer(kgrid, transducer);

   % print out transducer properties
   transducer.properties;

   % =========================================================================
   % 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';

   end

   % set the record mode such that only the rms and peak values are stored
   if USE_STATISTICS
   sensor.record_mode = 'statistics';
   end

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

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

   % input_args = {'DisplayMask', transducer.all_elements_mask, ...
   % 'PMLInside', false, 'PlotPML', false, ...
   % 'DataCast', DATA_CAST, '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, transducer, 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 [MPag]');
   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 [MPag]');
   axis image;

   % end the example
   return
   end

   % =========================================================================
   % 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, [Nx, Nj, length(kgrid.t_array)]);

   % compute the frequency axis
   freq = (0:ceil((kgrid.Nt + 1)/2) - 1) ./ (kgrid.dt*length(kgrid.t_array));

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

   % preallocate the beam pattern variables
   beam_pattern_f0 = zeros(Nx, Nj);
   beam_pattern_f1 = zeros(Nx, Nj);
   beam_pattern_total = zeros(Nx, Nj);

   % compute the amplitude spectrum of the time series recorded at each sensor
   % point, and then extract the corresponding amplitudes at the fundamental
   % frequency and second harmonic.
   for x_index = 1:Nx
   for j_index = 1:Nj

   % compute the amplitude spectrum
   amp_spect = spectrum(squeeze(sensor_data(x_index, j_index, :)), 1/kgrid.dt); %, 'Window', 'Hanning');

   % extract the amplitude at the source frequency and store
   beam_pattern_f0(x_index, j_index) = amp_spect(f0_index);

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

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

   end
   end

   figure;
   imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, beam_pattern_f0/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_f1/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;

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

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

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

   % 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_f0, 'k-', angles, directivity_f1, 'k--');
   axis tight;
   set(gca, 'FontSize', 12);
   xlabel('Angle [deg]');
   ylabel('Normalised Amplitude');
   legend('Fundamental', 'Second Harmonic', 'Location', 'NorthWest');

   _____________________________________________________________________

   However, after some time, in the animation charts, waves from the transducer disappear; but my signal is a continuos sinusoidal, so I don't understand why they stop to exit from transducer..

   Here is a picture:
   http://imageshack.us/photo/my-images/718/48364580.png

   Can you explain me please?

   Thank you,

   Daniele

   Posted 12 years ago #

3. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Daniele,

   Thanks for your feedback. To answer your previous questions:

   1.) Yes, the pressure is specified in units of Pascals, the [MPa] in the example comments is a typo.

   2.) When the optional input 'PMLInside' is set to false, the grid specified by the user is expanded by the values set for 'PMLSize' This region is used as an absorbing boundary layer. To force the total size of the computational grid to be a power of 2 (the FFT is fastest when the grid size is a power of two or has small prime factors), the PML size is subtracted from the desired size. This step is not required for the simulation to work (you can use any grid size).

   3.) The pressure is always acoustic (or excess) pressure. This is the dynamic perturbation to the ambient pressure in the medium caused by the acoustic wave.

   Regarding your simulation, I am also puzzled why your source stops transmitting. I can't see anything immediately wrong with your code. Leave it with me, and I'll get back to you. In the meantime, you could directly assign source.u_mask to represent your transducer, and then assign your signal to source.ux (i.e., try the same thing not using the transducer class).

   Brad.

   Posted 12 years ago #

4. ghost82
   

   Member
   Joined: Jan '12
   Posts: 4

   Thank you Bradley for yuor reply, I will try what you suggest and I will get back.

   Daniele

   Posted 12 years ago #

5. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Daniele,

   The problem is due to a bug in setting the data type for the delay mask used to apply the input signal to the grid points that form part of the ultrasound transducer. The data type is incorrectly based only on the delays, and not the delays + the length of the source.

   If you open kspaceFirstOrder_inputChecking.m in the private folder, and replace lines 441-448 with:

   % convert the data type depending on the maximum value of the delay
   % mask and the length of the source
   max_delay = max(delay_mask(:)) + length(transducer_input_signal) + 1;
   if max_delay < intmax('uint8');
       delay_mask = uint8(delay_mask);
   elseif max_delay < intmax('uint16');
       delay_mask = uint16(delay_mask);
   elseif max_delay < intmax('uint32');
       delay_mask = uint32(delay_mask);
   end

   That should fix the problem. Thanks for posting about the strange behaviour!

   Brad.

   Posted 12 years ago #

6. ghost82
   

   Member
   Joined: Jan '12
   Posts: 4

   Thank you Bradley, I think now it is fixed in the new release, I will try it.
   However I followed your suggestion and I was able to simulate with source.ux.

   Thank you,

   Daniele

   Posted 12 years ago #

7. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Daniele,

   Indeed, this issue should be fixed in the new release (B.0.5). If it's not, please let us know. Thanks for your feedback and for reporting this bug. It's very much appreciated!

   Brad.

   Posted 12 years ago #

RSS feed for this topic

Reply

You must log in to post.