kWave
A MATLAB toolbox for the timedomain
simulation of acoustic wave fields
 Getting Started
 Examples
 Initial Value Problems
 Example: Homogenous Propagation Medium
 Example: Using A Binary Sensor Mask
 Example: Defining A Sensor Mask By Opposing Corners
 Example: Loading External Image Maps
 Example: Heterogeneous Propagation Medium
 Example: Saving Movie Files
 Example: Recording The Particle Velocity
 Example: Defining A Gaussian Sensor Frequency Response
 Example: Comparison Of Modelling Functions
 Example: Setting An Initial Pressure Gradient
 Example: Simulations In One Dimension
 Example: Simulations In Three Dimensions
 Example: Photoacoustic Waveforms in 1D, 2D and 3D
 Time Varying Source Problems
 Example: Monopole Point Source In A Homogeneous Propagation Medium
 Example: Dipole Point Source In A Homogeneous Propagation Medium
 Example: Simulating Transducer Field Patterns
 Example: Steering A Linear Array
 Example: Snell's Law And Critical Angle Reflection
 Example: The Doppler Effect
 Example: Diffraction Through A Slit
 Example: Simulations In ThreeDimensions
 Sensor Directivity
 Example: Focussed Detector in 2D
 Example: Focussed Detector in 3D
 Example: Modelling Sensor Directivity in 2D
 Example: Modelling Sensor Directivity in 3D
 Example: Sensor Element Directivity in 2D
 Example: Focussed 2D Array with Directional Elements
 Photoacoustic Image Reconstruction
 Example: 2D FFT Reconstruction For A Line Sensor
 Example: 3D FFT Reconstruction For A Planar Sensor
 Example: 2D Time Reversal For A Line Sensor
 Example: 2D Time Reversal For A Circular Sensor
 Example: 3D Time Reversal For A Planar Sensor
 Example: 3D Time Reversal For A Spherical Sensor
 Example: Image Reconstruction With Directional Sensors
 Example: Image Reconstruction With Bandlimited Sensors
 Example: Iterative Image Improvement Using Time Reversal
 Example: Attenuation Compensation Using Time Reversal
 Example: Attenuation Compensation Using Time Variant Filtering
 Example: Automatic Sound Speed Selection
 Diagnostic Ultrasound Simulation
 Example: Defining An Ultrasound Transducer
 Example: Simulating Ultrasound Beam Patterns
 Example: Using An Ultrasound Transducer As A Sensor
 Example: Simulating Bmode Ultrasound Images
 Example: Simulating Bmode Images Using A Phased Array
 Numerical Analysis
 Example: Controlling The Absorbing Boundary Layer
 Example: Source Smoothing
 Example: Filtering A Delta Function Input Signal
 Example: Modelling Power Law Absorption
 Example: Modelling Nonlinear Wave Propagation
 Example: Optimising kWave Performance
 Using The C++ Code
 Elastic Wave Propagation
 Example: Explosive Source In A Layered Medium
 Example: Plane Wave Absorption
 Example: Shear Waves And Critical Angle Reflection
 Example: Simulations In Three Dimensions
 Functions  By Category
 Functions  Alphabetical List
 Release Notes
 License
kWave Toolbox 
kspaceFirstOrder1D
1D timedomain simulation of wave propagation
Syntax
sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor) sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, ...)
Description
kspaceFirstOrder1D
simulates the timedomain propagation of compressional waves through a onedimensional homogeneous or heterogeneous acoustic medium given four input structures: kgrid
, medium
, source
, and sensor
. The computation is based on a firstorder kspace model which accounts for power law absorption and a heterogeneous sound speed and density. If medium.BonA
is specified, cumulative nonlinear effects are also modelled. At each timestep (defined by kgrid.t_array
), the acoustic field parameters at the positions defined by sensor.mask
are recorded and stored. If kgrid.t_array
is set to 'auto'
, this array is automatically generated using makeTime
. An anisotropic absorbing boundary layer called a perfectly matched layer (PML) is implemented to prevent waves that leave one side of the domain being reintroduced from the opposite side (a consequence of using the FFT to compute the spatial derivatives in the wave equation). This allows infinite domain simulations to be computed using small computational grids.
For a homogeneous medium the formulation is exact and the timesteps are only limited by the effectiveness of the perfectly matched layer. For a heterogeneous medium, the solution represents a leapfrog pseudospectral method with a kspace correction that improves the accuracy of computing the temporal derivatives. This allows larger timesteps to be taken for the same level of accuracy compared to conventional pseudospectral timedomain methods. The computational grids are staggered both spatially and temporally.
An initial pressure distribution can be specified by assigning a matrix (the same size as the computational grid) of arbitrary numeric values to source.p0
. A time varying pressure source can similarly be specified by assigning a binary matrix (i.e., a matrix of 1's and 0's with the same dimensions as the computational grid) to source.p_mask
where the 1's represent the grid points that form part of the source. The time varying input signals are then assigned to source.p
. This can be a single time series (in which case it is applied to all source elements), or a matrix of time series following the source elements using MATLAB's standard columnwise linear matrix index ordering. A time varying velocity source can be specified in an analogous fashion, where the source location is specified by source.u_mask
, and the time varying input velocity is assigned to source.ux
.
The field values are returned as arrays of time series at the sensor locations defined by sensor.mask
. This can be defined in three different ways. (1) As a binary matrix (i.e., a matrix of 1's and 0's with the same dimensions as the computational grid) representing the grid points within the computational grid that will collect the data. (2) As the grid coordinates of two opposing ends of a line in the form [x1; x2]. This is equivalent to using a binary sensor mask covering the same region, however, the output is indexed differently as discussed below. (3) As a series of Cartesian coordinates within the grid which specify the location of the pressure values stored at each time step. If the Cartesian coordinates don't exactly match the coordinates of a grid point, the output values are calculated via interpolation. The Cartesian points must be given as a 1 by N matrix corresponding to the x positions, where the Cartesian origin is assumed to be in the center of the grid. If no output is required, the sensor
input can be replaced with an empty array []
.
If sensor.mask
is given as a set of Cartesian coordinates, the computed sensor_data
is returned in the same order. If sensor.mask
is given as a binary matrix, sensor_data
is returned using MATLAB's standard columnwise linear matrix index ordering. In both cases, the recorded data is indexed as sensor_data(sensor_point_index, time_index)
. For a binary sensor mask, the field values at a particular time can be restored to the sensor positions within the computation grid using unmaskSensorData
. If sensor.mask
is given as a list of opposing ends of a line, the recorded data is indexed as sensor_data(line_index).p(x_index, time_index)
, where x_index
corresponds to the grid index within the line, and line_index
corresponds to the number of lines if more than one is specified.
By default, the recorded acoustic pressure field is passed directly to the output sensor_data
. However, other acoustic parameters can also be recorded by setting sensor.record
to a cell array of the form {'p', 'u', 'p_max', ...}
. For example, both the particle velocity and the acoustic pressure can be returned by setting sensor.record = {'p', 'u'}
. If sensor.record
is given, the output sensor_data
is returned as a structure with the different outputs appended as structure fields. For example, if sensor.record = {'p', 'p_final', 'p_max', 'u'}
, the output would contain fields sensor_data.p
, sensor_data.p_final
, sensor_data.p_max
, and sensor_data.ux
. Most of the output parameters are recorded at the given sensor positions and are indexed as sensor_data.field(sensor_point_index, time_index)
or sensor_data(line_index).field(x_index, time_index)
if using a sensor mask defined as opposing ends of a line. The exceptions are the averaged quantities ('p_max'
, 'p_rms'
, 'u_max'
, 'p_rms'
, 'I_avg'
), the 'all' quantities ('p_max_all'
, 'p_min_all'
, 'u_max_all'
, 'u_min_all'
), and the final quantities ('p_final'
, 'u_final'
). The averaged quantities are indexed as sensor_data.p_max(sensor_point_index)
or sensor_data(line_index).p_max(x_index)
if using line ends, while the final and 'all' quantities are returned over the entire grid and are always indexed as sensor_data.p_final(nx)
, regardless of the type of sensor mask.
kspaceFirstOrder1D
may also be used for time reversal image reconstruction by assigning the time varying pressure recorded over an arbitrary sensor surface to the input field sensor.time_reversal_boundary_data
. This data is then enforced in time reversed order as a time varying Dirichlet boundary condition over the sensor surface given by sensor.mask
. The boundary data must be indexed as sensor.time_reversal_boundary_data(sensor_point_index, time_index)
. If sensor.mask
is given as a set of Cartesian coordinates, the boundary data must be given in the same order. An equivalent binary sensor mask (computed using nearest neighbour interpolation) is then used to place the pressure values into the computational grid at each time step. If sensor.mask
is given as a binary matrix of sensor points, the boundary data must be ordered using MATLAB's standard columnwise linear matrix indexing. If no additional inputs are required, the source
input can be replaced with an empty array []
.
Acoustic attenuation compensation can also be included during time reversal image reconstruction by assigning the absorption parameters medium.alpha_coeff
and medium.alpha_power
and reversing the sign of the absorption term by setting medium.alpha_sign = [1, 1]
. This forces the propagating waves to grow according to the absorption parameters instead of decay. The reconstruction should then be regularised by assigning a filter to medium.alpha_filter
(this can be created using getAlphaFilter
).
Note To run a simple photoacoustic image reconstruction example using time reversal (that commits the 'inverse crime' of using the same numerical parameters and model for data simulation and image reconstruction), the 
Inputs
The minimum fields that must be assigned to run an initial value problem (for example, a photoacoustic forward simulation) are marked with a *.

kWave grid structure returned by 

evenly spaced array of time values [s] (set to 



sound speed distribution within the acoustic medium [m/s] 

reference sound speed used within the kspace operator (phase correction term) [m/s] 

density distribution within the acoustic medium [kg/m^3] 

parameter of nonlinearity 

power law absorption exponent 

power law absorption coefficient [dB/(MHz^y cm)] 

optional input to force either the absorption or dispersion terms in the equation of state to be excluded; valid inputs are 

frequency domain filter applied to the absorption and dispersion terms in the equation of state 

two element array used to control the sign of absorption and dispersion terms in the equation of state 



initial pressure within the acoustic medium 

time varying pressure at each of the source positions given by 

binary matrix specifying the positions of the time varying pressure source distribution 

optional input to control whether the input pressure is injected as a mass source or enforced as a dirichlet boundary condition; valid inputs are 

time varying particle velocity in the xdirection at each of the source positions given by 

binary matrix specifying the positions of the time varying particle velocity distribution 

optional input to control whether the input velocity is applied as a force source or enforced as a dirichlet boundary condition; valid inputs are 



binary matrix or a set of Cartesian points where the pressure is recorded at each timestep 

cell array of the acoustic parameters to record in the form 

time index at which the sensor should start recording the data specified by 

time varying pressure enforced as a Dirichlet boundary condition over 

two element array specifying the center frequency and percentage bandwidth of a frequency domain Gaussian filter applied to the 
Note For a heterogeneous medium, 
Optional Inputs
Optional 'string', value pairs that may be used to modify the default computational settings.
Input  Valid Settings  Default  Description 




Interpolation mode used to extract the pressure when a Cartesian sensor mask is given. If set to 

(Boolean scalar) 

Boolean controlling whether the command line output is saved using the diary function with a date and time stamped filename. 

(string of data type) 

String input of the data type that variables are cast to before computation. For example, setting to 

(Boolean scalar) 

Boolean controlling whether the output data is cast back to double precision. If set to false, 

(binary matrix) or 

Binary matrix overlayed onto the animated simulation display. Elements set to 1 within the display mask are set to black within the display. 

(Boolean scalar) or 

Boolean controlling whether the pressure field is log compressed before display. The data is compressed by scaling both the positive and negative values between 0 and 1 (truncating the data to the given plot scale), adding a scalar value (compression factor) and then using the corresponding portion of a log10 plot for the compression (the negative parts are remapped to be negative thus the default color scale will appear unchanged). The amount of compression can be controlled by adjusting the compression factor which can be given in place of the Boolean input. The closer the compression factor is to zero, the steeper the corresponding part of the log10 plot used, and the greater the compression (the default compression factor is 0.02). 

(string cell array) 

Settings for 

(string) 

Name of the movie produced when 

(integer numeric scalar) 

The number of iterations which must pass before the simulation plot is updated. 

(Boolean scalar) 

Boolean controlling whether a four panel plot of the initial simulation layout is produced (initial pressure, sensor mask, sound speed, density). 

(Boolean scalar) 

Boolean controlling whether the perfectly matched layer is shown in the simulation plots. If set to 

(numeric two element vector) or 

[min, max] values used to control the scaling for 

(Boolean scalar) 

Boolean controlling whether the simulation iterations are progressively plotted. 

(numeric scalar or three element vector) 

Absorption within the perfectly matched layer in Nepers per grid point. 

(Boolean scalar) 

Boolean controlling whether the perfectly matched layer is inside or outside the grid. If set to 

(integer numeric scalar or three element vector) 

Size of the perfectly matched layer in grid points. To remove the PML, set the appropriate 

(Boolean scalar) 

Boolean controlling whether the displayed image frames are captured and stored as a movie using 

(Boolean scalar or three element vector) 

Boolean controlling whether 
Outputs
If sensor.record
is not defined by the user:

time varying pressure recorded at the sensor positions given by 
If sensor.record
is defined by the user:

time varying pressure recorded at the sensor positions given by 

maximum pressure recorded at the sensor positions given by 

minimum pressure recorded at the sensor positions given by 

rms of the time varying pressure recorded at the sensor positions given by 

final pressure field at all grid points within the domain (returned if 

maximum pressure recorded at all grid points within the domain (returned if 

minimum pressure recorded at all grid points within the domain (returned if 

time varying particle velocity in the xdirection recorded at the sensor positions given by 

maximum particle velocity in the xdirection recorded at the sensor positions given by 

minimum particle velocity in the xdirection recorded at the sensor positions given by 

rms of the time varying particle velocity in the xdirection recorded at the sensor positions given by 

final particle velocity field in the xdirection at all grid points within the domain (returned if 

maximum particle velocity in the xdirection recorded at all grid points within the domain (returned if 

minimum particle velocity in the xdirection recorded at all grid points within the domain (returned if 

time varying particle velocity in the xdirection recorded at the sensor positions given by 

time varying acoustic intensity in the xdirection recorded at the sensor positions given by 

average acoustic intensity in the xdirection recorded at the sensor positions given by 
Examples
 Simulations In One Dimension
 Photoacoustic Waveforms in 1D, 2D and 3D
 Filtering A Delta Function Input Signal
 Modelling Power Law Absorption
 Modelling Nonlinear Wave Propagation
See Also
fft
, ifft
, getframe
, kspaceFirstOrder2D
, kspaceFirstOrder3D
, makeGrid
, makeTime
, movie2avi
, smooth
, unmaskSensorData
interpftn  kspaceFirstOrder2D 
© 20092014 Bradley Treeby and Ben Cox.