© ESO, 2012

1. Introduction

Challenges in modern cosmology

The wealth of cosmological data in the last few decades (Larson et al. 2011; Schrabback et al. 2010; Percival et al. 2007b) has led to the establishment of a standard model of cosmology, which describes the Universe as composed today of approximately 4% baryons, 22% dark matter and 74% dark energy. The main challenges in modern cosmology are to understand the nature of both dark energy and dark matter, as well as the initial conditions of the Universe (Albrecht et al. 2006; Peacock et al. 2006). A thorough understanding of these three topics may lead to a revision of Einstein’s theory of General Relativity and our view of the early Universe.

New surveys are planned who aim to answer these important questions e.g. Planck for the CMB (Planck Collaboration 2011), DES (Dark Energy Survey, Annis et al. 2005), BOSS (Baryon Oscillation Spectroscopic Survey, Schlegel et al. 2007), LSST (Large Synoptic Survey Telescope, Tyson & LSST 2004) and Euclid (Laureijs et al. 2011) for weak lensing and the study of large scale structure with galaxy surveys.

The challenge with these upcoming large data-sets is to extract the cosmological information in the most suitable manner in order to test the cosmological paradigm. Depending on the signal one wishes to extract, and/or survey geometry, different bases may be more or less suitable (e.g., Fourier, spherical harmonic, configuration or wavelet space). Moreover, future surveys may be in 2D (e.g. Planck) or in 3D (e.g. galaxy or weak lensing surveys). Where 3D data is available, a tomographic analysis is possible (also known as 2D 1/2), or a full 3D analysis can be done. For data in spherical coordinates, this corresponds to a spherical Fourier-Bessel (SFB) decomposition (Heavens & Taylor 1995; Fisher et al. 1995; Castro et al. 2005; Erdoğdu et al. 2006a,b; Leistedt et al. 2012; Rassat & Refregier 2012).

Wavelet transform on the sphere

Wavelets are particularly well suited to the analysis of cosmological data (Martinez et al. 1993; Starck et al. 2006), since cosmological data can often be sparsely represented in wavelet space. 2D Wavelets have been used in many astrophysical studies (Starck & Murtagh 2006) for a broad range of applications such as denoising, deconvolution, detection, etc. CMB studies have motivated the development of 2D spherical wavelet decompositions. Continuous wavelet transforms on the sphere (Antoine 1999; Tenorio et al. 1999; Cayón et al. 2001; Holschneider 1996) have been proposed, mainly for non Gaussianity studies. In Starck et al. (2006), an invertible isotropic undecimated wavelet transform (UWT) on the sphere based on spherical harmonics was described, that can be also used for other applications such as deconvolution, component separation (Moudden et al. 2005; Bobin et al. 2008; Delabrouille et al. 2009), inpainting (Abrial et al. 2007, 2008), or Poisson denoising (Schmitt et al. 2010). A similar wavelet construction has been published in (Marinucci et al. 2008; Faÿ & Guilloux 2008; Faÿ et al. 2008) using so-called “needlet filters”, and in Wiaux et al. (2008), an algorithm was proposed which allows us to reconstruct an image from its steerable wavelet transform. Other multiscale transforms on the sphere such as ridgelets and curvelets have been developed (Starck et al. 2006), and are well adapted to detect anisotropic features. Other multiscale transforms on the sphere, such as ridgelets and curvelets, have been developed (Starck et al. 2006) and are well adapted to detect anisotropic features. An extension of this UWT has also been developed for polarised CMB data in Starck et al. (2009).

In this paper, we describe a new 3D isotropic spherical wavelet decomposition, which is reversible, and could therefore be useful for many different applications. It is based on the UWT proposed by Starck et al. (2006) and extended into 3D. The 3D-UWT proposed here can be used to analyse 3D data in spherical coordinates, such as a 3D galaxy or weak lensing survey with large (but not necessarily full) sky coverage.

2. Spherical 3D filtering using the spherical Fourier-Bessel transform

2.1. Spherical Fourier-Bessel transform

The SFB transform of a square integrable scalar field f(r,θ,φ) can be defined as: $\mathit{f̂}\mathit{lm}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\$\mathit{f}\mathrm{\left(}\mathit{r,\theta ,\phi }\mathrm{\right)}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}\overline{{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}}{\mathit{r}}^{\mathrm{2}}\sin \mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{d}\mathit{r}\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi ,}$(1)where ${\mathit{Y}}_{\mathit{l}}^{\mathit{m}}$ are spherical harmonics and j_l are spherical Bessel functions and $\overline{\mathit{Y}}$ represents the complex conjugate of Y. Note Eq. (1) uses the same convention as Heavens & Taylor (1995) which differs slightly from that of Castro et al. (2005), Leistedt et al. (2012) and Rassat & Refregier (2012), as explained below. This expression allows the expansion of a 3D field provided in spherical coordinates onto a set of orthogonal functions: ${\mathrm{\Psi }}_{\mathit{lmk}}\mathrm{\left(}\mathit{r,\theta ,\phi }\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathit{.}$(2)From the orthogonality of the Ψ_lmk’s, the original field f(r,θ,φ) can be reconstructed from its SFB transform using the following inversion formula: $\mathit{f}\mathrm{\left(}\mathit{r,\theta ,\phi }\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\sum _{\mathit{l}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathit{l}}\mathrm{\int }\mathit{f̂}\mathit{lm}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathit{k}}^{\mathrm{2}}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}\mathrm{d}\mathit{k}{\mathit{Y}}_{\mathit{lm}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathit{.}$(3)From this definition, the SFB transform can also be regarded as the commutative composition of two different transforms: a spherical harmonics transform (SHT) for the angular dimension and a spherical Bessel transform (SBT) for the radial dimension. In this work, the following convention is adopted to define the SHT of a function f(θ,φ) defined on the sphere: $\begin{array}{ccc}{\mathit{f}}_{\mathit{lm}}& \mathrm{=}& \\ \mathit{f}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}& \mathrm{=}& \end{array}$Our choice for Ψ_lmk allows us to give a symmetrical expression for the SBT and its inverse: $\begin{array}{ccc}\mathit{f̂}\mathit{l}\mathrm{\left(}\mathit{k}\mathrm{\right)}& & \mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\mathrm{\int }\mathit{f}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r,}\\ \mathit{f}\mathrm{\left(}\mathit{r}\mathrm{\right)}& & \mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\mathrm{\int }\mathit{f̂}\mathit{l}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}{\mathit{k}}^{\mathrm{2}}\mathrm{d}\mathit{k}\mathit{.}\end{array}$Although this symmetrical formulation for the SBT may differ from other works (e.g., Castro et al. 2005; Leistedt et al. 2012; Rassat & Refregier 2012), it will prove very convenient, especially to obtain a discretised transform (see Sect. 4).

2.2. Frequency filtering using the SFB transform

Spherical 3D filtering can be defined as the 3D convolution product of a 3D field in spherical coordinates with a 3D filter also provided in spherical coordinates. Using the relations presented in Sects. A.1 and A.2, it is possible to express such a product in terms of SFB coefficients and to relate those coefficients to regular Fourier coefficients.

To illustrate spherical 3D filtering on a simple case, let us consider a simple isotropic low-pass filter u whose 3D Fourier transform is U(k,θ_k,φ_k). The 3D Fourier transform of such a filter has a spherical symmetry and takes a simple form in spherical coordinates. Indeed, U(k,θ_k,φ_k) is a function only of k because of its symmetry, therefore U(k,θ_k,φ_k) = U(k). Furthermore, since U(k,θ_k,φ_k) is constant on every sphere centred on the origin in Fourier space, its SHT for a given k verifies ∀(l,m) ≠ (0,0), U_lm(k) = 0. As a result, the SFB transform of u has the following simple expression: $\mathit{û}\mathit{lm}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{ifl}\mathrm{=}\mathrm{m}\mathrm{=}\mathrm{0}\mathrm{,}\\ & \mathrm{otherwise}\mathrm{.}\end{array}$(6)Let us now consider a 3D field f defined in spherical coordinates. The filtered field v resulting from applying u to f is obtained using Eq. (A.8) simplified by the properties of g: $\begin{array}{ccc}\mathit{v̂}\mathit{lm}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& {\begin{array}{}[\\ \mathrm{(}\mathit{f}\mathrm{\ast }\mathit{u}\mathrm{)}\end{array}}_{\mathit{lm}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}\\ & \mathrm{=}& \mathrm{\left(}\mathit{i}{\mathrm{)}}^{\mathit{l}}\sqrt{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\sum _{{\mathit{l}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{{\mathit{m}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}{\mathit{l}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\right(}\mathrm{-}\mathit{i}{\mathrm{)}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathit{f̂}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\\ & & \mathrm{\times } \sum _{{\mathit{l}}^{\mathrm{\prime \prime }}\mathrm{=}\mathrm{|}\mathit{l}\mathrm{-}{\mathit{l}}^{\mathrm{\prime }}\mathrm{|}}^{\mathit{l}\mathrm{+}{\mathit{l}}^{\mathrm{\prime }}}{\mathit{c}}^{{\mathit{l}}^{\mathrm{\prime \prime }}}\mathrm{\left(}\mathit{l,m,}{\mathit{l}}^{\mathrm{\prime }}\mathit{,}{\mathit{m}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{(}\mathrm{-}\mathit{i}{\mathrm{)}}^{{\mathit{l}}^{\mathrm{\prime \prime }}}\mathit{û}{\mathit{l}}^{\mathrm{\prime \prime }}\mathit{m}\mathrm{-}{\mathit{m}}^{\mathrm{\prime }}\mathrm{(}\mathit{k}\mathrm{)}{\mathit{\delta }}_{{\mathit{l}}^{\mathrm{\prime \prime }}\mathrm{0}}{\mathit{\delta }}_{{\mathit{m}}^{\mathrm{\prime \prime }}\mathrm{0}}\mathit{,}\\ & \mathrm{=}& \mathrm{(}\mathrm{-}{\mathit{i}}^{\mathrm{2}}{\mathrm{)}}^{\mathit{l}}\sqrt{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\mathit{f̂}\mathit{lm}\mathrm{(}\mathit{k}\mathrm{\left)}{\mathit{c}}^{\mathrm{0}}\mathrm{\right(}\mathit{l,m,l,m}\mathrm{\left)}\mathrm{\right(}\mathrm{-}\mathit{i}{\mathrm{)}}^{\mathrm{0}}\mathit{û}\mathrm{00}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}\end{array}$(7)Knowing that ${\mathit{c}}^{\mathrm{0}}\mathrm{\left(}\mathit{l,m,l,m}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi }}$, the following expression is finally obtained: $\mathit{v̂}\mathit{lm}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{2}}\mathit{\pi }\mathit{û}\mathrm{00}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{f̂}\mathit{lm}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(8)In the special case of a 3D isotropic filter, frequency filtering is therefore easily obtained using the SFB transform and the filtered coefficients are simply the original coefficients multiplied by a function of k.

Although this filter seems overly simplistic, it will be shown in the following sections that such a low-pass filter is at the heart of the isotropic undecimated spherical 3D wavelet transform which makes direct use of Eq. (8).

3. Isotropic undecimated spherical 3D wavelet transform

An isotropic undecimated spherical wavelet transform defined on the sphere and based on the SHT was proposed in Starck et al. (2006). Now the aim is to transpose the ideas behind this transform to the case of data in 3D spherical coordinates. Indeed, the isotropic wavelet transform can be defined using only an isotropic low-pass filter. In the last section the necessary relations to apply such a filter have been obtained and the Isotropic Undecimated Spherical 3D Wavelet transform can now be defined.

3.1. Wavelet decomposition

Using the formalism introduced in the previous section, a Wavelet transform can be defined with the SFB transform. This isotropic transform is based on a scaling function φ^k_c(r,θ_r,φ_r) with cut-off frequency k_c and spherical symmetry. The symmetry of this function is preserved in the Fourier space and therefore, its SFB transform verifies $\mathit{\phi ̂}\begin{array}{}{\mathit{k}}_{\mathrm{c}}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathrm{0}$ whenever (l,m) ≠ (0,0). Furthermore, due to its cut-off frequency, the scaling function verifies $\mathit{\phi ̂}\begin{array}{}{\mathit{k}}_{\mathrm{c}}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathrm{0}$ for all k ≥ k_c. In other terms, the scaling function verifies: ${\mathrm{\Phi }}^{{\mathit{k}}_{\mathrm{c}}}\mathrm{\left(}\mathit{r,}{\mathit{\theta }}_{\mathit{r}}\mathit{,}{\mathit{\phi }}_{\mathit{r}}\mathrm{\right)}\mathrm{=}{\mathrm{\Phi }}^{{\mathit{k}}_{\mathrm{c}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{k}}_{\mathrm{c}}}\mathrm{\Phi ̂}\begin{array}{}{\mathit{k}}_{\mathrm{c}}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathit{k}}^{\mathrm{2}}{\mathit{j}}_{\mathrm{0}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}\mathrm{d}\mathit{k}{\mathit{Y}}_{\mathrm{0}}^{\mathrm{0}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{r}}\mathit{,}{\mathit{\phi }}_{\mathit{r}}\mathrm{\right)}\mathit{.}$(9)Using relation (8) the convolution of the original data f(r,θ,φ) with Φ^k_c becomes very simple: $\mathit{ĉ}\begin{array}{}\mathrm{0}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\begin{array}{}[\\ \left[{\mathrm{\Phi }}_{{\mathit{k}}_{\mathrm{c}}}\mathrm{\ast }\mathit{f}\right]\end{array}}_{\mathit{lm}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{2}}\mathit{\pi }\mathrm{\Phi ̂}\begin{array}{}{\mathit{k}}_{\mathrm{c}}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{f̂}\mathit{lm}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(10)Using this scaling function, it is possible to define a sequence of smoother approximations c^j(r,θ_r,φ_r) of a function f(r,θ_r,φ_r) on a dyadic resolution scale. Let Φ^2 − jk_c be a rescaled version of Φ^k_c with cut-off frequency 2 ^− jk_c. Then c^j(r,θ_r,φ_r) is obtained by convolving f(r,θ_r,φ_r) with Φ^2 − jk_c: $\begin{array}{ccc}{\mathit{c}}^{\mathrm{0}}& \mathrm{=}& {\mathrm{\Phi }}^{{\mathit{k}}_{\mathrm{c}}}\mathrm{\ast }\mathit{f,}\\ {\mathit{c}}^{\mathrm{1}}& \mathrm{=}& {\mathrm{\Phi }}^{{\mathrm{2}}^{-1}{\mathit{k}}_{\mathrm{c}}}\mathrm{\ast }\mathit{f,}\\ & \mathrm{·}\mathrm{·}\mathrm{·}& \\ {\mathit{c}}^{\mathit{j}}& \mathrm{=}& {\mathrm{\Phi }}^{{\mathrm{2}}^{\mathrm{-}\mathit{j}}{\mathit{k}}_{\mathrm{c}}}\mathrm{\ast }\mathit{f}\mathit{.}\end{array}$(11)Applying the SFB transform to the last relation yields: $\mathit{ĉ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{2}}\mathit{\pi }\mathrm{\Phi ̂}\begin{array}{}{\mathrm{2}}^{\mathrm{-}\mathit{j}}{\mathit{k}}_{\mathrm{c}}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{f̂}\mathit{lm}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(12)This leads to the following recurrence formula: $\mathrm{\forall }\mathit{k}\mathit{<}\frac{{\mathit{k}}_{\mathrm{c}}}{{\mathrm{2}}^{\mathit{j}}}\mathit{,} \mathit{ĉ}\begin{array}{}\mathit{j}\mathrm{+}\mathrm{1}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Phi ̂}\begin{array}{}{\mathrm{2}}^{\mathrm{-}\mathrm{(}\mathit{j}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathit{k}}_{\mathrm{c}}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}}{\mathrm{\Phi ̂}\begin{array}{}{\mathrm{2}}^{\mathrm{-}\mathit{j}}{\mathit{k}}_{\mathrm{c}}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathit{ĉ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(13)Just like with the à trous algorithm by Starck et al. (2010), the wavelet coefficients  { w^j }  can now be defined as the difference between two consecutive resolutions: ${\mathit{w}}^{\mathit{j}\mathrm{+}\mathrm{1}}\mathrm{\left(}\mathit{r,}{\mathit{\theta }}_{\mathit{r}}\mathit{,}{\mathit{\phi }}_{\mathit{r}}\mathrm{\right)}\mathrm{=}{\mathit{c}}^{\mathit{j}}\mathrm{\left(}\mathit{r,\theta ,\phi }\mathrm{\right)}\mathrm{-}{\mathit{c}}^{\mathit{j}\mathrm{+}\mathrm{1}}\mathrm{\left(}\mathit{r,\theta ,\phi }\mathrm{\right)}\mathit{.}$(14)This choice for the wavelet coefficients is equivalent to the following definition for the wavelet function Ψ^k_c: $\mathrm{\Psi ̂}\begin{array}{}{\mathrm{2}}^{\mathrm{-}\mathit{j}}{\mathit{k}}_{\mathrm{c}}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathrm{\Phi ̂}\begin{array}{}{\mathrm{2}}^{\mathrm{-}\mathrm{(}\mathit{j}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{k}}_{\mathrm{c}}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{-}\mathrm{\Phi ̂}\begin{array}{}{\mathrm{2}}^{\mathrm{-}\mathit{j}}{\mathit{k}}_{\mathrm{c}}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(15)so that: $\begin{array}{ccc}{\mathit{w}}^{\mathrm{0}}& \mathrm{=}& {\mathrm{\Psi }}^{{\mathit{k}}_{\mathrm{c}}}\mathrm{\ast }\mathit{f,}\\ {\mathit{w}}^{\mathrm{1}}& \mathrm{=}& {\mathrm{\Psi }}^{{\mathrm{2}}^{-1}{\mathit{k}}_{\mathrm{c}}}\mathrm{\ast }\mathit{f,}\\ & \mathrm{·}\mathrm{·}\mathrm{·}& \\ {\mathit{w}}^{\mathit{j}}& \mathrm{=}& {\mathrm{\Psi }}^{{\mathrm{2}}^{\mathrm{-}\mathit{j}}{\mathit{k}}_{\mathrm{c}}}\mathrm{\ast }\mathit{f}\mathit{.}\end{array}$(16)By applying the SFB transform to the definition of the wavelet coefficients and using the recurrence formula verified by the c^js yields: $\mathrm{\forall }\mathit{k}\mathit{<}\frac{{\mathit{k}}_{\mathrm{c}}}{{\mathrm{2}}^{\mathit{j}}}\mathit{,} \mathit{ŵ}\begin{array}{}\mathit{j}\mathrm{+}\mathrm{1}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{\Phi ̂}\begin{array}{}{\mathrm{2}}^{\mathrm{-}\mathrm{(}\mathit{j}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathit{k}}_{\mathrm{c}}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}}{\mathrm{\Phi ̂}\begin{array}{}{\mathrm{2}}^{\mathrm{-}\mathit{j}}{\mathit{k}}_{\mathrm{c}}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}}\right)\mathit{ĉ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(17)Equations (17) et (13) which define the wavelet decomposition are in fact equivalent to convolving the resolution at a given scale j with a low-pass and a high-pass filter in order to obtain respectively the resolution and the wavelet coefficients at scale j + 1.

The low-pass filter, denoted here by h^j, can be defined for each scale j by: $\mathit{ĥ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \frac{\mathrm{\Phi ̂}\begin{array}{}{\mathrm{2}}^{\mathrm{-}\mathrm{(}\mathit{j}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathit{k}}_{\mathrm{c}}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}}{\mathrm{\Phi ̂}\begin{array}{}{\mathrm{2}}^{\mathrm{-}\mathit{j}}{\mathit{k}}_{\mathrm{c}}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}}& \mathrm{if}\mathit{k}\mathit{<}\frac{{\mathit{k}}_{\mathrm{c}}}{{\mathrm{2}}^{\mathit{j}\mathrm{+}\mathrm{1}}} \mathrm{and} \mathit{l}\mathrm{=}\mathit{m}\mathrm{=}\mathrm{0}\mathit{,}\\ \mathrm{0}& \mathrm{otherwise}\mathrm{.}\end{array}$(18)Then the approximation at scale j + 1 is given by the convolution of scale j with h^j: ${\mathit{c}}^{\mathit{j}\mathrm{+}\mathrm{1}}\mathrm{=}{\mathit{c}}^{\mathit{j}}\mathrm{\ast }\frac{\mathrm{1}}{\sqrt{\mathrm{2}}\mathit{\pi }}{\mathit{h}}^{\mathit{j}}\mathit{.}$(19)In the same way, a high pass filter g^j can be defined on each scale j by: $\begin{array}{ccc}\mathit{ĝ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ & \\ & \end{array}& & \end{array}$(20)Given the definition of Ψ, g^j can also be expressed in the simple form: $\mathit{ĝ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\mathit{ĥ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(21)The wavelet coefficients at scale j are obtained by convolving the resolution at scale j − 1 with g^j − 1: ${\mathit{w}}^{\mathit{j}}\mathrm{=}{\mathit{c}}^{\mathit{j}\mathrm{-}\mathrm{1}}\mathrm{\ast }\frac{\mathrm{1}}{\sqrt{\mathrm{2}}\mathit{\pi }}{\mathit{g}}^{\mathit{j}\mathrm{-}\mathrm{1}}\mathit{.}$(22)In summary, the two relations necessary to recursively define the wavelet transform are: $\begin{array}{ccc}\mathit{ĉ}\begin{array}{}\mathit{j}\mathrm{+}\mathrm{1}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \mathit{ĥ}\begin{array}{}\mathit{j}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{ĉ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}\\ \mathit{ŵ}\begin{array}{}\mathit{j}\mathrm{+}\mathrm{1}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \mathit{ĝ}\begin{array}{}\mathit{j}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{ĉ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}\end{array}$

3.2. Reconstruction

Since the wavelet coefficients are defined as the difference between two resolutions, the reconstruction from the wavelet decomposition is straightforward and corresponds to the reconstruction formula of the à trous algorithm: ${\mathit{c}}^{\mathrm{0}}\mathrm{\left(}\mathit{r,}{\mathit{\theta }}_{\mathit{r}}\mathit{,}{\mathit{\phi }}_{\mathit{r}}\mathrm{\right)}\mathrm{=}{\mathit{c}}^{\mathit{J}}\mathrm{\left(}\mathit{r,}{\mathit{\theta }}_{\mathit{r}}\mathit{,}{\mathit{\phi }}_{\mathit{r}}\mathrm{\right)}\mathrm{+}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{J}}{\mathit{w}}^{\mathit{j}}\mathrm{\left(}\mathit{r,}{\mathit{\theta }}_{\mathit{r}}\mathit{,}{\mathit{\phi }}_{\mathit{r}}\mathrm{\right)}\mathit{.}$(25)However, given the redundancy of the transform, the reconstruction is not unique. It is possible to take advantage of this redundancy to reconstruct c^j from c^j + 1 and w^j + 1 by using a least squares estimate.

Fig. 1 Scaling function and wavelet function for kc = 1. The y-axis represents the amplitude and the x-axis the frequency associated with the scaling function. The units of the amplitude must be those of the SFB transform.

From the previous paragraph, the wavelet decomposition can be then recursively defined by: $\begin{array}{ccc}\mathit{ĉ}\begin{array}{}\mathit{j}\mathrm{+}\mathrm{1}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \mathit{ĥ}\begin{array}{}\mathit{j}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{ĉ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}\\ \mathit{ŵ}\begin{array}{}\mathit{j}\mathrm{+}\mathrm{1}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \mathit{ĝ}\begin{array}{}\mathit{j}\\ \mathrm{00}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{ĉ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}\end{array}$If these relations are respectively multiplied by $\overline{\mathit{ĥ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ and $\overline{\mathit{ĝ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ and then added together, the following expression is obtained for the least squares estimate of c^j from c^j + 1 and w^j + 1: $\mathit{ĉ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathit{ĉ}\begin{array}{}\mathit{j}\mathrm{+}\mathrm{1}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\widehat{\begin{array}{}\mathrm{˜}\\ \mathit{h}\end{array}}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}\mathit{ŵ}\begin{array}{}\mathit{j}\mathrm{+}\mathrm{1}\\ \mathit{lm}\end{array}\widehat{\begin{array}{}\mathrm{˜}\\ \mathit{g}\end{array}}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(28)where $\widehat{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{h}\end{array}}{}^{\mathit{j}}$ and $\widehat{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{g}\end{array}}{}^{\mathit{j}}$ are defined as follows: $\begin{array}{ccc}\widehat{\begin{array}{}\mathrm{˜}\\ \mathit{h}\end{array}}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \frac{{\overline{\mathit{ĥ}\mathit{j}}}_{\mathit{lm}}\mathrm{\left(}\mathit{k}\mathrm{\right)}}{\mathrm{\left|}\mathit{ĥ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\right(}\mathit{k}\mathrm{)}{\mathrm{|}}^{\mathrm{2}}\mathrm{+}\mathrm{|}\mathit{ĝ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathrm{|}}^{\mathrm{2}}}\mathit{,}\\ \widehat{\begin{array}{}\mathrm{˜}\\ \mathit{g}\end{array}}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \frac{{\overline{\mathit{ĝ}\mathit{j}}}_{\mathit{lm}}\mathrm{\left(}\mathit{k}\mathrm{\right)}}{\mathrm{\left|}\mathit{ĥ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\right(}\mathit{k}\mathrm{)}{\mathrm{|}}^{\mathrm{2}}\mathrm{+}\mathrm{|}\mathit{ĝ}\begin{array}{}\mathit{j}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathrm{|}}^{\mathrm{2}}}\mathrm{·}\end{array}$Among the advantages of using this reconstruction formula instead of the raw sum over the wavelet coefficients is that there is no need to perform an inverse and then direct SFB transform to reconstruct the coefficients of the original data. Indeed, both the wavelet decomposition and reconstruction procedures only require access to SFB coefficients and there is no need to revert back to the direct space (although this will be required if one wishes to observe the power of each coefficient in real space as we do in Sect. 5 for denoising purposes).

3.3. Choice of a scaling function

Any function with spherical symmetry and a cut-off frequency k_c would do as a scaling function but in this work we choose to use a B-spline function of order 3 to define our scaling function: $\mathrm{\Phi ̂}\begin{array}{}{\mathit{k}}_{\mathrm{c}}\\ \mathit{lm}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{B}}_{\mathrm{3}}\left(\frac{\mathrm{2}\mathit{k}}{{\mathit{k}}_{\mathrm{c}}}\right){\mathit{\delta }}_{\mathit{l}\mathrm{0}}{\mathit{\delta }}_{\mathit{m}\mathrm{0}}\mathit{,}$(31)where ${\mathit{B}}_{\mathrm{3}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{12}}\mathrm{\left(}\mathrm{\right|}\mathit{x}\mathrm{-}\mathrm{2}{\mathrm{|}}^{\mathrm{3}}\mathrm{-}\mathrm{4}\mathrm{|}\mathit{x}\mathrm{-}\mathrm{1}{\mathrm{|}}^{\mathrm{3}}\mathrm{+}\mathrm{6}\mathrm{|}\mathit{x}{\mathrm{|}}^{\mathrm{3}}\mathrm{-}\mathrm{4}\mathrm{|}\mathit{x}\mathrm{+}\mathrm{1}{\mathrm{|}}^{\mathrm{3}}\mathrm{+}\mathrm{|}\mathit{x}\mathrm{+}\mathrm{2}{\mathrm{|}}^{\mathrm{3}}\mathrm{)}\mathit{.}$(32)The scaling function and its corresponding wavelet function are plotted in SFB space for different values of j in Fig. 1.

Fig. 2 Comparaison between spline and needlet wavelet functions on the sphere. The y-axis corresponds to amplitude and the x-axis to position.

Other functions such as Meyer wavelets or the needlet function (Marinucci et al. 2008) can be used as well. Needlet wavelet functions have a much better frequency localisation than the wavelet function derived from the B₃-spline but present more oscillations in the direct space. To illustrate this, we show in Fig. 2 two different wavelet functions. Figure 2 left shows 1D profile of the spline (continuous line) and the needlet wavelet function (dotted line) at a given scale. Figure 2 right shows the same function, but we have plotted the absolute value in order to better visualise their respective ringing. As it can be seen, for wavelet functions with the same main lob, the needlet wavelet oscillates much more than the spline wavelet. Hence, the best wavelet choice certainly depends on the desired applications. For statistical analysis, detection or restoration applications, we may prefer to use a wavelet which does not oscillate too much and with a smaller support, in this case the spline wavelet is clearly the correct choice. For spectral or bispectral analysis, where the frequency localisation is fundamental, then needlets should be preferred to the spline wavelet.

4. The discrete spherical Fourier-Bessel transform (DSFBT)

The spherical 3D wavelet formalism derived in the previous section is based on a continuous SFB transform. However, a continuous transform is not practical numerically and a discrete equivalent is required. In the case of galaxy surveys, where the underlying density field is sampled at discrete points, this problem can be addressed by combining a boundary condition with the discrete nature of the survey (see for e.g. Eqs. (5) and (7) in Leistedt et al. 2012, even though they use a slightly different SFB expansion as in this paper) which yields in this specific case a practical way to estimate discrete SFB coefficients. However this method cannot be used for the purpose of performing wavelet treatments since this requires a continuous field.

Indeed, to perform operations on the wavelet scales (e.g. for denoising purposes), one first needs to recover the wavelet coefficients in physical space, apply a treatment to those coefficients (e.g. thresholding) and then convert them back into SFB space (as required by the wavelet transform algorithm which operates in SFB space). This last operation requires a SFB transform which takes as input fields and not discrete galaxy surveys.

In this section, we show how a discretisation of the k spectrum in combination with a Healpix type angular pixelisation, can be used to build a fast Discrete SFB Transform (DSFBT) which allows for a novel and practical two-way conversion between the discrete SFB coefficients and discrete values of the field in physical space sampled on a spherical 3D grid.

4.1. The discrete spherical Bessel transform

The discretisation of the SFB transform along the radial component uses the well known orthogonality property of the spherical Bessel functions on the interval  [0,R] . If f is a continuous function defined on  [0,R]  which verifies the boundary condition f(R) = 0 then the SBT defined by Eq. (5) can be expressed using SFB series: $\begin{array}{ccc}\mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{\mathit{ln}}\mathrm{\right)}& \mathrm{=}& \\ \mathit{f}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \end{array}$In this expression, ${\mathit{k}}_{\mathit{ln}}\mathrm{=}\frac{{\mathit{q}}_{\mathit{ln}}}{\mathit{R}}$ where q_ln is the nth zero of the Bessel function of order l and the weights ρ_ln are defined as: ${\mathit{\rho }}_{\mathit{ln}}\mathrm{=}\frac{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{R}}^{-3}}{{\mathit{j}}_{\mathit{l}\mathrm{+}\mathrm{1}}^{\mathrm{2}}\mathrm{\left(}{\mathit{q}}_{\mathit{ln}}\mathrm{\right)}}\mathrm{·}$(35)Although this formulation provides a discretisation of the inverse SBT and of the k spectrum, the direct transform is still continuous and another discretisation step is required.

If a boundary condition of the same kind is applied to so that , then by using the same result, the SFB expansion of is obtained by: $\begin{array}{ccc}\mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{r}}_{\mathit{ln}}\mathrm{\right)}& \mathrm{=}& \sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{K}}\mathit{f̂}\mathit{l}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}{\mathit{r}}_{\mathit{ln}}\mathit{k}\mathrm{\right)}{\mathit{k}}^{\mathrm{2}}\mathrm{d}\mathit{k,}\\ \mathit{f̂}\mathit{l}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{r}}_{\mathit{ln}}\mathrm{\right)}{\mathit{\kappa }}_{\mathit{ln}}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}{\mathit{r}}_{\mathit{ln}}\mathit{k}\mathrm{\right)}\mathit{,}\end{array}$where ${\mathit{r}}_{\mathit{ln}}\mathrm{=}\frac{{\mathit{q}}_{\mathit{ln}}}{{\mathit{K}}_{\mathit{l}}}$ and where the weights κ_ln are defined as: ${\mathit{\kappa }}_{\mathit{ln}}\mathrm{=}\frac{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{K}}_{\mathit{l}}^{-3}}{{\mathit{j}}_{\mathit{l}\mathrm{+}\mathrm{1}}^{\mathrm{2}}\mathrm{\left(}{\mathit{q}}_{\mathit{ln}}\mathrm{\right)}}\mathrm{·}$(37)The SBT being an involution, so that . Much like the previous set of equations had introduced a discrete k_ln grid, a discrete r_ln grid is obtained for the radial component. Since Eqs. (34) and (36b) can be used to compute f and for any value of r and k, they can in particular be used to compute f(r_ln) and where l′ does not have to match l. The SBT and its inverse can then be expressed only in terms of series: $\begin{array}{ccc}\mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{n}}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{p}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\mathit{f}\mathrm{\left(}{\mathit{r}}_{\mathit{lp}}\mathrm{\right)}{\mathit{\kappa }}_{\mathit{lp}}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}{\mathit{r}}_{\mathit{lp}}{\mathit{k}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{n}}\mathrm{\right)}\mathit{,}\\ \mathit{f}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{n}}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{p}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{\mathit{lp}}\mathrm{\right)}{\mathit{\rho }}_{\mathit{lp}}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{n}}{\mathit{k}}_{\mathit{lp}}\mathrm{\right)}\mathit{.}\end{array}$With these equations one can easily compute the SBT and its inverse without the need of evaluating any integral. Furthermore only discrete values of f and respectively sampled on r_ln and k_ln are required.

In practical applications, for a given value of l only a limited number N_max of and f(r_ln) coefficients can be stored so that r_lNmax = R and k_lNmax = K_l. Since r_ln is defined by ${\mathit{r}}_{\mathit{ln}}\mathrm{=}\frac{{\mathit{q}}_{\mathit{ln}}}{{\mathit{K}}_{\mathit{l}}}$, for n = N_maxR and K_l are bound by the following relationship: ${\mathit{q}}_{\mathit{lN}}\mathrm{=}{\mathit{K}}_{\mathit{l}}\mathit{R}\mathit{.}$(39)Therefore, from the value of R one can easily determine the appropriate value for K_l. Furthermore, because of the boundary conditions on and f, the series (38a) and (38b) can now be truncated to N_max terms.

This truncation allows the use of a very convenient matrix formalism to represent the transform. In their explicit form, the two truncated sums can be rewritten as: $\begin{array}{ccc}\mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{n}}\mathrm{\right)}& \mathrm{=}& {\mathit{K}}_{\mathit{l}}^{-3}\sum _{\mathit{p}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\max }}\mathit{f}\mathrm{\left(}{\mathit{r}}_{\mathit{lp}}\mathrm{\right)}\frac{\sqrt{\mathrm{2}\mathit{\pi }}}{{\mathit{j}}_{\mathit{l}\mathrm{+}\mathrm{1}}^{\mathrm{2}}\mathrm{\left(}{\mathit{q}}_{\mathit{lp}}\mathrm{\right)}}{\mathit{j}}_{\mathit{l}}\left(\frac{{\mathit{q}}_{\mathit{lp}}{\mathit{q}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{n}}}{{\mathit{q}}_{\mathit{l}{\mathit{N}}_{\max }}}\right)\mathit{,}\\ \mathit{f}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{n}}\mathrm{\right)}& \mathrm{=}& {\mathit{R}}^{-3}\sum _{\mathit{p}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\max }}\mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{\mathit{lp}}\mathrm{\right)}\frac{\sqrt{\mathrm{2}\mathit{\pi }}}{{\mathit{j}}_{\mathit{l}\mathrm{+}\mathrm{1}}^{\mathrm{2}}\mathrm{\left(}{\mathit{q}}_{\mathit{lp}}\mathrm{\right)}}{\mathit{j}}_{\mathit{l}}\left(\frac{{\mathit{q}}_{\mathit{lp}}{\mathit{q}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{n}}}{{\mathit{q}}_{\mathit{l}{\mathit{N}}_{\max }}}\right)\mathrm{·}\end{array}$From these equations, a transform matrix T^ll′ can be defined as: ${\mathit{T}}_{\mathit{pq}}^{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}}\mathrm{=}{\left(\frac{\sqrt{\mathrm{2}\mathit{\pi }}}{{\mathit{j}}_{\mathit{l}\mathrm{+}\mathrm{1}}^{\mathrm{2}}\mathrm{\left(}{\mathit{q}}_{\mathit{lq}}\mathrm{\right)}}{\mathit{j}}_{\mathit{l}}\left(\frac{{\mathit{q}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{p}}{\mathit{q}}_{\mathit{lq}}}{{\mathit{q}}_{\mathit{l}{\mathit{N}}_{\max }}}\right)\right)}_{\mathit{pq}}\mathrm{·}$(42)This matrix allows the computation of on any grid k_l′n from the values of f sampled on r_ln: $\left[\begin{array}{c}\\ \mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}^{\mathrm{\prime }}\mathrm{1}}\mathrm{\right)}\\ \mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}^{\mathrm{\prime }}\mathrm{2}}\mathrm{\right)}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}^{\mathrm{\prime }}{\mathit{N}}_{\max }}\mathrm{\right)}\end{array}\right]\mathrm{=}\frac{\mathrm{1}}{{\mathit{K}}_{\mathit{l}}^{\mathrm{3}}}{\mathit{T}}^{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}}\left[\begin{array}{c}\\ \mathit{f}\mathrm{\left(}{\mathit{r}}_{\mathit{l}\mathrm{1}}\mathrm{\right)}\\ \mathit{f}\mathrm{\left(}{\mathit{r}}_{\mathit{l}\mathrm{2}}\mathrm{\right)}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \mathit{f}\mathrm{\left(}{\mathit{r}}_{\mathit{l}{\mathit{N}}_{\max }}\mathrm{\right)}\end{array}\right]\mathit{.}$(43)Reciprocally, the inverse the values of f can be computed on any r_l′n grid from sampled on k_l′n using the exact same matrix: $\left[\begin{array}{c}\\ \mathit{f}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}^{\mathrm{\prime }}\mathrm{1}}\mathrm{\right)}\\ \mathit{f}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}^{\mathrm{\prime }}\mathrm{2}}\mathrm{\right)}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \mathit{f}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{N}}\mathrm{\right)}\end{array}\right]\mathrm{=}\frac{\mathrm{1}}{{\mathit{R}}^{\mathrm{3}}}{\mathit{T}}^{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}}\left[\begin{array}{c}\\ \mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{\mathit{l}\mathrm{1}}\mathrm{\right)}\\ \mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{\mathit{l}\mathrm{2}}\mathrm{\right)}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{\mathit{lN}}\mathrm{\right)}\end{array}\right]\mathit{.}$(44)The discrete SBT is defined by the set of Eqs. (43) and (44).

A simplified form of the transform could have been defined only for l = l′ so that T^ll′ = T^l. However, keeping the distinction between the order of the transform l and the order of the grid on which the results are provided l′ will be crucial to the implementation of the DSFBT. Indeed, the order of the grid on which the function is sampled has to match the order of the transform but the resulting transform coefficients do not. Therefore, it will be possible to compute the result of the inverse SBT of any order l on a grid of order l₀ so that only one radial grid of order l₀ will be required. Nevertheless, for the direct transform, if the field is sampled on the radial grid of order l₀, only the transform of order l₀ can be computed. An additional result is required to be able to relate the SBT of different orders between them. This is achieved by combining relations (33) and (34):

$\begin{array}{ccc}\mathit{f̂}\mathit{l}\mathrm{\left(}{\mathit{k}}_{\mathit{ln}}\mathrm{\right)}& \mathrm{=}& \sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{R}}\left[\sum _{\mathit{m}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\mathit{f̂}{\mathit{l}}_{\mathrm{0}}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}_{\mathrm{0}}\mathit{m}}\mathrm{\right)}{\mathit{\rho }}_{{\mathit{l}}_{\mathrm{0}}\mathit{m}}{\mathit{j}}_{{\mathit{l}}_{\mathrm{0}}}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}_{\mathrm{0}}\mathit{m}}\mathit{r}\mathrm{\right)}\right]{\mathit{j}}_{\mathit{l}}\mathrm{\left(}{\mathit{k}}_{\mathit{ln}}\mathit{r}\mathrm{\right)}{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r,}\\ & \mathrm{=}& \sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\sum _{\mathit{m}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\mathit{f̂}{\mathit{l}}_{\mathrm{0}}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}_{\mathrm{0}}\mathit{m}}\mathrm{\right)}{\mathit{\rho }}_{{\mathit{l}}_{\mathrm{0}}\mathit{m}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{R}}{\mathit{j}}_{{\mathit{l}}_{\mathrm{0}}}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}_{\mathrm{0}}\mathit{m}}\mathit{r}\mathrm{\right)}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}{\mathit{k}}_{\mathit{ln}}\mathit{r}\mathrm{\right)}{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r,}\\ & \mathrm{=}& \sum _{\mathit{m}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\mathit{f̂}{\mathit{l}}_{\mathrm{0}}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}_{\mathrm{0}}\mathit{m}}\mathrm{\right)}\frac{\mathrm{2}}{{\mathit{j}}_{{\mathit{l}}_{\mathrm{0}}\mathrm{+}\mathrm{1}}^{\mathrm{2}}\mathrm{\left(}{\mathit{q}}_{{\mathit{l}}_{\mathrm{0}}\mathit{m}}\mathrm{\right)}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}{\mathit{j}}_{{\mathit{l}}_{\mathrm{0}}}\mathrm{\left(}{\mathit{q}}_{{\mathit{l}}_{\mathrm{0}}\mathit{m}}\mathit{x}\mathrm{\right)}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}{\mathit{q}}_{\mathit{ln}}\mathit{x}\mathrm{\right)}{\mathit{x}}^{\mathrm{2}}\mathrm{d}\mathit{x,}\\ & \mathrm{=}& \end{array}$(45)where the weights ${\mathit{W}}_{\mathit{nm}}^{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}}$ are defined as: ${\mathit{W}}_{\mathit{nm}}^{{\mathit{l}}_{\mathrm{0}}\mathit{l}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}{\mathit{j}}_{{\mathit{l}}_{\mathrm{0}}}\mathrm{\left(}{\mathit{q}}_{{\mathit{l}}_{\mathrm{0}}\mathit{m}}\mathit{x}\mathrm{\right)}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}{\mathit{q}}_{\mathit{ln}}\mathit{x}\mathrm{\right)}{\mathit{x}}^{\mathrm{2}}\mathrm{d}\mathit{x}\mathit{.}$(46)The final expression in Eq. (45) is an important result which shows that the SBT of a given order can be expressed as the sum of the coefficients obtained for a different order of the transform, with the appropriate weighting. This means we can convert the Spherical Bessel coefficients of order l₀ into coefficients of any other order l, which considerably speeds up calculations for the SBT. It is also worth noticing that the weights ${\mathit{W}}_{\mathit{nm}}^{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}}$ are simply geometric terms, i.e. independent of the field and can thus be tabulated.

We note that this approach is an extension of the discrete SBT presented in Lemoine (1994) using spherical Bessel functions and where the transform in Lemoine (1994) can be considered as a special case where l = l′.

4.2. A spherical 3D grid compatible with a discrete spherical Fourier-Bessel spectrum

As presented in Sect. 2.1, the SFB transform is the composition of a SHT for the angular component and a SBT for the radial component. Since these two transforms can commute, they can be treated independently and by combining discrete algorithms for both transforms, one can build a DSFBT. In this work, the angular part of the transform is implemented using the HEALPix pixelisation scheme while the radial component uses the discrete SBT algorithm presented in Sect. 4.1. The choice of these two algorithms introduces a discretisation of the SFB coefficients as well as a 3D gridding of the density in physical space.

The SFB coefficients are defined by Eq. (1) for continuous values of k. Assuming a boundary condition on the density field f, the discrete SBT can be used to discretise the values of k. The Discrete SFB coefficients are therefore defined as: ${\mathit{f}}_{\mathit{lmn}}\mathrm{\equiv }\mathit{f̂}\mathit{lm}\mathrm{\left(}{\mathit{k}}_{\mathit{ln}}\mathrm{\right)}\mathit{,}$(47)for 0 ≤ l ≤ L_max,  − l ≤ m ≤ l and 1 ≤ n ≤ N_max. These discrete coefficients are simply obtained by sampling the continuous coefficients on the k_ln grid introduced in the previous section.

To this discretised SFB spectrum corresponds a dual grid of the 3D space defined by combining the HEALPix pixelisation scheme and the discrete SBT.

Used for the computation of the SHT, the HEALPix scheme maps the sphere using curvilinear quadrilaterals of different shapes but equal area. Although other pixelisation schemes for the sphere such as IGLOO or GLESP could have been used, the choice of HEALPix is essentially motivated by the homogeneity of the pixelisation and by its comprehensive software package. For a given value of r, the field f(r,θ,φ) can therefore be sampled on a finite number of points using HEALPix.

The radial component of the transform is conveniently performed using the discrete SBT. Indeed, this algorithm introduces a radial grid compatible with the discretised k_ln spectrum. Although this radial grid r_ln depends on the order l of the SBT, it will be justified in the next section that only one grid r_l0n is required to sample the field along the radial dimension. The value of l₀ is set to 0 because in this case the properties of the zeros of the Bessel function ensure that r_0n will be regularly spaced between 0 and R: ${\mathit{r}}_{\mathrm{0}\mathit{n}}\mathrm{=}\frac{\mathit{n}}{{\mathit{N}}_{\max }}\mathit{R}\mathit{.}$(48)For given values of θ and φ, the field f(r,θ,φ) can now be sampled on discrete values of r = r_0n.

Combining angular and radial grids, the 3D spherical grid is defined as a set of N_max HEALPix maps equally spaced between 0 and R. An illustration of this grid is provided in Fig. 3 where only one quarter of the space is represented for clarity.

Fig. 3 Representation of the spherical 3D grid for the DSFBT (R = 1 and Nmax = 4) as described in Sect. 4 where the axes represent x,y,z and the grid points correspond to the centre of Healpix pixels and the radial grid is defined in Sect. 4.1.

With this 3D grid, it will be possible to compute back and forth the SFB transform between a density field and its SFB coefficients. The details of the actual algorithm are provided in the next section.

4.3. Algorithm for the discrete spherical Fourier-Bessel transform (DSFBT)

Although the SHT and the SBT formally commute, some practical considerations on the order between the two operations are to be taken into account for their actual implementation. Here, a detailed algorithm of both the direct and inverse DSFBT is provided.

4.3.1. Direct transform

Given a density field f sampled on the spherical 3D grid, the SFB coefficients f_lmn are computed in three steps:

• 1)

  For each n between 1 and N_max the SHT of the HEALPix map of radius r_l0n is computed. This yields f_lm(r_l0n) coefficients.

• 2)

  The next step is to compute the SBT of order l₀ from the f_lm(r_l0n) coefficients for every (l,m). This operation is a simple matrix product: $\{\begin{array}{c}\\ \mathrm{\forall }& \mathrm{0}& \mathrm{\le }& \mathit{l}& \mathrm{\le }& {\mathit{L}}_{\max }\\ \mathrm{\forall }& \mathrm{-}\mathit{l}& \mathrm{\le }& \mathit{m}& \mathrm{\le }& \mathit{l,}\end{array}\left[\begin{array}{c}\\ \mathit{f̂}\begin{array}{}{\mathit{l}}_{\mathrm{0}}\\ \mathit{lm}\end{array}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}_{\mathrm{0}}\mathrm{1}}\mathrm{\right)}\\ \mathit{f̂}\begin{array}{}{\mathit{l}}_{\mathrm{0}}\\ \mathit{lm}\end{array}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}_{\mathrm{0}}\mathrm{2}}\mathrm{\right)}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \mathit{f̂}\begin{array}{}{\mathit{l}}_{\mathrm{0}}\\ \mathit{lm}\end{array}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}_{\mathrm{0}}{\mathit{N}}_{\max }}\mathrm{\right)}\end{array}\right]\mathrm{=}\frac{{\mathit{T}}^{{\mathit{l}}_{\mathrm{0}}{\mathit{l}}_{\mathrm{0}}}}{{\mathit{K}}^{\mathrm{3}}}\left[\begin{array}{c}\\ {\mathit{f}}_{\mathit{lm}}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}_{\mathrm{0}}\mathrm{1}}\mathrm{\right)}\\ {\mathit{f}}_{\mathit{lm}}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}_{\mathrm{0}}\mathrm{2}}\mathrm{\right)}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{f}}_{\mathit{lm}}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}_{\mathrm{0}}{\mathit{N}}_{\max }}\mathrm{\right)}\end{array}\right]\mathit{.}$(49)This operation yields $\mathit{f̂}\begin{array}{}{\mathit{l}}_{\mathrm{0}}\\ \mathit{lm}\end{array}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}_{\mathrm{0}}\mathit{n}}\mathrm{\right)}$ coefficients which are not yet SFB coefficients because the order of the SBT l₀ does not match the order of the spherical harmonics coefficients l. An additional step is required.

• 3)

  The last step required to gain access to the SFB coefficients f_lmn is to convert the Spherical Bessel coefficients for order l₀ to the correct order l that matches the spherical harmonics order. This is done by using relation (45): $\begin{array}{ccc}\{\begin{array}{c}\\ & & & & & \\ & & & & & \\ & & & & & \end{array}\mathit{f̂}\mathit{lm}\mathrm{\left(}{\mathit{k}}_{\mathit{ln}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{p}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\max }}\mathit{f̂}\begin{array}{}{\mathit{l}}_{\mathrm{0}}\\ \mathit{lm}\end{array}\mathrm{\left(}{\mathit{k}}_{{\mathit{l}}_{\mathrm{0}}\mathit{p}}\mathrm{\right)}\frac{\mathrm{2}{\mathit{W}}_{\mathit{np}}^{{\mathit{l}}_{\mathrm{0}}\mathit{l}}}{{\mathit{j}}_{{\mathit{l}}_{\mathrm{0}}\mathrm{+}\mathrm{1}}^{\mathrm{2}}\mathrm{\left(}{\mathit{q}}_{\mathit{lp}}\mathrm{\right)}}\mathit{,}& & \end{array}$(50)where ${\mathit{W}}_{\mathit{np}}^{{\mathit{l}}_{\mathrm{0}}\mathit{l}}$ are defined by Eq. (46). This operation finally yields the coefficients.

4.3.2. Inverse transform

Let f_lmn be the coefficients of the SFB transform of the density field f, where the f_lmn coefficients can be calculated for a continuous density field using the direct transform described above or for a galaxy or halo catalogue using the public code 3DEX (Leistedt et al. 2012). Note that the 3DEX code can account for masked regions of missing data. The reconstruction of f on the spherical 3D grid requires two steps:

• 1)

  First, from thef_lmn, the inverse discrete SBT is computed for all l and m. Again, this transform can be easily evaluated using a matrix product: $\{\begin{array}{c}\\ \mathrm{\forall }& \mathrm{0}& \mathrm{\le }& \mathit{l}& \mathrm{\le }& {\mathit{L}}_{\max }\\ \mathrm{\forall }& \mathrm{-}\mathit{l}& \mathrm{\le }& \mathit{m}& \mathrm{\le }& \mathit{l,}\end{array} \left[\begin{array}{c}\\ {\mathit{f}}_{\mathit{lm}}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}_{\mathrm{0}}\mathrm{1}}\mathrm{\right)}\\ {\mathit{f}}_{\mathit{lm}}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}_{\mathrm{0}}\mathrm{2}}\mathrm{\right)}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{f}}_{\mathit{lm}}\mathrm{\left(}{\mathit{r}}_{{\mathit{l}}_{\mathrm{0}}{\mathit{N}}_{\max }}\mathrm{\right)}\end{array}\right]\mathrm{=}\frac{{\mathit{T}}^{\mathit{l}{\mathit{l}}_{\mathrm{0}}}}{{\mathit{R}}^{\mathrm{3}}}\left[\begin{array}{c}\\ {\mathit{f}}_{\mathit{lm}\mathrm{1}}\\ {\mathit{f}}_{\mathit{lm}\mathrm{2}}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{f}}_{\mathit{lm}{\mathit{N}}_{\max }}\end{array}\right]\mathit{.}$(51)Here, it is worth noticing that the matrix T^ll₀ allows the evaluation of the SBT of order l and provides the results on the grid of order l₀.

• 2)

  From the spherical harmonics coefficients f_lm(r_l0n) given at specific radial distances r_l0n it is possible to compute the inverse SHT. For each n between 1 and N_max the HEALPix inverse SHT is performed on the set of coefficients  { f_lm(r_l0n) } _l,m. This yields N_max HEALPix maps which constitute the sampling of the reconstructed density field on the 3D spherical grid.

5. Wavelet decomposition of a test density field

To illustrate the wavelet transform described in Sect. 3, a set of SFB Coefficients was extracted from a 3D density field using the DSFBT described in in Sect. 4. The test density field was provided by a cosmological N-body simulation which was carried out by the Virgo Supercomputing Consortium using computers based at Computing Centre of the Max-Planck Society in Garching and at the Edinburgh Parallel Computing Centre. We use their large box simulations¹.

We first compute the SFB coefficients of the test density field by sampling the Virgo density field on the 3D grid illustrated in Fig. 3, for nside = 2048, l_max = 1023 and n_max = 512. In order to perform the SFB decomposition, we place ourselves in the middle of the box, and calculate the SFB coefficients out to r = 479/2   h^-1   Mpc, setting the over-density field to zero outside of this spherical volume. In order to better visualise the data, we present in Fig. 4a the data in a cube spanning half the size of the original box (i.e. 239.5   h^-1   Mpc), where the field has been reconstructed from the SFB coefficients, using the discrete inverse transform we discuss in Sect. 4.

From the SFB coefficients we calculate the wavelet decomposition, which yields the SFB coefficients of the various wavelet scales and the smoothed density field. Using the inverse DSFBT, the actual wavelet coefficients can be retrieved in the form of 3D density fields. These density fields corresponding to different wavelet scales are shown in Figs. 4b to e. The smoothed density field which arises from the wavelet decomposition is also shown in Fig. 4f, which is simply given by: $\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{-}\mathrm{\left(}\mathit{b}\mathrm{\right)}\mathrm{-}\mathrm{\left(}\mathit{c}\mathrm{\right)}\mathrm{-}\mathrm{\left(}\mathit{d}\mathrm{\right)}\mathrm{-}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathit{.}$(52)

Fig. 4 Spherical 3D wavelet decomposition of a density field for a box with side 239.5   h-1   Mpc. Panel a) shows the reconstructed density field directly using the initial SFB coefficients. Panels b–e) show the first four wavelet scales, while panel f) shows the smoothed density field.

6. Application to wavelet hard thresholding (denoising)

In this section, a noise removal application based on wavelet thresholding is presented as an example of a potential use for the Isotropic Undecimated Spherical 3D Wavelet transform.

Fig. 5 Wavelet Hard thresholding applied to a test density field.

6.1. Thresholding

Many wavelet filtering methods have been proposed in the last twenty years. Hard thresholding consists of setting to 0 all wavelet coefficients which have an absolute value lower than a threshold T_j (non-significant wavelet coefficient): $\begin{array}{ccc}\begin{array}{}\mathrm{˜}\\ {\mathit{w}}_{\mathit{j,k}}\end{array}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$where w_j,k is a wavelet coefficient at scale j and at spatial position k.

Soft thresholding consists of replacing each wavelet coefficient by the value $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{w}\end{array}$ where $\begin{array}{ccc}\begin{array}{}\mathrm{˜}\\ {\mathit{w}}_{\mathit{j,k}}\end{array}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$This operation is generally written as: $\begin{array}{ccc}\begin{array}{}\mathrm{˜}\\ {\mathit{w}}_{\mathit{j,k}}\end{array}\mathrm{=}\mathrm{soft}\mathrm{\left(}{\mathit{w}}_{\mathit{j,k}}\mathrm{\right)}\mathrm{=}\mathit{sgn}\mathrm{\left(}{\mathit{w}}_{\mathit{j,k}}\mathrm{\right)}\mathrm{\left(}\mathrm{\right|}{\mathit{w}}_{\mathit{j,k}}\mathrm{|}\mathrm{-}{\mathit{T}}_{\mathit{j}}{\mathrm{)}}_{\mathrm{+}}\mathit{,}& & \end{array}$(53)where (x) ₊  = max(0,x).

In the case of Gaussian noise with a given standard deviation σ, a reasonable choice for the threshold T_j is T_j = Kσ_j, where j is the scale of the wavelet coefficient, σ_j is the noise standard deviation at the scale j, and K is a constant generally chosen between to 3 and 5. The standard deviation σ_j depends only on the chosen wavelet function and the noise level σ. More details can be found in Starck et al. (2010).

Other threshold methods have been proposed, like the universal threshold (Donoho & Johnstone 1994; Donoho 1993), or the SURE (Stein unbiased risk estimate) method (Coifman & Donoho 1995), but they generally do not yield as good results as the hard thresholding method based on the significant coefficients. Concerning the threshold level, the universal threshold corresponds to a minimum risk. The larger the number of pixels, the larger the risk is, and it is normal that the threshold T depends on the number of pixels ($\mathit{T}\mathrm{=}{\mathit{\sigma }}_{\mathit{j}}\sqrt{\mathrm{2}\log \mathit{n}}$, n being the number of pixels). The Kσ threshold corresponds to a false detection probability, the probability to detect a coefficient as significant when it is due to the noise. The 3σ value corresponds to 0.27% false detection.

6.2. Denoising experiment

We test the hard thresholding algorithm by using it to remove artificially added Gaussian noise on the Virgo density field (the simulation is described in Sect. 5). To conduct this experiment, the original density field taken from the n-body simulation has been used to compute an initial set of SFB coefficients (with l_max = 1023 and N_max = 512) and out to r = 479   h^-1   Mpc. Figure 5a shows a slice of the 3D density field reconstructed from these coefficients. A Gaussian noise was then added to the SFB coefficients, the reconstruction of this noisy density field is shown in Fig. 5b. The hard thresholding algorithm was subsequently applied to the noisy SFB coefficients using 5 wavelet scales. The resulting density field is reconstructed in Fig. 5c. The residuals are shown in Fig. 5d. The wavelet analysis means we can successfully remove the noise we artificially added on entry, without much loss to the large scale structure, though some of the smaller structures are removed.

7. Conclusion

Modern cosmology requires the analysis of 3D fields on large areas of the sky, i.e. where the field is best viewed in spherical coordinates. In this configuration, a spherical Fourier Bessel (SFB) transform is the most natural way to statistically analyse the field. Wavelet transforms have been shown to be ideally suited for cosmological fields, which tend to be sparse in wavelet space. The wavelet transform can be used e.g. for denoising, but there is yet no 3D wavelet transform in spherical coordinates.

We present in this paper a new 3D spherical wavelet transform, based on the undecimated wavelet transform (UWT) described in (Starck et al. 2006). In order to perform operations on the wavelet transforms (such as denoising), we require a discrete version of the SFB transform for both the direct and inverse transforms. We show a novel way to perform such a fast discrete spherical Fourier-Bessel transform (DSFBT) based on both a discrete Bessel transform and the HEALPIX angular pixelisation scheme.

Using the 3D wavelet transform and the DSFBT, both introduced in this paper, we denoise a test large scale structure data set, taken from the Virgo large box simulations². We find we can satisfactorily remove artificially added Gaussian noise without much loss to the large scale structure. All the algorithms presented in this paper are available for download as a public code called MRS3D at http://jstarck.free.fr/mrs3d.html.

Acknowledgments

The authors are grateful to Boris Leistedt, Pirin Erdoğdu, Ofer Lahav and Alexandre Réfrégier for useful discussions regarding the Spherical Fourier-Bessel decomposition and discretisation scheme. The 3D wavelet library uses Healpix software (Górski et al. 2002, 2005). This research is in part supported by the Swiss National Science Foundation (SNSF).

References

Appendix A: Spherical Fourier-Bessel transform and 3D convolution products

A.1. Spherical Fourier-Bessel transform and relation to the 3D Fourier transform

In order to have a better understanding of the SFB coefficients and of how to use them to perform filtering, the SFB transform can be related to the 3D Fourier transform. We follow a similar definition as the one presented in Baddour (2010), but using our conventions for the different transforms.The following convention will be used for the Fourier transform: $\mathit{F}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}}\mathrm{\int }\mathit{f}\mathrm{\left(}{r}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{k}\mathrm{·}{r}}\mathrm{d}{r}\mathit{,}$(A.1)where F denotes the Fourier transform of f. This formulation does not assume any coordinate system. However, to relate this transform to the SFB transform, it is possible to express this equation in spherical coordinates using the following expansion for the Fourier kernel: ${\mathrm{e}}^{\mathrm{-}\mathrm{i}{k}\mathrm{·}{r}}\mathrm{=}\mathrm{4}\mathit{\pi }\sum _{\mathit{l}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathit{l}}\mathrm{(}\mathrm{-}\mathit{i}{\mathrm{)}}^{\mathit{l}}{\mathit{j}}_{\mathit{l}}\mathrm{(}\mathit{kr}\mathrm{\left)}\overline{{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{r}}\mathit{,}{\mathit{\phi }}_{\mathit{r}}\mathrm{\right)}}{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\right(}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{)}\mathit{,}$(A.2)where (k,θ_k,φ_k) and (r,θ_r,φ_r) are respectively the spherical coordinates of vectors k and r.

Substituting this expression for the kernel in the definition of the 3D Fourier transform yields: $\begin{array}{ccc}\mathit{F}\mathrm{\left(}\mathit{k,}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{4}\mathit{\pi }}{\sqrt{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathrm{\int }}_{\mathrm{\Omega }}\sum _{\mathit{l}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathit{l}}\mathrm{(}\mathrm{-}\mathit{i}{\mathrm{)}}^{\mathit{l}}\mathit{f}\mathrm{(}\mathit{r,}{\mathit{\theta }}_{\mathit{r}}\mathit{,}{\mathit{\phi }}_{\mathit{r}}\mathrm{)}\\ & & \mathrm{\times } {\mathit{j}}_{\mathit{l}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}\overline{{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{r}}\mathit{,}{\mathit{\phi }}_{\mathit{r}}\mathrm{\right)}}{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}\mathrm{d\Omega }{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r,}\\ & \mathrm{=}& \sum _{\mathit{l}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathit{l}}\mathrm{(}\mathrm{-}\mathit{i}{\mathrm{)}}^{\mathit{l}}[\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathrm{\int }}_{\mathrm{\Omega }}\mathit{f}\mathrm{\left(}\mathit{r,}{\mathit{\theta }}_{\mathit{r}}\mathit{,}{\mathit{\phi }}_{\mathit{r}}\mathrm{\right)}\\ & & \mathrm{\times } {\mathit{j}}_{\mathit{l}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}\overline{{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{r}}\mathit{,}{\mathit{\phi }}_{\mathit{r}}\mathrm{\right)}}\mathrm{d\Omega }{\mathit{r}}^{\mathrm{2}}\mathrm{d}{\mathit{r}}^{\mathrm{]}}{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}\mathit{,}\\ & \mathrm{=}& \end{array}$(A.3)In the last equation, the expression of the Spherical Harmonics Expansion of F(k,θ_k,φ_k) for a given value of k can be recognised from Eq. (4b). In the Fourier space, the are the spherical harmonics coefficients of F on a sphere of given radius k. In other words, the spherical harmonics coefficients F_lm(k) of the 3D Fourier transform F(k,θ_k,φ_k) on a sphere of given radius k in Fourier space are the SFB coefficients for the same value k but multiplied by factor ( − i)^l.

The relationship between 3D Fourier transform and SFB transform is therefore very simple. The SFB transform can be sought of as a mere Fourier transform in spherical coordinates. In the next sections, this relationship will be used to derive convolution and filtering relations for the SFB transform using the well known relations verified by the Fourier transform.

A.2. Expression of a 3D convolution product using the SFB transform

A prerequisite to the establishment of filtering relations is the expression of a 3D convolution product in terms of SFB coefficients. Let v(r,θ_r,φ_r) be the 3D convolution of f(r,θ_r,φ_r) and u(r,θ_r,φ_r). Then the 3D Fourier transform of v verifies: $\begin{array}{ccc}\mathit{V}\mathrm{\left(}\mathit{k,}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}& \mathrm{=}& \mathrm{ℱ}\mathrm{\{}\mathit{f}\mathrm{\ast }\mathit{u}\mathrm{\}}\mathrm{\left(}\mathit{k,}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}\mathit{,}\\ & \mathrm{=}& \sqrt{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\mathit{F}\mathrm{\left(}\mathit{k,}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}\mathit{U}\mathrm{\left(}\mathit{k,}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}\mathit{,}\end{array}$(A.4)where ℱ denotes the 3D Fourier transform. From Eq. (A.3) the expression of the 3D Fourier transform in spherical coordinates is known in terms of SFB coefficients. Applying this relationship to V(k,θ_k,φ_k) in the last equation yields: $\begin{array}{ccc}\mathrm{(}\mathrm{-}\mathit{i}{\mathrm{)}}^{\mathit{l}}\mathit{v̂}\mathit{lm}\mathrm{(}\mathit{k}\mathrm{)}& \mathrm{=}& \int \begin{array}{}\mathrm{2}\mathit{\pi }\\ \mathrm{0}\end{array}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\sqrt{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\mathit{F}\mathrm{\left(}\mathit{k,}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}\mathit{U}\mathrm{\left(}\mathit{k,}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}\\ & & \mathrm{\times } \overline{{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}}\sin \mathrm{\left(}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}\mathrm{d}{\mathit{\phi }}_{\mathit{k}}\mathrm{d}{\mathit{\theta }}_{\mathit{k}}\mathit{.}\end{array}$(A.5)Then, by applying Eq. (A.3) to F and U one gets: $\begin{array}{ccc}\mathit{v̂}\mathit{lm}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \mathrm{\left(}\mathit{i}{\mathrm{)}}^{\mathit{l}}\sqrt{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}”\sum _{{\mathit{l}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{{\mathit{m}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}{\mathit{l}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\right(}\mathrm{-}\mathit{i}{\mathrm{)}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathit{f̂}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathit{Y}}_{{\mathit{l}}^{\mathrm{\prime }}}^{{\mathit{m}}^{\mathrm{\prime }}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}\\ & & \mathrm{\times } \sum _{{\mathit{l}}^{\mathrm{\prime \prime }}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{{\mathit{m}}^{\mathrm{\prime \prime }}\mathrm{=}\mathrm{-}{\mathit{l}}^{\mathrm{\prime \prime }}}^{{\mathit{l}}^{\mathrm{\prime \prime }}}\mathrm{(}\mathrm{-}\mathit{i}{\mathrm{)}}^{{\mathit{l}}^{\mathrm{\prime \prime }}}\mathit{û}{\mathit{l}}^{\mathrm{\prime \prime }}{\mathit{m}}^{\mathrm{\prime \prime }}\mathrm{(}\mathit{k}\mathrm{\left)}{\mathit{Y}}_{{\mathit{l}}^{\mathrm{\prime \prime }}}^{{\mathit{m}}^{\mathrm{\prime \prime }}}\mathrm{\right(}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{)}\\ & & \mathrm{\times } \overline{{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}}\sin \mathrm{\left(}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}\mathrm{d}{\mathit{\phi }}_{\mathit{k}}\mathrm{d}{\mathit{\theta }}_{\mathit{k}}\mathit{,}\\ & \mathrm{=}& \mathrm{\left(}\mathit{i}{\mathrm{)}}^{\mathit{l}}\sqrt{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\sum _{{\mathit{l}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{{\mathit{m}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}{\mathit{l}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\right(}\mathrm{-}\mathit{i}{\mathrm{)}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathit{f̂}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\\ & & \mathrm{\times } \sum _{{\mathit{l}}^{\mathrm{\prime \prime }}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{{\mathit{m}}^{\mathrm{\prime \prime }}\mathrm{=}\mathrm{-}{\mathit{l}}^{\mathrm{\prime \prime }}}^{{\mathit{l}}^{\mathrm{\prime \prime }}}\mathrm{(}\mathrm{-}\mathit{i}{\mathrm{)}}^{{\mathit{l}}^{\mathrm{\prime \prime }}}\mathit{û}{\mathit{l}}^{\mathrm{\prime \prime }}{\mathit{m}}^{\mathrm{\prime \prime }}\mathrm{(}\mathit{k}\mathrm{)}\\ & & \mathrm{\times } ”{\mathit{Y}}_{{\mathit{l}}^{\mathrm{\prime }}}^{{\mathit{m}}^{\mathrm{\prime }}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}{\mathit{Y}}_{{\mathit{l}}^{\mathrm{\prime \prime }}}^{{\mathit{m}}^{\mathrm{\prime \prime }}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}\overline{{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{k}}\mathit{,}{\mathit{\phi }}_{\mathit{k}}\mathrm{\right)}}\mathrm{d}{\mathrm{\Omega }}_{\mathit{k}}\mathit{.}\end{array}$(A.6)The last integral over the two angular variables can be expressed as a Slater integral (which is a special case of the Gaunt integral) defined as: ${\mathit{c}}^{{\mathit{l}}^{\mathrm{\prime \prime }}}\mathrm{\left(}\mathit{l,m,}{\mathit{l}}^{\mathrm{\prime }}\mathit{,}{\mathit{m}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}”\overline{{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}}{\mathit{Y}}_{{\mathit{l}}^{\mathrm{\prime }}}^{{\mathit{m}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{\mathit{Y}}_{{\mathit{l}}^{\mathrm{\prime \prime }}}^{\mathit{m}\mathrm{-}{\mathit{m}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{d\Omega }\mathit{.}$(A.7)The Slater integrals are only nonzero for  | l − l′ |  ≤ l′′ ≤ l + l′ which simplifies the expression of .

The SFB transform of the 3D convolution product is therefore: $\begin{array}{ccc}{\begin{array}{}[\\ \mathrm{(}\mathit{f}\mathrm{\ast }\mathit{u}\mathrm{)}\end{array}}_{\mathit{lm}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \mathrm{\left(}\mathit{i}{\mathrm{)}}^{\mathit{l}}\sqrt{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\sum _{{\mathit{l}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{{\mathit{m}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}{\mathit{l}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\right(}\mathrm{-}\mathit{i}{\mathrm{)}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathit{f̂}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\\ & & \end{array}$(A.8)

All Figures

	Fig. 1 Scaling function and wavelet function for kc = 1. The y-axis represents the amplitude and the x-axis the frequency associated with the scaling function. The units of the amplitude must be those of the SFB transform.
In the text

	Fig. 2 Comparaison between spline and needlet wavelet functions on the sphere. The y-axis corresponds to amplitude and the x-axis to position.
In the text

	Fig. 3 Representation of the spherical 3D grid for the DSFBT (R = 1 and Nmax = 4) as described in Sect. 4 where the axes represent x,y,z and the grid points correspond to the centre of Healpix pixels and the radial grid is defined in Sect. 4.1.
In the text

	Fig. 4 Spherical 3D wavelet decomposition of a density field for a box with side 239.5   h-1   Mpc. Panel a) shows the reconstructed density field directly using the initial SFB coefficients. Panels b–e) show the first four wavelet scales, while panel f) shows the smoothed density field.
In the text

	Fig. 5 Wavelet Hard thresholding applied to a test density field.
In the text