2 Analysis

2.1 The wavelet decomposition

 Wavelet transforms have been widely investigated recently because of their suitability for a large number of signal and image processing tasks. Wavelet analysis involves a convolution of the signal with a convolving function or a kernel (the wavelet) of specific mathematical properties. When satisfied, these properties define an orthogonal basis, which conserves energy, and guarantee the existence of an inverse to the wavelet transform.

The principle behind the wavelet transform, as described by Grossman & Morlet (1984), Daubechies (1988) and Mallat (1989), is to hierarchically decompose an input signal into a series of successively lower resolution reference signals and their associated detail signals. At each decomposition level, L, the reference signal has a resolution reduced by a factor of 2^L with respect to the original signal. Together with its respective detail signal, each scale contains the information needed to reconstruct the reference signal at the next higher resolution level. Wavelet analysis can therefore be considered as a series of bandpass filters and be viewed as the decomposition of the signal into a set of independent, spatially oriented frequency channels. Using the orthogonality properties, a function in this decomposition can be completely characterised by the wavelet basis and the wavelet coefficients of the decomposition.

The multi-level wavelet transform (analysis stage) decomposes the signal into sets of different frequency localisations. It is performed by iterative application of a pair of Quadrature Mirror Filters (QMF). A scaling function and a wavelet function are associated with this analysis filter bank. The continuous scaling function $\phi_{\rm A}(x)$ satisfies the following two-scale equation:

$$\phi_{\rm A}(x)= \sqrt{2} \sum_n h_{0}(n)\phi_{\rm A}(2x-n),$$ (1)

where h₀ is the low-pass QMF. The continuous wavelet $\psi_{\rm A}(x)$ is defined in terms of the scaling function and the high-pass QMF h₁ through:

$$\psi_{\rm A}(x)=\sqrt{2} \sum_n h_{1}(n)\phi_{\rm A}(2x-n).$$ (2)

The same relations apply for the inverse transform (synthesis stage) but, generally, different scaling function and wavelet ($\phi_{\rm S}(x)$ and $\psi_{\rm S}(x)$) are associated with this stage:

$$\phi_{\rm S}(x)= \sqrt{2}\sum_n g_{0}(n)\phi_{\rm S}(2x-n),$$ (3)

$$\psi_{\rm S}(x)= \sqrt{2}\sum_n g_{1}(n)\phi_{\rm S}(2x-n).$$ (4)

Equations (1) and (3) converge to compactly supported basis functions when

$$\sum_n h_{0}(n)=\sum_n g_{0}(n)=\sqrt{2}.$$ (5)

The system is said to be bi-orthogonal if the following conditions are satisfied:

$$\begin{eqnarray} & & \int_{\mathbb{R}}\phi_{\rm A}(x)\phi_{\rm S}(x-k){\rm d}x=\... ... & & \int_{\mathbb{R}}\phi_{\rm S}(x)\psi_{\rm A}(x-k){\rm d}x=0 .\end{eqnarray}$$ (6) (7) (8)

Cohen et al. (1990) and Vetterli & Herley (1992) give a complete treatment of the relationship between the filter coefficients and the scaling functions.

  

Figure 1: In the upper panel, a 4 level "dyadic'' decomposition is used, in the lower, a 4 level "pyramidal'' decomposition is used. In both cases, the diagonal squares and the last decomposition level are identical and the grey levels indicate the amplitude of the wavelet coefficients

When applied to bi-dimensional data (typically images), three main types of decomposition can be considered: dyadic (or octave band), pyramidal and uniform.

1.

A "dyadic'' decomposition refers to a transform in which only the reference sub-band (low-pass part of the signal) is decomposed at each level. In this case, the analysis stage is applied in both directions of the image at each decomposition level. The total number of sub-bands after L levels of decomposition is then 3L+1 (Fig. 1, upper panel).

2.

A "pyramidal'' decomposition is similar to a "dyadic'' decomposition in the sense that only the reference sub-band is decomposed at each level, but it refers here to a transform that is performed separately in the two directions of the image. The total number of sub-bands after L levels of decomposition is then (L+1)² (Fig. 1, lower panel).

3.

A "uniform'' decomposition refers to one in which all sub-bands are transformed at each level. The total number of sub-bands after L levels of decomposition is then 4^L.

The wavelet functions are localised in space and, simultaneously, they are also localised in frequency. Therefore, this approach is an elegant and powerful tool for image analysis because the features of interest in an image are present at different characteristic scales. Moreover, if the input field is Gaussian distributed, the output is distributed the same way, regardless of the power spectrum. This arises from the linear transformation properties of Gaussian variables. The distribution of the wavelet coefficients of a Gaussian process is thus a Gaussian. Conversely, we expect that any non-Gaussian signal will exhibit a non-Gaussian distribution of its wavelet coefficients.

In our study, we have used bi-orthogonal wavelets, which are mainly used in data compression, because of their better performance than orthogonal wavelets, in compacting the energy into fewer significant coefficients. There exist bi-orthogonal wavelet bases of compact support that are symmetric or antisymmetric. Antisymmetric wavelets are proportional to, or almost proportional to, a first derivative operator (e.g. the 2/6 tap filter (filter #5) of Villasenor et al. (1995), or the famous Haar transform which is an orthogonal wavelet). Symmetric wavelets are proportional to, or almost proportional to, a second derivative operator (e.g. the 9/3 tap filter of Antonini et al. 1992). In the frame of detecting non-Gaussian signatures, the choice of the wavelet basis is critical because non-Gaussian features exhibit point sources or step edges. The wavelet must have a very good impulse response and a low shift variance, i.e. it better preserve the amplitude and the location of the details. Villasenor et al. (1995) have tested a set of bi-orthogonal filter banks, within this context, to determine the best ones for image compression. They conclude that even length filters have significantly less shift variance than odd length filters, and that their performance in term of impulse response is superior. In these filters, the high pass QMF is antisymmetric which is also a desirable property in the sense that we will also be interested in the statistical properties of the multi-scale gradients. Consequently, in our study, we have chosen the 6/10 tap filter (filter #3) of Villasenor et al. (1995) (Fig. 2) which represents the best compromise between all the criteria and energy compaction. Using this filter, we have chosen to perform a four level dyadic decomposition of our data. This particular wavelet and decomposition method have already been used for source detection by Aghanim et al. (1998).

  

Figure 2: Scaling function (top) and wavelet function (bottom) corresponding to the filter #3 of Villasenor et al. (1995). Note that the wavelet function is antisymmetric

With such a decomposition, we also benefit from correlations between the two directions at each level, which is not the case with the pyramidal decomposition that treats both directions as if they were independent. Another advantage of this transform is that, for each level of decomposition, or scale, we benefit from the maximum number of coefficients possible, which is crucial for the statistics.

2.2 The test maps

We use a test case which consists of sets of 100 maps of a Gaussian signal and 200 maps of a non-Gaussian signal, all having the same bell-shaped power spectrum and $512\times512$ pixels. One of the non-Gaussian sets (100 maps) consists of a distribution of disks with uniform amplitude (top-hat profiles), generating step edges. The disks have different sizes and amplitudes and are randomly distributed in the map. The signal is weakly skewed because the average skewness, which is the third moment of the distribution, computed over the 100 non-Gaussian realisations is $\mu_3=-0.10\pm0.05$ (one sigma error for one realisation). Whereas the excess of kurtosis, the fourth moment ($\mu_4$) of the distribution, is $k=\mu_4/\mu_2^2 -3=1.09\pm0.17$. The second set of 100 non-Gaussian maps consists of a distribution of Gaussian profiles with different sizes and amplitudes. The skewness and excess of kurtosis of the 100 statistical realisations of the non-Gaussian distribution of Gaussian profiles are respectively $-0.06\pm 0.04$ and $1.19\pm 0.13$. The same quantities computed over the 100 Gaussian maps are respectively $\mu_3=0.14\,10^{-2}\pm 2.15\,10^{-2}$ and $0.03\,10^{-2}\pm3.87\,10^{-2}$. These numbers should be equal to zero for a Gaussian distribution. They are not because there is a statistical dispersion over a finite set of realisations. The purpose of this paper being to develop suitable statistical tests for non-gaussianity, we will not study other effects such as noise or beam dilution. These effects will be considered in an application of the method to the CMB signal in a second paper (Aghanim & Forni 1999).