The continuous wavelet derived by smoothing function and its application in cosmology
本站小编 Free考研考试/2022-01-02
Yun Wang1, Ping He,1,21College of Physics, Jilin University, Changchun 130012, China 2Center for High Energy Physics, Peking University, Beijing 100871, China
Abstract The wavelet analysis technique is a powerful tool and is widely used in broad disciplines of engineering, technology, and sciences. In this work, we present a novel scheme of constructing continuous wavelet functions, in which the wavelet functions are obtained by taking the first derivative of smoothing functions with respect to the scale parameter. Due to this wavelet constructing scheme, the inverse transforms are only one-dimensional integrations with respect to the scale parameter, and hence the continuous wavelet transforms (CWTs) constructed in this way are more ready to use than the usual scheme. We then apply the Gaussian-derived wavelet constructed by our scheme to computations of the density power spectrum for dark matter, the velocity power spectrum and the kinetic energy spectrum for baryonic fluid. These computations exhibit the convenience and strength of the CWTs. The transforms are very easy to perform, and we believe that the simplicity of our wavelet scheme will make CWTs very useful in practice. Keywords：wavelet analysis;intergalactic medium;large-scale structure of Universe
PDF (771KB)MetadataMetricsRelated articlesExportEndNote|Ris|BibtexFavorite Cite this article Yun Wang, Ping He. The continuous wavelet derived by smoothing function and its application in cosmology. Communications in Theoretical Physics, 2021, 73(9): 095402- doi:10.1088/1572-9494/ac10be
1. Introduction
It is well known that the Fourier transform has many shortcomings, of which the most is that it cannot simultaneously provide information about the scale and position of the signal [1]. To overcome the drawbacks of the Fourier transform, the wavelet analysis technique was invented and has proven a powerful tool, and is now widely used in areas such as signal processing, image and data compactness [2, 3]. The wavelet analysis has also been extensively applied to astrophysics and cosmology for more than three decades. The continuous wavelet transform (CWT) is applied to, to name a few, the analysis of the large-scale structures of the Universe [4–8], the detection of point sources or patterns of astronomical objects [9–14], the detection of the non-Gaussianity in the CMB maps [15–17], the foreground or noise subtraction in the CMB maps or sky surveys [18–20], the cosmological N-body simulation with high performances [21, 22]. Besides the CWT, the discrete wavelet transform (DWT) is also used in studies of cosmological large-scale structures [23–26].
As mentioned above, there are two kinds of wavelet transforms, namely, the CWT and the DWT. The CWT provides an overcomplete representation of a signal by varying continuously the scale and translation parameter of the wavelets. With a and b as scale and translation parameter respectively, the CWT of function f(x) is defined as$\begin{eqnarray}{{WT}}_{f}(a,b)=\displaystyle \frac{1}{\sqrt{| a| }}{\int }_{-\infty }^{\infty }f(x)\bar{\psi }\left(\displaystyle \frac{x-b}{a}\right){\rm{d}}x,\end{eqnarray}$where ψ(x) is called the mother wavelet, and the over-bar indicates the complex conjugate. ψ(x) should meet the square-integrable condition$\begin{eqnarray}{\int }_{-\infty }^{\infty }| \psi (x){| }^{2}{\rm{d}}x=\displaystyle \frac{1}{2\pi }{\int }_{-\infty }^{\infty }| \hat{\psi }(k){| }^{2}{\rm{d}}k\lt \infty ,\end{eqnarray}$in which $\hat{\psi }(k)$ is the Fourier transform of ψ(x), and we use Parseval’s theorem. The inverse transform of equation (1) is$\begin{eqnarray}f(x)={C}_{\psi }^{-1}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{{WT}}_{f}(a,b)\displaystyle \frac{1}{\sqrt{| a| }}\psi \left(\displaystyle \frac{x-b}{a}\right)\displaystyle \frac{{\rm{d}}a{\rm{d}}b}{{a}^{2}},\end{eqnarray}$with$\begin{eqnarray}{C}_{\psi }={\int }_{-\infty }^{\infty }\displaystyle \frac{| \hat{\psi }(k){| }^{2}}{| k| }{\rm{d}}k,\end{eqnarray}$where Cψ should satisfy the so-called admissibility condition 0 < Cψ < ∞, which suggests that $\hat{\psi }(0)=0$.
The advantage of the CWTs is that they can provide arbitrary localization in the scale and position due to their continuity. On the contrary, the DWTs do not vary continuously in the scale and translation parameter, and for a function f(x), its DWT usually takes the dyadic form as$\begin{eqnarray}{\tilde{\epsilon }}_{j,l}=\sqrt{\displaystyle \frac{{2}^{j}}{L}}{\int }_{-\infty }^{\infty }f(x)\bar{\psi }({2}^{j}x/L-l){\rm{d}}x,\end{eqnarray}$and f(x) can be expressed by the wavelet coefficients ${\tilde{\epsilon }}_{j,l}$ from the inverse transform as$\begin{eqnarray}f(x)=\sqrt{\displaystyle \frac{{2}^{j}}{L}}\sum _{j=0}^{\infty }\sum _{l=-\infty }^{\infty }{\tilde{\epsilon }}_{j,l}\psi ({2}^{j}x/L-l).\end{eqnarray}$For a comprehensive understanding of the CWT and the DWT, we refer the reader to [2, 3].
The advantage of the DWTs constructed in this way is that they can provide a set of complete and orthonormal bases ${\psi }_{j,l}(x)={\left({2}^{j}/L\right)}^{1/2}\psi ({2}^{j}x/L-l)$ for the wavelet expansion.
The disadvantage of the usual CWTs is that the inverse transform is a two-dimensional integration as equation (3), which may be computationally very complex and difficult; while DWTs are awkward to use since both spatial translation and scale dilation have to be dyadically adopted, that is, the dilation or translation cannot be arbitrarily performed. As a result, the expressions of the DWT power spectrum are complicated and cumbersome [26–28].
In this work, we present a novel scheme of constructing continuous wavelet functions, in which the wavelet functions are obtained by taking the first derivative of smoothing functions with respect to the scale parameter, and the inverse transforms are only one-dimensional (1D) integrations, and hence the CWTs constructed in this way are more ready to use than the usual forms. The paper is organized as follows. We outline the basic theoretical framework in section 2, and give the fast algorithm of our wavelet transform in section 3. We present some simple applications in cosmology of the wavelet technique in section 4. In section 5, we give some discussions of the results, and present the summary and conclusions in section 6.
2. Basic theoretical framework
We describe the process of constructing the CWTs in the following. As a comparison, the detailed derivation of the traditional CWT and the inverse transform are given explicitly in appendix A.
Generally, a smoothing function sw(x) with the scale parameter w is an even function of x, i.e. sw(x) = sw(−x), and should satisfy the normalization condition$\begin{eqnarray}{\int }_{-\infty }^{+\infty }{s}_{w}(x){\rm{d}}x=1.\end{eqnarray}$Instead of a general presentation, however, we proceed with a concrete example, i.e. the Gaussian function ${g}_{\sigma }(x)\,={{\rm{e}}}^{-{x}^{2}/2{\sigma }^{2}}/(\sqrt{2\pi }\sigma )$, with σ > 0. Letting $\sigma =\sqrt{2}/w$, we have$\begin{eqnarray}{g}_{w}(x)=\displaystyle \frac{w}{2\sqrt{\pi }}{{\rm{e}}}^{-\tfrac{1}{4}{w}^{2}{x}^{2}},\end{eqnarray}$which satisfies the normalization condition equation (7) and is surely a smoothing function. In the following, we call w the scale parameter. When w → ∞, gw(x) degenerates to Dirac δ function, i.e. gw(x) → δD(x). The Fourier transform of gw(x) is$\begin{eqnarray}{\hat{g}}_{w}(k)={{\rm{e}}}^{-\tfrac{{k}^{2}}{{w}^{2}}}.\end{eqnarray}$We see that ${\hat{g}}_{w}(0)=1$, consistent with the normalization condition equation (7). Additionally, the reason for choosing w as $w=\sqrt{2}/\sigma $ is that w can one-to-one correspond to k, as can be seen from equation (9). In figure 1, we show the curves of gw(x) and ${\hat{g}}_{w}(k)$.
Figure 1.
New window|Download| PPT slide Figure 1.The upper panel is for the Gaussian and the wavelet function; the lower panel is for their corresponding Fourier transforms. We see that the Gaussian is a low-pass filter, while the wavelet is a band-pass filter.
For a random field, such as the density contrast δ(x) of the cosmological density, the smoothed field δg(w, x) under the scale parameter w can be obtained via the convolution with the smoothing function gw(x) as$\begin{eqnarray}{\delta }_{g}(w,x)={\int }_{-\infty }^{+\infty }\delta (u){g}_{w}(x-u){\rm{d}}u.\end{eqnarray}$When w → ∞, then gw(x) → δD(x), and hence δg(w, x) → δ(x). Taking the first derivative of both sides of equation (10) with respect to w, we have$\begin{eqnarray}\begin{array}{rcl}{W}_{\delta }(w,x) & \equiv & \displaystyle \frac{\partial {\delta }_{g}(w,x)}{\partial w}={\displaystyle \int }_{-\infty }^{+\infty }\delta (u)\displaystyle \frac{\partial {g}_{w}(x-u)}{\partial w}{\rm{d}}u\\ & = & {\displaystyle \int }_{-\infty }^{+\infty }\delta (u){{\rm{\Psi }}}_{g}(w,x-u){\rm{d}}u,\end{array}\end{eqnarray}$in which we define $\Psi$g as$\begin{eqnarray}{{\rm{\Psi }}}_{g}(w,x)\equiv \displaystyle \frac{\partial {g}_{w}(x)}{\partial w}=\displaystyle \frac{1}{4\sqrt{\pi }}(2-{w}^{2}{x}^{2}){{\rm{e}}}^{-\tfrac{1}{4}{w}^{2}{x}^{2}},\end{eqnarray}$which is the Gaussian-derived wavelet that we call, and is nothing but the 1D Mexican hat wavelet function. The Fourier transform of equation (12) is$\begin{eqnarray}{\hat{{\rm{\Psi }}}}_{g}(w,k)=\displaystyle \frac{2{k}^{2}}{{w}^{3}}{{\rm{e}}}^{-\tfrac{{k}^{2}}{{w}^{2}}}.\end{eqnarray}$We also plot $\Psi$g(w, x) and ${\hat{{\rm{\Psi }}}}_{g}(w,k)$ in figure 1, from which we see that ${\hat{g}}_{w}(k)$ is a low-pass filter, while ${\hat{{\rm{\Psi }}}}_{g}(w,k)$ is a band-pass filter. Note that ${\hat{g}}_{w}(k)$ and ${\hat{{\rm{\Psi }}}}_{g}(w,k)$ are related by$\begin{eqnarray}{\hat{{\rm{\Psi }}}}_{g}(w,k)=\displaystyle \frac{\partial {\hat{g}}_{w}(k)}{\partial w},\end{eqnarray}$which is consistent with the definition of $\Psi$g(w, x) in equation (12). Integrating equation (11) with respect to w, we have$\begin{eqnarray}{\delta }_{g}(w,x)={\delta }_{g}(0,x)+{\int }_{0}^{w}{W}_{\delta }(w^{\prime} ,x){\rm{d}}w^{\prime} ,\end{eqnarray}$where δg(0, x) is an integration constant. For cosmologically interesting objects, however, such as the density contrast field or the peculiar velocity field, the whole-cosmos averaged quantities should be usually vanishing, i.e. δg(0, x) = 0, or at least a constant that is independent of spatial positions, and hence we can safely neglect this term in the future. Notice that δ(x) = δg(∞, x), we obtain$\begin{eqnarray}\delta (x)={\int }_{0}^{\infty }{W}_{\delta }(w,x){\rm{d}}w,\end{eqnarray}$which is just the inverse transform of the wavelet transform equation (11). Compared with the usual inverse wavelet transform equation (3), we see that our inverse transform is only 1D integration, which is much easier to manipulate than equation (3).
3. Fast algorithm of wavelet transform
Due to the wavelet transform pair equations (11) and (16), Wδ(w, x) is equivalent to the original field δ(x), and hence we can use Wδ(w, x) for further studies instead of δ(x). Given a density contrast (or other) field δ(x), we can calculate its wavelet transform directly using equation (11). However, there exists a fast algorithm based on the technique of fast Fourier transform (FFT). From equation (11), we use the convolution theorem and obtain$\begin{eqnarray}{\hat{W}}_{\delta }(w,k)=\hat{\delta }(k){\hat{{\rm{\Psi }}}}_{g}(w,k),\end{eqnarray}$where quantities with a hat are Fourier transforms of corresponding quantities. ${\hat{{\rm{\Psi }}}}_{g}(w,k)$ is given by equation (13), $\hat{\delta }(k)$ can be obtained by FFT from δ(x), and hence Wδ(w, x) can be obtained by the inverse FFT from equation (17).
Taking the Fourier transform of both sides of equation (16), we obtain$\begin{eqnarray}\hat{\delta }(k)={\int }_{0}^{\infty }{\hat{W}}_{\delta }(w,k){\rm{d}}w.\end{eqnarray}$Notice that$\begin{eqnarray}{\int }_{0}^{\infty }{\hat{{\rm{\Psi }}}}_{g}(w,k){\rm{d}}w=1.\end{eqnarray}$Equation (18) can also be derived by integrating both sides of equation (17) with respect to w.
We see that the Fourier transform pair equations (17) and (18) are very simple and concise, and they compose the fast algorithm of wavelet transform.
4. Applications in cosmology
We apply our Gaussian-derived wavelet to analyze the IllustrisTNG simulation data [29–34], from which we select the sample IllustrisTNG100-1, whose simulation box is 75Mpc/h long. We use ‘cloud-in-cell’ scheme [35] to assign all the particle (dark matter and baryonic) mass or velocity into a 10243 mesh to acquire mass density or velocity at mesh points. Then we randomly select 100 000 lines, all vertical to the x–y plane of the simulation box, and record all the relevant data at each line. In this way, we have 100 000 1D data. In figure 2, we show the dark matter density contrast of one such line data, together with its scalogram Wδ(w, x).
Figure 2.
New window|Download| PPT slide Figure 2.Upper panel: the 1D density contrast field of dark matter. Lower panel: the corresponding wavelet scalogram of the 1D density contrast field. The data is taken from the IllustrisTNG simulation.
With the data, we compute the wavelet power spectrum of the density contrast field or peculiar velocity field. For one line data, the wavelet power spectrum can be defined as$\begin{eqnarray}{S}_{i}(w)\equiv \displaystyle \frac{1}{{L}_{b}^{2}}\int | {W}_{i}(w,x){| }^{2}{\rm{d}}x=\displaystyle \frac{1}{{N}_{p}{L}_{b}}\sum _{j=1}^{{N}_{p}}| {W}_{i}(w,{x}_{j}){| }^{2},\end{eqnarray}$where i is the ID of the ith line, Wi(w, x) is the wavelet transform of the line data, Lb = 75 Mpc/h is the length of the simulation box, and Np = 1024 is the point number of the line data. Averaging over all the 100 000 lines, we obtain the total power spectrum as$\begin{eqnarray}S(w)=\displaystyle \frac{1}{{N}_{l}}\sum _{i=1}^{{N}_{l}}{S}_{i}(w),\end{eqnarray}$in which Nl = 100 000.
We use equation (21) to compute the 1D wavelet power spectrum of the line data of dark matter density contrast, and show the power spectrum Sd(w) in figure 3, together with the 1D Fourier power spectrum Pd(k).
Figure 3.
New window|Download| PPT slide Figure 3.The Fourier and wavelet power spectrum of the 1D density contrast field for dark matter. The left vertical axis is for the Fourier power spectrum, and the right axis is for the wavelet power spectrum. The data is taken from the IllustrisTNG simulation.
According to Parseval’s theorem, we have$\begin{eqnarray}\begin{array}{rcl}\displaystyle \int | {W}_{\delta }(w,x){| }^{2}{\rm{d}}x & = & \displaystyle \frac{1}{2\pi }\displaystyle \int | {\hat{W}}_{\delta }(w,k){| }^{2}{\rm{d}}k\\ & = & \displaystyle \frac{1}{2\pi }\displaystyle \int | \hat{\delta }(k){| }^{2}| {\hat{{\rm{\Psi }}}}_{g}(w,k){| }^{2}{\rm{d}}k,\end{array}\end{eqnarray}$in which for the second equality we use equation (17). Equation (22) can provide an alternative method to compute the wavelet power spectrum in equation (20).
We can use equation (22) to derive an interesting relationship between the Fourier and wavelet power spectrum. We assume the Fourier power spectrum Pd(k) scales as a power law of k, i.e. ${P}_{{\rm{d}}}(k)\sim | \hat{\delta }(k){| }^{2}\sim {k}^{\alpha }$, in which α is the power index. If α > −5, then from equation (22) we have$\begin{eqnarray}\begin{array}{rcl}{\tilde{S}}_{{\rm{d}}}(w) & \sim & \displaystyle \int | {\hat{W}}_{\delta }(w,x){| }^{2}{\rm{d}}x\\ & \sim & \displaystyle \int {P}_{{\rm{d}}}(k){\hat{{\rm{\Psi }}}}_{g}{\left(w,k\right)}^{2}{\rm{d}}k\sim \displaystyle \frac{{P}_{{\rm{d}}}(w)}{w}.\end{array}\end{eqnarray}$From figure 3, we see that with a proper amplitude, the relationship by equation (23) is quite accurate when k < 10h/Mpc.
With the wavelet transform Wv(w, x) of the z-directional 1D velocity field for baryonic fluid, the wavelet power spectrum of the velocity field Sv(w), can be obtained in the same way as Sd(w) of the density contrast field. We show Sv(w) in figure 4, together with the Fourier power spectrum Pv(k) of the 1D velocity.
Figure 4.
New window|Download| PPT slide Figure 4.The Fourier and wavelet power spectrum of the 1D velocity field for baryonic fluid. The left vertical axis is for the Fourier power spectrum, and the right axis is for the wavelet power spectrum. The wavelet energy spectrum is also shown for comparison. The data is taken from the IllustrisTNG simulation.
Turbulent flows in the baryonic fluid of the Universe are an important area in cosmological studies [36–43]. The wavelet analysis technique is particularly suitable to investigate turbulent flows [44]. In the scale (or Fourier) space, the Kolmogorov theory assumes the existence of an energy cascade between the different excited wavenumbers of the turbulent flow, while in the physical (or real) space, turbulent flows are characterized by complex multi-scale and chaotic motions, which can be classified into more elementary components, namely the coherent structures3(3https://en.wikipedia.org/wiki/Coherent_turbulent_structure.)
To study turbulence in the baryonic fluid of the Universe, which is distributed greatly inhomogeneously in space, the kinetic energy spectrum of the fluid is more desirable than its velocity spectrum. However, one cannot define the usual power spectrum for kinetic energy [45]. The reason can be seen from below$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{1}{2}\langle \rho (x){v}^{2}(x)\rangle =\displaystyle \frac{1}{2}\displaystyle \int \rho (x){v}^{2}(x){\rm{d}}x\\ =\,\displaystyle \frac{1}{8{\pi }^{2}}\displaystyle \int {\hat{\rho }}^{* }({k}_{1}+{k}_{2})\hat{v}({k}_{1})\hat{v}({k}_{2}){\rm{d}}{k}_{1}{\rm{d}}{k}_{2}\\ =\,\displaystyle \frac{1}{8{\pi }^{2}}\displaystyle \int {B}_{\rho v}({k}_{1},{k}_{2}){\rm{d}}{k}_{1}{\rm{d}}{k}_{2},\end{array}\end{eqnarray}$where we do not have a power spectrum and we have to define a mixed bispectrum as ${B}_{\rho v}({k}_{1},{k}_{2})\equiv {\hat{\rho }}^{* }({k}_{1}\,+{k}_{2})\hat{v}({k}_{1})\hat{v}({k}_{2})$ [46]. With the wavelet transform, however, we can define a kinetic energy power spectrum for a 1D data readily as below$\begin{eqnarray}{S}_{{\rm{E}},i}(w)=\displaystyle \frac{1}{2{L}_{b}^{2}}\int {{\rm{\Delta }}}_{i}(x)| {W}_{v,i}(w,x){| }^{2}{\rm{d}}x,\end{eqnarray}$in which ${\rm{\Delta }}(x)={\rho }_{b}(x)/{\bar{\rho }}_{b}$ is the dimensionless baryonic density. Averaging over the total 100 000 lines by using equation (21), we can obtain the wavelet kinetic energy power spectrum SE(w), which is also shown in figure 4. From another viewpoint, one can consider SE(w) as the density-weighted velocity power spectrum. Moreover, due to the localization property of the wavelet transform, SE(w) can be generalized to the spectrum with the space-restricted (e.g. x1 < x < x2) or the density-restricted (e.g. ρ1 < ρ < ρ2) integration, which may be more suitable for investigations of inhomogeneous turbulent flows.
5. Discussions
Our wavelet scheme is a general method, which can be used to construct a large category of CWTs, but we emphasize that equations (12) and (14) are not sufficient conditions to construct a continuous wavelet. We show this with a counter-example in the following. The top-hat function in real space with the scale parameter w,$\begin{eqnarray}{h}_{w}(x)=\left\{\begin{array}{ll}\displaystyle \frac{w}{2}, & -\displaystyle \frac{1}{w}\leqslant x\leqslant \displaystyle \frac{1}{w},\\ 0, & \mathrm{otherwise},\end{array}\right.\end{eqnarray}$is a smoothing function, whose Fourier transform is ${\hat{h}}_{w}(k)=(w/k)\sin (k/w)$. According to our scheme, however, its derivative with respect to w,$\begin{eqnarray}{\hat{{\rm{\Psi }}}}_{h}(w,k)\equiv \displaystyle \frac{\partial {\hat{h}}_{w}(k)}{\partial w}=\displaystyle \frac{1}{k}\sin \left(\displaystyle \frac{k}{w}\right)-\displaystyle \frac{1}{w}\cos \left(\displaystyle \frac{k}{w}\right),\end{eqnarray}$is NOT a continuous wavelet, since ${\hat{{\rm{\Psi }}}}_{h}(w,k)$ is oscillatory when k → ∞, and hence does not satisfy the square-integrable condition equation (2). The real space couterpart of equation (27) is$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Psi }}}_{h}(w,x) & \equiv & \displaystyle \frac{\partial {h}_{w}(x)}{\partial w}=-\displaystyle \frac{1}{2w}\left[{\delta }_{D}\left(x+\displaystyle \frac{1}{w}\right)+{\delta }_{D}\left(x-\displaystyle \frac{1}{w}\right)\right]\\ & & +\left\{\begin{array}{ll}\displaystyle \frac{1}{2}, & -\displaystyle \frac{1}{w}\lt x\lt \displaystyle \frac{1}{w},\\ 0, & \mathrm{otherwise}.\end{array}\right.\end{array}\end{eqnarray}$To facilitate understanding of this counter-example, we show $\Psi$h(w, x) and ${\hat{{\rm{\Psi }}}}_{h}(w,k)$ in figure 5.
Figure 5.
New window|Download| PPT slide Figure 5.The functional form of the counter-example (upper panel), and its Fourier transform (lower panel) in section 5, shown as blue lines. The Gaussian-derived wavelet and its Fourier transform are also shown for comparison (red lines).
Generally, since the smoothing function sw(x) is an even function, its Fourier transform can be formally expressed as a Taylor expansion as follows,$\begin{eqnarray}{\hat{s}}_{w}(k)=1+\displaystyle \frac{{\hat{s}}^{(2)}(0)}{2!}{\left(\displaystyle \frac{k}{w}\right)}^{2}+\displaystyle \frac{{\hat{s}}^{(4)}(0)}{4!}{\left(\displaystyle \frac{k}{w}\right)}^{4}+...,\end{eqnarray}$from which we know ${\hat{s}}_{w}(0)=1$, consistent with the normalization condition equation (7). From equation (29), we can define a function $\Psi$s(w, x) by its Fourier transform as$\begin{eqnarray}{\hat{{\rm{\Psi }}}}_{s}(w,k)\equiv \displaystyle \frac{\partial {\hat{s}}_{w}(k)}{\partial w}=-\displaystyle \frac{k}{w}\displaystyle \frac{\partial {\hat{s}}_{w}(k)}{\partial k}.\end{eqnarray}$It is easy to see that ${\hat{{\rm{\Psi }}}}_{s}(w,k)$ satisfies the admissibility condition of wavelet, i.e. ${\hat{{\rm{\Psi }}}}_{s}(w,0)=0$. If it also satisfies the square-integrable condition equation (2), then $\Psi$s(w, x) should be a continuous wavelet.
Additionally, it is not difficult to generalize our scheme to the three-dimensional (3D) case. For example, the 3D anisotropic Gaussian function is$\begin{eqnarray}{g}_{{\boldsymbol{w}}}({\boldsymbol{x}})=\prod _{i=1}^{3}{g}_{{w}_{i}}({x}_{i})=\displaystyle \frac{{w}_{1}{w}_{2}{w}_{3}}{8{\pi }^{3/2}}{{\rm{e}}}^{-\tfrac{1}{4}({w}_{1}^{2}{x}_{1}^{2}+{w}_{2}^{2}{x}_{2}^{2}+{w}_{3}^{2}{x}_{3}^{2})},\end{eqnarray}$where w = (w1, w2, w3), x = (x1, x2, x3). Hence the 3D anisotropic Gaussian-derived wavelet can be obtained as$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Psi }}}_{g}({\boldsymbol{w}},{\boldsymbol{x}}) & \equiv & \displaystyle \frac{{\partial }^{3}{g}_{{\boldsymbol{w}}}({\boldsymbol{x}})}{\partial {w}_{1}\partial {w}_{2}\partial {w}_{3}}\\ & = & \displaystyle \prod _{i=1}^{3}\displaystyle \frac{\partial {g}_{{w}_{i}}({x}_{i})}{\partial {w}_{i}}=\displaystyle \prod _{i=1}^{3}{{\rm{\Psi }}}_{g}({w}_{i},{x}_{i}).\end{array}\end{eqnarray}$We see that the 3D anisotropic wavelet is not a 3D Mexican hat wavelet. In the future, we will explore the possible applications of this 3D wavelet in cosmology.
6. Summary and conclusions
The DWTs are constructed by dilation and translation both dyadically in the scale and position. They can provide a set of complete and orthonormal bases, based on which fast algorithms can be designed, and hence DWTs are very suitable for technical applications, such as data compactness, or image processing. While in some areas, such as cosmological investigations, the localization of continuous wavelets in both scale and position is much more desirable than the orthogonality of discrete wavelets. Nevertheless, the usual CWTs are not convenient to use since the inverse wavelet transforms are two-dimensional integrations, which may be computationally very cumbersome and time-consuming.
In this work, we present a novel scheme of constructing continuous wavelet functions, in which the wavelet functions are obtained by taking the first derivative of smoothing functions with respect to the scale parameter. Due to this wavelet constructing scheme, the inverse transforms are only 1D integrations with respect to the scale parameter, and hence CWTs are more ready to use than the usual scheme.
We then apply the Gaussian-derived wavelet constructed by our scheme to computations of the density power spectrum for dark matter, the velocity power spectrum and the kinetic energy spectrum for baryonic fluid. These computations exhibit the convenience and strength of the CWTs. From the transform pairs equations (11) and (16) in real space, and equations (17) and (18) in Fourier space, we see that our wavelet transform scheme is very simple, and we believe that the simplicity of our scheme will make CWTs very useful in cosmology.
Appendix A. Derivation of the traditional CWT
Given a wavelet function ψ(x), the traditional CWT of a function f(x) is defined as the convolution of f(x) with a scaled version of ψ(x). Or to put it another way, the traditional CWT can be viewed as the projection of f(x) on the dilation and translation of $\Psi$(x) [47]. The traditional CWT is shown in equation (1), and again as follows$\begin{eqnarray}{{WT}}_{f}(a,b)=\displaystyle \frac{1}{\sqrt{| a| }}{\int }_{-\infty }^{\infty }f(x)\bar{\psi }(\displaystyle \frac{x-b}{a}){\rm{d}}x.\end{eqnarray}$With Parseval’s theorem for traditional CWT (see theorem 3.10 in [3] for its proof)$\begin{eqnarray}\begin{array}{l}{C}_{\psi }{\displaystyle \int }_{-\infty }^{+\infty }f(x)\bar{g}(x){\rm{d}}x\\ \quad ={\displaystyle \int }_{-\infty }^{+\infty }{\displaystyle \int }_{-\infty }^{+\infty }{{WT}}_{f}(a,b){\overline{{WT}}}_{g}(a,b)\displaystyle \frac{{\rm{d}}a{\rm{d}}b}{{a}^{2}},\end{array}\end{eqnarray}$where Cψ is given in equation (4). Let g(x) be the Dirac delta function, i.e. $g(x)={\delta }_{D}(x-x^{\prime} )$, and from equation (A.2) we obtain$\begin{eqnarray}f(x)={C}_{\psi }^{-1}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{{WT}}_{f}(a,b)\displaystyle \frac{1}{\sqrt{| a| }}\psi (\displaystyle \frac{x-b}{a})\displaystyle \frac{{\rm{d}}a{\rm{d}}b}{{a}^{2}},\end{eqnarray}$which is the reconstruction formula or inverse transform of the conventional method. Notice that, although both traditional CWT equation (A.1) and our CWT equation (11) are the convolution of a general function with a wavelet function, the traditional reconstruction formula equation (A.3) is far more complex than our method equation (16).
It should be pointed out that the traditional CWT does not tell us how to construct a specific wavelet function. In fact, each existing continuous wavelet is constructed in some unique way. For example, the well-known Mexican hat wavelet is defined as the negative normalized second derivative of the Gaussian function, i.e.$\begin{eqnarray}{{\rm{\Psi }}}_{\mathrm{MH}}(x)\equiv C\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{x}^{2}}\left[\displaystyle \frac{1}{\sigma \sqrt{2\pi }}{{\rm{e}}}^{-\tfrac{{x}^{2}}{2{\sigma }^{2}}}\right],\end{eqnarray}$where C is a normalization constant, such that ${\int }_{-\infty }^{+\infty }| \psi (x){| }^{2}{\rm{d}}x=1$. If the standard deviation σ is set to $\sqrt{2}$, then we have$\begin{eqnarray}{{\rm{\Psi }}}_{\mathrm{MH}}(x)=\displaystyle \frac{1}{{\left(2\pi \right)}^{1/4}\sqrt{3}}(2-{x}^{2}){{\rm{e}}}^{-\tfrac{{x}^{2}}{4}},\end{eqnarray}$which is the same as our wavelet equation (12), except for the prefactor. However, the original function cannot be reconstructed by integrating WTf(a, b) with respect to the scale parameter a directly. In order to better understand this point, we perform the Fourier transform of equation (A.1) as$\begin{eqnarray}{\widehat{{WT}}}_{f}(a,k)=\hat{f}(k)\displaystyle \frac{a}{\sqrt{| a| }}{\hat{{\rm{\Psi }}}}_{\mathrm{MH}}({ak}),\end{eqnarray}$where ${\hat{{\rm{\Psi }}}}_{\mathrm{MH}}(k)=4{\left(8\pi /9\right)}^{1/4}{k}^{2}{{\rm{e}}}^{-{k}^{2}}$. Obviously, due to ${\int }_{-\infty }^{+\infty }(a/\sqrt{| a| }){\hat{{\rm{\Psi }}}}_{\mathrm{MH}}({ak}){\rm{d}}a=0$, we cannot recover the original function $\hat{f}(k)$ with ${\int }_{-\infty }^{+\infty }{\widehat{{WT}}}_{f}(a,k){\rm{d}}a$.
As another example, we describe briefly how to design the Meyer wavelet (see [47] for details). The Meyer wavelet is constructed by means of the multiresolution analysis in the Fourier space, which is much more tedious than constructing the Mexican hat wavelet. The procedure consists of two steps, i.e. constructing the scaling function with the conjugate mirror filter and then constructing the wavelet function with the wavelet equation. The conjugate mirror filter $\hat{h}(k)$ used here is given by$\begin{eqnarray}\hat{h}(k)=\left\{\begin{array}{ll}\sqrt{2}, & | k| \leqslant \displaystyle \frac{\pi }{3},\\ \sqrt{2}\cos \left[\displaystyle \frac{\pi }{2}\beta \left(\displaystyle \frac{3| k| }{\pi }-1\right)\right], & \displaystyle \frac{\pi }{3}\lt | k| \leqslant \displaystyle \frac{2\pi }{3},\\ 0, & \displaystyle \frac{2\pi }{3}\lt | k| \leqslant \pi ,\end{array}\right.\end{eqnarray}$where β(t) is an auxiliary function such that $| \hat{h}(k){| }^{2}+| \hat{h}(k\,+\pi ){| }^{2}\,=\,2$, which indicates the scaling functions are orthonormal under integer translations. The simplest case of β(t) [48] is$\begin{eqnarray}\beta (x)=\left\{\begin{array}{ll}0, & x\leqslant 0,\\ x, & 0\lt x\leqslant 1,\\ 1, & x\gt 1.\end{array}\right.\end{eqnarray}$Cutting off the terms of p > 1 in the scaling function $\hat{\phi }(k)={\prod }_{p=1}^{+\infty }\hat{h}({2}^{-p}k)/\sqrt{2}$ and letting $\hat{\phi }(k)=0$ for ∣k∣ >4π/3, we have$\begin{eqnarray}\hat{\phi }(k)=\left\{\begin{array}{ll}1, & | k| \leqslant \displaystyle \frac{2\pi }{3},\\ \sin \left(\displaystyle \frac{3| k| }{4}\right), & \displaystyle \frac{2\pi }{3}\lt | k| \leqslant \displaystyle \frac{4\pi }{3},\\ 0, & | k| \gt \displaystyle \frac{4\pi }{3},\end{array}\right.\end{eqnarray}$which is the Meyer scaling function. Then substituting equations (A.7) and (A.9) into the wavelet equation, $\hat{\psi }(k)=\hat{g}(k/2)\hat{\phi }(k/2)/\sqrt{2}$, where $\hat{g}(k)={{\rm{e}}}^{-{\rm{i}}{k}}\bar{\hat{h}}(k+\pi )$, we obtain the Meyer wavelet function as below,$\begin{eqnarray}{\hat{{\rm{\Psi }}}}_{\mathrm{MY}}\left(k\right)\equiv \left\{\begin{array}{ll}-\cos \left(\displaystyle \frac{3| k| }{4}\right){{\rm{e}}}^{-{\rm{i}}{k}/2}, & 2\pi /3\lt | k| \leqslant 4\pi /3,\\ \sin \left(\displaystyle \frac{3| k| }{8}\right){{\rm{e}}}^{-{\rm{i}}{k}/2}, & 4\pi /3\lt | k| \leqslant 8\pi /3,\\ 0, & \mathrm{otherwise},\end{array}\right.\end{eqnarray}$which is a complex function and incompatible with our scheme.
Unlike conventional continuous wavelet functions designed in various ways, our scheme offers a uniform approach to construct wavelets by differentiating a smoothing function with respect to its scale parameter w. Because of this wavelet-constructing method, the original signal can be reconstructed easily by integrating the CWT W(w, x) with respect to w, and hence, our method is much simpler than the traditional method.
Appendix B. Comparison with wavelets derived from different smoothing functions
In this work, we use the Gaussian function as the smoothing function to derive the wavelet function, due to its simplicity and wide applications. The Gaussian-derived wavelet is localized in both real and Fourier space. In other words, both $\Psi$g(w, x) and ${\hat{{\rm{\Psi }}}}_{g}(w,k)$ are approximately compactly supported, as shown in figure 1. Besides the Gaussian function, we also consider other bell-shaped functions, such as the Meyer scaling function [48] and the Epanechnikov function [49], as the smoothing function. The expressions of Meyer scaling function and Epanechnikov function are$\begin{eqnarray}{m}_{w}\left(x\right)=\left\{\begin{array}{ll}\left(\displaystyle \frac{1}{2\pi }+\displaystyle \frac{1}{{\pi }^{2}}\right)w, & x=0,\\ \displaystyle \frac{\pi \sin \left({wx}/2\right)+{wx}\cos \left({wx}\right)}{{\pi }^{2}x-{w}^{2}{x}^{3}}, & x\ne 0,\end{array}\right.\end{eqnarray}$and$\begin{eqnarray}{\varepsilon }_{w}\left(x\right)=\left\{\begin{array}{ll}\displaystyle \frac{3w}{4A}\left(1-{\left(\displaystyle \frac{{wx}}{A}\right)}^{2}\right), & -\displaystyle \frac{A}{w}\leqslant x\leqslant \displaystyle \frac{A}{w},\\ 0, & \mathrm{otherwise},\end{array}\right.\end{eqnarray}$respectively, where the constant A is approximately equal to 3.342. We will explain the reason of such a choice below. The Fourier transforms of equations (B.1) and (B.2) are given by$\begin{eqnarray}{\hat{m}}_{w}\left(k\right)=\left\{\begin{array}{ll}1, & | k| \leqslant \displaystyle \frac{w}{2},\\ \sin \left(\displaystyle \frac{\pi | k| }{w}\right), & \displaystyle \frac{w}{2}\lt | k| \leqslant w,\\ 0, & \mathrm{otherwise},\end{array}\right.\end{eqnarray}$and$\begin{eqnarray}{\hat{\varepsilon }}_{w}\left(k\right)=\left\{\begin{array}{ll}1, & k=0,\\ \displaystyle \frac{-(3{Ak}/w)\cos ({Ak}/w)+3\sin ({Ak}/w)}{{\left({Ak}/w\right)}^{3}}, & k\ne 0.\end{array}\right.\end{eqnarray}$The corresponding shapes of these smoothing functions in real space are shown in the upper left panel of figure B1, and the lower left panel is for their Fourier transforms.
Figure B1.
New window|Download| PPT slide Figure B1.Three different smoothing functions and their corresponding wavelets. Upper left panel: Gaussian (red line), Meyer scaling (blue line) and Epanechnikov (green line) function in real space. Upper right panel: the wavelet functions derived from these smoothing functions. Lower left panel: Fourier transforms of the three smoothing functions. Lower right panel: Fourier transforms of the wavelets derived from the corresponding smoothing functions.
Next, according to our scheme, we define the wavelets as the first derivative of smoothing functions with respect to w. Hence the Meyer-derived wavelet in real space is$\begin{eqnarray}{{\rm{\Psi }}}_{m}\left(w,x\right)\equiv \displaystyle \frac{\partial {m}_{w}(x)}{\partial w}=\left\{\begin{array}{ll}-\left(\displaystyle \frac{3}{16}+\displaystyle \frac{1}{4{\pi }^{2}}\right), & x=\pm \displaystyle \frac{\pi }{w},\\ {{\rm{\Psi }}}_{0}(w,x), & x\ne \pm \displaystyle \frac{\pi }{w},\end{array}\right.\end{eqnarray}$where the form of $\Psi$0(w, x) is$\begin{eqnarray}\begin{array}{l}{{\rm{\Psi }}}_{0}(w,x)=\displaystyle \frac{\cos \left({wx}\right)-\sin \left(\tfrac{{wx}}{2}\right)}{2{\left({wx}+\pi \right)}^{2}}+\displaystyle \frac{\sin \left(\tfrac{{wx}}{2}\right)+\cos \left({wx}\right)}{2{\left({wx}-\pi \right)}^{2}}\\ \quad +\displaystyle \frac{2\sin \left({wx}\right)-\cos \left(\tfrac{{wx}}{2}\right)}{4\left({wx}-\pi \right)}+\displaystyle \frac{2\sin \left({wx}\right)+\cos \left(\tfrac{{wx}}{2}\right)}{4\left({wx}+\pi \right)}.\end{array}\end{eqnarray}$The Epanechnikov-derived wavelet function is$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Psi }}}_{\varepsilon }\left(w,x\right) & \equiv & \displaystyle \frac{\partial {\varepsilon }_{w}(x)}{\partial w}\\ & = & \left\{\begin{array}{ll}\displaystyle \frac{3}{4A}\left[1-3{\left(\displaystyle \frac{{wx}}{A}\right)}^{2}\right], & -\displaystyle \frac{A}{w}\leqslant x\leqslant \displaystyle \frac{A}{w},\\ 0, & \mathrm{otherwise}.\end{array}\right.\end{array}\end{eqnarray}$The Fourier transforms of equations (B.5) and (B.7) are given as$\begin{eqnarray}{\hat{{\rm{\Psi }}}}_{m}\left(w,k\right)=\left\{\begin{array}{ll}-\displaystyle \frac{\pi \left|k\right|\cos \left(\tfrac{\pi \left|k\right|}{w}\right)}{{w}^{2}}, & \displaystyle \frac{w}{2}\lt | k| \leqslant w,\\ 0, & \mathrm{otherwise},\end{array}\right.\end{eqnarray}$and$\begin{eqnarray}{\hat{{\rm{\Psi }}}}_{\varepsilon }\left(w,k\right)=\left\{\begin{array}{ll}0, & k=0,\\ -\displaystyle \frac{3\left[{\left(\tfrac{{Ak}}{w}\right)}^{2}-3\right]\sin \left(\tfrac{{Ak}}{w}\right)+9\left(\tfrac{{Ak}}{w}\right)\cos \left(\tfrac{{Ak}}{w}\right)}{w{\left(\tfrac{{Ak}}{w}\right)}^{3}}, & k\ne 0,\end{array}\right.\end{eqnarray}$respectively. It is easy to verify that both the Meyer- and the Epanechnikov-derived wavelet satisfy the admissibility condition and the square-integrable condition. Notice that the scale parameter w from equations (B.1) to (B.9) is defined as the frequency at which the first peak of $\hat{{\rm{\Psi }}}(w,k)$ is located when k/w = ±1, i.e. the frequency at which $\hat{{\rm{\Psi }}}(w,k)$ takes the maximum value when k/w = ±1, as shown in the lower right panel in figure B1. That is the reason we choose A = 3.342 for Epanechnikov-derived wavelet.
With these two wavelets, we can now compare their properties with those of the Gaussian-derived wavelet. By visual inspection of the upper right panel of figure B1, we find that when x is farther away from the origin in real space, the Meyer-derived wavelet is still oscillatory obviously while other wavelets are already close to zero. On the other hand, as shown in the lower right panel of figure B1, the Epanechnikov-derived wavelet is more extended in the Fourier space than others. Hence, compared with the Gaussian-derived wavelet, the Meyer- and the Epanechnikov-derived wavelet are not well localized in both real and Fourier space simultaneously, which degrades their performance and applicability in practice.
Appendix C. Comparing the power spectrum computed with the traditional CWT
In order to fully demonstrate the strength and usefulness of our scheme, we compute the 1D power spectrum for the same data set of the dark matter density field at z = 0, with the traditional CWT and our CWT, respectively. Due to its similarity to the Gaussian-derived wavelet, we use the Mexican hat wavelet (shown in equation (A.5)) as the basis function of the traditional wavelet transform for the computation. Here, the wavelet power spectrum averaged for 100 000 lines is defined as$\begin{eqnarray}{S}_{{\rm{d}}}^{t}(a)=\displaystyle \frac{1}{{N}_{l}}\sum _{i=1}^{{N}_{l}}{\left\{\displaystyle \frac{1}{{L}_{b}}\int {{WT}}_{\delta }{\left(a,b\right)}^{2}{\rm{d}}b\right\}}_{i},\end{eqnarray}$where Nl = 100 000, and the superscript ‘t’ denotes ‘traditional’ to distinguish from our CWT power spectrum. Letting w = 1/a, we put ${S}_{{\rm{d}}}^{t}(a)$ and Sd(w) together in figure C1, and find that the magnitude of ${S}_{{\rm{d}}}^{t}(w)$ is greater than that of Sd(w) in the entire scale range.
Figure C1.
New window|Download| PPT slide Figure C1.The relationship between our wavelet power spectrum and the traditional CWT power spectrum. Upper panel: the Gaussian-derived wavelet power spectrum Sd(w) (red line), and the traditional wavelet power spectrum ${S}_{{\rm{d}}}^{t}(a=1/w)$ (blue line). Lower panel: the ratio of ${S}_{{\rm{d}}}^{t}(w)$ to LbwSd(w).
We reveal the relationship between these two wavelet power spectra in the following. Combining equations (A.1), (A.5) and (C.1), and with w = 1/a, we obtain$\begin{eqnarray}{S}_{{\rm{d}}}^{t}(w)=\displaystyle \frac{w}{3\sqrt{2\pi }{L}_{b}}\displaystyle \frac{1}{{N}_{l}}\sum _{i=1}^{{N}_{l}}{\left\{\int I(w,b){\rm{d}}b\right\}}_{i},\end{eqnarray}$in which I(w, b) is defined as$\begin{eqnarray}I(w,b)\equiv \iint {\rm{d}}x{\rm{d}}x^{\prime} \delta (x)\delta (x^{\prime} ){{\rm{\Psi }}}_{I}[w(x-b)]{{\rm{\Psi }}}_{I}[w(x^{\prime} -b)],\end{eqnarray}$where ${{\rm{\Psi }}}_{I}(x)=(2-{x}^{2}){{\rm{e}}}^{-{x}^{2}/4}$. Substituting equation (11) into (21), we get$\begin{eqnarray}{S}_{{\rm{d}}}(w)=\displaystyle \frac{1}{16\pi {L}_{b}^{2}}\displaystyle \frac{1}{{N}_{l}}\sum _{i=1}^{{N}_{l}}{\left\{\int I(w,b){\rm{d}}b\right\}}_{i}.\end{eqnarray}$From equations (C.2) and (C.4), one can see that ${S}_{{\rm{d}}}^{t}(w)$ and Sd(w) satisfy$\begin{eqnarray}\displaystyle \frac{{S}_{{\rm{d}}}^{t}(w)}{{L}_{b}{{wS}}_{{\rm{d}}}(w)}=\displaystyle \frac{16\pi }{3\sqrt{2\pi }}\approx 6.6843,\end{eqnarray}$which is reproduced by our numerical result shown in the lower panel of figure C1. Equation (C.5) indicates that our CWT and the traditional CWT are actually equivalent to each other in applicability.
Acknowledgments
We acknowledge the support by the National Science Foundation of China (No. 11947415, 12047569), and by the Natural Science Foundation of Jilin Province, China (No. 20180101228JC). In this work, we used the data from IllustrisTNG simulations. The IllustrisTNG simulations were undertaken with compute time awarded by the Gauss Centre for Supercomputing (GCS) under GCS Large-Scale Projects GCS-ILLU and GCS-DWAR on the GCS share of the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart (HLRS), as well as on the machines of the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany.
,1,21College of Physics, 
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
