On the resolution power of Fourier extensions
for oscillatory functions

The purpose of this paper is to discuss an alternative to polynomial approximation for nonperiodic functions; the so-called Fourier extension. Our main result is to show that Fourier extensions have excellent resolution properties. In particular, there is a user-determined parameter $T\in(1,\infty)$ that allows for continuous variation of the resolution constant from the value $2$ (in the limit $T\rightarrow 1$), the figure corresponding to Fourier series, to $\pi$ ($T\rightarrow\infty$), the value obtained by polynomial approximations. Thus, Fourier extensions are highly suitable tools for problems with oscillations at moderate to high frequencies.

Let us now describe Fourier extensions in more detail. As discussed, Fourier series are eminently suitable for approximating smooth and periodic oscillatory functions. It follows that a potential means to recover a highly accurate approximation of a nonperiodic function $f:[-1,1]\rightarrow\R$ is to seek to extend $f$ to a periodic function $g$ defined on a larger domain $[-T,T]$ and compute the Fourier series of $g$. Unfortunately, unless $f$ is itself periodic, no periodic extension $g$ will be analytic, and hence only spectral convergence can be expected. Preferably, for analytic $f$, we seek an approximation that converges exponentially fast.

To remedy this situation, rather than computing an extension $g$ explicitly, it was proposed in [7, 11] to directly compute a Fourier representation of $f$ on the domain $[-T,T]$ via the so-called Fourier extension problem:

Note that there are infinitely many ways in which to define a smooth and periodic extension of a function $f:[-1,1]\rightarrow\mathbb{C}$, of which (3) is but one. We say an extension $g$ is a Fourier extension if $g\in G_{n}$ is a finite Fourier series on $[-T,T]$. Note, however, that even within this stipulation, there are still infinitely many ways to define $g$. We refer to the extension $g_{n}$ defined by (3) as the continuous Fourier extension of $f$ (in the sense that it minimizes the continuous $L^{2}$-norm).

As numerically observed in [7, 11], the convergence of $g_{n}$ to $f$ is exponential, provided $f$ is analytic. When $T=2$, this was confirmed by the analysis presented in [25]. A pivotal role is played by the map

The importance of this map is that it transforms the trigonometric basis functions that span the space $G_{n}$ into polynomials in $y$. In this setting, Problem 1.1 reduces simply to the computation of expansions in a suitable basis of nonclassical orthogonal polynomials, since the least squares criterion corresponds exactly to an orthogonal projection in a particular weighted norm. Well-known results on orthogonal polynomials can then be used to establish convergence properties of the Fourier extension. We recall this analysis in greater detail in §2, along with its generalization to arbitrary $T$.

The map (4) demonstrates the close relationship that exists between Fourier series and polynomials. Note the similarity to the classical Chebyshev map $x=\cos\theta$ that takes Chebyshev polynomials to trigonometric basis functions. Compared to Chebyshev expansions, however, it is important to note that the roles of polynomials and trigonometric functions are interchanged. The Chebyshev expansion of a given function $f$ is a polynomial approximation, which is equivalent to the Fourier series of a related function. The Fourier extension of a function on the other hand is a trigonometric approximation, which is equivalent to the polynomial expansion of a related function. In that sense, Fourier extensions and Chebyshev expansions are dual to each other.

The main result of this paper is that the resolution constant $r=r(T)$ of the continuous Fourier extension satisfies

Accordingly, the Fourier extension $g_{n}$ of the function (1) will begin to converge once $n$ exceeds $\frac{1}{2}r(T)\omega$ (recall that $g_{n}$ involves $2n+1$ degrees of freedom). In particular, we find that $r(T)\sim 2T$ for $T\approx 1$ and $r(T)\sim\pi$ as $T\rightarrow\infty$. Note that the resolution of the function $f(x)=\exp(i\omega\pi x)$ using classical Fourier series on $[-T,T]$ would require a minimum of $2T\omega$ degrees of freedom (whenever $f$ is periodic). Thus, Fourier extensions exhibit comparable performance for $T$ close to $1$. However, as $T$ increases, $f$ can be resolved more efficiently via its Fourier extension than if it had been directly expanded in a Fourier series on $[-T,T]$. In particular, the resolution power is bounded above for all $T$ by $\pi$, which is precisely the resolution constant for polynomial approximations.

Aside from establishing when asymptotic convergence will occur, we also show that the convergence in this regime is exponential, and that for all sufficiently large $n$ the rate is precisely $\rho=E(T)$, where

(this result actually holds for all sufficiently analytic functions $f$, not just (1)). Here, we note that $E(1)=1$, implying no exponential convergence for $T=1$. This is of course a consequence of the Gibbs phenomenon of Fourier series on $[-1,1]$. Convergence is exponential for all $T$ greater than $1$, however, and the rate increases for larger $T$.

In summary, the main conclusion we draw in this paper is the following. Smaller $T$ yields better resolution power of the continuous Fourier extension, at a cost of slower, but still exponential, convergence. Conversely, larger $T$ yields faster exponential convergence, at the expense of reduced resolution power. Formally, one may also obtain a resolution constant of $2$ in the limit $n\rightarrow\infty$ by allowing $T\rightarrow 1$ as $n\rightarrow\infty$, and towards the end of the paper we shall discuss different strategies for doing this, some of which are quite effective in practice.

Unfortunately, the linear system of equations to be solved when computing (3) is severely ill-conditioned. As a result, the best attainable accuracy with a continuous Fourier extension is of order $\sqrt{\epsilon}$, where $\epsilon$ is the machine precision used. To overcome this, we present a new Fourier extension (referred to as the discrete Fourier extension), based on a judicious choice of interpolation nodes, which leads to far less severe condition numbers. Numerical examples demonstrate much higher attainable accuracies with this approach, typically of order $\epsilon$, whilst retaining both the same rates of convergence and resolution constant.

Inherent to Fourier extensions is the concept of redundancy: namely, there are infinitely many possible extensions of a given function. As we explain, this not only leads to the ill-conditioning mentioned above, it also means that the result of a numerical computation may differ somewhat from the ‘theoretical’ (i.e. infinite precision) continuous or discrete Fourier extension. Later in the paper, we discuss this observation and its consequences in some detail.

The question of resolution power of Fourier series and Chebyshev polynomial expansions was first developed rigorously by Gottlieb and Orszag [20, p.35] (see also [33]), who introduced and popularized the concept of points per wavelength. Therein the figures of $2$ and $\pi$ respectively were derived, the latter being generalized to arbitrary Gegenbauer polynomial expansions in [21] (in §4 we provide an alternative proof, valid for almost any orthogonal polynomial system (OPS)). Resolution was also discussed in detail in [6, chpt. 2], where, aside from Fourier and Chebyshev expansions, the extremely poor performance of finite difference schemes was noted.

One-dimensional Fourier extensions were, arguably, first studied in detail in [7, 11], where they were employed to overcome the Gibbs phenomenon in standard Fourier expansions, as well as the Runge phenomenon in equispaced polynomial interpolation (see also [9]). Application to surface parametrizations was also explored in [11].

Similar ideas (typically referred to as Fourier embeddings), were previously proposed to solve PDE’s in complex geometries. These work by embedding the domain in a larger bounding box, and computing a Fourier series approximation on this domain. See [4, 5, 8, 12, 13, 14, 15, 19, 26, 29, 34, 36]. Such methods can be shown to work well, in principle even for arbitrary regions, but they are typically prohibitively computationally expensive.

More recently, a very effective method was developed in [10, 32] to solve time-dependent PDE’s in complex geometries. This approach is based on a technique for obtaining one-dimensional Fourier extensions, known as the FE–Gram method, which is then combined with an alternating direction technique, as well as standard FFTs, to solve the PDE efficiently and accurately. This approach has also been applied to Navier–Stokes equations [3]. Interestingly, it is shown and emphasized in [32] (see also [3]) that this method leads to an absence of dispersion errors (or pollution errors) – another beneficial property for wave simulations shared with classical Fourier series, and very much related to resolution power.

Having said this, we note at this point that the FE–Gram approach is quite different to the Fourier extensions we consider in this paper. The FE–Gram approximation is not exponentially convergent, and in practice, the order of convergence is usually limited by the user to ensure stability in the overall PDE solver. However, the FE–Gram approximation is designed specifically to be computed efficiently, in $\mathcal{O}(n\log n)$ time, using only function evaluations on an equispaced grid. Although there has been some recent progress in the rapid computation of exponentially-convergent Fourier extensions of the form (2) from equispaced data [30, 31], this is not an issue we shall dwell on in this paper. Thus, the main contribution of this paper is approximation-theoretic: we show that one can construct exponentially-convergent Fourier extensions with a resolution constant arbitrarily close to optimal. In §6 we discuss ongoing and future work pertaining to these issues.

This aside, our use of Fourier extensions in this paper to resolve oscillatory functions is superficially quite similar to the Kosloff–Tal–Ezer mapping (see [27], [6, chpt 16.9] and references therein, and more recently, [24]) for improving the severe time-step restriction inherent in Chebyshev spectral methods. Such an approach also improves the very much related property of resolution power. In this technique, one replaces the standard Chebyshev interpolation nodes with a sequence of mapped points, and expands in a nonpolynomial basis defined via the particular mapping. Roughly speaking, with Fourier extensions, the situation is reversed. Rather than fixing interpolation nodes, one specifies a particular basis (i.e. the Fourier basis for a particular $T$) that gives rapidly convergent expansions and good resolution properties, and chooses an appropriate means for computing the extension (for example, a particular configuration of collocation nodes) via the corresponding mapping.

Finally, we remark in passing that there are a number of commonly used alternatives to Fourier extensions. Such methods typically arise from the desire to reconstruct a function directly from its Fourier coefficients (or pointwise values on an equispaced grid), whether via re-expanding in a sequence Gegenbauer polynomials [22, 23], or by smoothing the function by implicitly matching its derivatives at the domain boundary [18]. Whilst the latter retains a resolution constant of $2$, it only yields algebraic convergence of a finite order, and suffers from severe ill-conditioning. Conversely, the Gegenbauer reconstruction procedure offers exponential convergence, but with a significant deterioration in resolution power [21].

The outline of the remainder of this paper is as follows. In §2 we detail the convergence of Fourier extensions for arbitrary $T>1$, and in §3 we address numerical computation. In §4 we consider the resolution power of polynomial expansions, and in §5 we derive the resolution constant for Fourier extensions.

Standard approximations based on orthogonal polynomials (or Fourier series in the periodic case) converge exponentially fast provided $f$ is analytic, and spectrally fast if $f$ is only smooth [16, chpts 2,5]. Whenever $f$ has only finite regularity, convergence is algebraic at a rate determined by the degree of smoothness: specifically, if $f$ is $(k-1)$-times continuously differentiable in $[-1,1]$ and $f^{(k)}$ exists almost everywhere and is square-integrable (equivalently, $f\in H^{k}[-1,1]$ – the $k$th standard Sobolev space of functions defined on $[-1,1]$), then $f_{n}$ converges algebraically fast at rate $k$ . Our first result regarding Fourier extensions illustrates identical convergence in this setting:

The reason for ill-conditioning in the exact Fourier extension is simple: although the functions $\exp(i\frac{k\pi}{T}x)$ form an orthogonal basis on $[-T,T]$, they only constitute a frame when restricted to the interval $[-1,1]$.

Lemma 3.1 has several implications. Recall that frames are usually redundant: any given function $f$ typically has infinitely many representations in a particular frame. This implies that there will be approximate linear dependencies amongst the columns of the Gram matrix $A$ for all large $n$. This makes $A$ ill-conditioned, and in the case of the frame (15), leads to the aforementioned exponentially large condition numbers.

Let us consider several particular representations of a smooth function $f$ in the frame (15). Recall that associated to any frame $\{f_{k}\}$ is the so-called canonical dual frame $\{g_{k}\}$ [17]. A representation of $f$, the frame decomposition, is given in terms of this frame by

Since the frame (15) is tight, it coincides with its canonical dual up to a constant factor that is equal to the frame bound, in this case $T$. Therefore its frame decomposition is precisely

where the $a_{k}$’s are given by (16). However, as illustrated in the proof of the previous lemma, this is precisely the Fourier series of the discontinuous function $\tilde{f}$, and thus this infinite series converges only slowly and suffers from a Gibbs phenomenon.

On the other hand, as shown in the proof of Theorem 2.1, by extending $f$ more smoothly to $[-T,T]$, one obtain representations of $f$ in the frame (15) that converge algebraically at arbitrarily fast rates. Thus, we conclude the following: different frame expansions of the same function may give rise to wildly different approximations.

Clearly, given its exponential rate of convergence, the continuous Fourier extension $g_{n}$ will not coincide with any of the aforementioned representations (recall also that $g_{n}$ does not typically converge outside $[-1,1]$ [25], and therefore its limit cannot be a frame expansion in general). More importantly, however, it is not certain that the solution of the linear system (14), when computed in finite precision, will actually coincide with the continuous Fourier extension $g_{n}$ itself. The reason is straightforward: for large $n$, $A$ is approximately underdetermined, and therefore there will be many approximate solutions to (14) corresponding to different frame expansions of $f$. A typical linear solver will usually select an approximate solution with bounded coefficient vector $a$, and there is no guarantee that this need coincide with $g_{n}$.

This phenomenon has a significant potential consequence: since the numerically computed extension need not coincide with the continuous Fourier extension $g_{n}$, theoretical results for the behaviour of $g_{n}$ are not guaranteed to be inherited by the numerical solution. Having said this, in numerical experiments, one typically witnesses exponential convergence, exactly as predicted by Theorem 2.3. However, a difference can arise when investigating other properties, such as resolution power, as we discuss and explain in §5.

The above discussion is not intended to be rigorous. A full analysis of the behaviour of ‘numerical’ Fourier extensions is outside the scope of this paper. It transpires that this can be done, and a complete analysis will appear in a future paper [2]. Instead, in the remainder of this paper, we focus mainly on approximation-theoretic properties of theoretical extensions; resolution power, in particular. Nonetheless, in §5.3 we revisit the issue of numerical extensions insomuch as it relates to resolution, and give an explanation of the phenomenon of differing resolution power mentioned above.

The continuous Fourier extension is defined by a least squares criterion, and as such it suffers from the drawback of requiring computation of the integrals $\langle f,\phi_{k}\rangle$. To avoid this issue one may instead define a Fourier extension by a collocation condition. In other words, if $\{x_{i}\}^{2n+1}_{i=1}$ is a set of nodes in $[-1,1]$ we replace the matrix $A$ and the vector $B$ in (14) by $\tilde{A}$ and $\tilde{B}$ with entries $\phi_{k}(x_{j})$ and $f(x_{j})$ respectively. We then define $\tilde{A}a=\tilde{B}$ once more.

The key question is how to determine good nodes. Recall that $g_{n}$, as defined by (3), is essentially a sum of two polynomial approximations in the variable $y\in[c(T),1]$. The polynomial interpolant of an analytic function in $y$ at Chebyshev nodes (appropriately scaled to the interval $[c(T),1]$) converges exponentially fast at the same rate $\rho$ as its orthogonal polynomial expansion. Therefore, such nodes, when mapped to the original domain via $x=\frac{T}{\pi}\cos^{-1}y$, will ensure exponential convergence at rate $E(T)$ of the resulting Fourier extension.

It is now slightly easier to redefine $G_{n}$ to be the space of dimension $2n+2$ spanned by the functions $\{\cos\frac{\pi}{T}kx\}^{n}_{k=0}$ and $\{\sin\frac{\pi}{T}kx\}^{n+1}_{k=1}$. The $2n+2$ collocation nodes in $[-1,1]$ therefore take the form $\{x_{i}\}^{n}_{i=0}\cup\{-x_{i}\}^{n}_{i=0}$, where

Recall that $c(T)=\cos\frac{\pi}{T}$. We refer to such nodes as symmetric mapped Chebyshev nodes, and the corresponding extension $g_{n}$ as the discrete Fourier extension of the function $f$.

Besides removing the requirement to compute the integrals $\langle f,\phi_{n}\rangle$, this choice of nodes also has the significant effect of improving conditioning, as we now explain. We first note the following:

This lemma demonstrates that the normal equations for the discrete Fourier extension are the equations of a continuous Fourier extension based on the weighted $L^{2}_{W}$ optimization problem. It is well-known that forming normal equations leads to worse numerical results [28, Ch.19], and exactly the same is true in this case. Indeed, since $W$ is an integrable weight function, we may expect that $\kappa(A_{W})\approx E(T)^{2n}$, i.e. exactly as in the unweighted case. Thus collocation leads to the significant reduction in condition number, with $\kappa(\tilde{A})\approx E(T)^{n}$ as opposed to $E(T)^{2n}$. As a result, one typically observes a much higher accuracy, ${\mathcal{O}}(\epsilon)$ as opposed to ${\mathcal{O}}(\sqrt{\epsilon})$, with this approach.