Exponential shapelets: basis functions for data analysis of isolated features

Bergé, Joel; Massey, Richard; Baghi, Quentin; Touboul, Pierre

doi:10.1093/mnras/stz787

Abstract

We introduce one- and two-dimensional ‘exponential shapelets’: orthonormal basis functions that efficiently model isolated features in data. They are built from eigenfunctions of the quantum mechanical hydrogen atom, and inherit mathematics with elegant properties under Fourier transform, and hence (de)convolution. For a wide variety of data, exponential shapelets compress information better than Gauss–Hermite/Gauss–Laguerre (‘shapelet’) decomposition, and generalize previous attempts that were limited to 1D or circularly symmetric basis functions. We discuss example applications in astronomy, fundamental physics, and space geodesy.

methods: data analysis

1 INTRODUCTION

A frequent task in data analysis is to categorize and quantify the shapes of localized objects – such as transient events in a (one-dimensional) time-series, or regions of interest in a (two-dimensional) image. It is such a universal challenge that methods developed for one field frequently turn out to be useful in others. For example, astrophysicists measure the shapes of distant galaxies by decomposing them into orthogonal basis functions, such as CHEFs (Jiménez-Teja & Benítez 2012) or (Gaussian) shapelets (Bernstein & Jarvis 2002; Refregier 2003; Refregier & Bacon 2003; Massey & Refregier 2005). Shapelets have been used to analyse data in other branches of astrophysics, modelling extrasolar planets (Hoekstra et al. 2005; Amara & Quanz 2012), the distribution of dark matter (Birrer, Amara & Refregier 2015; Tagore & Jackson 2016), or flashes of pulsars (Ellis & Cornish 2015; Lentati et al. 2015; Desvignes et al. 2016). They have also been used in medical imaging (Weissman, Hancewicz & Kaplan 2004; Apostolopoulos et al. 2017), pattern recognition in human vision (Sharpee & Victor 2009) or artificial vision (Sabzmeydani & Mori 2007), data compression (Holbrey 2006), and the manufacture of nanoscale thin films (Suderman & Lizotte 2015; Akdeniz, Lizotte & Mohieddin Abukhdeir 2018).

Gaussian shapelets are based on eigenfunctions of the quantum mechanics harmonic oscillator. In one-dimensional form, they are Gauss–Hermite functions, which seem to be well adapted in several time-series analyses, especially when transients (whose shape is close to damped oscillations) must be detected, characterized, and/or corrected for. This is the case for instance in fundamental physics for MICROSCOPE ‘crackles’ (e.g. Baghi et al. 2015; Bergé et al. 2015), space geodesy for GRACE ‘twangs’ (Flury, Bettadpur & Tapley 2008; Peterseim, Jakob & Schlicht 2010; Peterseim et al. 2014), or ‘glitches’ in gravitational wave searches with LIGO (Cornish & Littenberg 2015; Powell et al. 2015, 2017; Principe & Pinto 2017 and references therein). In two-dimensional form, they can be expressed equivalently as either Gauss–Hermite (Cartesian; Refregier 2003) or Gauss–Laguerre (polar; Bernstein & Jarvis 2002; Massey & Refregier 2005) functions. Owing to their quantum mechanical origin, they have elegant properties (they make a complete set) under Fourier transform, and are hence efficient for operations involving convolution or deconvolution of two images.

The main limitation of shapelets is that they are perturbations around a Gaussian, which is flat near its peak. They inefficiently parametrize peaky features, including the distant galaxies for which they were originally suggested (Melchior et al. 2010). Galaxies have a two-dimensional surface density that decreases approximately exponentially with distance from the centre. Capturing the steep gradient near the centre requires a weighted sum of many Gaussian shapelet basis functions, which then overfit noise in the extended wings. Attempting to overcome this limitation, Ngan et al. (2009) developed ‘Sersiclets’, although they were forced to be circularly symmetric, and have not seen wider applications.

In this paper, we extend the quantum mechanical framework underpinning Gaussian shapelets, and define 1D and 2D exponential shapelets based on wavefunctions of the hydrogen atom. These functions are perturbations around a decreasing exponential, and should efficiently model any arbitrarily shaped but centrally peaked regions of interest. The perturbations themselves are Laguerre polynomials.

We introduce 1D exponential shapelets in Section 2, and 2D exponential shapelets in Section 3. In both cases, we describe their main properties, compare them to Gaussian shapelets, and provide example uses. We conclude in Section 4.

2 1D EXPONENTIAL SHAPELETS

2.1 Definition

The 1D hydrogen atom is the solution to the motion of a particle in a 1D Coulomb potential 1/|x|. In this paper, we will neither dwell on its rich history, nor the debates about the finitude of its ground state, and about the existence of even wavefunctions (see Nieto 1979; Palma & Raff 2006; Nú nez-Yépez, Salas-Brito & Solis 2011, 2014 and references therein). Instead, we simply exploit the normalized 1D hydrogen atom wavefunctions as given by Palma & Raff (2006) to define the 1D exponential shapelet basis functions

\begin{eqnarray*} \Psi _n^+(x; \beta) \! = \! \frac{(-1)^{n-1}}{\sqrt{n^3\beta }}\, \frac{2x}{n\beta }\, L_{n-1}^1\left(\frac{2x}{n\beta }\right) \exp \left(-\frac{x}{n\beta }\right)~~ \forall x \geqslant 0, \end{eqnarray*}

(1)

for n ≥ 1, where |$L_{n-1}^1$| are the generalized Laguerre polynomials (see e.g. Massey & Refregier 2005). The characteristic scale size β corresponds to the Bohr radius, and n is the eigenfunction energy level.

Note that, contrary to normal procedure in quantum physics, we restrict 1D exponential shapelets on positive x, and refrain from defining the negative-x part |$\Psi _n^-(x; \beta) = -\Psi _n^+(-x; \beta)$| (for x < 0). Hence, we do not follow Palma & Raff (2006) and do not define the even and odd wavefunctions |$\Psi _n^+(x; \beta) \pm \Psi _n^-(x; \beta)$|⁠. Events in time-series can often be adequately described by their behaviour after a certain moment (in this case, x = 0). Henceforth, we shall therefore drop the ⁺ superscript, to write more simply |$\Psi _n(x; \beta) \equiv \Psi _n^+(x; \beta)$|⁠. These functions are continuously differentiable, smoothly departing from zero at x = 0 and tending back to zero as x → ∞.

Because the basis functions (1) are eigenfunctions of the 1D hydrogen atom’s Hamiltonian, they form an orthogonal basis of the square integrable functions L²([0, ∞[, 〈 ·, ·〉) Hilbert space equipped with the inner product |$\langle f(x),g(x) \rangle =\int _0^\infty f(x)g^*(x) {\rm d}x$| and an asterisk denotes complex conjugation. They are orthonormal, in the sense that |$\int _0^\infty \Psi _n(x; \beta)\Psi _m(x; \beta)\, {\rm d}x = \delta _{nm}$|⁠, where δ_nm is the Kronecker symbol, and complete, in the sense that |$\sum _{n=1}^{\infty } \Psi _n(x; \beta) \Psi _n(x^{\prime }; \beta) = \delta (x - x^{\prime })$|⁠.

They thus form a basis on which we can uniquely decompose a function as

\begin{eqnarray*} f(x) = \sum _{n=1}^\infty f_n \Psi _n(x; \beta), \end{eqnarray*}

(2)

with coefficients f_n given by an overlap integral

\begin{eqnarray*} f_n = \int _0^\infty f(x) \Psi _n(x; \beta)\, {\rm d}x. \end{eqnarray*}

(3)

Bessel’s inequality then assures us that for any function f ∈ L²([0, ∞[, 〈 ·, ·〉),

\begin{eqnarray*} \sum _{n=1}^{\infty } |f_n|^2 \leqslant ||f||^2, \end{eqnarray*}

(4)

where ||.|| is the L² norm, so that the series (2) converges and coefficients f_n must vanish as n increases. The series can therefore be truncated for suitably localized functions, at some value n ≤ n_max.

2.2 Properties

2.2.1 Maximum and minimum effective scales

The first few 1D exponential shapelet basis functions Ψ_n(x; β) are shown in Fig. 1. As the order n increases (with constant β), the largest-scale oscillation dominates, and rapidly acquires a larger extent. However, the smallest oscillations always remain roughly the same size.

Figure 1.

The first few 1D exponential shapelets Ψ_n(x), for β = 1. The inset shows a zoom on the smallest x.

Open in new tab Download slide

To act as a convenient figure of merit, we follow Refregier (2003) in defining the largest scale that can be described by an exponential shapelet model as

\begin{eqnarray*} \theta _{\rm max} \approx \sqrt{\frac{\int x^2 \Psi _{n_{\rm max}}(x;\beta)\, {\rm d}x}{\int \Psi _{n_{\rm max}}(x;\beta)\, {\rm d}x}} = \beta \sqrt{2{n_{\rm max}}(2{n_{\rm max}}+1)}\, . \end{eqnarray*}

(5)

Empirically, the smallest scale that can be described is

\begin{eqnarray*} \theta _{\rm min} \approx \beta . \end{eqnarray*}

(6)

As n_max → ∞, the range of scales modelled by a 1D exponential shapelet decomposition, θ_max/θ_min ∝ n_max = n_coeffs. However, the resolution of an exponential shapelet model is greatest near the origin, and decreases away from it. This behaviour is very different from Gaussian shapelets, where the resolution is more spatially uniform: Increasing n increases the scale of the model slowly, while simultaneously adding smaller-scale oscillations (see fig. 1 of Refregier 2003). The wide effective range of scales for exponential shapelets is what will allow them to efficiently describe spiky but long-duration events.

2.2.2 Fourier and Laplace transforms

To compute the Fourier transform of the 1D exponential shapelet basis functions, we adopt a convention where the Fourier transform is defined as

\begin{eqnarray*} \begin{array}{ccc}\tilde{f}(k) & = & \displaystyle(2\pi)^{-\frac{1}{2}} \int _{-\infty }^\infty f(x) {\rm e}^{ikx}\, {\rm d}x, \\ f(x) & = &\displaystyle (2\pi)^{-\frac{1}{2}} \int _{-\infty }^\infty \tilde{f}(k) {\rm e}^{-ikx}\, {\rm d}k \end{array} \end{eqnarray*}

(7)

for any function f(x), where i² = −1. The Fourier transform of the basis function Ψ_n(x; β) is then (setting f(x) = 0 for x < 0)

\begin{eqnarray*} \tilde{\Psi }_n(k; \beta) = (-1)^n \sqrt{\frac{2n\beta }{\pi }} \frac{(n\beta k -i)^{2n}}{[(n\beta k)^2 + 1]^{n+1}}\, . \end{eqnarray*}

(8)

Note that the real part of the Fourier transform is even for all n, while its imaginary part is odd for all n. Its modulus is a Lorentzian function centred on k = 0 whose width parameter Γ ≡ 2/(nβ) is inversely proportional to β (as typical for exponentially suppressed oscillations)

\begin{eqnarray*} \left|\tilde{\Psi }_n(k; \beta)\right| = \sqrt{\frac{2\pi }{n\beta }} \frac{1}{\pi } \frac{\Gamma /2}{k^2 + (\Gamma /2)^2}\, . \end{eqnarray*}

(9)

With conventions similar to those for the Fourier transform, the Laplace transform is given by

\begin{eqnarray*} \mathcal {L}(\Psi _n(x; \beta))(s) = (-1)^{n-1} 2\sqrt{n\beta } \frac{(n\beta s-1)^{n-1}}{(n\beta s+1)^{n+1}}\, . \end{eqnarray*}

(10)

Although equations (8) and (10) are not as simple as for Gaussian shapelets (the Fourier transform of a Gaussian shapelet is itself a Gaussian shapelet), closed forms exist for exponential shapelets that are easy to implement. The Fourier and Laplace transforms of 1D exponential shapelets are also well localized in frequency space (Fig. 2).

Figure 2.

Exponential shapelets are well localized in frequency space. This shows the modulus of the Fourier transform of the first few 1D exponential shapelets Ψ_n(x; β), for β = 1.

Open in new tab Download slide

2.2.3 Convolution

Convolution is an inevitable operation during signal acquisition, whereby an input (‘real’) signal undergoes the measurement apparatus’ transfer function. What is ultimately measured is the convolution of the input signal with the transfer function.

Let us consider the convolution of two functions f(x) and g(x), whose convolution product is

\begin{eqnarray*} h(x) \equiv (f \star g)(x) \equiv \int _{-\infty }^\infty {\rm d}x^{\prime } f(x-x^{\prime })g(x^{\prime }). \end{eqnarray*}

(11)

Following the same arguments as Refregier (2003), the functions can all be decomposed into exponential shapelets, perhaps with different scale sizes α, β, and γ. We can then relate the 1D exponential shapelet space coefficients h_{n; γ} to f_{m; β} and g_{l; α} via

\begin{eqnarray*} h_n = \sum _{m,l=1}^\infty C_{nml}f_m g_l, \end{eqnarray*}

(12)

where the convolution tensor is given by

\begin{eqnarray*} C_{nml}(\gamma , \alpha , \beta) &=& \sqrt{\frac{8\, \alpha \beta \gamma \, nml}{\pi } } \sum _{u=0}^{2m}\sum _{v=0}^{2l}\sum _{w=0}^{2n} i^{u+v+w} \binom{2m}{u}\nonumber \\ && \times \,\binom{2l}{v}\binom{2n}{w} (m\alpha)^u (l\beta)^v (n\gamma)^n \mathcal {I}_{u,v,w}^{m,l,n} \end{eqnarray*}

(13)

and |$\binom{n}{m}$| is the binomial coefficient. A proof of equation (13) is provided in Appendix A1.

The integral

\begin{eqnarray*} \mathcal {I}_{u,v,w}^{m,l,n} \!\! \equiv \!\!\! \int _{-\infty }^\infty \!\!\frac{ (-1)^w ~ k^{u+v+w}~~{\rm d}k}{\left[(m\alpha k)^2 + 1\right]^{m+1} \left[(l\beta k)^2 + 1\right]^{l+1} \left[(n\gamma k)^2 + 1\right]^{n+1}} \nonumber \\ \end{eqnarray*}

(14)

appearing in equation (13) is zero if u + v + w is odd (a proof is provided in Appendix A2). Hence, C_nml(γ, α, β) is always wholly real. Under some conditions, it can be estimated analytically and expressed as a sum of converging series (see equations A18, A22, and A29 in Appendix A3). Those conditions are often not met, but the integrand is a well-behaved function, so the integration can also be performed numerically. We find that numerical integration is most convenient if the function is cut into three segments (see equations A14, A20, and A27 in Appendix A3).

2.2.4 Rescaling

One-dimensional exponential shapelets obey the integral property

\begin{eqnarray*} \int _0^\infty \Psi _n(x;\beta) {\rm d}x = 2\sqrt{n\beta }. \end{eqnarray*}

(15)

Using this, it can be shown that under a rescaling x → x′ = ax, the coefficients f_n of a model (3) transform as

\begin{eqnarray*} f^{\prime }_{n;\beta } = a^{-1/2} f_{n;a\beta }\, , \end{eqnarray*}

(16)

and under f(x) → f′(x) = kf(x), the coefficients are themselves multiplied by k.

2.3 Shape characterization

Let f(x) = ∑_nf_nΨ_n(x; β) be a 1D object decomposed into 1D exponential shapelets. In terms of its shapelet coefficients, its integral (total ‘flux’) is

\begin{eqnarray*} F\equiv \int _{-\infty }^\infty f(x){\rm d}x = 2\sqrt{\beta } \sum _{n=1}^\infty \sqrt{n}\, f_n\, , \end{eqnarray*}

(17)

which can be readily shown from the integral property (15).

Provided F ≠ 0, its barycentre (centroid) is

\begin{eqnarray*} x_c \equiv \frac{1}{F}\int _{-\infty }^\infty x f(x){\rm d}x = \frac{4\sqrt{2}\beta ^{3/2}}{F} \sum _{n=1}^\infty n^2 f_n\, , \end{eqnarray*}

(18)

and it has a characteristic size

\begin{eqnarray*} r_{\rm c}^2 \equiv \frac{1}{F} \int _{-\infty }^\infty x^2 f(x){\rm d}x = \frac{4 \beta ^{5/2}}{F} \sum _{n=1}^\infty n^3 (2n^2+1)\, f_n\, . \end{eqnarray*}

(19)

This size can be used to determine the exponential decay rate of the object, for instance when modelling damped oscillations.

Note that the series (19) converges only if the amplitudes of the 1D exponential shapelet coefficients decrease faster than n⁻⁶, which may not always be the case. Care must therefore be taken to check for convergence when characterizing the shape of a feature using this technique.

2.4 Exponential shapelet modelling in practice

As shown above (equation 2), a 1D feature is straightforward to model for a given couple (n_max, β), where n_max is the maximum order of the truncated sum (2). For example, a linear least-square method can be efficiently used for this purpose. Then the model depends non-linearly on the two parameters n_max and β that can be optimized by iteratively minimizing the residuals between the observed feature and its model. This procedure was described at length, in the 2D case, by Massey & Refregier (2005).

2.5 Example applications

This section demonstrates three possible applications of 1D exponential shapelets. We start by modelling exponentially suppressed oscillations, which are measurements throughout experimental physics, including the response of accelerometers¹ onboard the space missions MICROSCOPE (Touboul, Métris & Rodrigues 2017), GRACE (Tapley et al. 2004), and GOCE (Rummel, Yi & Stummer 2011). We then discuss a potential application to unmodelled bursts in the analysis of gravitational waves. Exponential shapelets may also be convenient to model charge transfer inefficiency trailing due to radiation damage in spacebourne imaging detectors (Massey et al. 2010), although we do not explore that further here.

2.5.1 Cleaning accelerometer time series data, and modelling an experiment’s response function

MICROSCOPE tested the weak equivalence principle (WEP) by precisely measuring the differential acceleration experienced by two concentric cylindrical test masses onboard a drag-free satellite in Low Earth Orbit.² In theory, any non-zero difference at a well-defined frequency f_EP (which depends on the orbit and attitude characteristics of the satellite) is a signature of violation of the WEP.

In practice, many transients are apparent in MICROSCOPE data (the upper panel of Fig. 3 shows a typical high-signal-to-noise example). These transients are generally caused by crackles of the satellite’s coating (because of temperature variations), crackles of the satellite’s gas tanks (as their pressure decreases as the gas is consumed), or impacts with micro-meteorites. Such transients occupy frequencies higher than f_EP, so they do not directly impact a possible WEP violation signal. However, it is necessary to detect and mask them in measured time-series, and then either reconstruct the corresponding ‘missing data’ (Bergé et al. 2015; Pires et al. 2016) to allow for a least-square fit of the expected WEP violation signal in the frequency domain (Touboul et al. 2017), or adapt the maximum likelihood technique to take missing data into account (Baghi et al. 2015, 2016). Existing techniques are suboptimal, as they may affect the noise characteristics. Moreover, transients could be considered as conveying useful information. They are created by an external impulse; if this is assumed to be instantaneous (Dirac), the shape and relaxation time of the observed signal are a measurement of MICROSCOPE’s transfer function.

Figure 3.

Test data for 1D exponential shapelets. Upper panel: a transient in MICROSCOPE accelerometer time series, as observed (black) and reconstructed from a 1D exponential model (red), with model residuals (green). Lower panel: convergence of the model as we increase the order of the model n_max.

Open in new tab Download slide

1D exponential shapelets are well-matched to these transient signals. The upper panel of Fig. 3 shows a transient in the time domain: the observed data are shown in black, a model (fitted between t = 10 s and t = 25 s) and its 68 per cent confidence interval are shown in red, and residuals consistent with noise are shown in green. This model uses n_max = 15, meaning that it has compressed 60 data points in the time domain description (15 s sampled at 4 Hz) into 15 shapelet coefficients. None the less, the model has rich, empirical freedom to capture multiple response modes of the instrument’s complex structure. The lower panel of Fig. 3 shows the convergence of the model as we increase n_max, quantified as the square residuals between the observed transient and the model. Most of the information is contained in coefficients 3 ≤ n ≤ 13.

An extension to this process will be presented in a future paper. By fitting 1D exponential shapelet coefficients to many transients, it is possible to model temporal variations in relaxation time of MICROSCOPE’s instrument via the parameter β or the exponential envelope scale r_rms. Interpolating these models then yields a model of the transfer function at any time, to either gain insight into the instrument’s performance or correct (deconvolve) signals with an ‘inverse transfer function’ transform. Note that this procedure is similar to techniques developed for 2D astronomical image processing, where a point spread function (PSF) is determined from the shape of stars then interpolated across the data (see e.g. Bergé et al. 2012; Gentile, Courbin & Meylan 2013; Kilbinger 2015).

2.5.2 Characterization of perturbing signals in space-borne geodesy missions

The GRACE mission (Tapley et al. 2004) revolutionized geodesy by measuring the Earth’s gravitational field with unprecedented precision. Two satellites followed each other on the same orbit, monitoring the distance between them via microwave ranging. In theory, any variation in their relative speed or distance can be ascribed to variations in the Earth’s geopotential.

In practice, the satellites were also subject to non-gravitational forces. These were measured by ultrasensitive accelerometers, for removal during post-processing (Flury et al. 2008; Peterseim et al. 2010). Peterseim (2014) and Peterseim et al. (2014) modelled transients (known as ‘spikes’) in accelerometer data using a piecewise function made from the derivative of a Gaussian, a third-order polynomial, and a damped oscillation. Some were successfully classified as either due to the satellite’s heaters or due to activation of its magnetic torquers – but no physical origin could be assigned to others, known as ‘twangs’.

Tentative correlations of twangs with the position of the satellite along its orbit hint at a possible geophysical origin. Both categories of twangs are compactly represented as 1D exponential shapelets, so we will investigate in a future paper whether these provide a cleaner set of shape parameters to categorize and understand their origin.

2.5.3 Unmodelled bursts and glitches in gravitational wave data analysis

A wealth of methods have been developed to search for, characterize, and classify unmodelled bursts and instrumental glitches in searches for gravitational waves (see e.g. Cornish & Littenberg 2015; Powell et al. 2015, 2017; Principe & Pinto 2017, and references therein). Glitches are often modelled using ‘sine-Gaussian waveforms’ (Principe & Pinto 2017). These have similar properties to Gaussian shapelets, although shapelets can encode more details. One-dimensional exponential shapelets could further optimize the data compression of complex glitch shapes that often exhibit damped oscillations.

Shapelets might therefore improve the detection, characterization, and classification of instrumental glitches in gravitational wave detectors. Indeed, bursts and glitches are usually detected in the time–frequency domain, which is easily reproduced in Gaussian shapelet-time space. If 1D exponential shapelets more efficiently model the information in a glitch, exponential shapelet-time space would be even better. We will investigate in a future paper whether glitches can be detected by scanning a matched filter and correlating the measured signal with a 1D exponential shapelet model while leaving β (and possibly n_max) free.

2.6 Comparison with Gaussian shapelets

The MICROSCOPE transient modelled as exponential shapelets (n_max = 15) in Fig. 3 is modelled as Gaussian shapelets (n_max = 18) in Fig. 4. Achieving a similar precision of fit requires more coefficients, and creates more small-scale artefacts: The data are overfitted near the centre of the model, and underfitted at the extremes. We also find that its coefficients are highly covariant, with the good fit produced by large positive and negative basis functions almost precisely cancelling each other. This is much less apparent for the 1D exponential shapelet coefficients. Indeed, we find that estimating the model again on points different from those observed data fails to reproduce the overall shape of the transient, whereas the 1D exponential shapelet model is itself predictive.

Figure 4.

Gaussian shapelet model of the MICROSCOPE transient whose 1D exponential shapelet model is shown in Fig. 3.

Open in new tab Download slide

In this example, it appears that 1D exponential shapelets outperform Gaussian shapelets. This is because their perturbations around an exponential decay are better matched to the underlying signal, and because their wider extent spreads information density more evenly. In general, the type of shapelet to choose should depend on the decay rate of the target function.

3 2D EXPONENTIAL SHAPELETS

3.1 Definition

The quantum mechanics of a hydrogen atom restricted to two dimensions is well studied (see e.g. Zaslow & Zandler 1967; Yang et al. 1991; Chaos-Cador & Ley-Koo 2007). The bound states of the electron (i.e. eigenfunctions of the Hamiltonian) provide a natural set of basis functions to represent a bounded distribution. After renormalizing these, we define the 2D exponential shapelet basis functions

\begin{eqnarray*} \Psi _{n,m}(r,\phi ;\beta) &=& (-1)^n \sqrt{\frac{2}{\beta \pi (2n+1)^{3}}\, \frac{(n-|m|)!}{(n+|m|)!}} \nonumber\\ && \times \left(\frac{2r}{\beta (2n+1)} \right)^{|m|} L_{n-|m|}^{2|m|}\left(\frac{2r}{\beta (2n+1)} \right) \nonumber\\ && \times \exp \left(-\frac{r}{\beta (2n+1)} \right) \exp (-im\phi), \end{eqnarray*}

(20)

for n ≥ 0, where |$L_i^j\left(x\right)$| are generalized Laguerre polynomials. The (− 1)ⁿ term is used to ensure that the integral of each basis function is positive, but this is otherwise the form used in quantum theory. In terms of quantum mechanics, n is the principal quantum number (the eigenfunction energy level) and m the magnetic quantum number, which takes integer values between −n and n. The characteristic scale β is linked to the Bohr radius.

These functions form an orthonormal set of basis functions in the L²([0, ∞[ × [0, 2π[, 〈 ·, ·〉) space equipped with inner product |$\langle \Psi _{n,m}(r,\phi), \Psi _{n^{\prime },m^{\prime }}(r,\phi)\rangle = \int _0^{2\pi } \int _0^\infty \Psi _{n,m}(r,\phi)\, \Psi ^*_{n^{\prime },m^{\prime }}(r,\phi) \, r\, {\rm d}r\, {\rm d}\phi = \delta _{nn^{\prime }}\delta _{mm^{\prime }}$|⁠.

Any localized function f(r, ϕ) can be uniquely decomposed into a weighted sum of these basis functions

\begin{eqnarray*} f(r,\phi)=\sum _{n=0}^\infty \sum _{m=-n}^n f_{n,m}\, \Psi _{n,m}(r,\phi ;\beta), \end{eqnarray*}

(21)

where the 2D exponential shapelet coefficients are given by

\begin{eqnarray*} f_{n,m}=\int _0^{2\pi } \int _0^\infty f(r,\phi)\, \Psi _{n,m}(r,\phi ;\beta)\, r\, {\rm d}r\, {\rm d}\phi . \end{eqnarray*}

(22)

Using Bessel’s inequality like in the 1D case, we find that series (21) converges, and the coefficients f_n,m must vanish when n and m increase. For a real function |$f(r,\phi)\in \mathbb {R}$|⁠, |$f_{n,m}=f_{n-m}^*$|⁠. In this case, truncating the series at n ≤ n_max results in n_coeffs = (n_max + 1)² independent coefficients.

It may also be possible to define elliptical 2D exponential shapelets by applying a shear transformation to the coordinate system, as Nakajima & Bernstein (2007) suggested for Gaussian shapelets. This preserves their orthonormality, and can increase their rate of data compression, at the cost of two additional non-linear parameters to specify the shear.

3.2 Properties

3.2.1 Maximum and minimum effective scales

The first few 2D exponential shapelet basis functions are illustrated in Fig. 5, and their radial component is shown in Fig. 6, defined via Ψ_n,m(r, ϕ; β) ≡ R_n,m(r; β)exp (−imϕ). Their resemblance with Gaussian polar shapelets is striking (see fig. 2 of Massey & Refregier 2005). However, due to their exponential kernel, exponential shapelets are both more peaky and further spread out than Gaussian shapelets. The size difference between the lowest-order basis function Ψ₀₀ and even a low-order one like Ψ₄₀ is striking. It is this behaviour that will help to describe spatially extended features.

Figure 5.

The first few 2D exponential shapelet basis functions: real part (left) and imaginary part (right). Red is positive; blue is negative. The normalization of both the colour scale and the spatial scale is arbitrary, but is the same in every box.

Open in new tab Download slide

Figure 6.

Radial component of the first few 2D exponential shapelet basis functions, R_n,m(r). Functions with the same n are of the same colour (blue for n = 0, green for n = 1, red for n = 2, black for n = 3). Functions with m ≠ 0 are shown as dashed lines.

Open in new tab Download slide

Generalizing the 1D case (Section 2.2.1), we define the largest effective scale in a 2D exponential shapelet model as

\begin{eqnarray*} \theta _{\rm max} &\approx & \sqrt{\frac{\textstyle\iint _{\mathbb {R}^2} r^2\, \Psi _{n_{\rm max},0}(r,\phi ;\beta)\, r{\rm d}r{\rm d}\phi }{\textstyle\iint _{\mathbb {R}^2} \Psi _{n_{\rm max},0}(r,\phi ;\beta)\, r{\rm d}r{\rm d}\phi }} \nonumber \\ &=& \beta (2n_{\rm max}+1)\sqrt{2(2n_{\rm max}^2+2n_{\rm max}+3)}\, . \end{eqnarray*}

(23)

Again empirically, the smallest resolved scales are roughly constant,

\begin{eqnarray*} \theta _{\rm min} \approx \beta . \end{eqnarray*}

(24)

The range of scales included in a 2D exponential shapelet model as n_max → ∞ is |$(\theta _{\rm max}/\theta _{\rm min})^2\propto n_{\rm max}^4\propto n_{\rm coeffs}^2$|⁠. The resolution of an exponential shapelet model is greatest near the origin, and the information density is concentrated there. This behaviour is different from the Gaussian shapelets, where resolution is more spatially uniform, and the information density is constant, with |$(\theta _{\rm max}/\theta _{\rm min})^2=n_{\rm max}^2+1\propto n_{\rm coeffs}$|⁠.

3.2.2 Fourier transformation and convolution

Using the same convention as in the 1D case, it can be shown that the Fourier transform of 2D exponential shapelets is given by

\begin{eqnarray*} \tilde{\Psi }_{n,m}(\vec{k};\beta) = 2\pi i^m {\rm e}^{-im\phi _k} \left[ F_m(k)\right]_{n,\beta }, \end{eqnarray*}

(25)

where |$k=|\vec{k}|$|⁠, ϕ_k is the angle between the |$\vec{k}$| direction and the ϕ = 0 direction in polar space, and [F_m(k)]_n,β is the Hankel transform of the R_nm(r, β) radial function. Consequently, the convolution tensor involved in the convolution of two objects modelled in 2D exponential shapelets is the integral of products of Hankel transforms of R_n,m functions. We were not able to find an analytic form for this Hankel transform, so it must instead be computed numerically.

3.2.3 Coordinate transforms

It can be useful to know the response of the basis functions to 2D linear coordinate transforms: either to know how to mix the coefficients of a shapelet decomposition in order to perform that transform, or to form combinations of coefficients that are invariant under some transforms. This was used to construct estimators of the shear distortion applied to images of galaxies by the effect of weak gravitational lensing (Kuijken 2006; Massey et al. 2007).

A convenient shortcut to calculating those transforms for Gaussian shapelets was provided by the ladder operators associated with the quantum mechanical harmonic oscillator (Refregier & Bacon 2003). Unfortunately, we have not yet found a useful form of the ladder operators for the 2D hydrogen atom. If it becomes necessary to perform linear transformations on 2D exponential shapelets, it will be necessary to apply the transforms manually to the basis functions (a long and arduous task, but one that is guaranteed to yield a closed form, because the basis is complete).

3.3 Object shape measurement

Once a feature has been decomposed into 2D exponential shapelets, its coefficients can be used to construct characteristic measurements of its shape. In this section, we derive expressions for an object’s (azimuthally averaged) radial profile, flux, centroid, unweighted quadrupoles, size, and ellipticity, in terms of its shapelet coefficients. We have not attempted it here, but it should also be possible to develop an exponential-shapelet classification of galaxy morphologies via e.g. their concentration, asymmetry, and symmetry (Conselice, Bershady & Jangren 2000).

3.3.1 Radial profile

Azimuthally averaging an object’s signal yields its mean radial profile

\begin{eqnarray*} \bar{f}(r)\equiv \frac{1}{2\pi }\int _0^{2\pi }f(r,\phi)\, {\rm d}\phi . \end{eqnarray*}

(26)

Noting that all m ≠ 0 basis functions average to zero, the radial profile reduces to

\begin{eqnarray*} \bar{f}(r)=\sum _{n=0}^\infty f_{n,0} \Psi _{n,0}(r;\beta), \end{eqnarray*}

(27)

where the (rotationally invariant) m = 0 basis functions are

\begin{eqnarray*} \Psi _{n,0}(r;\beta)&=& (-1)^n \frac{2}{\beta \sqrt{2\pi }} (2n+1)^{-3/2} \nonumber\\ &&\times \exp \left(-\frac{r}{\beta (2n+1)}\right) L_n\left(\frac{r}{\beta (n+\frac{1}{2})} \right). \end{eqnarray*}

(28)

Equation (26) is identical to the equivalent derivation for Gaussian shapelets (see equation 16 of Massey & Refregier 2005), up to the fact that m = 0 basis functions with odd n do not exist in the Gaussian case.

3.3.2 Flux

Integrating the signal inside a circular aperture of radius R yields its ‘flux’

\begin{eqnarray*} F_R \equiv \int _0^{2\pi } \int _0^R f(r,\phi)\, r\, {\rm d}r\, {\rm d}\phi . \end{eqnarray*}

(29)

To evaluate this integral, it is useful to note that, once again, all m ≠ 0 basis functions cancel out to zero during integration over ϕ, and also a closed form for the generalized Laguerre polynomials,

\begin{eqnarray*} L_n^\alpha (x) = \sum _{k=0}^n (-1)^k \binom{n+\alpha }{n-k} \frac{x^k}{k!}. \end{eqnarray*}

(30)

Using this, it can be shown that

\begin{eqnarray*} F_R &=& 2 \sqrt{2\pi } \beta \sum _{n=0}^\infty f_{n,0}\, (2n+1)^{1/2} \nonumber \\ && \times \,\sum _{k=0}^n \frac{2^k (-1)^{n+k}}{k!} \binom{n}{k} ~ \gamma \left(k+2, \frac{R}{\beta (2n+1)}\right), \end{eqnarray*}

(31)

where γ(y, x) is the lower incomplete gamma function. Extrapolating F_R to a large radius, and taking the limit R → ∞, we obtain

\begin{eqnarray*} F \equiv \iint _{\mathbb {R}^2} f(r,\phi)\, r{\rm d}r{\rm d}\phi ~=~ 2\sqrt{2\pi } \beta \sum _{n=0}^\infty (2n+1)^{3/2}f_{n,0}\, . \end{eqnarray*}

(32)

3.3.3 Centroid

Similarly, it can be shown that the unweighted centroid (x_c, y_c), defined by

\begin{eqnarray*} x_\mathrm{ c}+i y_\mathrm{ c} \equiv \frac{\textstyle\iint _{\mathbb {R}^2}(x+i y) f(x,y)\, {\rm d}x{\rm d}y}{\textstyle\iint _{\mathbb {R}^2} f(x,y)\, {\rm d}x{\rm d}y}, \end{eqnarray*}

(33)

is, in terms of 2D exponential shapelet coefficients,

\begin{eqnarray*} x_\mathrm{ c}+i y_\mathrm{ c} = -\frac{4\sqrt{2\pi }\beta ^2}{F} \sum _{n=1}^\infty \sqrt{n(n+1)(2n+1)^{5}}\, f_{n,1}\, . \end{eqnarray*}

(34)

3.3.4 Quadrupole moments

The unweighted quadrupole moments of an object

\begin{eqnarray*} F_{ij} = \iint _{\mathbb {R}^2}x_ix_j f(x,y)\, {\rm d}x{\rm d}y \end{eqnarray*}

(35)

can be used to define quantities such as its rms size and ellipticity (see below). They are, in terms of 2D exponential shapelet coefficients,

\begin{eqnarray*} F_{11} &=& 2\sqrt{2\pi } \beta ^3 \sum _{n=0}^\infty (2n+1)^{7/2} \nonumber\\ && \times \left[(2n^2+2n+3)f_{n,0} + 2\sqrt{\frac{(n+2)!}{(n-2)!}}\, f_{n,2} \right] \end{eqnarray*}

(36)

\begin{eqnarray*} F_{22} &=& 2\sqrt{2\pi } \beta ^3 \sum _{n=0}^\infty (2n+1)^{7/2} \nonumber\\ && \times \left[(2n^2+2n+3)f_{n,0} - 2\sqrt{\frac{(n+2)!}{(n-2)!}}\, f_{n,2} \right] \end{eqnarray*}

(37)

\begin{eqnarray*} F_{12}\! =\! F_{21}\! =\! -4\sqrt{2\pi } \beta ^3 i \! \sum _{n=0}^\infty (2n+1)^{7/2}\! \sqrt{\frac{(n+2)!}{(n-2)!}}\, f_{n,2} \, . \end{eqnarray*}

(38)

3.3.5 Size and ellipticity

Using equations (35)–(37), expressions can be obtained to quantify a feature’s size

\begin{eqnarray*} R^2 & \equiv & \frac{F_{11}+F_{22}}{F} \nonumber \\ & = & \frac{4\sqrt{2\pi } \beta ^3}{F} \sum _{n=0}^\infty (2n+1)^{7/2}\, (2n^2+2n+3)\, f_{n,0} \end{eqnarray*}

(39)

and ellipticity

\begin{eqnarray*} \varepsilon & \equiv & \frac{F_{11}-F_{22}+2iF_{12}}{F_{11}+F_{22}} \nonumber \\ & = & \frac{8\sqrt{2\pi }\beta ^3}{FR^2} \sum _{n=2}^\infty (2n+1)^{7/2} \sqrt{\frac{(n+2)!}{(n-2)!}}\, f_{n,2}\, . \end{eqnarray*}

(40)

In this complex notation with ε ≡ |ε|cos (2ϕ) + i|ε|sin (2ϕ), ε = 0 denotes rotational invariance, and a positive real (imaginary) component denotes elongation along (at 45° to) the x-axis.

3.3.6 Convergence

Expressions for the shape estimators resemble those obtained for polar Gaussian shapelets in section 6 of Massey & Refregier (2005). In particular, a feature’s radial profile, flux, and size are encoded within the n = 0 coefficients, while the centroid is described by the n = 1 coefficients and ellipticity by the n = 2 coefficients. The expressions also contain the same power of β as those for Gaussian shapelets, because that encodes information about the units in which the data is expressed.

However, shape measures for 2D exponential shapelets contain higher powers of n than those for 2D Gaussian shapelets, and will converge more slowly. Some of this is due to the normalization of the basis functions against an inner product. In the Gaussian shapelet case, this happens to yield an integral of basis functions (equation 31) that is independent of n, while in exponential shapelets it is ∝n^3/2. That power of n could have been included in the basis functions and removed from the coefficients – although doing so would merely make them look more stable rather than actually altering convergence. In the limit that R → ∞, convergence of the flux estimator requires |f_n,0| to decrease faster than n^−5/2; convergence of centroid requires |f_n,1| to decrease faster than n^−11/2; and convergence of size and ellipticity requires |f_n,0| and |f_n,2| to decrease faster than n^−13/2. This may not be a problem, since we expect 2D exponential shapelets to converge faster than shapelets (i.e. to have a lower n_max for a given galaxy). However, if this does not hold true, e.g. due to measurement noise, care should be taken to restrict the decomposition and measurement of photometry inside a finite aperture, or to truncate the 2D exponential shapelet model at finite n_max.

3.4 Example application

This section demonstrates the use of 2D exponential shapelets to model the shapes of galaxies seen outside our own Milky Way (they may also be convenient to model the shapes of gamma-ray events; S. Pires, private communication). We examine the convergence speed of 2D exponential shapelets, and compare their performance to that of 2D Gaussian shapelets.

3.4.1 Distant galaxies in astronomical imaging

Galaxies are collections of (hundreds of billions of) stars, with a combined surface brightness that peaks sharply near the centre, and decreases to large radii in a way that is often exponential. Four galaxies in the COSMOS survey, (Scoville et al. 2007a,b) observed by the Hubble Space Telescope (Koekemoer et al. 2007), are shown in Fig. 7. As mentioned in Section 2.4 for the 1D case, we follow the algorithm described by Massey & Refregier (2005) to choose the best non-linear parameters of the models –n_max, β, and centroid). Note that galaxies in this figure were modelled into shapelets with the ‘diamond’ truncation scheme introduced in Massey et al. (2007), which is intended to reduce small-scale oscillations by ignoring the higher m terms of the decomposition. In this case, the number of coefficients is divided by 2, allowing also a better compression of the information.

Figure 7.

Test images of four galaxies from the Hubble Space Telescope COSMOS survey (Scoville et al. 2007a,b). For each galaxy, we show the observed data (upper left); the 2D exponential shapelet model (upper middle) and its residuals (upper right); the Gaussian shapelet model (lower middle) and its residuals (lower right). Also shown for each galaxy are the maximum order of decomposition n_max and the total number of coefficients n_coeffs in each model. For a given galaxy, the grey-scales are the same in all panels.

Open in new tab Download slide

In all four cases, 2D exponential shapelets require fewer coefficients (lower n_max; here n_coeffs = 85) to provide a model whose residuals are consistent with the noise than Gaussian shapelets (here n_coeffs = 121). Even more interestingly, exponential shapelets tend to provide a (visually) cleaner model of the central region of the galaxies. The top left galaxy has an irregular morphology, and the 2D exponential shapelet model shows fewer artefacts than the polar shapelet model. For an elliptical (top right) or edge-on spiral galaxy (bottom left), the ringing common in Gaussian shapelets disappears with 2D exponential shapelets: They do not possess high-frequency ripples at large radii that need large coefficients alternating in sign for them to cancel. Even for highly eccentric galaxies (bottom right), 2D exponential shapelets provide cleaner models than polar shapelets, especially in the outskirts of the galaxies, as the exponential wings are fundamentally a better match to the underlying signal.

3.4.2 Convergence speed

To assess the ability of 2D exponential shapelets to model galaxies, we require noise-free images whose true properties are known. We simulate elliptical galaxies with a Sersic radial profile

\begin{eqnarray*} \ln \left[\frac{I(r)}{I(r_\mathrm{ e})}\right] = -b_\mathrm{ n}\left[\left(\frac{r}{r_\mathrm{ e}}\right)^{1/n_\mathrm{ s}}-1\right], \end{eqnarray*}

(41)

where I(r) is the galaxy’s s profile, n_s is the Sersic index, b_n is a normalization constant, and r_e is the effective radius (which contains 50 per cent of the galaxy’s flux). We simulate three types of galaxies: exponential (n_s = 1), intermediate (n_s = 2), and de Vaucouleurs (n_s = 4).

Fig. 8 compares the radial profiles fitted with 2D exponential shapelets (upper row) or Gaussian shapelets (lower row), for different maximum orders of decomposition n_max, and coarsely optimised β (this could be further optimized at low n_max in the presence of noise). The normalization of the ordinate depends on b_n and pixelization, and so is arbitrary; only its relative scaling is informative. For all three galaxy types, 2D exponential shapelets greatly outperform Gaussian shapelets. A decomposition to n_max = 2 or 4 is sufficient to model an exponential galaxy. Depending on the signal to noise of real data, a decomposition to n_max = 6 may be sufficient for an intermediate n_s = 2 galaxy. It is more challenging to model a de Vaucouleurs galaxy, and even exponential shapelets require n_max = 12 to capture its extended tails three orders of magnitude fainter than the peak. However, they do so without the oscillatory ‘ringing’ introduced by the cancellation of positive and negative basis functions in Gaussian shapelets.

Figure 8.

Exponential (upper panels) and Gaussian (lower panels) shapelet models’ speed to converge on the profile of different types of (noiseless) elliptical galaxies: exponential (n_s = 1, left), intermediate n_s = 2 (centre), and de Vaucouleurs galaxy (n_s = 4, right). In each panel, squares represent the profile of the galaxy estimated in pixel space.

Open in new tab Download slide

To illustrate the convergence speed, Fig. 9 shows the mean square residual of the models of the same galaxies as in Fig. 8, as a function of the maximum shapelet order n_max (upper row), or the total number of coefficients n_coeffs. The normalization of the ordinate is arbitrary, and only its relative scaling is informative. Once again, it is clear that 2D exponential shapelets (solid lines) outperform Gaussian shapelets (dashed lines). For exponential and intermediate galaxies (left-hand and middle columns), the residual of the 2D exponential shapelet model decreases consistently, and ends up several orders of magnitude better than that of Gaussian shapelets. The enhanced performance of exponential shapelets is even more prominent for a de Vaucouleurs galaxy. The limitation for exponential shapelets in this case is numerical precision of our algorithm to pixelate the basis functions. At high n, this breaks down, causing spurious residuals at the very centre of the model, shown as an artificial upturn in the bottom right-hand panel of Fig. 9. This could be circumvented by more accurate pixelization, or by imposing limits on θ_min.

Figure 9.

Convergence speed as a function of shapelets’ n_max (upper panels) and number of coefficients (lower panels), for exponential (solid lines) and Gaussian (dashed lines) shapelets, for the same Sersic galaxies as in Fig. 8 (n_s = 1, 2, 4 from left to right).

Open in new tab Download slide

4 CONCLUSION

In this paper, we introduced exponential shapelets, a family of orthonormal basis functions that can be efficiently used to model 1D and 2D objects. They borrow many concepts from the original Gaussian shapelets, but are perturbations around a decaying exponential rather than a Gaussian, and more efficiently compress the information in damped 1D oscillations or centrally concentrated 2D features. In particular, they can simultaneously capture the central peak and the extended wings of galaxies in astronomical imaging, thereby solving the main criticism levelled at Gaussian shapelets in the field for which they were originally intended.

Modelling a feature using exponential shapelets first requires a choice of characteristic scale size β and expansion order n_max; these can be selected using an non-linear optimization method developed for Gaussian shapelets (Massey & Refregier 2005). The function is then decomposed into a weighted sum of exponential shapelet basis functions using linear regression. A least-squares fit is often sufficient, although in some cases modelling a transient may be stabilized by regularization techniques such as sparsity constraints.

Once a feature is described as a weighted sum of exponential shapelets, simple combinations of the weights yield expressions for its total flux, centroid, size, and ellipticity. These are similar to those for Gaussian shapelets. However, convolution and deconvolution in exponential shapelets are significantly more difficult than in the Gaussian case, with the convolution tensor being time-consuming to calculate via numerical integration.

We described possible applications of exponential shapelets in several fields of experimental and observational science. One-dimensional exponential shapelets can be used to model (and subtract or understand) spurious transients in time-series, such as measurements by the MICROSCOPE or GRACE satellites, or the LIGO search for gravitational waves. Thanks to their exponential decay, 2D exponential shapelets overcome the main criticism aimed at Gaussian shapelets (though at the cost of losing some simplicity) and can be used to measure the brightness, shape, and size of distant galaxies in astronomical imaging. The characteristics of exponential shapelets should make them well suited to data compression and analysis in a wide range of fields; their convenient mathematical properties should see them frequently adopted.

ACKNOWLEDGEMENTS

The concepts in this paper matured over many years, after a preliminary study to generalize Gaussian polar shapelets for Hubble Space Telescope data, for which we acknowledge support from NASA through grant HST-AR-11747 – and which was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. We are indebted to Barnaby Rowe for researching Kummer functions and the literature of the 2D hydrogen atom. We also acknowledge helpful and motivating discussions with Jason Rhodes, Richard Ellis, Alexandre Refregier, Sandrine Pires, James Nightingale, Bernard Foulon, Bruno Christophe, Gilles Métris, Manuel Rodrigues, Alain Robert, and Nicolas Touquoy. This work makes use of technical data from the CNES-ESA-ONERA-CNRS-OCA Microscope mission, and has received financial support from ONERA and CNES. JB acknowledges the financial support of the UnivEarthS Labex programme at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02) and of CNES through the APR programme (LISA project). RM is supported by a Royal Society University Research Fellowship and through STFC grant ST/N001494/1.

Footnotes

1

http://www.onera.fr/en/dphy

2

The mission ended on 2018 October 16; data analysis is still underway.

REFERENCES

Akdeniz

T. J.

,

Lizotte

D.

,

Mohieddin Abukhdeir

N.

,

2018

,

IOP Nanotechnology

,

30

,

075703

Crossref

Search ADS

Amara

A.

,

Quanz

S. P.

,

2012

,

MNRAS

,

427

,

948

10.1111/j.1365-2966.2012.21918.x

Crossref

Search ADS

Apostolopoulos

G.

,

Koutras

A.

,

Christoyianni

I.

,

Dermatas

E.

,

2017

,

IEEE 25th European Signal Processing Conference (EUSIPCO)

.

IEEE

, p.

56

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Baghi

Q.

,

Métris

G.

,

Bergé

J.

,

Christophe

B.

,

Touboul

P.

,

Rodrigues

M.

,

2015

,

Phys. Rev. D

,

91

,

062003

10.1103/PhysRevD.91.062003

Crossref

Search ADS

Baghi

Q.

,

Métris

G.

,

Bergé

J.

,

Christophe

B.

,

Touboul

P.

,

Rodrigues

M.

,

2016

,

Phys. Rev. D

,

93

,

122007

10.1103/PhysRevD.93.122007

Crossref

Search ADS

Bergé

J.

,

Price

S.

,

Amara

A.

,

Rhodes

J.

,

2012

,

MNRAS

,

419

,

2356

10.1111/j.1365-2966.2011.19888.x

Crossref

Search ADS

Bergé

J.

,

Pires

S.

,

Baghi

Q.

,

Touboul

P.

,

Métris

G.

,

2015

,

Phys. Rev. D

,

92

,

112006

10.1103/PhysRevD.92.112006

Crossref

Search ADS

Bernstein

G. M.

,

Jarvis

M.

,

2002

,

AJ

,

123

,

583

10.1086/338085

Crossref

Search ADS

Bezrodnykh

S. I.

,

2016

,

Math. Notes

,

100

,

318

10.1134/S0001434616070282

Crossref

Search ADS

Birrer

S.

,

Amara

A.

,

Refregier

A.

,

2015

,

ApJ

,

813

,

102

10.1088/0004-637X/813/2/102

Crossref

Search ADS

Chaos-Cador

L.

,

Ley-Koo

E.

,

2007

,

Int. J. Quantum Chem.

,

107

,

12

10.1002/qua.21052

Crossref

Search ADS

Conselice

C. J.

,

Bershady

M. A.

,

Jangren

A.

,

2000

,

ApJ

,

529

,

886

10.1086/308300

Crossref

Search ADS

Cornish

N. J.

,

Littenberg

T. B.

,

2015

,

Class. Quantum Gravity

,

32

,

135012

10.1088/0264-9381/32/13/135012

Crossref

Search ADS

Desvignes

G.

et al. .,

2016

,

MNRAS

,

458

,

3341

10.1093/mnras/stw483

Crossref

Search ADS

Ellis

J. A.

,

Cornish

N. J.

,

2016

,

Phys. Rev. D

,

93

,

084048

10.1103/PhysRevD.93.084048

Crossref

Search ADS

Flury

J.

,

Bettadpur

S.

,

Tapley

B. D.

,

2008

,

Adv. Space Res.

,

42

,

1414

10.1016/j.asr.2008.05.004

Crossref

Search ADS

Gentile

M.

,

Courbin

F.

,

Meylan

G.

,

2013

,

A&A

,

549

,

A1

10.1051/0004-6361/201219739

Crossref

Search ADS

Hasanov

A.

,

Srivastava

H. M.

,

2007

,

Comput. Math. Appl.

,

53

,

1119

10.1016/j.camwa.2006.07.007

Crossref

Search ADS

Hoekstra

H.

,

Wu

Y.

,

Udalski

A.

,

2005

,

ApJ

,

626

,

1070

10.1086/430299

Crossref

Search ADS

Holbrey

R.

,

2006

,

Dimension Reduction Algorithms for Data Mining and Visualization

.

University of Leeds

,

Edinburgh

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Jiménez-Teja

Y.

,

Benítez

N.

,

2012

,

ApJ

,

745

,

150

10.1088/0004-637X/745/2/150

Crossref

Search ADS

Kilbinger

M.

,

2015

,

Rep. Prog. Phys.

,

78

,

086901

10.1088/0034-4885/78/8/086901

Crossref

Search ADS

PubMed

Koekemoer

A. M.

et al. .,

2007

,

ApJS

,

172

,

196

10.1086/520086

Crossref

Search ADS

Kuijken

K.

,

2006

,

A&A

,

456

,

827

10.1051/0004-6361:20054794

Crossref

Search ADS

Lentati

L.

,

Alexander

P.

,

Hobson

M. P.

,

2015

,

MNRAS

,

447

,

2159

10.1093/mnras/stu2611

Crossref

Search ADS

Massey

R.

,

Refregier

A.

,

2005

,

MNRAS

,

363

,

197

10.1111/j.1365-2966.2005.09453.x

Crossref

Search ADS

Massey

R.

,

Rowe

B.

,

Refregier

A.

,

Bacon

D. J.

,

Bergé

J.

,

2007

,

MNRAS

,

380

,

229

10.1111/j.1365-2966.2007.12072.x

Crossref

Search ADS

Massey

R.

,

Stoughton

C.

,

Leauthaud

A.

,

Rhodes

J.

,

Koekemoer

A.

,

Ellis

R.

,

Shaghoulian

E.

,

2010

,

MNRAS

,

401

,

371

10.1111/j.1365-2966.2009.15638.x

Crossref

Search ADS

Melchior

P.

,

Böhnert

A.

,

Lombardi

M.

,

Bartelmann

M.

,

2010

,

A&A

,

510

,

A75

Crossref

Search ADS

Nakajima

R.

,

Bernstein

G.

,

2007

,

AJ

,

133

,

1763

10.1086/511957

Crossref

Search ADS

Ngan

W.

,

van Waerbeke

L.

,

Mahdavi

A.

,

Heymans

C.

,

Hoekstra

H.

,

2009

,

MNRAS

,

396

,

1211

10.1111/j.1365-2966.2009.14787.x

Crossref

Search ADS

Nieto

M. M.

,

1979

,

Am. J. Phys.

,

47

,

1067

10.1119/1.11976

Crossref

Search ADS

Nú nez-Yépez

H. N.

,

Salas-Brito

A. L.

,

Solis

D. A.

,

2011

,

Phys. Rev. A

,

83

,

064101

10.1103/PhysRevA.83.064101

Crossref

Search ADS

Nú nez-Yépez

H. N.

,

Salas-Brito

A. L.

,

Solis

D. A.

,

2014

,

Phys. Rev. A

,

89

,

049908

10.1103/PhysRevA.89.049908

Crossref

Search ADS

Palma

G.

,

Raff

U.

,

2006

,

Can. J. Phys.

,

84

,

787

10.1139/p06-072

Crossref

Search ADS

Peterseim

N.

,

2014

,

PhD Thesis

.

Technische Universität

,

München

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Peterseim

N.

,

Jakob

F.

,

Schlicht

A.

,

2010

,

38th COSPAR Scientific Assembly

,

38

,

4

Peterseim

N.

,

Schlicht

A.

,

Flury

J.

,

Dahle

C.

,

2014

, in

Flechtner

F.

,

Sneeuw

N.

,

Schuh

W.

, eds,

Identification and Reduction of Satellite-Induced Signals in GRACE Accelerometer Data; Observation of the System Earth from Space - CHAMP, GRACE, GOCE and Future Missions. GEOTECHNOLOGIEN Science Report No. 20 , Vol. 2014

.

Springer

, p.

53

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Pires

S.

,

Bergé

J.

,

Baghi

Q.

,

Touboul

P.

,

Métris

G.

,

2016

,

Phys. Rev. D

,

94

,

123015

10.1103/PhysRevD.94.123015

Crossref

Search ADS

Powell

J.

,

Trifirò

D.

,

Cuoco

E.

,

Heng

I. S.

,

Cavaglià

M.

,

2015

,

Class. Quantum Gravity

,

32

,

215012

10.1088/0264-9381/32/21/215012

Crossref

Search ADS

Powell

J.

,

Torres-Forné

A.

,

Lynch

R.

,

Trifirò

D.

,

Cuoco

E.

,

Cavaglià

M.

,

Heng

I. S.

,

Font

J. A.

,

2017

,

Class. Quantum Gravity

,

34

,

034002

10.1088/1361-6382/34/3/034002

Crossref

Search ADS

Principe

M.

,

Pinto

I. M.

,

2017

,

Phys. Rev. D

,

95

,

082006

10.1103/PhysRevD.95.082006

Crossref

Search ADS

Refregier

A.

,

2003

,

MNRAS

,

338

,

35

10.1046/j.1365-8711.2003.05901.x

Crossref

Search ADS

Refregier

A.

,

Bacon

D.

,

2003

,

MNRAS

,

338

,

48

10.1046/j.1365-8711.2003.05902.x

Crossref

Search ADS

Rummel

R.

,

Yi

W.

,

Stummer

C.

,

2011

,

J. Geod.

,

85

,

777

10.1007/s00190-011-0500-0

Crossref

Search ADS

Sabzmeydani

P.

,

Mori

G.

,

2007

,

IEEE Conference on Computer Vision and Pattern Recognition

, Vol.

1

, p.

90

Scoville

N.

et al. .,

2007a

,

ApJS

,

172

,

1

10.1086/516585

Crossref

Search ADS

Scoville

N.

et al. .,

2007

b,

ApJS

,

172

,

38

10.1086/516580

Crossref

Search ADS

Sharpee

T. O.

,

Victor

J. D.

,

2009

,

J. Comput Neurosci

,

26

,

203

10.1007/s10827-008-0107-5

Crossref

Search ADS

PubMed

Suderman

R.

,

Lizotte

D. J.

,

Abukhdeir

N. M.

,

2015

,

Phys. Rev. E

,

91

,

033307

10.1103/PhysRevE.91.033307

Crossref

Search ADS

Tagore

A. S.

,

Jackson

N.

,

2016

,

MNRAS

,

457

,

3066

10.1093/mnras/stw057

Crossref

Search ADS

Tapley

B. D.

,

Bettadpur

S.

,

Watkins

M.

,

Reigber

C.

,

2004

,

Geophys. Res. Lett.

,

31

,

L09607

10.1029/2004GL019920

Crossref

Search ADS

Touboul

P.

et al. .,

2017

,

Phys. Rev. Lett.

,

119

,

231101

10.1103/PhysRevLett.119.231101

Crossref

Search ADS

PubMed

Weissman

J.

,

Hancewicz

T.

,

Kaplan

P.

,

2004

,

Opt. Express

,

12

,

5760

10.1364/OPEX.12.005760

Crossref

Search ADS

PubMed

Yang

X. L.

,

Guo

S. H.

,

Chan

F. T.

,

Wong

K. W.

,

Ching

W. Y.

,

1991

,

Phys. Rev. A

,

43

,

1186

10.1103/PhysRevA.43.1186

Crossref

Search ADS

PubMed

Zaslow

B.

,

Zandler

M. E.

,

1967

,

Am. J. Phys.

,

35

,

1118

10.1119/1.1973790

Crossref

Search ADS

APPENDIX A: DERIVATION OF CONVOLUTION KERNEL TENSOR C_nml FOR 1D EXPONENTIAL SHAPELETS

A1 Proof of equation (13)

Let f(x) and g(x) be two one-dimensional functions, whose convolution product is h(x) = f(x)⋆g(x). In Fourier space, |$\tilde{h}(k) = \tilde{f}(k) \tilde{g}(k)$|⁠, where

\begin{eqnarray*} \tilde{f}(k) = \frac{1}{\sqrt{2\pi }} \int f(x) \mathrm{ e}^{{\rm i}kx} {\rm d}x, \end{eqnarray*}

(A1)

and similarly for |$\tilde{g}(k)$|⁠.

Decomposing f(x) and g(x) into exponential shapelets (⁠|$f(x) = \sum _{m=1}^\infty f_m \Psi _m(x,\alpha)$|⁠), substituting in equation (A1) and rearranging terms, we get

\begin{eqnarray*} \tilde{h}(k) = \sum _{m=1}^\infty \sum _{l=1}^\infty f_m g_l \tilde{\Psi }_m(k,\alpha) \tilde{\Psi }_l(k,\beta), \end{eqnarray*}

(A2)

where we recall that

\begin{eqnarray*} \tilde{\Psi _l}(k,\beta) = (-1)^l \sqrt{\frac{2l\beta }{\pi }} \frac{(l\beta k -{\rm i})^{2l}}{\left[ (l\beta k)^2 +1 \right]^{l+1}}. \end{eqnarray*}

(A3)

Using the Parseval–Plancherel theorem, we get

\begin{eqnarray*} h_n = \int {\rm d}k \tilde{h}(k) \overline{\tilde{\Psi }_n(k,\gamma)}, \end{eqnarray*}

(A4)

where the bar denotes the complex conjugate.

Substituting equation (A2) in equation (A4), and rearranging terms, we find

\begin{eqnarray*} h_n = \sum _{m=1}^\infty \sum _{l=1}^\infty f_m g_l \int _{-\infty }^\infty {\rm d}k (-1)^{n+m+l} \sqrt{\alpha \beta \gamma } ~ \sqrt{mnl} ~ \frac{(m\alpha k - {\rm i})^{2m}}{[(m\alpha k)^2 + 1]^{m+1}} \frac{(l\beta k - {\rm i})^{2l}}{[l\beta k)^2 + 1]^{l+1}} \frac{(n\gamma k - {\rm i})^{2n}}{[(n\gamma k)^2 + 1]^{n+1}}, \end{eqnarray*}

(A5)

which defines C_nml such that

\begin{eqnarray*} h_n = \sum _{m,l=1}^\infty f_m g_l C_{nml}. \end{eqnarray*}

(A6)

The final expression for C_nml (equation 13) is then found using the binomial decompositions of (mαk − i)^2m, (lβk − i)^2l, and (nγk − i)²ⁿ.

It can be noted that C_nml is a complex number, unless u + v + w is even, such that i^{u + v + w} = −1. If u + v + w is odd, then the integral is 0 (see equation A2). Hence, we do not need to impose any restriction on u + v + w for C_nml to be real. We can also note that since u ≤ 2m, v ≤ 2l, and w ≤ 2n, the integral in equation (13) converges.

We will then aim to look for an analytic expression for the integral in equation (13). We introduce

\begin{eqnarray*} \mathcal {I}_{u,v,w}^{m,l,n} \equiv \int _{-\infty }^\infty {\rm d}k (-1)^w \frac{k^{u+v+w}}{\left[(m\alpha k)^2 + 1\right]^{m+1} \left[(l\beta k)^2 + 1\right]^{l+1} \left[(n\gamma k)^2 + 1\right]^{n+1}}. \end{eqnarray*}

(A7)

A2 Properties of |$\mathcal {I}_{u,v,w}^{m,l,n}$|’s integrand

In this section, we derive some properties of the integrand that appears in the definition of |$\mathcal {I}_{u,v,w}^{m, l, n}$|⁠. Let us name it f(x) (rigorously, we should write |$f_{u,v,w}^{m,l,n}(x)$|⁠, but we drop the u, v, w, n, m, l indices to simplify the notation), such that

\begin{eqnarray*} f(x) = (-1)^w \frac{x^{u+v+w}}{\left[(m\alpha x)^2 + 1\right]^{m+1} \left[(l\beta x)^2 + 1\right]^{l+1} \left[(n\gamma x)^2 + 1\right]^{n+1}}, \end{eqnarray*}

(A8)

where 0 ≤ u ≤ 2m, 0 ≤ v ≤ 2l, and 0 ≤ w ≤ 2n (see equation 13).

A2.1 Alternative definition

When developing binomial expressions in f, we get another form for the function, which appears as the inverse of a series of x

\begin{eqnarray*} 1/f(x) = (-1)^w \sum _{\nu _m=0}^{m+1} \sum _{\nu _l=0}^{l+1} \sum _{\nu _n=0}^{n+1} \binom{m+1}{\nu _m} \binom{l+1}{\nu _l} \binom{n+1}{\nu _n} (m\alpha)^{2\nu _m} (l\beta)^{2\nu _l} (n\gamma)^{2\nu _n} x^{2(\nu _m + \nu _l + \nu _n) -u-v-w} . \end{eqnarray*}

(A9)

A2.2 Parity

It is obvious that f(− x) = (− 1)^{u + v + w}f(x), so that f is even (resp. odd) when u + v + w is even (resp. odd). Hence, |$\mathcal {I}_{u,v,w}^{m,l,n}$| vanishes when u + v + w is odd, in which case C_nml = 0. As mentioned in the main text, C_nml is real if u + v + w is even (then, C_nml is real for all combinations of u, v, w, n, m, l). In the remainder of this appendix, we restrict ourselves to x ≥ 0.

A2.3 Value in 0

Evidently, f(0) = 0 if u + v + w ≠ 0. If u + v + w = 0, then f(0) = 1.

A2.4 Limit at x → ∞

Since 0 ≤ u ≤ 2m, 0 ≤ v ≤ 2l, and 0 ≤ w ≤ 2n, it is easy to see that lim_{x → ∞} = 0.

A2.5 Dependence on n, m, l

For a given x, f quickly decreases as m^{2(m + 1)} (and similarly for n and l), showing that only the contributions of low n, m, and l are significant in C_nml. Then, the series h_n (equation A6) converges quickly.

The left-hand panel of Fig. A1 shows f(x) for some realistic combinations (α, β, γ, n, ml, u, v, w). It can be seen that f(x) is well behaved, and tends quickly towards 0 for large x. Its extent depends on α, β, and γ. As mentioned above, for a given x, it quickly decreases when n, m, l increase. The right-hand panel of Fig. A1 shows f(x) when γ = α = β = 0.1, n = 2, m = 1 l = 1, w = 2, u = 1 and v = 1, together with the values x = min ((mα)⁻¹, (lβ)⁻¹, (nγ)⁻¹) and x = max ((mα)⁻¹, (lβ)⁻¹, (nγ)⁻¹) (borders of the grey area), which are key values in computing |$\mathcal {I}_{u,v,w}^{m,l,n}$| (see below).

f(x) function. Left: f(x) for different values of α, β, γ, n, m, l, u, v, w. Right: f(x) for γ = α = β) = 1, n = 2, m = 1 l = 1, w = 2, u = 1, v = 1; the grey region shows the domain [min ((mα)−1, (lβ)−1, (nγ)−1), max ((mα)−1, (lβ)−1, (nγ)−1)].

Figure A1.

f(x) function. Left: f(x) for different values of α, β, γ, n, m, l, u, v, w. Right: f(x) for γ = α = β) = 1, n = 2, m = 1 l = 1, w = 2, u = 1, v = 1; the grey region shows the domain [min ((mα)⁻¹, (lβ)⁻¹, (nγ)⁻¹), max ((mα)⁻¹, (lβ)⁻¹, (nγ)⁻¹)].

Open in new tab Download slide

A3 Computation of the integral in C_nml

We now turn to look for an analytic expression for |$I_{u,v,w}^{m,l,n}$|⁠. We first note that, due to the parity property of |$f_{u,v,w}^{m,l,n}(x)$|⁠,

\begin{eqnarray*} \mathcal {I}_{u,v,w}^{m,l,n} \equiv 2 \mathcal {I}_{u,v,w,\gt }^{m,l,n}= 2 \int _0^\infty f_{u,v,w}^{m,l,n}(k) {\rm d}k \end{eqnarray*}

(A10)

if u + v + w is even, or |$I_{u,v,w}^{m,l,n}=0$| otherwise.

We decompose |$I_{u,v,w,\gt }^{m,l,n}$| as

\begin{eqnarray*} \mathcal {I}_{u,v,w,\gt }^{m,l,n} = \int _0^\mu f_{u,v,w}^{m,l,n}(k) {\rm d}k + \int _\mu ^M f_{u,v,w}^{m,l,n}(k) {\rm d}k + \int _M^\infty f_{u,v,w}^{m,l,n}(k) {\rm d}k, \end{eqnarray*}

(A11)

where, as we show below, we impose μ < min ((mα)⁻¹, (lβ)⁻¹, (nγ)⁻¹) and M > max ((mα)⁻¹, (lβ)⁻¹, (nγ)⁻¹).

A3.1 Computation of first integral on r.h.s. of equation (A11)

We first introduce

\begin{eqnarray*} I_{u,v,w}^{m,l,n}(\mu) \equiv \int _0^\mu f_{u,v,w}^{m,l,n}(k) {\rm d}k, \end{eqnarray*}

(A12)

where, for clarity, we wrote f with all its indices.

Changing variable such that k = yμ, we get

\begin{eqnarray*} I_{u,v,w, 1}^{m,l,n}(\mu) = (-1)^w \mu ^{u+v+w+1} \int _0^1 {\rm d}y [(m\alpha \mu)^2y^2+1]^{-m-1} [(l\beta \mu)^2y^2+1]^{-l-1} [(n\gamma \mu)^2y^2+1]^{-n-1} y^{u+v+w}. \end{eqnarray*}

(A13)

Another change of variable, t = y², gives

\begin{eqnarray*} I_{u,v,w}^{m,l,n}(\mu) = \frac{(-1)^w}{2} \mu ^{u+v+w+1} \int _0^1 t^{\frac{u+v+w-1}{2}} [1+ (m\alpha \mu)^2 t]^{-m-1} [1+(l\beta \mu)^2 t]^{-l-1} [1+(n\gamma \mu)^2 t]^{-n-1} {\rm d}t, \end{eqnarray*}

(A14)

where we can recognize, if μ < min ((mα)⁻¹, (lβ)⁻¹, (nγ)⁻¹), Lauricella’s function of the fourth kind (e.g. Hasanov & Srivastava 2007; Bezrodnykh 2016)

\begin{eqnarray*} F_D^{(r)}(a, b1, \dots , b_r, c ; x_1, \dots , x_r) = \frac{\Gamma (c)}{\Gamma (a) \Gamma (c-a)} \int _0^1 t^{a-1} (1-t)^{c-a-1} \prod _{i=1}^r (1-x_it)^{-b_i} {\rm d}t \end{eqnarray*}

(A15)

with r = 3, |$a=\frac{u+v+w+1}{2}$|⁠, |$c = \frac{u+v+w+3}{2}$|⁠, b₁ = m + 1, b₂ = l + 1, c₂ = n + 1, x₁ = −(mαμ)², x₂ = −(lβμ)², and x₃ = −(nγμ)², such that

\begin{eqnarray*} I_{u,v,w}^{m,l,n}(\mu) &=& \frac{(-1)^w}{2} \mu ^{u+v+w+1} \frac{\Gamma \left((u+v+w+1)/2\right)}{\Gamma \left((u+v+w+3)/2\right)} \nonumber\\ && \times F_D^{(3)}\left(\frac{u+v+w+1}{2}, m+1, l+1, n+1, \frac{u+v+w+3}{2}; -(m\alpha L)^2, -(l\beta \mu)^2, -(n\gamma L)^2 \right). \end{eqnarray*}

(A16)

|$F_D^{(3)}$| is defined (as long as max (|x_i|) < 1) as a converging series (equation 1.4 in Hasanov & Srivastava 2007), such that

\begin{eqnarray*} F_D^{(r)}(a, b1, \dots , b_r, c ; x_1, \dots , x_r) = \sum _{\nu 1, \dots , \nu _r=0}^\infty \frac{(a)_{\nu _1 + \dots + \nu _r} (b_1)_{\nu _1} \dots (b_r)_{\nu _r}}{(c)_{\nu _1+\dots +\nu _r}} \frac{x_1^{\nu _1}}{\nu _1!}\dots \frac{x_r^{\nu _r}}{\nu _r!}, \end{eqnarray*}

(A17)

where (a)_ν is the Pochhammer symbol, such that (after simplifying Pochhammer symbols):

\begin{eqnarray*} I_{u,v,w}^{m,l,n}(\mu) &=& (-1)^w \sum _{\nu _1,\nu _2,\nu _3=0}^\infty \frac{(-1)^{\nu _1+\nu _2+\nu _3}}{u+v+w+2(\nu _1+\nu _2+\nu _3) + 1} \frac{(m+1)_{\nu _1} (l+1)_{\nu _2} (n+1)_{\nu _3}}{\nu _1!\nu _2!\nu _3!} \nonumber\\ && \times (m\alpha)^{2\nu _1} (l\beta)^{2\nu _2} (n\gamma)^{2\nu _3} \mu ^{u+v+w+2(\nu _1+\nu _2+\nu _3)}. \end{eqnarray*}

(A18)

A3.2 Computation of second integral on r.h.s. of equation (A11)

We then introduce

\begin{eqnarray*} J_{u,v,w}^{m,l,n}(\mu ,M) \equiv \int _\mu ^M f_{u,v,w}^{m,l,n}(k) {\rm d}k, \end{eqnarray*}

(A19)

which we will compute in a similar fashion to the previous integral.

With successive changes of variables y = k², z = y − μ², and t = z/(M² − μ²) (under our assumptions, M ≠ μ; if M = μ, |$J_{u,v,w}^{m,l,n}(\mu ,M)=0$|⁠), we obtain

\begin{eqnarray*} J_{u,v,w}^{m,l,n}(\mu ,M) &=& \frac{(-1)^w (M^2-\mu ^2) \mu ^{u+v+w-1}}{2 \left[ 1 + (\mu m\alpha)^2\right]^{m+1} \left[ 1 + (\mu l\beta)^2\right]^{l+1} \left[ 1 + (\mu n\gamma)^2\right]^{n+1}}\int _0^1 \left[ 1 - \frac{\mu ^2-M^2}{\mu ^2}t\right]^{\frac{u+v+w-1}{2}} \nonumber\\ &&\times \, \left[1 - \frac{(\mu ^2-M^2)(m\alpha)^2}{1+(\mu m\alpha)}t \right]^{-m-1} \left[1 - \frac{(\mu ^2-M^2)(l\beta)^2}{1+(\mu l\beta)}t \right]^{-l-1}\left[1 - \frac{(\mu ^2-M^2)(n\gamma)^2}{1+(\mu n\gamma)}t \right]^{-n-1} {\rm d}t \end{eqnarray*}

(A20)

If |$M \geqslant \sqrt{2}\mu$|⁠, then |$J_{u,v,w}^{m,l,n}(\mu ,M)$| must be calculated numerically. This is the case in the example of Fig. A1. However, if |$M\lt \sqrt{2}\mu$| (i.e. max {(mα)⁻¹, (lβ)⁻¹, (nγ)⁻¹} < |$\sqrt{2} \min \left\lbrace (m\alpha)^{-1}, (l\beta)^{-1}, (n\gamma)^{-1}\right\rbrace$|⁠, then we recognize a Lauricella |$F_D^{(4)}$| function, such that

\begin{eqnarray*} J_{u,v,w}^{m,l,n}(\mu ,M) &=& \frac{(-1)^w (M^2-\mu ^2) \mu ^{u+v+w-1}}{2 \Gamma (2) \left[ 1 + (\mu m\alpha)^2\right]^{m+1} \left[ 1 + (\mu l\beta)^2\right]^{l+1} \left[ 1 + (\mu n\gamma)^2\right]^{n+1}} \nonumber\\ &&\times \, F_D^{(4)} \left(1, \frac{1-u-v-w}{2}, m+1, l+1, n+1, 2; \frac{\mu ^2-M^2}{\mu _2}, \frac{(\mu ^2-M^2)(m\alpha)^2}{1+(\mu m\alpha)^2}, \frac{(\mu ^2-M^2)(l\beta)^2}{1+(\mu l\beta)^2} \frac{(\mu ^2-M^2)(n\gamma)^2}{1+(\mu n\gamma)^2}\right).\nonumber\\ \end{eqnarray*}

(A21)

If |$M\lt \sqrt{2}\mu$|⁠, |$J_{u,v,w}^{m,l,n}(\mu ,M)$| can finally be expressed as a converging series

\begin{eqnarray*} J_{u,v,w}^{m,l,n}(\mu ,M) &=& \frac{(-1)^{w+1}}{2} \sum _{\nu _1,\nu _2,\nu _3, \nu _4=0}^\infty \frac{\left(\frac{1-u-v-w}{2}\right)_{\nu _1} (m+1)_{\nu _2} (l+1)_{\nu _3} (n+1)_{\nu _4}}{(\nu _1+\nu _2+\nu _3+\nu _4) \nu _1! \nu _2! \nu _3! \nu _4!} \nonumber\\ &&\times \, \frac{\mu ^{-2\nu _1} (\mu ^2-M2)^{\nu _1+\nu _2+\nu _3+\nu _4} (m\alpha)^{2\nu _2} (l\beta)^{2\nu _3} (n\gamma)^{2\nu _4}}{\left[ 1+ (\mu m\alpha)^2\right]^{m+1+\nu _2} \left[ 1+ (\mu l\beta)^2\right]^{l+1+\nu _3} \left[ 1+ (\mu n\gamma)^2\right]^{n+1+\nu _4}}. \end{eqnarray*}

(A22)

A3.3 Computation of third integral on r.h.s of equation (A11)

We now introduce

\begin{eqnarray*} K_{u,v,w}^{m,l,n}(M) \equiv \int _M^\infty f_{u,v,w}^{m,l,n}(k) {\rm d}k, \end{eqnarray*}

(A23)

which we will compute in a similar fashion to |$I_{u,v,w}^{m,l,n}$| and |$J_{u,v,w}^{m,l,n}$|⁠.

We first note that |$\int _M^\infty f_{u,v,w}^{m,l,n}(k) {\rm d}k = \lim _{\xi \rightarrow \infty } I_\xi (M)$|⁠, where, for simplicity, we dropped the indices from the K_ξ definition, and ξ is some positive cut-off, and

\begin{eqnarray*} K_\xi \equiv \int _M^\xi f_{u,v,w}^{m,l,n}(k) {\rm d}k . \end{eqnarray*}

(A24)

K_ξ can be computed using successive changes of variable. First, we set y = 1/k, so that (after rearranging some terms)

\begin{eqnarray*} K_\xi = \int _{1/\xi }^{1/M} \frac{(-1)^w y^{-2-(u+v+w)-2(3+n+m+l)}}{(m\alpha)^{2(m+1)}(l\beta)^{2(l+1)}(n\gamma)^{2(n+1)}} \left[ 1 + \frac{y^2}{(m\alpha)^2}\right]^{-m-1} \left[ 1 + \frac{y^2}{(l\beta)^2}\right]^{-l-1} \left[ 1 + \frac{y^2}{(n\gamma)^2}\right]^{-n-1} {\rm d}y. \end{eqnarray*}

(A25)

Additional changes of variables (z = y², t′ = z − 1/ξ², |$t= \frac{M^2\xi ^2}{\xi ^2-M^2}t$|⁠) then yield

\begin{eqnarray*} K_\xi &=& \frac{(-1^w)}{2(m\alpha)^{2(m+1)}(l\beta)^{2(l+1)}(n\gamma)^{2(n+1)}} \frac{\xi ^2-M^2}{\xi ^2M^2} \left(\frac{\xi ^2 (m\alpha)^2+1}{\xi ^2(m\alpha)^2}\right)^{-m-1} \left(\frac{\xi ^2 (l\beta)^2+1}{\xi ^2(l\beta)^2}\right)^{-l-1} \left(\frac{\xi ^2 (n\gamma)^2+1}{\xi ^2(n\gamma)^2}\right)^{-n-1} \nonumber\\ &&\times \, \int _0^1 \left(\frac{1}{\xi ^2} + \frac{\xi ^2-M^2}{\xi ^2M^2}t\right)^g \left(1 + \frac{\xi ^2-M^2}{M^2[(\xi m\alpha)^2-1]}t\right)^{-m-1} \left(1 + \frac{\xi ^2-M^2}{M^2[(\xi l\beta)^2-1]}t\right)^{-l-1} \left(1 + \frac{\xi ^2-M^2}{M^2[(\xi n\gamma)^2-1]}t\right)^{-n-1} {\rm d}t \end{eqnarray*}

(A26)

where g = 3 + n + m + l − (3 + u + v + w)/2. Taking the limit ξ → ∞, we obtain

\begin{eqnarray*} K_{u,v,w}^{m,l,n}(M) &=& \frac{(-1)^w M^{-2(3+n+m+l)+1+u+v+w}}{2(m\alpha)^{2(m+1)}(l\beta)^{2(l+1)}(n\gamma)^{2(n+1)}} \nonumber\\ &&\times \, \int _0^1 t^{3+n+m+l-(3+u+v+w)/2} \left[1+\frac{t}{(Mm\alpha)^2}\right]^{-m-1} \left[1+\frac{t}{(Ml\beta)^2}\right]^{-l-1} \left[1+\frac{t}{(Mn\gamma)^2}\right]^{-n-1} {\rm d}t. \end{eqnarray*}

(A27)

Under the assumption that M > max ((mα)⁻¹, (lβ)⁻¹, (nγ)⁻¹), we recognize a Lauricella |$F_D^{(3)}$| function, such that

\begin{eqnarray*} K_{u,v,w}^{m,l,n}(M) &=& \frac{(-1)^w M^{-2(3+n+m+l)+1+u+v+w}}{2(m\alpha)^{2(m+1)}(l\beta)^{2(l+1)}(n\gamma)^{2(n+1)}} \frac{\Gamma \Big(n+m+l-\frac{u+v+w+5}{2}\Big)}{\Gamma \Big(n+m+l-\frac{u+v+w+7}{2}\Big)} \nonumber\\ &&\times \, F_D^{(3)} \bigg(n+m+l-\frac{u+v+w+5}{2}, m+1, l+1, n+1, n+m+l-\frac{u+v+w+7}{2};\nonumber\\ &&\quad \quad\quad \quad -\,(Mm\alpha)^{-2}, -(Ml\beta)^{-2}, -(Mn\gamma)^{-2}\bigg). \end{eqnarray*}

(A28)

and we can finally express |$K_{u,v,w}^{m,l,n}(M)$| as a converging series

\begin{eqnarray*} K_{u,v,w}^{m,l,n}(M) &=& (-1)^w \sum _{\nu _1,\nu _2,\nu _3=0}^{\infty } \frac{(-1)^{\nu _1+\nu _2+\nu _3}}{5+2(n+m+l)-(u+v+w)} \frac{(m+1)_{\nu _1} (l+1)_{\nu _2} (n+1)_{\nu _3}}{\nu _1!\nu _2!\nu _3!} \nonumber\\ &&\times \, (m\alpha)^{-2(m+1-\nu _1)} (l\beta)^{-2(l+1-\nu _2)} (n\gamma)^{-2(n+1-\nu _3)} M^{-2(3+n+m+l)+u+v+w+1-2(\nu _1+\nu _2+\nu _3)}. \end{eqnarray*}

(A29)

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Download all slides

Month:	Total Views:
March 2019	5
April 2019	9
May 2019	6
June 2019	23
July 2019	9
August 2019	7
September 2019	14
October 2019	7
November 2019	4
December 2019	6
January 2020	6
February 2020	2
March 2020	13
April 2020	14
May 2020	6
June 2020	12
July 2020	7
August 2020	15
September 2020	13
October 2020	5
November 2020	7
December 2020	28
January 2021	11
February 2021	15
March 2021	23
April 2021	19
May 2021	11
June 2021	10
July 2021	14
August 2021	13
September 2021	6
October 2021	11
November 2021	25
December 2021	14
January 2022	17
February 2022	21
March 2022	22
April 2022	21
May 2022	25
June 2022	15
July 2022	35
August 2022	16
September 2022	27
October 2022	35
November 2022	21
December 2022	21
January 2023	16
February 2023	31
March 2023	35
April 2023	38
May 2023	22
June 2023	46
July 2023	20
August 2023	25
September 2023	16
October 2023	22
November 2023	13
December 2023	19
January 2024	21
February 2024	14
March 2024	27
April 2024	24
May 2024	19
June 2024	11
July 2024	25
August 2024	22
September 2024	17
October 2024	23
November 2024	14
December 2024	15
January 2025	8
February 2025	19
March 2025	16
April 2025	10
May 2025	3

Article Contents

Exponential shapelets: basis functions for data analysis of isolated features

Abstract

1 INTRODUCTION

2 1D EXPONENTIAL SHAPELETS

2.1 Definition

2.2 Properties

2.2.1 Maximum and minimum effective scales

2.2.2 Fourier and Laplace transforms

2.2.3 Convolution

2.2.4 Rescaling

2.3 Shape characterization

2.4 Exponential shapelet modelling in practice

2.5 Example applications

2.5.1 Cleaning accelerometer time series data, and modelling an experiment’s response function

2.5.2 Characterization of perturbing signals in space-borne geodesy missions

2.5.3 Unmodelled bursts and glitches in gravitational wave data analysis

2.6 Comparison with Gaussian shapelets

3 2D EXPONENTIAL SHAPELETS

3.1 Definition

3.2 Properties

3.2.1 Maximum and minimum effective scales

3.2.2 Fourier transformation and convolution

3.2.3 Coordinate transforms

3.3 Object shape measurement

3.3.1 Radial profile

3.3.2 Flux

3.3.3 Centroid

3.3.4 Quadrupole moments

3.3.5 Size and ellipticity

3.3.6 Convergence

3.4 Example application

3.4.1 Distant galaxies in astronomical imaging

3.4.2 Convergence speed

4 CONCLUSION

ACKNOWLEDGEMENTS

Footnotes

REFERENCES

APPENDIX A: DERIVATION OF CONVOLUTION KERNEL TENSOR Cnml FOR 1D EXPONENTIAL SHAPELETS

A1 Proof of equation (13)

A2 Properties of |$\mathcal {I}_{u,v,w}^{m,l,n}$|’s integrand

A2.1 Alternative definition

A2.2 Parity

A2.3 Value in 0

A2.4 Limit at x → ∞

A2.5 Dependence on n, m, l

A3 Computation of the integral in Cnml

A3.1 Computation of first integral on r.h.s. of equation (A11)

A3.2 Computation of second integral on r.h.s. of equation (A11)

A3.3 Computation of third integral on r.h.s of equation (A11)

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

APPENDIX A: DERIVATION OF CONVOLUTION KERNEL TENSOR C_nml FOR 1D EXPONENTIAL SHAPELETS

A3 Computation of the integral in C_nml