Representing arbitrary acoustic source and sensor distributions in Fourier collocation methods^a)

One of the fundamental tasks in numerical acoustics is to compute solutions p(x,t) to the wave equation

Here, p(x,t) is the acoustic pressure at position $x∈ℝd$ (in d-dimensions) and time $t∈ℝ+$⁠, and the acoustic medium is specified by a constant scalar sound-speed c₀. A source term S(x,t) is also included. In many cases, the source can be separated into a spatial distribution s(x) and a temporal waveform f(t), where

Several different classes of numerical techniques are available for solving the wave equation. For time-domain modeling, collocation methods are popular. These are numerical methods that find approximate solutions, which satisfy the model equations at a finite number of grid points.¹ This paper is concerned in particular with Fourier spectral collocation methods,² which use an equispaced, orthogonal, grid of collocation points {x_j} covering the domain. For example, in one-dimension, the set of discrete grid points can be written as x_j = jΔx, where Δx is the grid spacing, j = 0, 1,…, N − 1, and N is the number of grid points.

Fourier collocation methods seek an approximate solution to the wave equation whose spatial part can be written as a Fourier series. As well as leading to efficient numerical algorithms based on the fast Fourier transform,^3–5 Fourier schemes have the advantage that the nature of the approximation has a clear physical interpretation: the spatial part of any solution must be in the set $B$ of functions that are supported by the set of wavenumbers bounded from above and below by ±π/Δx (where the grid spacing has been assumed to be equal in each dimension).² 

Despite their widespread use, one question that has not been widely studied is how sources in Fourier schemes should be modelled. In particular, how can a source be incorporated if it is not spatially band limited, i.e., $s(x)∉B$? For example, consider a source distribution corresponding to the surface of a physical ultrasound transducer, which may be a bowl or a planar disk. In this case, the support of the source—the region in which s(x) is non-zero—is a two-dimensional surface embedded in $ℝ3$ and is therefore not band limited. Furthermore, it is likely that few, if any, grid points coincide exactly with this surface. This highlights a related question: as the source in a collocation method can only be defined by assigning values at the grid points, which of the grid points should be used as source points, and what should the source grid weights s_j be to best approximate the source distribution s(x)? In typical practical applications, grid points close to the required source surface are defined as source points,⁶ but this naive approach can result in “staircasing” and serious errors in the acoustic field.⁷ 

Progress can be made by realising that the closest approximation to a source that can be made within a Fourier collocation method is the projection of the source distribution onto the set of band limited functions $B$⁠. Thus, a band limited source distribution can be defined, and the source grid weights can be generated by sampling this distribution at the grid points. This paper describes a method for performing this band-limiting operation, and examines the process using a range of examples relevant to problems in ultrasound. Section II describes how the band-limiting is performed via convolution of the source shape with the band limited delta function. Section III then provides several numerical experiments which demonstrate the accuracy and convergence of this approach for different shaped ultrasound sources.

Any source can be written as a convolution of the source distribution function with a point source. Let $δ(x;ξ)$ denote a delta function centred on a point $ξ∈ℝd$⁠. The spatial distribution of the source is then, trivially

where Γ = supp(s) is the support of s(x). To band-limit this source, the delta function $δ∉B$ should be replaced with its band limited version $b∈B$⁠, which is the projection of δ onto $B$⁠. This gives

where $s̃(x)$ is the band limited source distribution. In the general case this convolution cannot be solved exactly, and so it must be numerically approximated. This can be done by replacing the integral with a discrete sum,

where s_j is the source grid weight at the jth grid point x_j, i are a set of M indices (where M is the number of integration points), C_i are quadrature weights, and $ξi$ are the integration points.

Computing the discrete convolution given in Eq. (5) involves two main tasks. First, an analytical expression is required for the band limited approximant $b(x;ξ)$ to a Dirac delta function centred at an arbitrary point $ξ$⁠. Second, a strategy is needed for effectively and efficiently discretising the convolution. This involves choosing discrete integration points $ξi$ covering Γ (see Fig. 1), and selecting corresponding quadrature weights C_i.

There have been a number of past works which have attempted to represent delta functions in the context of particular numerical methods. For example, Waldén,⁸ approximated one-dimensional delta functions in space within finite-difference and finite-element methods using compactly supported functions that satisfied some number of moment conditions. Similarly, Tornberg and Engquist⁹ did this for the finite-difference method in the multidimensional case. Petersson et al.¹⁰ note that discretisations that only satisfy moment conditions will introduce spurious oscillations, and so add a number of smoothness conditions to their delta function approximations. For Fourier collocation methods, the band-limiting of delta functions can be accomplished analytically.

The one-dimensional delta function (positioned at the origin) has a Fourier transform which is equal to one for all wavenumbers. However, in Fourier collocation methods, the wavenumbers are restricted to a finite set. Specifically, the use of discrete Fourier transforms leads to the set of spatial wavenumbers

It follows then that the band limited approximation to a delta function is the function whose Fourier transform is equal to one for all supported wavenumbers, and equal to zero for all others.

To translate the band limited delta function to an arbitrary position, a Fourier-space shifting operator can be applied. To begin with, consider a one-dimensional grid in which the number of grid points N is odd. In Fourier-space, the jth component of the band limited delta function is then simply given by the shift operator

To get an expression for b(x;ξ), the Fourier series is evaluated including the shift operator

The imaginary components in this sum cancel because of conjugate symmetry about j = 0, which gives

and the series simplifies to yield

A scaled version of this function is often referred to as the Dirichlet kernel, appearing in many texts including Hesthaven et al.¹¹ 

If the number of grid points N is even, further steps are required. As a real-valued point source should have a real-valued representation on the grid, the Fourier coefficients of the band limited delta function must be conjugate symmetric. A full derivation follows in the  Appendix, but in short, the imaginary part of the Nyquist coefficient must be dropped so that

The Fourier series then sums to

This band limited delta function expression contains an imaginary sinusoid, but this is zero at the grid points for any shift ξ. It also contains a real-valued Nyquist sinusoid that is zero at the grid points. Note that $b(x;ξ)≠b(x−ξ;0)$ in the even case, except when ξ is a multiple of Δx. Figure 2 depicts the odd and even (real component only) band limited delta functions in one dimension.

The band limited delta functions given in Eqs. (10) and (12) can be extended into higher dimensions via the product of one-dimensional band limited delta functions

where d is the number of dimensions, and x^(p) denotes the pth component of the vector x.

Discretisation of the band-limiting convolution in Eq. (5) requires a number of steps. First is the selection of a finite number of integration points $ξi$⁠. In general, these should be placed according to a given quadrature rule. In this paper, a uniform sampling strategy is used, with the spacing between integration points being equal (or nearly equal if this is not feasible), and the outermost integration points being offset from the source boundary by half the inter-point spacing. The latter ensures that outer and inner integration points cover equal portions of the source region. To avoid staircasing effects, integration points are placed such that they also conform to the source boundary, rather than to the computational grid.

Second, the quadrature weights C_i must be chosen. These account for any difference in the spacing of the integration points relative to the grid spacing, and are all equal to the ratio of these spacings with a uniform sampling strategy. As an example, for a two-dimensional source embedded in a three-dimensional domain, the quadrature weights are given by

where A is the area of the source's support, M_grid is the area of the source in units of grid squares, and M_integration is the number of integration points that has been used. Note that the number of integration points will be measured relative to the grid in this paper. The phrase “integration point density” will refer to the ratio M_integration/M_grid, and a source will be referred to as upsampled if M_integration/M_grid > 1 and undersampled if M_integration/M_grid < 1.

A set of grid points ${xl}⊆{xj}$ must be chosen over which the band limited source $s̃$ is evaluated to give the source grid weights $sl=s̃(xl)$⁠. In principle, every grid point should be used because band-limiting means that $supp(s̃)=ℝd$ (a function with compact support in the wavenumber domain will have infinite support in the spatial domain, and vice versa). However, it can be computationally beneficial to restrict the discretised source to grid points which lie near the support of the true source. This limits the number of band limited delta function evaluations and the subsequent memory requirements for the source grid weights. To see why the number of grid points used to discretise the source can be considerably reduced without introducing substantial error, note that for a large domain size, the band limited delta functions given in Eqs. (10) and (12) can be approximated by a sinc function

In an analogous manner to Eqs. (10) and (12), this function derives from the continuous inverse Fourier transform of a boxcar function whose limits are the minimum and maximum supported wavenumbers. Figure 3 illustrates the accuracy of this approximation as N increases for a shift distance of ξ = Δx/2. It can be seen that the error drops below 1% when the grid size reaches approximately 100 grid points, which is a relatively modest size in the context of ultrasound simulation.⁷ 

The envelope of the sinc approximation decays at a rate of approximately $Δx/π|x−ξ|$⁠, thus a magnitude threshold can be defined beyond which contributions from a given integration point can be ignored. Denoting this threshold ε, a given band limited delta function thus only needs to be evaluated to

grid points on either side of the centre of the integration point ξ_i. For example, with ϵ = 10% each band limited delta function needs to be evaluated to a distance of only m = 4 grid points, and for ϵ = 1% this becomes m = 32 grid points. Note that in multiple dimensions the benefit of truncation compounds, as diagonal decay rates are higher than those along the grid axes. The benefit of truncation also increases with domain size, as the truncation distance m is independent of this. Figure 4 depicts the distance at which the truncation thresholds lie in two dimensions, illustrating the reduction in extent that truncation provides.

To demonstrate the accuracy and utility of using source grid weights calculated using the approach described in Sec. II, a series of numerical experiments was conducted using spatial source distributions relevant to ultrasound. For each numerical experiment, the wave equation was solved using one of two Fourier collocation methods. The first was the acoustic field propagator (AFP),⁵ which uses a Green's function method to solve Eq. (1) when the source is time-harmonic. Free-space is approximated by evaluating the field at a time when it has propagated over the whole domain, but with computations performed using an extended domain to prevent periodic wrapping effects. The second was the open-source k-Wave toolbox,¹² which is not restricted to time-harmonic problems and uses a dispersion-corrected finite-difference scheme for time-stepping.⁴ Here, free-space is approximated using a perfectly matched layer. Note, for time-harmonic problems, the AFP and k-Wave give solutions that match to a high degree of accuracy.

To distinguish the band limited source distributions from the conventional approach of using the nearest available grid points to represent the source shape, the following terminology is introduced. When the integration points are allowed to lie anywhere within the support of the true source, the resulting set of source grid weights will be called an off-grid source. When the integration points are instead restricted to the grid points (e.g., nearby to the true source), the resulting set of source grid weights will be called an on-grid source. On-grid sources will often be subject to staircasing effects and the errors that result from this. Figure 5 depicts distributions of integration points for both on- and off-grid sources with different geometries:

1. (Top-left) An on-grid approximation to an arc source. Here, the integration points are restricted to the grid points and are thus misaligned with the true source, resulting in staircasing.

2. (Top-right) A staircase-free, off-grid arc source. The integration points are spread equally over the arc, with the end points offset from the ends of the source by half the inter-point spacing. The integration point spacing is approximately half that of the grid spacing.

3. (Bottom-left) A disk source. The integration points are chosen as concentric circles whose number increases linearly with radius. This ensures all points are approximately equidistant from their neighbours. The outermost points are offset from the edge of the source, and the integration point spacing is approximately half that of the grid spacing.

4. (Bottom-right) A square source that is not aligned with the grid. Here the integration points form a regular grid, but one which is aligned with the source boundaries rather than the computational grid. Once again the outermost points are offset from the edge of the source, and the integration point spacing is approximately half that of the grid spacing.

To illustrate the elimination of staircasing errors when off-grid sources are used, the field from a 5 mm line source was simulated in two dimensions using the AFP. A 12.6 mm square domain was discretised using a grid spacing of 98 μm giving a grid size of 128 × 128 grid points. The source was placed 5 mm from the centre of the grid and emitted a continuous sinusoidal pressure waveform at 3 MHz (corresponding to 5 spatial points per wavelength or PPW). The source was then rotated around the centre of the grid and the acoustic field was computed for each rotation angle. Both on- and off-grid sources were used. For the off-grid approach, source grid weights were calculated using an integration point spacing half that of the grid spacing, and were based on the exact, untruncated band limited delta functions.

Figure 6 depicts the time-harmonic amplitude of the field generated by the on- and off-grid line sources at an angle of 30°. The on-grid source produces considerable staircasing errors. Of particular note is the irregular interference pattern in the near-field. This is caused by phase errors resulting from forcing the integration points to lie on nearby grid points. In contrast, the field generated with an off-grid source shows no evidence of staircasing errors and appears symmetric about the beam axis.

Figure 7 depicts the amplitude and phase at the centre of the grid for each angle of rotation of the source. Significant phase errors are evident for the on-grid source, and an angular dependence is seen in the amplitude due to the larger spacing between on-grid integration points as the source becomes diagonal to the grid. In contrast, the amplitude and phase for the off-grid source remains constant regardless of the orientation of the source relative to the grid.

As described in Sec. II C, the integration points used to discretise the band-limiting convolution need not be chosen with the same spacing as the grid points. To investigate how the integration point density affects accuracy, a circular piston was simulated in three-dimensions. The piston diameter was 20 mm and the driving waveform was a 1 MHz sinusoid. The sound speed was 1500 m s⁻¹. The wavefield was computed using the AFP to a distance of 50 mm, at spatial resolutions of 3, 5, and 7 PPW. The source grid weights were computed using exact, untruncated band limited delta functions centred on integration points distributed like those of the disk source in Fig. 5. The source was positioned such that it aligned with the grid along the axial direction. An analytical reference solution for the axial pressure given by Pierce¹³ was used to compute the errors.

Figure 8 depicts the relative L^∞ error in the axial pressure amplitude with a varying integration point density for the off-grid sources. The error can be seen to converge algebraically as the integration point density increases, with less than 2% error achieved for all three PPWs at an upsampling rate of approximately M_integration/M_grid = 4. This corresponds to an integration point spacing that is half that of the grid spacing, since the source is a two-dimensional surface.

To demonstrate the convergence of the band limited source distribution on the true source distribution as the grid resolution increases, a focused bowl source was simulated. This source geometry is especially prone to staircasing errors, as it is impossible to align any portion of it with an orthogonal grid.⁶ The bowl had an aperture diameter of 20 mm, a radius of curvature of 20 mm, and was driven by a 1 MHz sinusoid. The sound speed was set to 1500 m s⁻¹. The wavefield was computed with the AFP to a distance of 47 mm, using a varying number of PPW for the grid spacing. The off-grid sources used integration points that were generated with a spiral phyllotaxis pattern (this can produce uniform samples covering any surface of revolution). This pattern is depicted in Fig. 9, along with illustrative source grid weight slices. A reference solution given by O'Neil¹⁴ was used for the axial pressure (ignoring the first two wavelengths, as these fall behind the bowl's aperture plane where the reference solution is inaccurate). This reference is valid when the transducer diameter is large compared to both the transducer height and acoustic wavelength, as is the case here.

Figure 10 depicts the convergence of the relative L^∞ error in the axial pressure amplitude with an increasing number of PPW for three different source discretisations. The on-grid source converges slowly and produces considerably higher errors than either off-grid source. Indeed, it is known that this error will not converge to zero, because diagonally-aligned portions of the source will be undersampled and hence produce lower amplitudes than they should.⁶ The off-grid sources used exact, untruncated band limited delta functions. They differed in their density of integration points relative to the grid points, with one being undersampled (0.25×) and one being upsampled (4×). The errors resulting from an undersampled off-grid source are considerably worse than those produced by an upsampled off-grid source, but are nonetheless much better than those resulting from an on-grid source. The errors resulting from both off-grid sources also decrease steadily as the PPW increases (for a fixed integration point density), unlike those resulting from an on-grid source. The upsampled off-grid source converges much more quickly than the others, dropping below 0.3% relative error with only 3 PPW.

Figure 11(a) shows the axial pressure amplitude for the on-grid and undersampled off-grid sources at 3 PPW. For the on-grid source, the amplitude is substantially underestimated at the focus due to the undersampling of portions of the source which are diagonally-aligned with the grid.⁶ In the near-field, some of the local pressure maxima are also overestimated, likely due to phase errors, and there is a misalignment of the zero-amplitude points that occur due to destructive interference. In contrast, the undersampled off-grid source produces pressures which are visually indistinguishable from the reference solution in the far-field, and produces very small errors in the near-field. To show this in more detail, Fig. 11(b) depicts the relative error for the undersampled and upsampled off-grid sources. It can be seen that the errors for both off-grid sources are greatest in the near-field, and that the pressure resulting from the upsampled off-grid source oscillates about the reference solution. The error arising from the undersampled off-grid source also oscillates but with an offset from zero, indicating a misalignment of the two solutions.

To examine the use of the sinc approximation to the band limited delta functions given in Eq. (15), and in particular to determine an appropriate truncation threshold ϵ, the focused bowl simulations above were repeated using a number of sinc-based off-grid sources. The error convergence for these simulations is shown in Fig. 12. The exact band limited delta functions and untruncated sinc approximation can be seen to produce nearly identical levels of error. This is expected, as the domain size is large enough that the sinc approximation is accurate (all dimensions had more than 100 grid points). With a truncation threshold, the error resulting from sinc approximations is considerably higher than produced by exact band limited delta functions. However, at ϵ = 1% the error drops below 1% by 3 PPW and reaches 0.2% by 7 PPW. To give a sense of the difference in compute times between these source discretisations, the time taken to compute source grid weights was recorded. At 2 PPW, it took approximately the same time to generate a source using an untruncated sinc approximation as it did using the exact band limited delta function. However, with a sinc truncation threshold of ϵ = 1% the source generation was around 9× faster, and with ϵ = 10% it was around 130× faster. At 7 PPW these advantages become 326× and 4760×, respectively, demonstrating that the computational benefit improves with domain size, as expected.

Thus far, all of the examples in this paper have included time-harmonic source terms. However, the discretisation procedure described here can be used to band-limit any spatial distribution. Hence, it can also be applied to initial value problems such as those arising in photoacoustics. Indeed, an initial value condition is equivalent to a source term of the form $S(x,t)=s(x)(∂/∂t)δ(t)$⁠. To demonstrate this, a disk-shaped initial pressure distribution was simulated in two dimensions using k-Wave. The pressure over the disk was 1 Pa, and the disk radius was 0.4 mm, or 4 times the grid spacing. A sensor was placed 3.2 mm from the disk to record the waveform that was generated. Both on- and off-grid discretisations of the initial condition were used, with the off-grid initial condition using the sinc approximation to the band limited delta functions with a truncation threshold of ϵ = 1%. A high-resolution reference solution was generated using an on-grid initial condition with the grid spacing reduced by a factor of 64×. Note that the initial condition in the reference simulation was filtered in the wavenumber domain such that only the wavenumbers which were supported in the low-resolution simulations were present. This procedure ensures that no higher-frequency waves propagate.

Figure 13 depicts the recorded waveforms for all three initial conditions. The off-grid initial condition can be seen to produce a waveform which is nearly identical to the reference. In contrast, there are significant deviations in the waveform generated by the on-grid initial condition. The oscillations that are present in all three waveforms correspond to the band limited nature of the simulation, as expected.

The band-limiting operation introduced in Sec. II to generate source grid weights can equally be applied to modeling acoustic sensors. In this case, the source (now sensor) grid weights s_l should be used as quadrature weights. For example, if the sensor is a two-dimensional surface with area A, then the average pressure over the sensor is given by

where p_l are the pressure values on the grid. Note that the sensor distribution s(x) in Eq. (5) encodes the sensitivity of the sensor, and can be used to model spatially varying sensitivities, or to convert between units such as pressure and voltage. If a simple average of the field variable is desired, then s(x) should be made equal to one over its region of support.

To illustrate the elimination of staircasing errors using off-grid sensors, a 5 mm line sensor was simulated in two dimensions using the AFP. This sensor was placed in the path of a 1.2 MHz time-harmonic plane wave, and the average pressure over the sensor was computed for a number of orientation angles. The sound speed was 1500 m s⁻¹, making the sensor length 4× the source wavelength. This experiment measures the sensor's directivity, or directional sensitivity, and an exact reference solution is given by Blackstock.¹⁵ Both on- and off-grid sensors were used. For the off-grid sensor, grid weights were calculated using a uniform integration point spacing that was half that of the grid spacing, and were based on the exact, untruncated band limited delta functions. The centre of the sensor was placed on a grid point, so that the on-grid sensor was aligned with the grid whenever the orientation angle was a multiple of π/2. The domain was discretised at 8.3 PPW so that the sensor length was an odd multiple of the grid spacing.

Figure 14 depicts the pressure recorded at multiple orientation angles, normalised relative to values recorded with a point sensor. The directivity of the on-grid sensor exhibits considerable staircasing errors, as multiple orientations produce the same sensor grid weights. In contrast, the directivity of the off-grid sensor shows no evidence of staircasing errors, and matches the reference solution.

An alternative viewpoint on the source discretisation procedure described in Sec. II can be had by reinterpreting the discrete source convolution given in Eq. (5) in terms of the Huygens–Fresnel principle. Instead of considering each point $ξi$ to be an integration point, it can be considered as the location of a point source emitting the desired waveform. The task of source discretisation is then to define a weighted collection of off-grid point sources that cover the true source. This interpretation also gives some insight into the required number of integration points since, as the number of point sources increases, the Huygens–Fresnel principle becomes better satisfied. It may also explain why errors are greatest in the near field: when a finite number of integration points are used, it will take some distance before spherical spreading causes the individual point source wavefronts to merge.

One aspect of the proposed source discretisation method that warrants discussion is its computational expense. Specifically, off-grid source discretisations have a wider region of support than on-grid discretisations, which means memory requirements may be increased. While it has already been said that the use of a truncated sinc approximation to the band limited delta functions provides substantial computational savings by reducing both the number of evaluation points and the subsequent number of source grid weights, there are additional implementations that may complement this in certain circumstances. The most obvious is that source grid weights need not be computed in advance, and instead can be added to the pressure field on-the-fly (and in parallel if the computer hardware supports this). In addition, a middle-ground between on-the-fly computation and full precomputation is possible. Each band limited delta function can be decomposed into a separable product across the spatial dimensions. This means the corresponding source grid weights can be stored as a set of vectors, one for each dimension. For N_p integration points in a three-dimensional simulation, the overall memory required is then O(N_p(N_x + N_y + N_z)), which will typically be less than the O(N_xN_yN_z) required for a fully precomputed set of source grid weights. To reconstruct the full set of source grid weights from the set of vectors, a tensor product must be performed for each integration point, followed by a sum over all of the integration points. This approach also allows one-off, precomputation of expensive trigonometric operations (to compute the vector set for each integration point), and subsequent element-wise elementary arithmetic at each time-step to reassemble the full set of source grid weights.

The derivation given in Sec. II only considered individual source distributions. However, the techniques described straightforwardly extend to the case of multiple sources (each with their own waveform) due to linearity

Thus, each source distribution can be individually discretised, and the resulting source terms can be summed at each time-step in the simulation.

In addition, consideration has only been given to acoustic sources which are fixed in space, and which can be separated into a product S(x,t) = s(x)f(t). However, moving sources can also be accommodated in time-stepping models by considering the moving source as a series of stationary off-grid sources separately defined at each time-step.

The examples given in Sec. III all relate to pressure sources in a homogeneous and lossless medium. However, the proposed method for defining the source grid weights is equally applicable to more complex wave equations, to particle-velocity sources (and the staggered spatial grids they are often implemented on), to problems beyond ultrasound, and indeed any problem to which a Fourier collocation method is applied. In addition, this paper has also only considered sources for which s(x) is constant over its support, but the discrete convolution in Eq. (5) allows for a source distribution to take an arbitrarily complex form. One such example would be an apodised ultrasound transducer.

A method for band-limiting arbitrary source distributions has been derived. The process is based on a discrete convolution between the source distribution and a band limited delta function. This allows for the accurate discretisation of sources with regions of support that do not conform to the equispaced, orthogonal grids that are used with Fourier collocation methods. When applied to a range of source geometries, simulated acoustic fields converge much more quickly than those resulting from staircased source discretisations. The technique was also applied to initial value problems and acoustic sensors, with similarly good results observed for both. A number of codes that implement the ideas in this paper will be released as part of the open-source k-Wave toolbox.¹² 

The authors would like to thank Dr. James L. B. Robertson for useful discussions. This work was supported in part by the Engineering and Physical Sciences Research Council, UK, Grant Nos. EP/L020262/1 and EP/P008860/1, and in part by the European Union's Horizon 2020 research and innovation programme H2020 ICT 2016-2017 under Grant Agreement No. 732411 (as an initiative of the Photonics Public Private Partnership).

Let $b′(x;ξ)$ be a band limited delta function that is derived assuming the number of grid points N is even. Without conjugate symmetry, this is defined in Fourier-space by

for which the Fourier series sums to

To ensure that grid samples of this expression are real-valued, an additional term f(x;ξ) is required that provides conjugate symmetry in the Fourier-domain, i.e.,

where $b=b′+f$⁠. Comparing the Nyquist terms for $b̂′$ and $b̂$⁠, the additional term can be seen to be

for which

Starting with the real components, adding f to $b′$ yields

Now, the imaginary component of $b′$ is

which can be expanded using the trigonometric product-to-sum identities

and rearranged to yield

Noting that $cos(−k−nx)=cos(k−nx)$⁠, the second term in this expression can be seen to be negated by the imaginary component of f

and hence it is clear that

Finally, combining the real and imaginary components and substituting the wavenumbers in Eq. (6) yields

The last two terms in this expression are zero at all grid points regardless of the shift ξ. For a shift that is a multiple of the grid node spacing Δx, they are also zero between the grid points. In this case, the expression matches that derived by Trefethen¹⁶ using a modified Fourier series that treats wavenumbers symmetrically.