Computation of the Fourier transform for a continuous integrable function via NFFT

1 Introduction

The non-equispaced fast Fourier transform (NFFT) is a well-established tool for the efficient evaluation of trigonometric polynomials at arbitrary nodes [2–7]. Traditionally, the NFFT is employed once the Fourier coefficients of a trigonometric polynomial are already known. In many applications, however, the primary quantity of interest is not the trigonometric polynomial itself but the continuous Fourier transform

f^(v)=∫ℝf(x)e-2πixvdx, v∈ℝ,(1)

This naturally raises the question whether the NFFT can also be used as a numerical method for computing values of the Fourier transform.

The present paper provides a positive answer to this question. Building upon ideas developed by Ehler et al. [1], we interpret trigonometric sampling polynomials as approximations of the Fourier transform and derive explicit uniform error estimates. These results show that the NFFT can be employed not only for the evaluation of trigonometric polynomials but also for the efficient computation of values of the continuous Fourier transform. Consequently, the NFFT becomes more than a fast evaluation scheme for trigonometric polynomials. Combined with the approximation results developed in this paper, it yields a practical and theoretically justified method for the computation of continuous Fourier transforms.

More precisely, we consider the computation of the Fourier transform of f^, see Equation 1, under certain decay conditions on f and f^. We restrict f to the spatial domain [−L/2, L/2] and f^ to the frequency domain [−P/2, P/2], where L∈2ℕ and P∈2ℕ are sufficiently large positive integers. We sample f at the equispaced points j/P, j∈[LP], with spatial step size 1/P, where

[LP]:={j∈ℤ:-LP2<j≤LP2}

denotes the centered index set. The P-periodic trigonometric sampling polynomial of f with degree LP/2 is defined by

sP,L(v):=1P∑n∈[LP]f(nP)e-2πinv/P, v∈ℝ,(2)

which approximates f^(v) on the compact frequency interval [−P/2, P/2]. We determine approximate values of f^(kL), k∈[LP], by evaluating the trigonometric sampling polynomial s_{*P, L*}(v) at v = k/L, giving

sP,L(kL)=1P∑j∈[LP]f(jP)e-2πijk/(LP), k∈[LP].(3)

This standard approach to computing the Fourier transform is well-documented in the literature (see, e.g., [8, pp. 17–23], [5, pp. 507–508], or [1]). Note that the values (Equation 3) can be computed efficiently via the fast Fourier transform (FFT), provided L and P are products of powers of the small primes (see [9]). More generally, the trigonometric sampling polynomial s_{*P, L*}(v) can be evaluated at arbitrary points v with |v| ≤ P/2 by means of the non-equispaced fast Fourier transform (NFFT). Recall that the NFFT is a fast approximate algorithm to evaluate a 1-periodic trigonometric polynomial

∑k∈[M]cke2πikx, x∈ℝ,(4)

with given coefficients c_*k*∈ℂ, k∈[M], at given non-equispaced points x_*j*∈[0, 1], j = 1, …, N, where M∈2ℕ. The fundamental idea of the NFFT is to exploit the FFT to efficiently switch between the spatial and frequency domains, thereby enabling a computationally efficient approximation at arbitrary non-uniform points. Throughout this paper, we assume that L, P∈2ℕ are sufficiently large.

Ehler et al. [1] established error estimates for this computation of the values f^(kL) for k∈[LP]. Specifically, they studied how the scaled Euclidean approximation error

EP,L*(f):=1L(∑k∈[LP]|f^(kL)-sP,L(kL)|2)1/2

depends on P and L, given known decay behavior of f∈C(ℝ) and f^∈C(ℝ). These estimates are formulated in terms of norms of Wiener amalgam spaces [10].

In this paper, we address the analogous question for the scaled maximum approximation error

MP,L(f):=1Lmax|v|≤P/2|f^(v)-sP,L(v)|,(5)

again under decay conditions on f and f^. In particular, we study the discretized version

MP,L*(f):=1Lmaxk∈[LP]|f^(kL)-sP,L(kL)|.(6)

Obviously, it holds MP,L*(f)≤MP,L(f). Therefore, we are interested in practicable, explicit upper bounds of M_{*P, L*}(f). We discuss the following questions:

• Under which decay assumptions on f and f^ can the Fourier transform f^ be uniformly approximated by the trigonometric sampling polynomial s_{*P, L*} with explicitly controllable error?
• How can the NFFT be employed to compute values of the continuous Fourier transform f^(v) efficiently and with rigorous ℓ_∞ error bounds?

We introduce the so-called decay rate δ(f, L) of f∈W(C(ℝ), ℓ¹(ℤ)) with respect to step size L∈2ℕ, which is defined by Equation 8. If f∈C(ℝ) has polynomial or exponential decay (see Example 3.2), then the corresponding decay rate δ(f, L) can be easily estimated (see Lemma 3.9 and Lemma 3.15). Applying the Poisson summation Equation 10, we show in Theorem 3.5 that an upper bound of M_{*P, L*}(f) is a linear combination of the decay rates δ(f^,P) and δ(f, L), if f, f^∈W(C(ℝ),ℓ1(ℤ)) and P, L∈2ℕ are given.

Theorem 3.5 shows that the Fourier transform f^ can be approximated uniformly by the trigonometric sampling polynomial s_{*P, L*} with an error controlled by the decay rates of f and f^. Since s_{*P, L*} is a trigonometric polynomial, its values can be computed efficiently by FFT or NFFT techniques. Consequently, the theorem provides a rigorous foundation for the computation of continuous Fourier transforms by means of the NFFT.

The main contributions of this paper can be summarized as follows:

• (C1) We introduce the decay rate δ(f, L) in Equation 8 and prove in Theorem 3.5 the uniform error bound MP,L(f)≤L-1/2δ(f^,P)+L1/2δ(f,L), valid on the whole interval [−P/2, P/2] and not only on the grid {k/L}. This is the property which is needed if s_{*P, L*} is evaluated at arbitrary nodes by the NFFT.
• (C2) All constants in the resulting estimates (Equations 18, 21, and22) are explicit and are computed directly from the decay parameters of f and f^. Moreover, the exponential case is treated for all α, β>0, whereas Ehler et al. [1] is restricted to 0 < α, β ≤ 1.
• (C3) We prove in Theorem 3.20 a matching lower bound, which shows that the rate L^1/2−*a* in the truncation parameter cannot be improved, and we determine in Remark 3.13 the exact range [1,LP] of the quotient EP,L(f)/MP,L(f).
• (C4) We state the resulting method explicitly as Algorithm 1, including its arithmetic cost and its a priori error bound, and we compare direct evaluation, FFT, and NFFT numerically in Section 5.

Algorithm 1 Computation of the Fourier transform via NFFT. Algorithm description for function approximation using NFFT, presented in a fixed-width font. Steps include input requirements, bandwidth and truncation parameter choices, sampling procedure, NFFT parameter selection, evaluation, output guarantees, and computational cost analysis.

This paper is organized as follows: Section 2 reviews the relevant literature and states precisely which of our results are new. In Section 3, we recall the Wiener amalgam space W(C(ℝ), ℓ¹(ℤ)), introduce the decay rate δ(f, L), and derive uniform error estimates for the scaled maximum approximation error M_{*P, L*}(f) under polynomial, exponential, and mixed decay conditions on f and f^, together with a matching lower bound and the resulting Algorithm 1. Section 4 demonstrates how the theoretical error bounds translate into practical parameter choices for FFT- and NFFT-based approximations of the Fourier transform. In particular, the numerical experiments confirm that the NFFT yields highly accurate evaluations of f^(v) throughout the frequency interval while preserving the expected computational efficiency. Section 5 is devoted to the computational cost and compares direct evaluation, FFT, and NFFT. Section 6 concludes the paper.

2 Literature review

In this section, we relate the present paper to the state of the art. We group the discussion into three strands—the approximation of the Fourier transform by discrete sums, the function-space framework, and fast algorithms—and close with a precise account of what is new here.

2.1 Approximation of the Fourier transform by discrete sums

Replacing the Fourier integral (Equation 1) by a finite sum of the form (Equation 3) is the oldest and most widely used numerical device in Fourier analysis; classical accounts are given in Briggs and Henson [8, pp. 17–23] and Plonka et al. [5, pp. 507–8]. The two error sources—aliasing, caused by the discretization of the integral with step size 1/P, and truncation, caused by restricting the sum to the window [−L/2, L/2]—have long been understood qualitatively. Quantitative, non-asymptotic statements have been obtained only recently.

The paper of Ehler et al. [1] is the direct starting point of the present work. The authors measure the accuracy of Equation 3 by the scaled Euclidean error EP,L*(f) and prove, for f and f^ in Wiener amalgam spaces, bounds of the form C(L^−α+P^−β) with α<a-12 and β<b-12 in the case of polynomial decay, and of the form C(e^−*r*′(L/2)^α+e^−*s*′(P/2)^β) with r′ < r, s′ < s and 0 < α, β ≤ 1 in the case of exponential decay. Their results are asymptotically sharp in the exponents, and they give an asymptotically optimal recipe for balancing the number of samples, the sampling step and the grid size. The loss in the exponents is intrinsic to their argument.

Closely related questions arise in the numerical realization of Shannon's sampling theorem, where a band-limited function is reconstructed from finitely many samples. Sharp bounds for the norm of the Shannon sampling operator and regularized sampling formulas with sinh-type or continuous Kaiser–Bessel windows are studied in Kircheis et al. [11, 12]; see also the earlier work [13, 14]. These papers are concerned with the reconstruction of f from its samples, i.e., with the situation of Remark 3.6, whereas we are interested in the forward problem, the evaluation of f^.

2.2 The function-space framework

The natural setting for the Poisson summation formula (Equations 9, 10) is the Wiener amalgam space W(C(ℝ), ℓ¹(ℤ)) introduced by Wiener [15, p. 73] and developed systematically by Feichtinger [10]; see also [16, pp. 103–5]. The example of Katznelson [17] shows that the space C(ℝ)∩L¹(ℝ) is genuinely too large, and Gröchenig [18] clarified the role of the underlying uncertainty principle. We adopt this framework unchanged from Ehler et al. [1], since it is exactly the setting in which the two Poisson summation formulae hold pointwise and with absolute and uniform convergence, which is what a uniform error bound requires.

2.3 Fast algorithms

The NFFT goes back to Dutt and Rokhlin [2], Beylkin [3] and Steidl [4]; a survey is Greengard and Lee [6] and a textbook treatment is Plonka et al. [5, Chapter 7]. Widely used software libraries are NFFT3 [19, 20] and FINUFFT [7, 21]; the experiments in Section 5 use the Python interface pyNFFT3 [22] of NFFT3.

The accuracy of the NFFT itself—i.e., the error committed when the trigonometric polynomial (Equation 4) with given coefficients is evaluated approximately—is by now well-understood. Uniform error estimates for compactly supported window functions are given in Potts and Tasche [23], continuous window functions are analyzed in Potts and Tasche [24], and the case of non-equispaced data in both domains, together with fast sinc transforms, is treated in Kircheis et al. [25]. The FINUFFT kernel is analyzed in Barnett [26]. We stress that this line of work is complementary to ours: it bounds the difference between s_{*P, L*} and its NFFT approximation, whereas we bound the difference between f^ and s_{*P, L*}. Only the combination of both, as in Algorithm 1 and Remark 3.22, yields a rigorous bound for the quantity a user actually computes.

2.4 What is new in this paper

Since the relation to Ehler et al. [1] is close, we state explicitly which of our results are new, which cannot be obtained from Ehler et al. [1], and which assumptions have been weakened. A synopsis is given in Table 1.

• (N1) A different, strictly finer error measure. We estimate the maximum error M_{*P, L*}(f) over the whole interval [−P/2, P/2], while Ehler et al. [1] estimates the Euclidean error EP,L*(f) over the finite grid {k/L:k∈[LP]}. A maximum error does occur in Ehler et al. [1, Section 2.7, Equation 12], but only through the trivial estimate MP,L(f)≤EP,L(f), i.e., the discrete maximum error is bounded there by the Euclidean error and therefore inherits its rate. This is exactly the step that is lossy: the two quantities satisfy

MP,L(f)≤EP,L(f)≤LPMP,L*(f),

Table 1. Comparison of Ehler et al. [1] with the present paper.

Here a, b>1 denote the polynomial and r, α, s, β>0 the exponential decay parameters of f and f^; ε>0 is arbitrarily small.

and both inequalities are attained up to constants, see Remark 3.13 and Figure 1b. The two error measures are therefore not equivalent: their quotient is unbounded as LP → ∞. Consequently, Theorem 3.10 cannot be deduced from Ehler et al. [1, Theorem 2.1] via [1, (12)]; that route yields only the exponent P^{−*b*+1/2+ε} instead of P^−*b*, and the gap is clearly visible in Figure 1b. Conversely, an estimate on a finite grid carries no information about the values of f^-sP,L between the grid points, which is precisely the information needed when s_{*P, L*} is evaluated by an NFFT at arbitrary nodes.

• (N2) Explicit constants. All constants in Equations 18, 21, and22 are given in closed form in terms of the decay parameters a, b, c, d, r, s, α, β. This makes the estimates usable as an a priori stopping criterion, i.e., for the actual choice of L and P (see Algorithm 1).
• (N3) No loss in the exponents, and a weaker assumption in the exponential case. In the polynomial case we obtain the exponent b in P and the exponent a-12 in L, whereas Ehler et al. [1] must give away an arbitrarily small ε>0 in both exponents. In the exponential case, [1, Theorems 2.3 and 2.4] assume 0 < α, β ≤ 1 and lose an arbitrarily small amount in the decay parameters r and s; Theorem 3.16 and Theorem 3.18 hold for all α, β>0 and retain the original r and s. In particular, the Gaussian (α = β = 2) of Example 3.24 is not covered by Ehler et al. [1].
• (N4) Sharpness. Theorem 3.20 provides a lower bound which matches the upper bound of Theorem 3.10 in the truncation parameter L up to a constant factor. In Ehler et al. [1] the sharpness of the estimates is conjectured and confirmed by numerical experiments, but no lower bound is proved there.
• (N5) Algorithm and computational assessment. Algorithm 1 combines the estimate of Theorem 3.5 with the NFFT and states the total error and the arithmetic cost; Section 5 compares direct evaluation, FFT, and NFFT.

The general strategy of splitting the Poisson summation formula into a principal part and two remainders, and the use of Wiener amalgam spaces, are taken from Ehler et al. [1], as is the scaling L^−1/2 in the definition (Equation 5) of the error, which we have kept in order to make the two sets of results directly comparable.

3 Uniform approximation by trigonometric sampling polynomials

The Wiener amalgam space W(C(ℝ), ℓ¹(ℤ)) is the Banach space consisting of all functions f:ℝ → ℂ which are locally in C(ℝ) and have globally an ℓ¹(ℤ) behavior at infinity in the sense that the uniform norms of f over the intervals [k, k+1], k∈ℤ, form a sequence in ℓ¹(ℤ). The norm of W(C(ℝ), ℓ¹(ℤ)) is given by

The Wiener amalgam space was introduced by Wiener [15, p. 73] and was used as a convenient function space in Fourier analysis, cf. [10], [16, pp. 103–5], and [1]. The condition f∈W(C(ℝ), ℓ¹(ℤ)) is strong enough to exclude many pathological functions which play a role in Fourier analysis, but are of little practical interest. The Wiener amalgam space W(C(ℝ), ℓ¹(ℤ)) contains all continuous, compactly supported functions and is therefore a dense subspace of the Lebesgue space L¹(ℝ). This Wiener amalgam space has the following property:

Lemma 3.1. (see[16, Lemma 6.1.2]) If f∈W(C(ℝ), ℓ¹(ℤ)) and R>0 are given, then

supx∈ℝ∑k∈ℤ|f(x+Rk)|≤(1R+1)||f||W(C(ℝ),ℓ1(ℤ)).

From Lemma 3.1 with R=1L, L∈2ℕ, it follows immediately that

∑k∈ℤ|f(kL)|≤(1R+1)‖f‖W(C(ℝ),l1(ℤ))<∞(7)

such that (f(kL))k∈ℤ∈ℓ1(ℤ) for each L∈2ℕ.

For each f∈W(C(ℝ), ℓ¹(ℤ)) we introduce the decay rate of f with respect to the step size L∈2ℕ by

δ(f,L):=max|x|≤L/2∑k∈ℤ\{0}|f(x+kL)|.(8)

From Lemma 3.1 it follows immediately that δ(f, L) < ∞ for arbitrary L∈2ℕ.

Example 3.2. A function f belongs to the Wiener amalgam space W(C(ℝ), ℓ¹(ℤ)) if and only if f∈C(ℝ) can be majorized by the step function s with jumps at integers, given by s(x) = m_*k* for x∈[k, k+1), k∈ℤ, with mk:=maxx∈[k,k+1]|f(x)| and ∑k∈ℤmk<∞. Hence it holds

W(C(ℝ),ℓ1(ℤ))⊂C(ℝ)∩L1(ℝ).

In particular, each compactly supported function f∈C(ℝ) belongs to W(C(ℝ), ℓ¹(ℤ)). If f∈C(ℝ) has polynomial decay, i.e.,

supx∈ℝ|f(x)|(1+|x|)a≤c<∞, a>1,

or exponential decay, i.e.,

supx∈ℝ|f(x)|er|x|α≤d<∞, r,α>0,

then f belongs to W(C(ℝ), ℓ¹(ℤ)). For corresponding proofs, see Lemma 3.9 and Lemma 3.15.

3.1 Poisson summation formula

The key tool for estimating the uniform approximation error is the Poisson summation formula. It allows us to express the error between f^ and the trigonometric sampling polynomial s_{*P, L*} in terms of the tails of f and f^, which can then be controlled by decay conditions on these functions.

Lemma 3.3. see[18] Let P∈2ℕ be fixed. For each f∈W(C(ℝ), ℓ¹(ℤ)) with f^∈W(C(ℝ),ℓ1(ℤ)), the Poisson summation formulae

∑k∈ℤf(x+kP)=1P∑n∈ℤf^(nP)e2πinx/P,(9)

∑k∈ℤf^(v+kP)=1P∑n∈ℤf(nP)e-2πinv/P(10)

hold for all x, v∈ℝ. All series converge absolutely and uniformly on ℝ.

Proof . From Lemma 3.1 it follows that

fP(x):=∑k∈ℤf(x+kP), x∈ℝ,

is a continuous, P-periodic function. By the assumption f^∈W(C(ℝ),ℓ1(ℤ)), it follows from Equation 7 that

∑k∈ℤ|f^(nP)|<∞,

so that the P-periodic Fourier series 1P∑n∈ℤf^(nP)e2πinx/P converges absolutely and uniformly on ℝ. Both functions coincide, since f_*P* has the P-periodic Fourier coefficients

1P∫0PfP(x)e-2πinx/Pdx=1P∫ℝf(x)e-2πinx/Pdx=1Pf^(nP), n∈ℤ.

This completes the proof of Equation 9. Replacing f by f^ yields Equation 10.

Remark 3.4. By an example, Katznelson [17] has shown that the Poisson summation formula (Equation 9) does not hold for all f∈C(ℝ)∩L¹(ℝ) with f^∈C(ℝ)∩L1(ℝ), illustrating the essential difference between the spaces W(C(ℝ), ℓ¹(ℤ)) and C(ℝ)∩L¹(ℝ).

3.2 The trigonometric sampling polynomial and its fast evaluation

The trigonometric sampling polynomial provides a natural approximation of the Fourier transform f^ on a compact frequency interval, and can be evaluated efficiently at arbitrary points using the NFFT. Now we make this precise.

For f∈W(C(ℝ), ℓ¹(ℤ)) with f^∈W(C(ℝ),ℓ1(ℤ)) and sufficiently large L, P∈2ℕ, we define the P-periodic trigonometric sampling polynomial of f with sampling support [−L/2, L/2] and degree LP/2 by s_{*P, L*}, see Equation 2. Under suitable decay conditions on f and f^, the polynomial s_{*P, L*} provides a good uniform approximation of f^ on the compact interval [−P/2, P/2], see Theorem 3.5 and the following results.

The values s_{*P, L*}(k/L), k∈[LP], on the uniform frequency grid can be computed via the FFT in O(LPlogLP) operations, provided L and P are products of powers of small primes (see [9]). More generally, the trigonometric sampling polynomial s_{*P, L*} can be evaluated at K arbitrary points v₁, …, v_*K*∈[−P/2, P/2] by means of the non-equispaced fast Fourier transform (NFFT) in only O(LPlogLP+K) operations (see [5]). Hence, for each v with |v| ≤ P/2, the value s_{*P, L*}(v) is a computationally efficient approximation of f^(v).

3.3 The uniform approximation error

Let L, P∈2ℕ be given. By convenient decay conditions on f and f^, we can show that s_{*P, L*} is a good uniform approximation of f^ on the compact interval [-P2,P2]. To this end, we use the scaled uniform approximation error (Equation 5). If we form the maximum only on the grid {kL:k∈[LP]}, then we obtain the discretized version of the scaled uniform approximation error (Equation 6). We derive explicit estimates for M_{*P, L*}(f) and MP,L*(f) under various decay conditions on f and f^. The key decomposition is provided by Theorem 3.5 below, which splits the error into two contributions: one governed by the decay of f^ and one by the decay of f. Explicit bounds for each contribution are then given in Lemma 3.9 and Lemma 3.15, covering the cases of polynomial and exponential decay, respectively.

Theorem 3.5. Let L, P∈2ℕ be fixed. Let f∈W(C(ℝ), ℓ¹(ℤ)) with f^∈W(C(ℝ),ℓ1(ℤ)) be given.Then f^ can be uniformly approximated on [−P/2, P/2] by the trigonometric sampling polynomial s_{*P, L*}, and it holds the estimate

MP,L*(f)≤MP,L(f)≤1Lδ(f^,P)+Lδ(f,L),(11)

where δ(f^,P) and δ(f, L) denote the corresponding decay rates (Equation 8).

Proof . Splitting both series in the second Poisson summation formula (10) into the principal terms and the remainders, we obtain for all v∈ℝ

f^(v)+σ1(v)=sP,L(v)+σ2(v),

with

σ1(v):=∑k∈ℤ\{0}f^(v+kP),σ2(v):=1P∑n∈ℤ\[LP]f(nP)e-2πinv/P.

Hence, this provides that

MP,L*(f)≤MP,L(f)≤1Lmax|v|≤P/2|σ1(v)|+1Lmax|v|≤P/2|σ2(v)|.

Applying the triangle inequality and Equation 8, we obtain

max|v|≤P/2|σ1(v)|≤max|v|≤P/2∑k∈ℤ\{0}|f^(v+kP)|=δ(f^,P).

In the expression σ₂, the summation indices n∈ℤ\[LP] can be written in the form n = m+kLP with m∈[LP] and k∈ℤ\{0}, such that

σ2(v)=1P∑m∈[LP]∑k∈ℤ\{0}f(mP+kL)e-2πimv/Pe-2πikLv.

Using triangle inequality and Equation 8, this yields

max|v|≤P/2|σ2(v)|≤LPPmaxm∈[LP]∑k∈ℤ∖{0}|f​(mP+kL)|≤Lδ(f,L).

This completes the proof. ■

Remark 3.6. For given f,f^∈W(C(ℝ),ℓ1(ℤ)), we can analogously estimate the maximum approximation error

1Lmax|x|≤P/2|f(x)-1P∑n∈[LP]f^(nP)e2πinx/P|(12)

by applying the first Poisson summation Equation 9. The error (Equation 12) has the upper bound

1Lδ(f,P)+Lδ(f^,L),

which can be further estimated under concrete decay conditions on f and f^ by Lemma 3.9 and Lemma 3.15. For the discretized version of Equation 12 we obtain

1Lmaxℓ∈[LP]|f(ℓL)-1P∑n∈[LP]f^(nP)e2πinℓ/(LP)|≤1Lδ(f,P)+L δ(f^,L).

Remark 3.7. The reconstruction in Remark 3.6 presupposes the exact Fourier samples f^(n/P), n∈[LP]. In practice, however, only the computed trigonometric sampling polynomial s_{*P, L*}, see Equation 2, is available. We therefore sample s_{*P, L*} on the finer equispaced grid {j/M:j∈[MP]} of the frequency interval [−P/2, P/2), where M∈2ℕ with M≥2L denotes a frequency oversampling factor, and form the fully computable approximation

fM,P,L(x):=1M∑j∈[MP]sP,L(jM)e2πijx/M, x∈[-L/2,L/2).(13)

This is the M-point rectangle-rule discretization of the truncated inverse Fourier integral ∫-P/2P/2sP,L(v)e2πivxdv; on the grid x = ℓ/P, ℓ∈[LP], it can be evaluated by a single inverse FFT of length MP.

Introducing the truncated cardinal series

gP,L(x):=∫-P/2P/2sP,L(v)e2πivxdv=∑n∈[LP]f(nP)sinc(Px-n),

sinc(t):=sin(πt)πt,(14)

the approximation error splits into three contributions,

f(x)-fM,P,L(x)=∫|v|>P/2f^(v)e2πivxdv︸band truncation+∫-P/2P/2(f^-sP,L)(v)e2πivxdv︸sampling polynomial+(gP,L(x)-fM,P,L(x))︸frequency discretization.

The first term is bounded by the L¹-tail η(f^,P):=∫|v|>P/2|f^(v)|dv, the second by Pmax|v|≤P/2|f^(v)-sP,L(v)|=PLMP,L(f), see Equation 5. For the third term, a direct evaluation of the geometric sum in Equation 13 gives the closed form

fM,P,L(x)=1MP∑n∈[LP]f(n/P)eiπαn/Msin(πPαn)sin(παn/M)

with α_*n*: = x−n/P, hence

gP,L(x)−fM,P,L(x)=1P∑n∈[LP]f​(nP)sin(πPαn)[1παn

−eiπαn/MMsin(παn/M).]

Since |sin(πPα_*n*)| ≤ 1 and, by the elementary inequality 1−βcotβ ≤ β² for |β| ≤ π/2, the bracket is bounded in modulus by M-11+π2/4 whenever |α_*n*| ≤ M/2 (which holds for |x| ≤ L/2 as M≥2L), we obtain with Lemma 3.1 (for R = 1/P and |x| ≤ L/2)

≤(P+1)1+π2/4MP||f||W(C(ℝ),ℓ1(ℤ)).

Combining the three bounds with Theorem 3.5 and (P+1)/P ≤ 2 yields the uniform error estimate

max|x|≤L/2|f(x)-fM,P,L(x)|≤η(f^,P)+P(δ(f^,P)+Lδ(f,L))

+21+π2/4M||f||W(C(ℝ),ℓ1(ℤ)).(15)

The three terms reflect, respectively, the truncation of the inverse transform to the band [−P/2, P/2], the approximation of f^ by s_{*P, L*} on that band, and the rectangle-rule discretization in the frequency variable.

If f^ has polynomial decay supv|f^(v)|(1+|v|)b≤c<∞ with b>1, then η(f^,P)≤2cb-1(1+P/2)1-b, and with Lemma 3.9 the estimate Equation 15, becomes

max|x|≤L/2|f(x)-fM,P,L(x)|≤2cb-1(1+P2)1-b+2c(2b-1)ζ(b)P1-b+PLδ(f,L)

+21+π2/4M||f||W(C(ℝ),ℓ1(ℤ)).(16)

For compactly supported f, where δ(f, L) = 0 for sufficiently large L, the spatial error therefore decays as O(P1-b)+O(M-1), and balancing the two contributions suggests the choice M≳P^*b*−1. We remark that the leading M⁻¹-coefficient is in fact governed by |sP,L(P/2)|≈|f^(P/2)|, which is already small at the band edge, so that in practice, the frequency discretization is markedly less critical than the truncation with respect to P.

Remark 3.8. Roughly speaking, smoothness of f implies polynomial decay of f^, and smoothness of f^ implies polynomial decay of f. If each derivative f^(ℓ), ℓ = 0, …, b, belongs to L¹(ℝ), then the Fourier transform of f^(ℓ) equals (2πvi)ℓf^(v) and is bounded on ℝ, so that

supv∈ℝ|f^(v)|(1+|v|)b≤c<∞.

Conversely, if each derivative f^(ℓ), ℓ = 0, …, b, belongs to L¹(ℝ), then the inverse Fourier transform of f^(ℓ) equals (−2πix)^ℓf(x) and is bounded on ℝ, so that

supx∈ℝ|f(x)|(1+|x|)b≤c<∞.

Now we estimate the decay rate (Equation 8) of a function f∈C(ℝ) with polynomial decay.

Lemma 3.9. Let L∈2ℕ be given. Assume that f∈C(ℝ) has polynomial decay

supx∈ℝ|f(x)|(1+|x|)a≤c<∞, a>1.

Then the decay rate of f with respect step size L can be estimated by

δ(f,L)≤2c(2a-1)ζ(a)L-a,(17)

where ζ(a) denotes the Riemann zeta function

ζ(a):=∑n=0∞(n+1)-a<1+1a-1, a>1.

Further f belongs to W(C(ℝ), ℓ¹(ℤ)) with the norm

Proof . By the polynomial decay of f, for all k∈ℤ\{0} and |x| ≤ L/2 it holds that

Hence it follows that

δ(f,L)=max|x|≤L/2∑k∈ℤ\{0}|f(x+kL)|≤2c∑k=1∞(1+(k-12)L)-a=2cL-a∑k=1∞(k-12+1L)-a<2cL-a∑n=0∞(n+12)-a=2cζ(a,12)L-a,

where ζ(a,12):=∑n=0∞(n+12)-a denotes the Hurwitz zeta function. By ζ(a,12)=2a∑n=0∞(2n+1)-a we obtain

ζ(a)=∑n=0∞(2n+1)-a+∑n=0∞(2n+2)-a=∑n=0∞(2n+1)-a +2-a∑n=0∞(n+1)-a =2-aζ(a,12)+2-aζ(a)

such that ζ(a,12)=(2a-1)ζ(a), which gives Equation 17. For a∈2ℕ, explicit values of ζ(a) are known, namely ζ(a)=|Ba|2a!(2π)a with the Bernoulli number B_*a* (see [27, pp. 266–7]). The integral test for convergence of series provides that for a>1,

ζ(a)=∑n=0∞(n+1)-a≤1+∫0∞(x+1)-adx=1+1a-1.

Obviously, we have f∈W(C(ℝ), ℓ¹(ℤ)), since it holds by assumption that

‖f‖W(C(ℝ),l1(ℤ))≤c∑n=0∞(1+n)−a+c∑n=1∞(1+(n−1))−a=2cζ(a) <∞.

This completes the proof.

Using the corresponding decay rates of f and f^ which have both polynomial decays, we can estimate the scaled uniform maximum error M_{*P, L*}(f):

Theorem 3.10. Let L, P∈2ℕ be given. If both functions f,f^∈C(ℝ) have polynomial decays

supx∈ℝ|f(x)|(1+|x|)a≤c<∞, a>1,supv∈ℝ|f^(v)|(1+|v|)b≤d<∞, b>1,

then it holds the error estimate

MP,L*(f)≤MP,L(f)≤2d(2b-1)ζ(b)L-1/2P-b

+2c(2a-1)ζ(a)L-a+1/2.(18)

Proof . Using Theorem 3.5, it follows that

MP,L*(f)≤MP,L(f)≤1Lδ(f^,P)+Lδ(f,L).

Then Lemma 3.9 provides that

δ(f^,P)≤2d(2b-1)ζ(b)P-b,δ(f,L)≤2c(2a-1)ζ(a)L-a.

This completes the proof.

Remark 3.11. We compare Theorem 3.10 with the corresponding result of [1, Theorem 2.1]. From [1, Theorem 2.1] it follows that

MP,L(f)≤EP,L(f)≤C(L-α+P-β)

with arbitrary α∈(0,a-12), β∈(0,b-12), and some constant C>0. Our estimate (18) is explicit and more practicable, since all constants are described directly by the assumptions of Theorem 3.10, and it does not lose an arbitrarily small ε>0 in the exponents. Moreover, using Theorem 3.10, the scaled Euclidean approximation error can be bounded by

EP,L(f)≤LPMP,L(f)≤2d(2b-1)ζ(b)P-b+1/2+2c(2a-1)ζ(a)L1-aP1/2.

The estimate (Equation 18) is not symmetric in the two parameters: the decay of f^ enters with the full exponent b, whereas the decay of f enters with the reduced exponent a-12. Since the corresponding bound of [1, Theorem 2.1] is symmetric, we explain in the following two remarks where this difference comes from. It turns out that it is caused neither by the proof technique nor by the normalization, but by a genuine structural difference between the two error contributions σ₁ and σ₂ in the proof of Theorem 3.5.

Remark 3.12. Let f, f^ have polynomial decay with exponents a, b as in Theorem 3.10. Consider the two remainders

σ1(v)=∑k∈ℤ\{0}f^(v+kP), σ2(v)=1P∑n∈ℤ\[LP]f(nP)e-2πinv/P

separately, and compare their maximum norm on [−P/2, P/2] with their quadratic mean (1P∫-P/2P/2|σj(v)|2dv)1/2.

The aliasing remainder σ₁ is a periodization of f^. It is non-oscillatory, and both its maximum and its quadratic mean are of the exact order P^−*b*; the quotient of the two is bounded above and below by absolute constants.

The truncation remainder σ₂ behaves completely differently. By the orthogonality of the exponentials e^{−2π*inv*/*P*} on [−P/2, P/2] we have exactly

1P∫-P/2P/2|σ2(v)|2dv=1P2∑n∈ℤ\[LP]|f(nP)|2=Θ(P-1L1-2a),

so that the quadratic mean of σ₂ is of the order P^−1/2L^1/2−*a*, whereas for non-negative f at v = 0 all summands add up constructively and give

max|v|≤P/2|σ2(v)| ≥ |σ2(0)|=1P∑n∈ℤ\[LP]f(nP)=Θ(L1-a).

Here, the two-sided order estimates presuppose a matching two-sided polynomial decay c0(1+|x|)-a≤|f(x)|≤c(1+|x|)-a; the lower estimate for |σ₂(0)| is exactly the one established in Theorem 3.20. Hence, the maximum of σ₂ exceeds its quadratic mean by the factor LP, while for σ₁ the two agree up to constants.

This is the source of the asymmetry. Passing from an averaged error to a uniform error costs nothing in the aliasing part, but costs the full factor LP in the truncation part, and no symmetric bound for the uniform error can therefore exist. We add that after the normalization by L^−1/2 both error measures predict the same rate in the truncation parameter, namely L^1/2−*a*; the difference between Equation 18 and [1, Theorem 2.1] is confined to the aliasing parameter, where our bound is better by the factor LP. This is confirmed numerically in Example 3.27 and Figure 1b.

Remark 3.13. The two error measures MP,L(f) and EP,L(f) are not equivalent. Indeed, for any f we have the elementary two-sided estimate

1≤EP,L(f)MP,L(f)≤LP,(19)

since a maximum of LP numbers is at most their Euclidean norm and at least (LP)^−1/2 times that norm. Both bounds in Equation 19 are essentially attained:

• If the error is truncation-dominated, i.e., if δ(f^,P)≪Lδ(f,L), then by Remark 3.12 the quotient stays bounded. For the band-limited function g of Example 3.27 one finds numerically EP,L(g)/MP,L(g)→1.244 as L → ∞, independently of P.
• If the error is aliasing-dominated, i.e., if Lδ(f,L)≪δ(f^,P), the quotient grows like LP. For the compactly supported B-spline M₂ with L = 2 one measures EP,L/MP,L=22.18 at P = 1024, while LP=45.25; see Figure 1b.

Consequently, an estimate for EP,L(f) may overestimate MP,L(f) by a factor which grows without bound, and Theorem 3.10 cannot be recovered from [1, Theorem 2.1]. We stress in addition that M_{*P, L*}(f), which controls the error at every v with |v| ≤ P/2, is not comparable to EP,L(f) at all, since the latter contains no information about the values of f^-sP,L off the grid {k/L:k∈[LP*]}.

As is known, a non-trivial function cannot be band-limited and space-limited (see [5, pp. 103–104]). If f is space-limited or band-limited, then we obtain simpler error estimates:

Corollary 3.14. Let L, P∈2ℕ be given. If f,f^∈C(ℝ) satisfy

supx∈ℝ|f(x)|(1+|x|)a≤c<∞, a>1, suppf^⊆[-M2,M2], 0<M≤P,

then it holds the error estimate

MP,L*(f)≤MP,L(f)≤2c(2a-1)ζ(a)L1/2-a.

If f,f^∈C(ℝ) satisfy

suppf⊆[-M2,M2], 0<M≤L,supv∈ℝ|f^(v)|(1+|v|)b≤d<∞, b>1,

then it holds the error estimate

MP,L*(f)≤MP,L(f)≤2d(2b-1)ζ(b)L-1/2P-b.

Proof . This result follows immediately from Theorem 3.10. If suppf^⊆[-M2,M2] and M ≤ P, then f^(v+kP)=0 for all k∈ℤ\{0} and |v| ≤ P/2, so that δ(f^,P)=0. If suppf⊆[-M2,M2] and M ≤ L, then all samples f(nP) vanish for n∈ℤ\[LP], so that δ(f, L) = 0.

Now we estimate the decay rate (Equation 8) of a function with exponential decay.

Lemma 3.15. Let L∈2ℕ be given. Let f∈C(ℝ) be a function with exponential decay

supx∈ℝ|f(x)|er|x|α≤c<∞, r,α>0.

Assume that in the case 0 < α < 1 it holds r(L/2)α>1α-1.Then the decay rate of f with respect to step size L can be estimated by

δ(f,L)≤2cγ(r,α,L)e-r(L/2)α(20)

with

γ(r,α,L):={1+r1/αL(αr(L/2)α−1+α) 0<α<1,1+1rL α=1,1+2(r(L/2)α+bα)1/α−r1/αL4bα α>1.

For α>1, it holds bα:=Γ(1α+1)α/(1-α). Further f belongs to W(C(ℝ), ℓ¹(ℤ)) and has the norm

‖f‖W(C(ℝ),l1(ℤ))≤2c+2cαr−1/αΓ(1α)<∞.

Proof . By the exponential decay of f, for all k∈ℤ\{0} and |x| ≤ L/2 it holds that

Hence, it follows that

δ(f,L)=max|x|≤L/2∑k∈ℤ\{0}|f(x+kL)|≤2c∑k=1∞e-r(L/2)α(2k-1)α =2c∑n=0∞e-r(L/2)α(2n+1)α.

Setting u: = r(L/2)^α for shortness, the integral test for convergence of series yields

∑n=0∞e-u(2n+1)α≤e-u+∫0∞e-u(2x+1)αdx=e-u +12αu-1/αΓ(1α,u),

where

Γ(1α,u):=∫u∞t1/α−1e−tdt

denotes the upper incomplete gamma function. A simple check shows that

−1αu−1/αΓ(1α,tαu)

is a primitive of e^−*u**t*α and therefore

∫1∞e−u tαdt=1αu−1/αΓ(1α,u).

By [28, Theorem 1.1 and Proposition 2.7], an upper bound of the upper incomplete gamma function Γ(1α,u) is given by

Γ(1α,u)≤{u1/αu+1−1/αe−u 0<α<1,(1+12u)e−u α=1,α(u+bα)1/α−u1/αbαe−u α>1.

This implies the estimate (Equation 20). Obviously, we have f∈W(C(ℝ), ℓ¹(ℤ)), since it holds by assumption and [29, Formula 3.478] that

‖f‖W(C(ℝ),l1(ℤ))≤c∑n=0∞e−r nα+c∑n=1∞e−r (n−1)α=2c∑n=0∞e−r nα ≤2c+2c∫0∞e−r tαdt=2c+2cαr−1/αΓ(1α)<∞.

This completes the proof.

Using the decay rates of f and f^, which have both exponential decays, we can estimate the scaled uniform maximum error M_{*P, L*}(f):

Theorem 3.16. Let L, P∈2ℕ be given. Let f,f^∈C(ℝ) be functions with exponential decay

supx∈ℝ|f(x)|er|x|α≤c<∞, r,α>0,supv∈ℝ|f^(v)|es|v|β≤d<∞, s,β>0.

Assume that r(L/2)α>1α-1 in the case 0 < α < 1 and that s(P/2)β>1β-1 in the case 0 < β < 1.Then it holds the error estimate

MP,L(f)≤2dγ(s,β,P)L-1/2e-s(P/2)β+2cγ(r,α,L)L1/2e-r(L/2)α,(21)

where the constants γ(s, β, P) and γ(r, α, L) are defined in Lemma 3.15.

Proof . Using Theorem 3.5, it follows that

MP,L(f)≤1Lδ(f^,P)+Lδ(f,L).

Then Lemma 3.15 provides that

δ(f^,P)≤2dγ(s,β,P)e-s(P/2)β,δ(f,L)≤2cγ(r,α,L)e-r(L/2)α.

This completes the proof.

Remark 3.17. We compare Theorem 3.16 with the corresponding result in [1, Theorem 2.3]. Then from [1, Theorem 2.3] it follows immediately that

MP,L∗(f)≤EP,L∗(f)≤C(e−r′(L/2)α+e−s′(P/2)β)

with arbitrary r′∈(0, r), s′∈(0, s), and certain constant C>0. Note that in [1, Theorem 2.3] it is considered only the case 0 < α, β ≤ 1. Thus, we see that our estimate (Equation 21) is more general, explicit, and more practicable, since all constants in Equation 21 are known from the decay conditions of Theorem 3.16.

Using the decay rates of f and f^, which have mixed decays, we can estimate the scaled uniform maximum error M_{*P, L*}(f):

Theorem 3.18. Let L, P∈2ℕ be given. Let f,f^∈C(ℝ) be functions with mixed decay

supx∈ℝ|f(x)|(1+|x|)a≤c<∞, a>1,supv∈ℝ|f^(v)|es|v|β≤d<∞, s,β>0.

Assume that s(P/2)β>1β-1 in the case 0 < β < 1. Then it holds the error estimate

MP,L(f)≤2dγ(s,β,P)L-1/2e-s(P/2)β+2c(2a-1)ζ(a)L-a+1/2.

Alternatively, let f,f^∈C(ℝ) be functions with mixed decay

supx∈ℝ|f(x)|er|x|α≤c<∞, r,α>0,supv∈ℝ|f^(v)|(1+|v|)b≤d<∞, b>1.

Assume that r(L/2)α>1α-1 in the case 0 < α < 1. Then it holds the error estimate

MP,L(f)≤2d(2b-1)ζ(b)L-1/2P-b+2cγ(r,α,L)L1/2e-r(L/2)α.(22)

Proof . In both cases, Theorem 3.5 yields

MP,L(f)≤1Lδ(f^,P)+Lδ(f,L).

Then the assertion follows immediately from Lemma 3.9 and Lemma 3.15.

In practice, L, P∈2ℕ are sufficiently large, so that the exponential error terms in the estimates of Theorem 3.18 are negligible. Thus M_{*P, L*}(f) is dominated by the corresponding polynomial error term.

Remark 3.19. We compare the second case of Theorem 3.18 with the corresponding result in [1, Theorem 2.4]. Then from [1, Theorem 2.4 and (4.21)] it follows immediately that

MP,L∗(f)≤EP,L∗(f)≤C(e−r′(L/2)α+P−β)

with arbitrary r′∈(0, r), β∈(0,b-12), and certain constant C>0. Note that in [1, Theorem 2.4] it is considered only the case 0 < α ≤ 1. Thus, we see that our estimate (Equation 22) is more general, explicit, and more practicable, since all constants in Equation 22 are known from the decay conditions of Theorem 3.18.

3.4 A lower bound for the uniform approximation error

All estimates obtained so far are upper bounds. We now show that the rate in the truncation parameter L predicted by Theorem 3.10 is sharp. This is in particular relevant for Remark 3.12, where we claimed that the asymmetry of Equation 18 is a genuine effect.

Theorem 3.20. Let L, P∈2ℕ be given, and let f∈W(C(ℝ), ℓ¹(ℤ)) with f^∈W(C(ℝ),ℓ1(ℤ)). Assume that f is non-negative, that f is non-increasing on [x₀, ∞) for some x₀≥0, and that

f(x) ≥ c0(1+x)-a, x≥x0,

with constants c₀>0 and a>1. Then for all L≥2x₀ it holds

MP,L(f) ≥ MP,L*(f) ≥ 1L(c0a-1(1+L2+1P)1-a-δ(f^,P)).(23)

The bound (Equation 23) is informative in the truncation-dominated regime δ(f^,P)<c0a-1(1+L2+1P)1-a. In particular, if f^ is band-limited with suppf^⊆[-M2,M2] and M ≤ P, then δ(f^,P)=0 and

MP,L(f) ≥ c0a-1(1+L2+1P)1-aL-1/2 = Ω(L-a+1/2) (L→∞),

uniformly in P≥M. Together with the upper bound O(L-a+1/2) of Theorem 3.10 and Corollary 3.14, this yields the two-sided estimate MP,L(f)=Θ(L-a+1/2).

Proof . As in the proof of Theorem 3.5 we have f^(v)-sP,L(v)=σ2(v)-σ1(v) for all v∈ℝ. Evaluating at the grid point v = 0, which belongs to {k/L:k∈[LP]}, we obtain

LMP,L*(f) ≥ |f^(0)-sP,L(0)|=|σ2(0)-σ1(0)| ≥ |σ2(0)|-|σ1(0)|.

By the proof of Theorem 3.5 we have |σ1(0)|≤δ(f^,P). Since f≥0, the value

σ2(0)=1P∑n∈ℤ\[LP]f(nP)

is non-negative, so that |σ₂(0)| = σ₂(0), and discarding the indices n≤-LP2 only decreases the sum,

σ2(0) ≥ 1P∑n>LP/2f(nP)=1P∑n≥LP/2+1f(nP).

Since L≥2x₀, the function f is non-increasing on [L2,∞), and the right-hand side is a lower Riemann sum with step size 1P of the integral of f over [L2+1P,∞). Hence

σ2(0) ≥ ∫L/2+1/P∞f(x)dx ≥ c0∫L/2+1/P∞(1+x)-adx =c0a-1(1+L2+1P)1-a.

Together with MP,L*(f)≤MP,L(f) this gives Equation 23.

Remark 3.21. Theorem 3.20 shows that the exponent a-12 of L in Equation 18 is optimal, and that the maximum of |f^-sP,L| is, in the truncation-dominated regime, attained near v = 0 and is of the order of the L¹-tail ∫|x|>L/2|f(x)|dx of f. The assumptions of Theorem 3.20 are satisfied, for example, by the band-limited function g of Example 3.27 and, more generally, by every non-negative even function with a monotone tail of exact polynomial order. They are not satisfied by oscillating functions, for which cancellation in σ₂(0) may occur; in that case the upper bound of Theorem 3.10 need not be attained.

3.5 The resulting algorithm

We now summarize the method. The input consists of the function f, the decay data of f and f^, a target accuracy ε>0 and the nodes at which f^ is to be evaluated. The output consists of approximate values of f^ together with a rigorous a priori error bound.

Remark 3.22. The total error of Algorithm 1 splits into the approximation error analyzed in this paper and the error of the NFFT itself,

where s~P,L denotes the value actually computed by the NFFT. It is essential here that the first term is bounded uniformly in v, because the nodes v_*j* are arbitrary; a bound on the grid {k/L} only, as provided by [1], would not suffice. Since the NFFT error decays exponentially in the cut-off parameter m, the second term can always be made negligible against the first at a cost which is linear in m; see Table 4.

4 Numerical examples

In the previous section, we derived explicit error estimates for the uniform approximation error M_{*P, L*}(f) under polynomial, exponential, and mixed decay conditions on f and f^. We now illustrate these theoretical bounds by means of concrete examples. For each example, we specify the relevant decay parameters and state the resulting error estimate from Section 3 explicitly before comparing it numerically with the true approximation error M_{*P, L*}(f).

Example 3.23. Let M₂ be the centered linear B-spline given by

M2(x):={1−|x|x∈[−1,1],0x∈ℝ∖[−1,1].

Then f = M₂ has the Fourier transform

f^(v)=∫-11(1-|x|)cos(2πvx)dx=(sinc(πv))2, v∈ℝ.

Since suppf = [−1, 1], the function f∈C(ℝ) is space-limited with M = 2, and f^∈C(ℝ) has polynomial decay with b = 2. If we estimate the even function f^(v)(1+|v|)2 separately on the intervals [0,12] and [12,∞), we see that d=94 is a possible choice. By Corollary 3.14, for any L≥2 it holds

MP,L(f)≤2d(22-1)ζ(2)L-1/2P-2=π2LdP-2.

Hence the error decays like P⁻² as P → ∞, which is confirmed by the numerical results in Figure 2a. More generally, for the B-spline M_2*m* of even order 2m, the Fourier transform M2m^ has polynomial decay with b = 2m, and Corollary 3.14 yields the estimate

MP,L(f)≤2d(22m-1)ζ(2m)L-1/2P-2m.

This P^−2*m* decay is illustrated in Figure 2b (m = 2) and Figure 3a (m = 3).

Alternatively, consider the band-limited function

g(x):=M2(sinc(M2πx))2, x∈ℝ,

whose Fourier transform is ĝ(v)=M2(2Mv) with suppĝ⊆[-M2,M2]. Since g has polynomial decay with a = 2, Corollary 3.14 gives, for P≥M,

MP,L(g)≤2c(22-1)ζ(2)L1/2-2=π2cL-3/2.

Hence the error decays like L^−3/2 as L → ∞, which is confirmed numerically in Example 3.27 and Figure 1a.

Example 3.24. The Gaussian function

f(x)=e-πx2, x∈ℝ,

has exponential decay with r = π and α = 2, and its Fourier transform f^(v)=e-πv2 has the same exponential decay with s = π and β = 2. Provided π(L/2)²≥1 and π(P/2)²≥1, Theorem 3.16 yields

MP,L(f)≤1L(2(1+18)e-π(P/2)2+2(1+18)Le-π(L/2)2)(24)

=94L(e-πP2/4+Le-πL2/4).(25)

For fixed L = P we obtain

ML,L(f)≤9(L+1)4Le-πL2/4.

In Figure 3b we plot the error M_{*P, P*} as well as the analytical estimate for P = 1, …8.

Example 3.25. This example is taken from [1, Section 3.2], where it is used to illustrate the mixed-decay case for the scaled Euclidean error; we use it here in order to make the two sets of results directly comparable, and additionally in Section 5. Let

f(x):=e-2π|x-δ|, x∈ℝ,

with shift parameter δ∈ℝ. Its Fourier transform is

f^(v)=e-2πiδvπ(1+v2), v∈ℝ.

Hence f has exponential decay with r = 2π and α = 1, and f^ has polynomial decay with b = 2 and constant d = 1/π. This is the mixed-decay case of Theorem 3.18. Provided 2π(L/2)≥1, i.e., L≥1/π, we obtain

MP,L(f)≤1L(2π(22-1)ζ(2)P-2+2(1+12)Le-2π(L/2)) =1L(2πP2+3Le-πL).

For sufficiently large L, the exponential term is negligible, and the error is dominated by the polynomial term 2πLP2. A good parameter choice is L = 8 and P = 2^*t*, t≥6; for these values the exponential contribution is below machine precision, and the estimate reduces to MP,L(f)≈2π8P2, confirming the observed P⁻² decay in Figure 4a.

Example 3.26. For fixed parameter β∈ℕ we consider the function

f(x):={(1−x2)β−1/2 x∈[−1,1],0 x∈ℝ∖[−1,1].

Using [30, p. 8], we determine the corresponding Fourier transform

f^(v)=∫-11(1-x2)β-1/2e-2πivxdx=2∫01(1-x2)β-1/2cos(2πvx)dx[1ex]=π(2β)!4ββ!{(πv)-βJβ(2πv) v∈ℝ\{0},1β! v=0.(26)

and obtain the estimate [23, Section 5.3]

For β = 2 we observe by Corollary 3.14 with M = 1 and b = 2 the decay P^−2.5, see Figure 4b.

Example 3.27. So far, all experiments were carried out for a fixed truncation parameter L and increasing bandwidth P, i.e. they probe the aliasing term L-1/2δ(f^,P) of Equation 11. We now isolate the truncation term L^1/2δ(f, L) and keep P fixed while L increases. To this end, we use the band-limited function of Example 3.23,

g(x)=M2(sinc(M2x))2, ĝ(v)=M2(2vM) M=2

for which suppĝ⊆[−1, 1], hence δ(ĝ, P) = 0 for all P≥2, and g has polynomial decay with a = 2 and c=supx|g(x)|(1+|x|)2=1.267505. By Corollary 3.14 we obtain the estimate

MP,L(g) ≤ 2c(22-1)ζ(2)L-3/2=π2cL-3/2,

and by Theorem 3.20, the matching lower bound of the same order L^−3/2.

Figure 1a shows the measured errors for L = 2¹, …, 2⁹ and P = 8. The observed convergence order, computed from consecutive values, is −1.457, −1.487, −1.497, −1.499, −1.500, −1.500, −1.500, −1.500, i.e., it converges to the predicted exponent -32=12-a. Repeating the experiment with P = 32 reproduces the same values to five significant digits, which confirms that the truncation term of Equation 11 is indeed independent of P.

The same figure shows the scaled Euclidean error EP,L(g) of Ehler et al. [1]. In this truncation-dominated regime, both errors have the same order, and the quotient EP,L(g)/MP,L*(g) increases monotonically from 1.216 at L = 2 to the limit 1.244, in agreement with Remark 3.12. Hence, for this class of functions, our uniform estimate loses nothing compared with the Euclidean estimate, although it controls the error at every v and not only at the LP grid points.

5 Computational cost

The estimates of Section 3 answer the question how accurate the sampling polynomial s_{*P, L*} is. We now address the complementary question how expensive it is to evaluate it, and compare the three possibilities listed in Algorithm 1:

• Direct evaluation of Equation 2 at K nodes, which requires O(LPK) arithmetic operations,
• The FFT, which delivers the LP values s_{*P, L*}(k/L), k∈[LP], on the equispaced grid in O(LPlogLP) operations,
• The NFFT, which delivers the values s_{*P, L*}(v_*j*), j = 1, …, K, at arbitrary nodes in O(LPlogLP+K) operations.

All computations were performed in double precision on an Intel^® Core™ i7-6700 CPU at 3.40 GHz, using Python 3.13 with NumPy 2.4 and the NFFT3 library [19] through its Python interface pyNFFT3 1.0.2 [22]. As a test function, we use the two-sided exponential f(x) = e^−2π|*x*| of Example 3.25 with L = 8, for which f^(v)=(π(1+v2))-1 is known in closed form. Timings are the minimum over five runs after one warm-up run.

Table 2 compares the FFT with direct evaluation on the equispaced grid. Both compute the same quantity, and indeed they agree to machine precision; the FFT is faster by up to four orders of magnitude. The last column recalls the corresponding approximation error M_{*P, L*}(f), which is the quantity bounded by Theorem 3.18. It is worth noting that the difference between the two algorithms is purely a matter of speed: the accuracy of the computed Fourier transform is governed entirely by M_{*P, L*}(f) and not by the evaluation scheme.

Table 2. Equispaced output nodes v = k/L, k∈[LP]: direct evaluation vs. FFT, for f(x) = e^−2π|*x*| and L = 8.

The column ||·||_∞ contains the maximum deviation between the two computed results.

The situation of real interest, however, is the evaluation at arbitrary nodes, where the FFT is not applicable at all. Table 3 compares the NFFT with direct evaluation at K = 10, 000 nodes drawn uniformly at random from [−P/2, P/2). The NFFT is used with an oversampling factor σ = 2 and cut-off parameter m = 8; the resulting relative deviation e_NFFT from the directly computed values stays at the level of the round-off error, while the speed-up reaches a factor of 4, 262 at P = 4, 096. The column maxj|f^(vj)-sP,L(vj)| shows the error that is actually committed with respect to the true Fourier transform. It agrees with LMP,L(f) to three digits in every row, i.e., the random nodes do attain the uniform bound. This illustrates once more that it is the uniform error (Equation 5), and not its discrete version (Equation 6), which governs the accuracy of the NFFT-based computation.

Table 3. K = 10, 000 random non-equispaced nodes in [−P/2, P/2): direct evaluation vs. NFFT (σ = 2, m = 8), for f(x) = e^−2π|*x*| and L = 8.

Here e_NFFT is the relative maximum deviation of the NFFT from the direct evaluation.

Finally, Table 4 documents step 4 of Algorithm 1, i.e., the choice of the NFFT cut-off parameter m. As predicted by Potts and Tasche [23, 24], the NFFT error decays exponentially in m, while the runtime grows only mildly. For P = 256 and L = 8, the approximation error is LMP,L(f)=2.85·10-5, so that already m = 3 makes the NFFT error negligible; the frequently used default m = 8 reduces it to the level of round-off. In other words, the accuracy of the whole computation is determined by the parameters L and P through Theorem 3.5, and the NFFT contributes no relevant error at negligible additional cost.

Table 4. Influence of the NFFT cut-off parameter m for f(x) = e^−2π|*x*|, P = 256, L = 8, K = 10, 000 and σ = 2.

For comparison, the approximation error of the sampling polynomial is LMP,L(f)=2.85·10-5.

6 Conclusion

We have studied how well the continuous Fourier transform f^ of an integrable function f can be approximated by the P-periodic trigonometric sampling polynomial s_{*P, L*} of degree LP/2, measured in the uniform norm on the whole frequency interval [−P/2, P/2].

The central result is Theorem 3.5, which bounds the scaled maximum approximation error M_{*P, L*}(f) by the linear combination L-1/2δ(f^,P)+L1/2δ(f,L) of the decay rates (Equation 8) of f^ and f. Combined with the explicit estimates of the decay rates in Lemma 3.9 and Lemma 3.15, this yields fully explicit error bounds under polynomial, exponential, and mixed decay assumptions (Theorems 3.10, 3.16, and 3.18), in which every constant is determined by the decay parameters of f and f^. Theorem 3.20 shows that the resulting rate in the truncation parameter L is sharp.

These estimates have two consequences of practical importance. First, they can be used a priori: given a target accuracy, the parameters L and P can be determined before the computation is started, which is what Algorithm 1 does. Second, since the bound is uniform in v, it remains valid when s_{*P, L*} is evaluated at arbitrary nodes by an NFFT. The NFFT thereby becomes a numerical method for the Fourier transform itself and not merely a fast evaluation scheme for trigonometric polynomials.

Our results also clarify the relation to the work of Ehler et al. [1]. The scaled maximum error and the scaled Euclidean error are not equivalent; their quotient ranges over [1,LP], and both extremes occur (Remark 3.13). In the truncation-dominated regime, the two error measures agree up to a constant, whereas in the aliasing-dominated regime, the uniform estimate derived here is better by the factor LP; Remark 3.12 traces this back to the different oscillatory behavior of the aliasing and the truncation remainder, and explains why no symmetric bound can hold for the uniform error.

Correction note

This article has been corrected with minor changes. These changes do not impact the scientific content of the article.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Firdous A. Shah, University of Kashmir, India

Reviewed by: Amitava Ghosh, Institute of Engineering and Management (IEM), India

Yunfei Yang, Sun Yat-sen University, China

Kusum Sharma, National Institute of Technology Uttarakhand, India