Rapid calculation of acoustic fields from arbitrary continuous-wave sources

The calculation of the acoustic field from single element and phased array transducers in homogeneous media is typically performed using semi-analytical approaches based on the spatial impulse response^1–3 or Rayleigh integral.^4,5 These methods have been widely used and validated, particularly for the design of multi-element linear arrays. However, one restriction of these approaches is that the compute time is proportional to the number of grid points at which the pressure field is evaluated. For computing the three-dimensional (3D) field produced by a multi-element array, this can become computationally prohibitive.⁶ If the input field is known in a two-dimensional (2D) plane, alternative methods based on the angular spectrum approach can be much faster.⁷ However, this also has limitations if the waves are propagating in more than one direction (in this case, several simulations must be performed and the total field obtained by superposition).

Here, an alternative formulation is derived to calculate the 3D wave field from a phased array transducer (or other acoustic source) of arbitrary shape when driven by a single frequency continuous wave excitation with spatially varying amplitude and phase. The method belongs to the family of Fourier and k-space methods that use exact propagators, typically expressed in the Fourier domain, to map from an input field to an output field at a later time or position. For example, the angular spectrum method uses an exact propagator to map from one spatial 2D plane to the next,⁷ and k-space methods for initial value problems use an exact propagator to map from an initial impulsive pressure distribution to the pressure field at time t > 0.^8,9 In the current work, an exact propagator is derived, which maps from the spatially varying amplitude and phase at t = 0 to the pressure field at time t > 0 when the source is subject to single frequency continuous wave excitation. This allows the acoustic pressure field at all spatial positions to be calculated in a single step without numerical quadrature. The formulation is validated by comparison with the fast near-field method (FNM)⁵ for several transducer geometries.

The linear wave equation for a homogeneous medium subject to a time-varying source term S(x, t) is given by

where p(x, t) is the acoustic pressure as a function of position $x∈ℝn, n=1,2,3$⁠, and time $t∈ℝ+$⁠, and c₀ is the small signal sound speed. The free-space Green's function for the wave equation can be written in the spatial frequency domain (or k-space) in the following form:^9,10

where T_P is the time propagator

Here $k=|k|$ is the scalar wavenumber, and $k∈ℝn$ is the wavevector. Using the Green's function, it is possible to calculate the pressure field at time t by convolving the source term S(x, t) (which is assumed to be causal) with the Green's function¹⁰ 

The last two integrals account for the initial conditions at time t = 0.

If the acoustic source is defined as a conventional mass source M(x, t), which represents the time rate of input of mass per unit volume in units of kg m⁻³ s⁻¹, the source term S(x, t) in Eq. (1) is given by

Here, the mass source is assumed to be a single frequency continuous wave sinusoid with a spatially varying amplitude A(x) and phase ϕ(x) in the form

where ω₀ is the source frequency in rad s⁻¹, and M(x, t) = 0 for t < 0. The use of a complex exponential allows two continuous wave sources to be encoded onto the real and imaginary parts of M(x, t) with a relative phase offset of π/2. This allows the magnitude and phase of the resulting pressure field to be calculated in a single step. Alternatively, a real-valued sine or cosine source term could be used, and the magnitude and phase calculated from the real-valued pressure field at two different times or phases. Using Eq. (5), the source term then becomes

where the initial conditions for the acoustic pressure at time t = 0 are given by

The second initial condition is derived from the time derivative of the linear pressure density relation (equation of state), where the time rate of change of pressure is proportional to the time rate of mass injection.

Combining Eqs. (2), (3), (4), (7), and (8) and changing the order of integration then gives

where

The choice of a single frequency continuous wave source in Eq. (6) allows the time integral to be solved analytically, which gives

where the limits as the denominator goes to zero are

Note, these limits must be included when numerically calculating I(k, t).

Returning to Eq. (9), the integrals over $x′$ and k can be recognised as the forward and inverse Fourier transforms, and thus an exact expression for the complex acoustic pressure at time t can be written succinctly as

where $F$ and $F−1$ are the forward and inverse Fourier transforms, respectively. As mentioned above, the spatially varying wave fields for cosine and sine excitation are encoded on the real and imaginary parts of p(x, t), respectively. This allows the amplitude B(x) and phase θ(x) of the acoustic pressure to be extracted from the complex pressure field, where

Here the phase θ(x) varies between + π and −π.

Practically, Eq. (13) allows the acoustic pressure field at any time t to be calculated directly (i.e., without time stepping or numerical quadrature) from a map of the spatially varying amplitude and phase at time t = 0 using two Fourier transforms. There are no restrictions on the input amplitude A(x) and phase ϕ(x), and thus single or multi-element transducers (or other acoustic sources) with any shape, apodisation, and phase delay can be modeled easily.

The numerical implementation of Eq. (13) can be achieved as follows. First, given a source distribution within a computational domain $Ω⊂ℝn$⁠, the domain Ω is discretised using a uniform Cartesian grid. The Fourier transforms can then be computed efficiently using the fast Fourier transform (FFT). For a given distribution of amplitude A(x) and phase ϕ(x) at time t = 0 (each represented by an n-dimensional matrix), the forward n-dimensional Fourier transform of the source term is calculated. Next, this is multiplied by the propagator I(k, t) given in Eq. (11), where the components of the wavevector for a grid spacing of Δx are defined as

if the number of grid points N_x is even, and

if N_x is odd (similarly for the y and z directions). Finally, the complex pressure field is calculated by taking the inverse Fourier transform according to Eq. (13).

Note, in some cases, the band-limiting imposed on the numerical implementation through the use of a discrete Fourier transform can result in Gibbs' oscillations at the front edge of the wave.¹¹ These can be avoided by including a smooth ramp function into the source term. If a suitable choice is made for the ramp function, the time integration can still be performed analytically. Here, a half-cosine ramp is used. In this case, two additional integrals are added to the expression for I(k, t) given in Eq. (10) to account for the ramp. The modified propagator is derived in the  Appendix, and is used for the numerical experiments presented in Sec. IV.

The use of the FFT implicitly assumes that the source distribution is periodically repeated, the effect of which is that the waves leaving one side of the computational domain will reappear at the opposite side. This can be alleviated by zero padding the size of the domain appropriately. For calculating the steady state amplitude and phase in Ω, this can be achieved as follows. First, the time t at which the wave field is calculated is selected to be slightly larger (by a factor of d) than the time for the waves to propagate across the longest grid diagonal, i.e., in 3D,

Here Δx is the spacing between the grid points (which is assumed to be isotropic), and N_x, N_y, and N_z are the number of grid points in each Cartesian direction. Second, the amplitude and phase matrices A(x) and ϕ(x) are zero padded in each Cartesian direction by a propagation distance that is slightly larger than c₀t (by a factor of h) such that waves at the edge of the domain will not wrap back to the opposite edge in time t (the same padding is used in all directions).

In addition to avoiding wave wrapping, it is important to choose grid sizes with small prime factors to maintain the computational efficiency of the FFT. To achieve this, the size of the padding in each Cartesian direction should be selected as the value in the range [pad_min, pad_min + s], which gives the total grid size with the smallest prime factors, where s is the search range and

After grid expansion, the amplitude and phase are calculated following Eqs. (13) and (14). Finally, the matrices B(x) and θ(x) are truncated back to the original domain size Ω. For the examples shown in Sec. IV, d was set to 1.5, h was set to 1.1, and s was set to 50. For brevity, the calculation of B(x) and θ(x) using this method is referred to herein as the acoustic field propagator (AFP) model.

One constraint of the AFP model is that sources (i.e., the matrices of spatially varying amplitude and phase) must be defined on a regular Cartesian grid. This can lead to staircasing errors if the source geometry does not conform exactly to the grid. However, for spectral methods, the underlying band-limited interpolant (BLI) is known analytically.^11–13 Therefore, it is possible to represent off-grid sources by convolving the BLI with the source geometry in continuous space, and then evaluating the output at the grid nodes. This approach is described in more detail in Ref. 14. As the discretisation must already be sufficient to meet the Nyquist criterion for the source frequency ω₀, any source shape can be represented in this manner. As with other methods based on a Fourier collocation spectral method using non-impulsive sources (e.g., the k-Wave toolbox¹⁵), an additional scaling factor of 2c₀/Δx is also required to account for the spatial spread of the BLI.

The AFP model was implemented using the c++ programming language with calculations and data storage performed using the single precision (float) data type, with the exception of the propagator I(k, t) given in Eq. (A10), which was calculated in double precision due to the large exponents. The code was designed to be executed on a single shared-memory computer. FFTs were computed in-place using complex-to-complex transforms from the Intel Math Kernel Library (Intel Corp., Santa Clara, CA) library. Computational kernels were parallelised at the thread-level using OpenMP 4.0. The inner-most loops were vectorised using OpenMP single instruction, multiple data (SIMD) directives, as the operations (element-wise matrix multiplications and additions) do not carry any dependencies between the iterations. Input and output files were stored in hierarchical data format (HDF5). The creation of input files and the visualisation of output files was performed using matlab.

The memory usage of the AFP model can be estimated directly from the grid size. As the calculations are performed in-place, only one complex matrix is required, where

Here M_x, M_y, and M_z are the number of grid points in each Cartesian direction after the grid expansion to avoid wave wrapping. Space complexity is therefore of the order O(m), where m = M_xM_yM_z, and time complexity is of the order $O(m log m)$⁠. For representative grid sizes, ∼30% of the total computation time is spent on memory allocation and reading and writing the input and output files, ∼50% is spent performing the FFTs, and ∼20% is spent on other matrix operations, including calculation of the propagator I(k, t).

To illustrate that the AFP model correctly encodes the Green's function, a point source in one dimension (1D), 2D, and 3D was simulated. In each case, the point source was modeled as an amplitude matrix with a single grid point set to an amplitude of one (with zeros elsewhere), and the phase set to zero. The grid size along the longest dimension was set to 500 grid points, the grid spacing was Δx = 100 μm, the sound speed was c₀ = 1500 m s⁻¹, and the frequency was 1 MHz. The amplitude decay with distance is plotted in Fig. 1(a), along with the theoretical values for geometric spreading shown for comparison. To avoid the origin singularity for the analytical values, all amplitudes are normalised to 1 at a distance of 1 mm. There is no discernible difference between the curves.

To investigate the accuracy more quantitatively, and to study the effect of the time expansion factor d discussed in Sec. III B, the simulations were repeated for different values of d, and the L_∞ error between the numerical simulations and the theoretical values calculated. (No difference was observed for changes in the grid expansion factor h provided h > 1.) The error is plotted in Fig. 1(b), and slightly reduces as the time expansion factor is increased. This is due to very small Gibbs' oscillations present at the front edge of the wave (see the  Appendix) being moved further from the region of interest. The absolute error is also larger for simulations in higher dimensions, due to the more rapid amplitude decay. However, in all cases, the maximum error is below 10^–3, which is sufficient for most practical purposes.

In addition to simple sources, several numerical simulations with more complex transducer geometries were performed. First, the steady state response of a rectangular aperture was compared with the FNM as implemented in the FOCUS toolbox.^5,16 This approach is comparable to evaluating the Rayleigh–Sommerfeld integral, but converges more rapidly by using an equivalent integral expression that removes numerical singularities. A comparison between the two models for a 4.1 × 4.1 mm rectangular piston uniformly driven by a continuous wave sinusoid at 2 MHz is shown in Fig. 2. The sound speed was set to c₀ = 1500 m s⁻¹. The AFP model was discretised using a 192 × 64 × 64 grid with a grid point spacing of Δx = 100 μm, and an extended grid size of 576 × 432 × 432. The FNM model was evaluated on the central plane using 200 integration points. There is excellent agreement between the two models, both in the beam plots and the axial and lateral profiles.

To quantify the error, the simulations were repeated using the same source and discretisation, with the frequency varied from 2 points per wavelength (7.5 MHz) to 20 points per wavelength (750 kHz). The L_∞ error was then computed from the on-axis pressure, not including the source plane. The error is shown in Fig. 2(d). At the Nyquist limit, the maximum L_∞ error is significant (10% of the peak pressure). However, this falls rapidly, and the maximum error is on the order of 0.1% by 3 points per wavelength, which is sufficient for most practical modeling scenarios. Beyond 17 points per wavelength, the maximum error is reduced to below 0.01%. The error rates are similar to the point source comparison shown in Fig. 1(b).

A second comparison between the two models for a single-element focused bowl transducer uniformly driven by a continuous wave sinusoid at 1 MHz is shown in Fig. 3. The bowl was defined with an aperture diameter of 20 mm and a radius of curvature of 20 mm, and the sound speed was set to c₀ = 1500 m s⁻¹. The AFP model was discretised using a 320 × 96 × 96 grid with a grid point spacing of Δx = 250 μm, and an extended grid size of 900 × 675 × 675. The spatially varying amplitude for the AFP model was defined by taking the sum of BLIs uniformly distributed over the surface of the bowl at twice the resolution of the simulation grid.¹⁴ The bowl for the FNM model was defined as a spherical shell in an infinite rigid baffle,¹⁷ and the output was evaluated on the central plane using 200 integral points. Again, there is excellent agreement between the two models, both in the magnitude and phase of the beam plots, and the axial profile. Note, for a focused bowl transducer, simulations based on the FNM or Rayleigh integral will only be accurate when the transducer diameter is large compared to both the transducer height and the acoustic wavelength (this is satisfied in this example).

One of the key advantages of the AFP model is its computational speed when calculating the pressure field over a large 3D volume. To demonstrate this, a series of benchmark calculations were performed and compared against the FNM model as implemented in the FOCUS toolbox (version 0.922, using the fnm_cw function set to run using 8 threads and 20 integration points). The calculations were performed using a desktop computer running Windows 10 (Microsoft Corp., Redmond, WA) with an 8-core Xeon E5-1660 v3 3.0 GHz processor, and 64 GB of 2133 MHz DDR4 memory. The source was defined as a single 4.1 × 4.1 mm rectangular piston uniformly driven by a continuous wave sinusoid at 500 kHz. The output was evaluated on a cubic grid of 8 × 8 × 8 mm, with an increasing number of grid points within the 3D volume. The FNM was also used to compute the pressure on the central 2D plane to compare the compute time when the entire 3D volume may not be of interest.

The compute times for a range of grid sizes based on ten averages are shown in Fig. 4. The grid dimensions correspond to the actual domain size. An additional grid extension was used for the AFP model to avoid wave wrapping as described in Sec. III B, e.g., for a grid dimension of 256, the expanded grid size used by the AFP model was 1024. For larger grid sizes, the AFP model is almost an order of magnitude faster than the FNM model for computing the entire 3D volume. This is a significant gain when the wave field across a large domain is of interest. Moreover, the compute time for methods based on the Rayleigh integral (including the FNM model) scales with the number of sources, while the compute time for the AFP model is only dependent on the grid size, and not the complexity of the source. This makes the AFP model particularly efficient for modeling large multi-element arrays. However, if only a subset of the domain is required (e.g., a 2D plane), the FNM model is substantially faster. The key limitation of the AFP model is that the entire 3D domain must always be calculated, and the grid must have at least two points per acoustic wavelength. Conversely, models based on the Rayleigh integral do not have either restriction, and in the limit can be used to compute the pressure field at just a single point. Considering memory usage, if the whole domain is of interest, typically at least the complex pressure over the 3D volume will be stored, regardless of the model used. However, for the AFP model, the pressure over the expanded grid must be calculated (see Sec. III D), which in some cases may be an order of magnitude larger.

To demonstrate the wider capabilities of the AFP model for simulating the output from multi-element arrays, two numerical simulations were performed. First, the steady state response of a linear array transducer typical of those used for diagnostic ultrasound imaging was computed. The array was defined with 32 active rectangular elements, with an element pitch of 281.25 μm and an elevation height of 6 mm. The array was driven at 4 MHz, and the phase delays were defined to focus the beam at an axial distance of 20 mm at a steering angle of 20°. The AFP model was discretised using a 384 × 256 × 85 grid with a grid point spacing of Δx = 93.75 μm, and an extended grid size of 1200 × 1080 × 864. The amplitude and phase over the entire 3D domain was computed in 13 s using 8.4 GB of memory (using the hardware described in Sec. IV C). The amplitude of the response in the lateral plane is shown in Fig. 5, both with and without amplitude apodisation.

A second example for a multi-element hemispherical array with 512 bowl-shaped elements typical of those used for transcranial ultrasound therapy is shown in Fig. 6.¹⁸ The individual elements had a diameter of 5 mm, a frequency of 500 kHz, and were randomly positioned on the surface of a hemisphere with radius 150 mm. The AFP model was discretised using a 266 × 512 × 512 grid with a grid point spacing of Δx = 625 μm and an extended grid size of 1568 × 1800× 1800. The amplitude and phase over the entire 3D domain was computed in 71 s using 38 GB of memory.

A Green's function solution that maps from a continuous wave source with an arbitrary distribution of amplitude and phase at t = 0 to the pressure field at a later time t has been derived. The solution is based on the free-space Green's function expressed in the spatial Fourier domain or k-space, where the time convolution is performed analytically. This allows the steady state amplitude and phase of the pressure field at all spatial positions to be calculated without time-stepping or numerical quadrature using only two FFTs. Compared to semi-analytical approaches based on the spatial impulse response or Rayleigh integral, this allows the rapid calculation of the pressure field in large 3D domains. This is of particular relevance for modeling the output from non-planar transducer arrays, including hemispherical arrays used for transcranial ultrasound therapy. The codes will be released as part of the open-source k-Wave toolbox.¹⁵ 

This work was supported by the Engineering and Physical Sciences Research Council, United Kingdom, Grant Nos. EP/L020262/1, EP/M011119/1, and EP/P008860/1. This work was also supported by The Ministry of Education, Youth and Sports of the Czech Republic from the Large Infrastructures for Research, Experimental Development and Innovations project IT4Innovations National Supercomputing Center–LM2015070, and the National Programme of Sustainability (NPU II) project IT4Innovations Excellence in Science–LQ1602.

The band-limiting imposed by the use of the discrete Fourier transform can lead to Gibbs' oscillations in the calculated pressure field, particularly at the front edge of the wave where the pressure changes rapidly from zero. To avoid this, a time varying ramp function R(t) can be added to the mass term source, where

If a suitable choice is made for the ramp function, the time integration can still be performed analytically. Here, a half-cosine ramp is chosen, which varies smoothly from 0 to 1 over one period of the input signal (i.e., the ramp has a frequency of ω₀/2). This gives

For $0≤t≤2π/w0$⁠, the source term then becomes

where the initial conditions for the acoustic pressure at time t = 0 are given by

The complex pressure field is given by Eq. (9), where the integral is now

where

Again these integrals can be computed analytically. For I₁ this gives

and for I₂,

and for I₃,

Combining the integrals and simplifying then gives the propagator when a half-cosine ramp is included

The limits as the denominator go to zero are given by

An illustrative example of a pressure field in 1D computed using the propagator with and without the start-up ramp is given in Fig. 7. In this case, the time t was set such that the front of the wave was within the grid. The three rows show the real and imaginary parts of the complex pressure (corresponding to cosine and sine input signals, respectively), and the pressure magnitude. When no ramp is included, Gibbs' oscillations can be seen at the beginning of the signal, particularly for the cosine propagator, which has a sharp edge. When a ramp is included, these oscillations are damped.

A similar approach can be used to compute propagators for other start-up ramps, e.g., a half-cosine ramp that varies smoothly over two periods of the input signal. If the medium is absorbing with the absorption following a frequency power law, the lossless wave equation could also be replaced with the fractional Laplacian wave equation,^19,20 for which the time propagator for the k-space Green's function is also known.¹¹