Rapidly decaying Fourier-like bases

: In many applications it is natural to seek to extract a characteristic scale for a function’s variations by reference to a frequency spectrum. Although the moments of a spectrum appear to promise simple options to make such a connection, standard Fourier methods fail to yield ﬁnite moments when the function’s domain is itself ﬁnite. We investigate a family of Fourier-like bases with rapidly decaying spectra that yield well-deﬁned moments for such cases. These bases are derived by considering classes of functions for which a normalised mean square derivative is stationary. They are shown to provide precisely the type of spectrum needed to complete a recent investigation of mid-spatial frequency structure on optical surfaces [K. Liang, Opt. Express 27 , 3390-3408 (2019)].


Introduction
Fourier methods have a well-earned status in the toolbox of science and technology. The standard methods include Fourier transforms -both continuous and discrete-as well as Fourier series. The associated spectra, typically in terms of spatial or temporal frequencies, have widespread applications. Oftentimes, however, generalisations or variations of the standard methods give a vital boost to utility in certain contexts. These variations include such things as time-frequency analysis [1], wavelets [2], and linear canonical transformations [3]. Further, since the Fourier series is just one special case of Sturm-Liouville methods [4], other instances are close cousins with similar orthogonality-related properties and potential advantages.
It turns out that Fourier methods often encounter significant challenges in situations where a function is known over only a finite domain. Standard methods can sometimes be applied in these cases by taking arbitrary steps such as replicating the function to make it complete and periodic or by assuming the function to be zero outside the finite domain (for the Fourier series or transform, respectively). In both these cases, however, even when the finite domain can be replicated to effectively tesselate a complete function that is periodic in each independent variable, it is difficult to avoid spectral artefacts associated with the edges of the original domain. The methods described below can be considered as potentially useful tools in these frequently encountered situations.
A relevant context for this discussion is the study of the surface shape of as-built optical surfaces with complex nominal shapes including aspheric and freeform optics [5][6][7]. The tool marks left by the associated fabrication steps are generally referred to as mid-spatial frequency structure (MSF). In order to (i) characterise, (ii) model the optical impact of, and (iii) set tolerances for MSF, it is natural to seek the benefits that typically come from an intuitive spectral decomposition. For most of these purposes, an orthogonal polynomial basis has been shown to be well matched [8]. The second of these steps, however, is typically performed by using a perturbation approximation and we recently showed that the accuracy of that approximation can be estimated in terms of the mean squared Laplacian of the MSF [9]. In a simple one-dimensional model, it is the mean squared second derivative that is required, of course. Although we claimed that such entities could be related to a fourth moment of the associated Fourier spectrum, we did not explicitly show how such a spectrum can be computed in general. In fact, we did not appreciate the subtleties involved and, once they were understood, came to the view that the required methods and ideas are valuable in themselves and potentially of more general interest.
To have a finite fourth moment, a spectrum must be sufficiently well localised. Figure 1(a) shows for an analytic function that orthogonal (Legendre) polynomials yield spectra that decay exponentially [10], and hence possess finite spectral moments of any order. On the other hand, the standard Fourier basis decays slowly and is not localised enough for the calculation of a fourth -or even second -spectral moment. It is shown below that a class of rapidly decaying Fourier-like (RDF) decompositions has increasingly compact spectra, as demonstrated in Fig. 1, and can provide a useful connection between these two limits. In this way, we finally resolve the nature of the spectra referred to in [9] that gives access to an intuitive characteristic scale of a function's variations. Interestingly, some of these RDF bases are useful in a range of other applications [11][12][13][14][15].  1. Two different depictions of the power spectra of various bases as functions of (a) basis element index and (b) eigenvalue, for the antisymmetric sinusoid with 8.8 cycles over a finite region shown in the inset at top right of (a). Note that only the coefficients of the odd elements in each basis are shown since those of the even components vanish. The Fourier basis and the Legendre polynomials, shown only in (a), give the slowest and fastest decay rates, respectively. The artificial high-frequency components in the Fourier series are associated with the mismatched behavior at the interval's ends; in keeping with Parseval's theorem, the spectrum for the low frequency content is lowest for the Fourier basis. The RDF bases exhibit intermediate decay rates and admit progressively more finite moments.

RDF bases in one dimension
The derivation of the basis elements with the desired properties proceeds by requiring a certain normalized mean square derivative over the finite domain a = [−L/2, L/2] to take a stationary value. In particular, we look for real functions V(x) such that the ratio is stationary, where N ≥ 1 is a positive integer and V (N) (x) denotes the N-th derivative. This can be achieved by introducing α, a Lagrange multiplier, in order to construct the following functional: Upon slightly perturbing V, namely V(x) → V(x) + ε(x), there should be no change to L at first order [or equivalently to the value of the expression in Eq. (1)]. By making this substitution and expanding, we find Stationarity follows upon choosing V so that the quantity within square brackets in Eq. (3) vanishes for any ε(x). It is possible to rewrite this term by performing integration by parts N times on the first integral in order to shift all of the derivatives from ε to V. The resulting expression contains terms evaluated at the endpoints plus an integral in x whose integrand has a global factor of ε(x). Requiring the result to vanish for arbitrary ε(x) leads to the following 2N-th order Sturm-Liouville differential (eigenvalue) system: The solutions to Eq. (4) naturally split into even and odd forms given by where ( Although it is not immediately obvious, it turns out that the RDF elements in Eqs. (5) are real-valued. Because this is a Sturm-Liouville problem, it is guaranteed that there are an infinite number of eigenvalues (found to be roughly proportional to i).
In addition to the solutions in Eqs. (5), there are N solutions (degenerate for N>1) with eigenvalue equal to zero that, in keeping with orthonormality, can be taken to be where the Legendre polynomial P i (x) is taken to be defined over [−1, 1]. Note that these degenerate solutions are also alternately even and odd and automatically satisfy the boundary conditions of Eq. (4). In the limit N → ∞, only this degenerate polynomial subspace remains, with eigenvalue α = 0. This observation reveals an explicit connection of this family of bases between polynomials (N → ∞) and the Fourier basis (similar to the case of N = 1). For intermediate values of N, the non-degenerate subspace of solutions given in Eqs. (5) has properties that transition between orthogonal polynomials and Fourier series.
Since Eq. (4) is a Sturm-Liouville differential equation, it is well-known that the values of α e,i and α o,i alternate in ascending order. Hence, it is straightforward to assign a global ordering to RDF elements simply through the order of their associated eigenvalues. The degenerate orthogonal polynomials (which all have zero eigenvalue) can be ordered as the first N elements through their polynomial degree. Once ordered, the basis elements can simply be referred to through N V n (x) with eigenvalue α n , where n corresponds to the number of zeros. Figure 2 shows some of these functions for N = 1, 2, 3, and 4, along with the standard Fourier basis and the Legendre polynomials (N → ∞) for comparison. It is evident that larger N values give rise to non-degenerate basis elements that are steeper near the boundaries of the region. RDF functions are orthogonal under two different inner products: While the first of these relations is a natural result of the fact that these are the eigenfunctions of a Hermitian operator, namely (d/dx) 2N along with Eq. (4), the proof of the second requires N iterations of integration by parts as well as Eq. (4). It follows from the first of these orthogonality relations that the RDF expansion of any function µ(x) (e.g. the height of a surface with MSF structure) can be written as Additionally, the squared norm of µ(x), referred to here as M, can be expressed as the sum of the square of the coefficients following Parseval's relation as By using Eq. (4) and integration by parts it is possible to rewrite the integral in Eq. (9) as As can be seen from Eqs. (5), the first term inside the square brackets of Eq. (11) has a leading behavior that is proportional to α N−1 n . Therefore, in the limit of large n, the spectral coefficients | N c n | 2 generally decay as α −(2N+2) n , so the lines in the log-log plot in Fig. 1(a) approach a slope of −(2N + 2) for large n. Note that continuing to integrate Eq. (11) by parts leads to a closed-form expression for the RDF expansion whenever µ(x) is a polynomial.
The second orthogonality relation, Eq. (8b), gives access to a simple calculation for the normalized mean square of the N-th derivative of µ(x) akin to the form appearing in Eq. (1): where Eq. (8b) was used in the last step and M is defined in Eq. (10). This equation reveals that the mean square N-th derivative corresponds simply to a normalized spectral moment of order 2N. Furthermore, the n-th RDF element has n zeros across this region; this provides an intuitive understanding as to why the RDF spectra shown in Fig. 1(a) peak around the seventeenth basis element (the sinusoid in the inset of Fig. 1 has 17 zeros). In Fig. 1(b), which shows the same coefficients plotted against their eigenvalues, the peak is approximately at 2π 8.8 ≈ 55.3.
To give a clearer appreciation, the RDF basis for each of the cases of N = 1 and N = 2 is derived in the following subsections.

N = 1
For N = 1, Eq. (4) can be solved to give even and odd orthonormal solutions of the form where the allowed values of α e,i and α o,i are given by for i ≥ 1. The zero eigenvalue has the solution 1 V d,0 = 1/ √ L. After ordering the even, odd, and degenerate solutions in the manner described earlier, the eigenvalues can be written as α n = πn/L, with n ≥ 0. Note that the solutions given by Eqs. (13) satisfy Neumann boundary conditions, whereas the standard Fourier basis functions satisfy periodic boundary conditions (see Fig. 2).

N = 2
For N = 2, Eq. (4) can be solved to give orthonormal, non-degenerate solutions of the form The roots to Eq. (16) are well approximated by the simple expressions for i ≥ 1. The two degenerate solutions are given by 2 V d,0 (x) = 1/L and 2 V d,1 (x) = 12/L 3 x.
It should be noted that the RDF elements for N = 2 correspond directly to the vibrational modes of an unclamped one-dimensional bar [11][12][13]. As alluded to in Sec. 1, this basis is also a natural choice for the assessment of the validity of the perturbation model for the propagation of one-dimensional MSF structures in imaging systems and -as given in Eq. (12) -it has a valuable fourth spectral moment in that context [9].

RDF elements in two dimensions: flat space
The analysis in Secs. 2.1 and 2.2 is readily extended to find two-dimensional RDF bases for the case of a circular region a with radius R, where plane polar coordinates (r, φ) are a natural choice. For this case, there are multiple ways to generalize higher derivatives from one to two dimensions. In keeping with the motivation of minimizing RMS-derived quantities of the surface structure, consider the following functional that is analogous to Eq. (2): where we adopt the notation This means that the vertical bars in the first integral of Eq. (18) represent either a scalar (N even) or a vector (N odd) absolute value.
Similar variational calculations, described in Appendix A, can now be performed to ultimately arrive at the following 2N-th order Sturm-Liouville differential (eigenvalue) system: where the differential operatorD j is given bŷ The where m ∈ Z. It is always possible to label and order n such that n ∈ Z + . The remaining ordinary differential equation for N R m n (r), see Appendix B, has the solution where c 0 = 1 and N Γ nm is a normalization constant determined by The remaining values of c and the allowed values of α m n are determined by the boundary conditions in Eq. (20).
The differential (eigenvalue) system in Eq. (20), for a given value of m and N>1, has N/2 solutions (degenerate for N>3) with eigenvalue α = 0. These are given by where the radial Zernike polynomial Z m j (x) is defined as in [16]. In addition, for odd N and m = 0, there is one other solution for α = 0 of the form Once again, this degenerate subspace is all that remains as the solutions for Eq. (20) in the limit of N → ∞. All these solutions satisfy the following orthogonality relations: We again consider the cases of N = 1 and N = 2 explicitly.

N = 1
For N = 1, the two-dimensional version of Eq. (2) becomes The analogous variational calculation (see Appendix A) in two dimensions leads to the following second-order Sturm-Liouville partial differential system: Equation (29) can be solved to yield orthonormal solutions of the form 1 where J m is the m-th order regular Bessel function of the first kind. The allowed values of α m n are given by the roots of the following transcendental equation: For a given value of |m|, these eigenvalues are approximately given by for large n. The α = 0 solution, as indicated in Eq. (26), is the constant function

N = 2
For N = 2, the two-dimensional version of Eq. (2) becomes Once again, variational calculations can be performed to yield a fourth-order Sturm-Liouville partial differential equation: Again with Eq. (22), Eq. (35) can be solved to yield orthonormal solutions of the form for n ≥ 1, where I m is the m-th order modified Bessel function of the first kind. The allowed (discrete) values of α m n are the roots of the following transcendental equation: .

(37)
For a given value of |m|, these eigenvalues satisfy α m n ≈ (n + |m|)π/R, m even, (n + |m| − 1/2)π/R, m odd, for large n. The solutions with α = 0, as indicated by Eq. (25), are Figures 3 and 4 shows some of these basis elements. In analogy with the one-dimensional solutions in Eqs. (15), the solutions given by 2 V m n (r, φ) happen to describe the elastic modes of unclamped thin disks [11][12][13], and have been used to study deformable mirrors [14]. Furthermore, in two-dimensional paraxial imaging systems with MSF structures, these solutions are a natural basis for the assessment of the perturbation model [15]; the basis required for analogous studies of non-paraxial imaging systems is obtained by solving Eq. (35) on a sphere rather than on flat space, and this is considered in Sec. 4.

RDF bases in two dimensions: spherical space
Equation (35) can also be solved with the spherical polar coordinate system (ϑ, φ), where ϑ is the colatitude (polar) angle and φ is once again the azimuthal angle. We define the region a to be a spherical cap centered at the north pole, with colatitude boundary ϑ 0 . The derivation of the Sturm-Liouville differential system, analogous to that in Eq. (20), is shown in Appendix A. Furthermore, orthogonality relations like those in Eqs. (27) of the corresponding solutions to this differential system (the RDF basis elements for a spherical cap) are also stated in Appendix A.
For brevity, we analyze only the case N = 2, whose elements happen to correspond to the vibrational modes of an unclamped solid spherical cap. As mentioned at the end of Sec. 3.2, this basis arises naturally in the analysis of the perturbation model for MSF structures in non-paraxial imaging systems. For the coordinate system (ϑ, φ), the RDF elements are solutions to the following differential system: where here ∇ 2 denotes the Laplacian in the spherical space. In this case, we find that the RDF basis elements can be written in terms of the regular and modified spherical harmonics: where m ∈ Z, n ≥ 1, P m f ± is an associated Legendre function with The allowed values of the eigenvalue α m n correspond to the roots of the following transcendental equation: where the explicit dependence of f ± on α m n in Eq. (43) has been left off for convenience. For a given value of |m|, these eigenvalues satisfy α m n ≈ (n + |m|)π/ϑ 0 , m even, (n + |m| − 1/2)π/ϑ 0 , m odd, for large n. The α = 0 solutions are given by where the denominator is a normalization constant with 2 F 1 (a, b, c, z) being a hypergeometric function. Figure 5 shows some of these basis elements.

Concluding remarks
We were motivated to pursue the ideas reported above by the fact that the fourth moment of an RDF-spectrum gives intuitive access to a physically significant characteristic length in investigations of MSF structure on optical surfaces [9,15]. As indicated in Eqs. (12) and more generally from (27b), these RDF bases allow the determination of the mean square N-th derivative of a function in terms of a spectral moment of order 2N. For example, mean square slopes can be determined from the second moment of the spectra associated with the RDF bases presented in Secs. 2.1 and 3.1. Notice that this is quite different from a slope-orthogonal polynomial basis (akin to the Q polynomials [17], for example) where the mean-square slope is found by summing the squares of the spectral coefficients themselves. As pointed out in the Sec. 1, however, those polynomials are well matched to other aspects of the investigation of the characterization and modelling of MSF structure. By contrast, it follows from Eqs. (8a) and (27a), that the sum of squares of the RDF spectral coefficients is the mean square value of the decomposed function itself, see Eq. (10). In some contexts, it may also be of value to exploit the curious fact that the RDF bases are also bi-orthogonal in the more general sense of Eq. (55). Such properties are indicators of strong similarities to Fourier bases and give RDF bases potential value for a variety of applications that involve functions defined over finite domains. Indeed, certain subsets of RDF bases are already used in the study of elastic deformations of finite domains in different geometries [11][12][13][14][15].

Appendix A. Details of RDF derivation in two dimensions
In this appendix, we show the details of several steps in the derivation of the RDF bases in two dimensions. This will be done in general curvilinear coordinates (u, v), and can be specialized to treat both the plane polar and spherical polar coordinates used in Secs. 3 and 4, respectively. We begin with where dA is the infinitesimal area element of the corresponding metric associated with (u, v).
As in the one-dimensional case, one perturbs V(u, v) by ε(u, v) and requires L to be stationary (unchanged to first order), hence where the · operation is a scalar product between two scalar quantities or two vector quantities, depending on whether N is even or odd, respectively. Stationarity then requires the quantity within the square brackets of Eq. (47) to vanish for any ε. The first step in this process is to transfer the derivatives from ε over to V in the first integral within the square brackets. That integral can be written, using the product rule, as (Note that, depending on the parity of N, either the product within the integrand of the left-hand side or that in the second term of the right-hand side are a scalar product of two vectors while the other is a multiplication of two scalars.) We now use Gauss's theorem to express the first integral on the right side of Eq. (48) as a boundary term, hence At this point, one derivative has been transferred from ε to V, but we are left with a boundary term that must vanish for any ε. This condition is given bŷ wheren a unit vector that is parallel to d ì S ⊥ . Note that Eq. (51) is consistent with Eq. (21) when a is chosen to be a circular region of radius R and (u, v) = (r, φ), the plane polar coordinates.
Repeating the process from Eqs. (48) to (49) N times and forcing each boundary term to vanish, the quantity in the square brackets in Eq. (47) ultimately becomes By requiring Eq. (52) and the boundary terms to vanish for any ε, we end up with the following differential system: Equation (53)  We note for completeness that the solutions N V m n (u, v) also satisfy the following more general orthogonality relation: valid for any non-negative integer N 1 and N 2 for which N 1 + N 2 is an integer multiple of 2N. This relation includes Eqs. (54) as special cases. The one-dimensional analogue for these relations arises by replacing ∇ with a simple derivative. For completeness, we state the gradient and divergence for plane polar coordinates (r, φ) and the spherical polar coordinates (ϑ, φ). The former is given by ∇g(r, φ) = ∂ r g, ∇ · [f r (r, φ), f φ (r, φ)] = 1 r ∂ r (rf r ) and the latter by The corresponding area elements are given by dA = rdrdφ and dA = sin(ϑ)dϑdφ, respectively.

Appendix B. Separation of variables
Here, we show that the method of separation of variables may be used to solve the differential equation in Eq. (53) in the specific case of plane polar coordinates (r, φ); the procedure in spherical polar coordinates is similar since both coordinate systems have the same azimuthal coordinate φ.
It follows from Eqs. (56) and (57) that the differential equation in Eq. (53) is By assuming a solution of the form V(r, φ) = R(r) Φ m (φ), Eq. (60) can be written as Equation (61) is an ordinary differential equation for R(r) and its solutions are given by Eq. (23).