k-Wave Toolbox |
![]() ![]() |
On this page… |
---|
Introduction to the k-Wave Toolbox |
k-Wave is a freely available third-party toolbox for MATLAB (released under the GNU LGPLv3 license) developed by Bradley Treeby and Ben Cox. The toolbox is designed to make the time-domain simulation of acoustic wave propagation fast and easy. Some of the main features include:
The MATLAB functions included within the toolbox can be divided into four basic categories:
The simulation functions compute the time evolution of an acoustic wavefield in a homogeneous or heterogeneous medium in either one, two, or three dimensions. The computations are based on a k-space psuedospectral solution to coupled first-order acoustic equations defined on a regularly spaced Cartesian grid. The photoacoustic image reconstruction functions compute an estimate of initial photoacoustic pressure given a set of time-varying boundary measurements. The two algorithms included in the toolbox are (a) time reversal, which can accept data recorded from an arbitrarily-shaped measurement surface, and (b) FFT-based reconstruction, which can estimate the initial photoacoustic pressure from data recorded over a linear (2D) or planar (3D) measurement surface. The geometry creation functions allow many useful objects to be defined easily, including circles, arcs, discs, spheres, shells, and balls. They can return the shapes either fitted to a known spatial grid (i.e. as a binary matrix of 1's and 0's, with the 1's corresponding to the location of shape) or in exact geometric coordinates. The utility functions are used to perform additional tasks such as grid creation, matrix smoothing, matrix interpolation, file loading, etc.
The computational architecture of the k-Wave simulation functions is illustrated
below (see kspaceFirstOrder1D
,
kspaceFirstOrder2D
, and
kspaceFirstOrder3D
for further details).
The functions are given information about the discretisation of the propagation medium (the k-Wave grid), its acoustic properties (sound speed, absorption, etc), the initial (or time varying) pressure or particle velocity distribution, and the location and characteristics of the measurement surface which detects the acoustic wave-field. These properties are assigned as fields within four input structures:
kgrid
medium
,source
, andsensor
.The propagation of the wave-field through the medium is then computed step-by-step, with the pressure (and velocity) values over the sensor stored after each iteration. These values are returned when the time loop has completed.
A simple example of a 2D initial value problem in a heterogeneous medium is given below (1D and 3D simulations are performed in an analogous fashion).
% create the computational grid Nx = 128; % number of grid points in the x (row) direction Ny = 256; % number of grid points in the y (column) direction dx = 50e-6; % grid point spacing in the x direction [m] dy = 50e-6; % grid point spacing in the y direction [m] kgrid = makeGrid(Nx, dx, Ny, dy); % define the medium properties medium.sound_speed = 1500*ones(Nx, Ny); % [m/s] medium.sound_speed(1:50, :) = 1600; % [m/s] medium.density = 1000*ones(Nx, Ny); % [kg/m^3] medium.density(1:50, :) = 1040; % [kg/m^3] % define the initial pressure using makeDisc disc_x_pos = 75; % [grid points] disc_y_pos = 120; % [grid points] disc_radius = 8; % [grid points] disc_mag = 3; % [au] source.p0 = disc_mag*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius); % define a Cartesian sensor mask of a centered circle with 50 sensor elements sensor_radius = 2.5e-3; % [m] num_sensor_points = 50; sensor.mask = makeCartCircle(sensor_radius, num_sensor_points); % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);
The medium discretisation is created using the utility function
makeGrid
. The size (dx, dy) and
number (Nx, Ny) of grid points in each Cartesian direction are used to calculate
the Cartesian and k-space discretisations and an object of the kWaveGrid class (called kgrid
in this example) is returned. The discretisations are calculated to satisfy the requirements of the FFT-based spatial derivatives. This object is used extensively by both the simulation and utility functions within k-Wave. The time steps used in the simulation are defined by kgrid.t_array
which is set to 'auto'
by default. In this case, the time array is automatically calculated within the simulation functions using the utility function makeTime
based on the size and properties of the k-Wave grid and sensible stability criterion (a Courant-Friedrichs-Lewy stability value of 0.3).
For a homogeneous medium, medium.sound_speed
is given
simply as a scalar value (the density is not required unless velocity inputs or outputs are used). For a heterogeneous medium, medium.sound_speed
and medium.density
are given as Nx by Ny matrices with arbitrary numeric values. The initial acoustic pressure source.p0
is similarly defined as an Nx by Ny matrix. The measurement surface sensor.mask
is given either as a binary grid (i.e., an Nx by Ny matrix of 1's and 0's) representing the grid points within the computational grid at which the data will be recorded, or as a 2 by N matrix of Cartesian coordinates where the pressure values are calculated at each time step using interpolation. Note, in 3D the input matrices are instead Nx by Ny by Nz in size, and the Cartesian sensor points are given by a 3 by N matrix.
For the example given here, the geometric function makeDisc
is used to define an initial pressure distribution of a small filled disc, while makeCartCircle
is used to define
a Cartesian sensor mask with a set of evenly spaced points on a circle.
The simulation is invoked by calling kspaceFirstOrder2D
with the inputs described 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 through reds to black, zero pressures as white, and negative pressures as light to dark blue-greys. As the function runs, status updates and computational parameters are printed to the command line. When the time loop has completed, the function returns the time series recorded at the detector locations defined by the sensor mask. If the sensor mask is given as a set of Cartesian coordinates, the computed sensor_data
is returned in the same order. If the sensor mask is given as a binary grid, sensor_data
is returned using MATLAB's standard column-wise linear index ordering. In both cases, the recorded data is indexed as sensor_data(sensor_position, time)
. For a binary sensor mask, the pressure values at a particular time can be restored to the sensor positions within the computational grid using the utility function unmaskSensorData
. The final pressure within the computational domain can also be obtained by calling the simulation functions with two output arguments [sensor_data, field_data] = kspaceFirstOrder2D(...)
.
Additional properties of the medium, source, and sensor can be assigned using the remaining structure fields. For example, arbitrary power law absorption can be assigned using medium.alpha_power
and medium.alpha_coeff
, nonlinear governing equations can be used by setting medium.BonA
, a time varying pressure source can be used by defining source.p
and source.p_mask
, and a time varying velocity source can be used by defining source.ux
, source.uy
, source.uz
and source.u_mask
.
The default behaviour of the simulation functions within the toolbox can also be controlled through the use of optional input parameters. These are given as 'param'
, value
pairs following the primary function inputs. For example, the visualisation can be automatically saved as a movie by setting 'RecordMovie'
to true
; a plot of the initial pressure distribution can be automatically generated by setting 'PlotLayout'
to true
; the properties of the absorbing boundary condition can be adjusted; the interpolation method used to calculate the pressure values on a Cartesian sensor mask can be set; the smoothing of input data can be controlled; and the particle velocity recorded by the sensor mask can be returned.
For additional help using the k-Wave toolbox, see the list of examples and the individual function references. The examples are provided with source code that can both be viewed and executed from within the help browser.
![]() |
Troubleshooting | Modelling Wave Propagation | ![]() |
© 2009-2012 Bradley Treeby and Ben Cox.