k-Wave Toolbox Previous   Next

Homogeneous Propagation Medium Example

Overview

This example provides a simple demonstration of using k-Wave for the simulation and detection of the pressure field generated by an initial pressure distribution within a two-dimensional homogeneous propagation medium.

 Back to Top

Creating the k-space grid

The medium discretisation is performed by makeGrid. Both the size (dx, dz) and number (Nx, Nz) of pixels in each Cartesian direction are used to compute the discretisation, and a k-Wave grid structure is returned. As the k-space method is based heavily on the FFT, the simulations will be fastest for grids where the number of pixels is given by a power of two (see fft2).

% create the computational grid
Nx = 128;           % number of pixels in the x (column) direction
Nz = 128;           % number of pixels in the z (row) direction
dx = 0.1e-3;        % pixel width [m]
dz = 0.1e-3;        % pixel height [m]
kgrid = makeGrid(Nx, dx, Nz, dz);

The structure kgrid contains numerous fields, including matrices of the computational wavenumbers (kgrid.kx, kgrid.kz, and kgrid.k) and positions (kgrid.x, kgrid.z) used by many k-Wave functions (for a full list, see the help file for makeGrid).

 Back to Top

Avoiding the perfectly matched layer

When the acoustic waves reach the edge of the computational domain, they are absorbed by a special type of anisotropic absorbing boundary layer known as a perfectly matched layer. The effects of the layer can be seen by watching what happens to propagating waves as they get close to the edge of the computational domain. By default, this layer occupies a 20 pixel strip (10 voxels in 3D) around the edge of the domain inside the computational domain specified using makeGrid. To avoid strange effects, care must be taken not to place the source or sensor positions inside this layer. Alternatively, the perfectly matched layer can be set to be outside the computational domain set by the user. See Controlling The Absorbing Boundary Layer Example for more detailed instructions on how to modify the properties and position of the perfectly matched layer.

 Back to Top

Defining the medium properties

For a homogeneous medium, the sound speed is set as a scalar in SI units. Power law acoustic absorption can also be optionally set by assigning values to medium.alpha_power and medium.alpha_coeff. It is useful to note that for a homogeneous medium, the computation of the acoustic pressure is not dependent on the ambient density. Consequently, if no velocity inputs or outputs are required, medium.density does not need to be specified.

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]
medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;

The grid spacing and the sound speed govern the maximum frequency that the simulation grid is able to propagate. This frequency (for each Cartesian direction) can be calculated by f_max_x = medium.sound_speed/(2*dx).

 Back to Top

Defining the initial pressure distribution

The initial pressure distribution is set as a matrix of arbitrary numeric values over the medium discretisation given by the k-space grid. Several functions are included in the toolbox for the creation of simple geometric shapes. Here makeDisc is used to create an initial pressure distribution of two small discs with different diameters. This is assigned to the p0 field of the source structure. There are no restrictions on source.p0 except that it must be the same size as the computational grid, i.e., it must have Nx columns and Nz rows.

% create initial pressure distribution using makeDisc
disc_magnitude = 5;
disc_x_pos = 50;    % pixels
disc_z_pos = 50;    % pixels
disc_radius = 8;    % pixels
disc_1 = disc_magnitude*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius);

disc_magnitude = 3;
disc_x_pos = 60;    % pixels
disc_z_pos = 80;    % pixels
disc_radius = 5;    % pixels
disc_2 = disc_magnitude*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius);

source.p0 = disc_1 + disc_2;

 Back to Top

Defining the sensor mask

The sensor mask defines the locations where the pressure field is recorded at each time-step. The sensor mask can be given either as a binary grid of sensor positions (see the Using A Binary Sensor Mask Example), or as a set of arbitrary Cartesian coordinates (given in 2D as a 2 x N matrix) within the dimensions of the k-space grid (the Cartesian origin is assumed to be in the center). Here makeCartCircle is used to set the sensor mask to the Cartesian coordinates of a set of evenly spaced points on a circle.

% define a centered circular sensor
sensor_radius = 4e-3;   % [m]
num_sensor_points = 50;
sensor.mask = makeCartCircle(sensor_radius, num_sensor_points);

A plot of the initial pressure distribution and the sensor mask (using nearest neighbour interpolation) is given below. By default, the pressure at the Cartesian points in 2D is computed using linear interpolation.

 Back to Top

Running the simulation

The computation is invoked by calling kspaceFirstOrder2D with the four input structures defined above. By default, a visualisation of the propagating wave-field and a status bar are displayed, with frame updates every ten time-steps. The default k-Wave color map displays positive pressures as yellows to reds to black, and negative pressures as light to dark blue-greys (see getColorMap).

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

As the function runs, status updates and computational parameters are printed to the command line.

Running k-space simulation...
  start time: 09-Feb-2011 11:37:42
  dt: 20ns, t_end: 12.06us, time steps: 604
  input grid size: 128 by 128 pixels (12.8 by 12.8mm)
  maximum supported frequency: 7.5MHz
  smoothing p0 distribution...
  calculating Delaunay triangulation (TriScatteredInterp)...
  precomputation completed in 0.58273s
  starting time loop...
  estimated simulation time 2.626s...
  simulation completed in 4.2276s
  total computation time 4.84s

When the time loop has completed, the function returns the recorded time series at each of sensor locations defined by sensor.mask. This is indexed as sensor_data(sensor position, time step). For a Cartesian sensor mask, the time series are returned in the same order as the Cartesian coordinates specified in sensor.mask. A visualisation of the sensor data recorded in this example is given below.

% plot the simulated sensor data
figure;
imagesc(sensor_data, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;

Running the same example with num_sensor_points = 500; records the wave-field with a much higher spatial resolution.

 Back to Top


© 2009, 2010, 2011 Bradley Treeby and Ben Cox.