Noise reduction in astronomical spectra using wavelet packets

SUPPLEMENT SERIES Astron. Astrophys. Suppl. Ser. 124, 579-587 (1997) Noise reduction in astronomical spectra using wavelet packets M. Fligge and S.K. Solanki Institut für Astronomie, ETH-Zentrum, CH-8092 Zürich, Switzerland Received May 15; accepted December 6, 1996 Abstract. The wavelet representation of a signal offers greater flexibility in de-noising astronomical spectra than classical Fourier smoothing due to the additional wavelength resolution of the decomposed signal. We present here a new wavelet-based approach to noise reduction. It is similar to an application of the splitting algorithm of a wavelet packets analysis using non-orthogonal wavelets. It clearly separates the signal from the noise, in particular also at the noise-dominated highest frequencies. This allows a better suppression of the noise, so that the spectrum de-noised in this manner provides a closer approximation of the uncorrupted signal than in the case of a single wavelet transformation or a Fourier transform. We test this method on intensity and circularly polarized spectra of the sun and compare with Fourier and other wavelet-based de-noising algorithms. Our technique is found to give better results than any other tested denoising algorithm. It is shown to be particularly successful in recovering weak signals that are practically drowned by the noise. Key words: methods: data analysis — methods: numerical 1. Introduction Astronomical observations are often photon starved. Consequently many astronomical spectra have a poor signal-to-noise ratio (SNR) and are significantly corrupted by (white Gaussian) photon noise. This is even true for solar observations, since high resolution measurements of polarized light soon run out of photons (e.g., Stix 1991). The reduction of this noise is highly desirable for a number of reasons (cf. the papers in the volume edited by Cayrel De Strobel & Spite 1988, for examples of the merits of high SNR spectra). Send offprint requests to: M. Fligge Fourier smoothing has long been the method of choice to suppress noise (Brault & White 1971), but recently methods based on the wavelet transformation have become increasingly popular (Starck & Bijaoui 1994; Starck & Murtagh 1994; Murtagh et al. 1995). In principle they offer much greater flexibility for analyzing and processing data. The main advantage of wavelets lies in the additional “spatial” resolution of the transformed signal. In contrast to the Fourier transformation, the signal is decomposed into waves of finite length, i.e. into waves which are spatially localized – hence the name wavelets. The wavelet transform of a one-dimensional signal has two independent variables – a frequency and a spatial location variable. It leads to a decomposition of, say, a spectrum into a series of spectra at finer and coarser resolutions. Indeed, there is a close mathematical relationship between the wavelet transformation and the multi-resolution analysis of a signal (Mallat 1989). Hence the wavelet transform furnishes us with the frequency spectrum of a signal at every spatial location. This feature, besides others, opens new and fruitful ways of processing and analyzing data of various kinds. In particular, it allows the smoothing of astronomical spectra, typically composed of a continuum with interspersed spectral lines, to be optimized. Since highfrequency signals are present only at the wavelengths of the spectral lines, only at these wavelengths need they be kept in the de-noised spectrum. At the remaining positions, i.e. in the continuum, only the lowest frequencies are due to the source itself. Whereas Fourier filtering affects all data points in the same manner, wavelets allow different parts of spectra to be filtered individually, in principle promising a considerably refined and improved treatment. In the present paper, we compare different wavelet smoothing methods to each other and to Fourier smoothing for the specific case of astronomical spectra. We also present a wavelet-packets based smoothing scheme which we find to be superior in recovering the true signal from a combination of signal and noise, at least for the cases we have considered. 2. Some relevant properties of wavelets 3. Methods The Fourier transformation decomposes a signal into sines and cosines of different frequencies. The wavelet transformation acts similarly, but instead of non-local, strictly periodic sines and cosines, it uses a set of spatially localized functions ψa,b (x) called wavelets (Daubechies 1988; Meyer 1993; see Press et al. 1992 for a simple introduction to the subject). The wavelets are constructed by translating and dilating a mother wavelet ψ(x)   1 x−b ψa,b (x) = p ψ (a 6= 0), (1) a |a| 3.1. Filtering and de-noising methods based on wavelets where the scale parameter a plays the role of a frequency and b is the position parameter. By increasing a, the wavelet ψa,b (x) is broadened, while changing b translates it along the x-axis. The set of parameters (a, b) describes a point in the so called scale-space plane. The continuous wavelet transform of a function f (x), Wc f (a, b), is defined by Z ∞ ∗ Wc f (a, b) = hψa,b (x), f (x)i = f (x)ψa,b (x)dx. (2) −∞ It is invertible (Grossmann & Morlet 1984) and the function f (x) can be recovered by evaluating the double integral: Z Z dadb 1 ∞ ∞ f (t) = Wc f (a, b)ψa,b (x) 2 . (3) Λ −∞ −∞ a Note that, unlike the Fourier transform, the wavelet transform is not its own inverse. This implies that a signal may be transformed several times using wavelets and be further decomposed at each transformation. Such a sheme of repeated application of the wavelet transform leads to the splitting algorithm of a wavelet packets analysis (Wickerhauser 1991, 1994; Chui 1992b) which lies at the heart of the technique we propose. For practical applications, the continuous set of parameters (a, b), must be discretized (Daubechies 1988; Mallat 1989). The parameterization of the discrete (a, b) pairs is of crucial significance to the discrete wavelet transform and especially to the stability of the reconstruction algorithm (Daubechies 1990). For most of the parameterizations of (a, b) the set of {ψa,b (x)} is highly redundant, i.e. each subset of them can be generated by linear combinations of the others (Daubechies et al. 1986). Although in such cases the reconstruction is not exact anymore, such a decomposition has a remarkable advantage when considering de-noising applications (Daubechies 1990). The particular class of non-orthogonal 1-D wavelets and the corresponding discrete wavelet transform together with the numerical algorithm we have used was proposed by Mallat & Zhong (1992, see their Appendix A). The multi-level decomposition described in Sect. 3.2 is a direct application of wavelet-packets using this particular kind of base functions. Most of the wavelet coefficients of the transform of a noiseless signal are close to zero. Therefore the most obvious way of filtering in the wave (...truncated)