k-Wave Toolbox Previous   Next

Point Source In A Homogeneous Propagation Medium Example

Overview

This example provides a simple demonstration of using k-Wave for the simulation and detection of a time varying pressure source 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.

% 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 = 50e-3/Nx;      % pixel width [m]
dz = dx;            % pixel height [m]
kgrid = makeGrid(Nx, dx, Nz, dz);

As the k-space computation is based heavily on the FFT, the computation will be fastest for grids where the number of pixels is given by a power of two (see fft2).

 Back to Top

Defining the medium properties

For a homogeneous medium, the sound speed and density are set as scalars in SI units. Power law absorption can also be optionally set by assigning values to medium.alpha_power and medium.alpha_coeff.

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

It is useful to note that for a homogeneous medium, the wave propagation is not dependent on the ambient density. However, an arbitrary value for medium.density must still be given due to code's inherent ability to deal with heterogeneous media.

 Back to Top

Defining the time varying pressure source

A time varying pressure source is defined 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 pixels that form part of the source. The time varying input signal is then assigned to source.p. This must be the same length as kgrid.t_array and 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 column-wise linear matrix index ordering. Here a sinusoidal input is assigned to a single source element. To avoid numerical stabilities, the input should first be filtered using filterTimeSeries (see the Filtering A Delta Function Input Signal Example for more information).

% create the time array
[kgrid.t_array dt] = makeTime(kgrid, medium.sound_speed);

% define a single source element
source.p_mask = zeros(Nz, Nx);
source.p_mask(end - Nz/4, Nx/2) = 1;

% define a time varying sinusoidal source
source_freq = 0.25e6;
source_mag = 2;
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);

% smooth the source
source.p = filterTimeSeries(kgrid, medium, source.p);

 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 centered origin k-space grid (see the Homogeneous Propagation Medium Example). Here a single sensor position is defined within a binary sensor mask.

% define a single sensor element
sensor.mask = zeros(Nz, Nx);
sensor.mask(Nz/4, Nx/2) = 1;

 Back to Top

Running the simulation

The computation is invoked by calling kspaceFirstOrder2D with the inputs 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 p_final] = kspaceFirstOrder2D(kgrid, medium, source, sensor);

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

Running k-space simulation...
  dt: 78.125ns, t_end: 47.1094us, time steps: 604
  input grid size: 128 by 128 pixels (50 by 50mm)
  maximum supported frequency: 1.92MHz
  precomputation completed in 0.26407s
  starting time loop...
  computation completed in 7.1625s

When the time loop has completed, the function returns the recorded time series at each of sensor locations defined by sensor.mask. A plot of the time series driving the source and the time series recorded at the sensor position is given below.

The final pressure field within the computational domain is also returned.

 Back to Top


© 2009, 2010 Bradley Treeby and Ben Cox.