Fourier transform in combination with deconvolution tech- ... odic signals, the behavior of these algorithms highly depends on intrinsic ... of both t...

0 downloads 0 Views 467KB Size

No documents

Printed 24 October 2018

(MN LATEX style file v2.2)

Detection of superimposed periodic signals using wavelets X. Otazu1,2⋆, M. Rib´o1, M. Peracaula1, J. M. Paredes1, J. N´ un ˜ez1,3 1

Departament d’Astronomia i Meteorologia, Universitat de Barcelona, Av. Diagonal 647, E-08028, Barcelona, Spain de Visi´ o per Computador, Edifici O, Campus Universitat Aut` onoma de Barcelona (UAB), 08193 Cerdanyola (Barcelona), Spain 3 Observatori Fabra, Cam´ ı del Observatori sn, 08035 Barcelona, Spain

arXiv:astro-ph/0202107v1 5 Feb 2002

2 Centre

Accepted YYYY MMMMMMMM DD. Received YYYY MMMMMMMM DD

ABSTRACT

In this paper we present a wavelet based algorithm that is able to detect superimposed periodic signals in data with low signal-noise ratio. In this context, the results given by classical period determination methods highly depend on the intrinsic characteristics of each periodic signal, like amplitude or profile. It is then difficult to detect the different periods present in the data set. The results given by the wavelet based method for period determination we present here are independent of the characteristics of the signals. Key words: methods: data analysis – stars: variables: other – X-rays: general.

1

INTRODUCTION

The search for periodic signals is common in many areas of astronomy, and, often, simple Fourier based methods or epoch folding methods are able to detect the periods that are present in the data sets. However, this turns out to be difficult when the periodic signals are not sinusoidal or when the signal-noise ratio is very low. Specifically, in high energy astrophysics two problems are responsible for reduced signal-noise ratios in the data: the relatively small number of high energy photons when compared to the number of photons available at other wavelengths, and the low efficiency in photon counting of high energy detectors. Hence, the astronomical data supplied by Xray telescopes, especially flux monitorings, often suffer from poor statistics. This poses many problems when processing and analysing the data in order to determine the flux, detect periodic signals, etc. In addition, some X-ray sources present several periods, which arise from different physical processes such as pulses, eclipses in two-body systems or occultation by a precessing accretion disc. Some examples of this kind of sources are SMC X-1 (Wojdowski et al. 1998), LMC X-4 (Levine et al. 1991), and LS I +65◦ 010 (Corbet, Finley & Peele 1999). Classical period determination methods are based either on epoch folding techniques or on Fourier decomposition analysis. The first group of methods are based on the analysis of the dispersion of the diferent light curves produced by folding the data over a range of trial periods (see, for example, Lafler & Kinman 1965; Jurkevich 1971 and Stellingwerf 1978). The second type of methods use the Fourier transform in combination with deconvolution tech⋆

[email protected]

niques to deal with the data sampling function (see, for example, Lomb 1976; Scargle 1982, 1989; Roberts, Leh´ ar & Dreher 1987 and Press & Rybicky 1989). Epoch folding methods and Fourier based methods work well when applied to data sets that present a unique period. However, when the analysed data set contains several periodic signals, the behavior of these algorithms highly depends on intrinsic signal characteristics. Every individual periodic signal may present different spectral behavior (amplitude and profile) and depending on them, the algorithm will detect some kinds of signal better than others. Fourier based techniques will, in general, successfully detect two combined sinusoidal signals, even if one of them has a low signal-noise ratio. However, they have difficulties in detecting the nonsinusoidal signals that might be present in a data set (as in the astronomical case of pulsed emission superimposed on an orbital period). On the contrary, epoch folding methods are well suited for detecting non-sinusoidal signals but they tend to fail when two or more periods are present in the data set, especially when the signal-noise ratio is low. In general, we can say that the greater the difference between the spectral characteristics of each signal, the more difficult it is to detect each period with classical methods, and the less statistically significant are the results given by them. Therefore, a new algorithm becomes immediately necessary for the detection of superimposed periodic signals. In this paper, we show how the methodology of wavelet theory is very well suited to this problem, since it is completely oriented to decompose functions into their several spectral characteristics. This allows us to isolate every signal present in our data and to analyse them separately, avoiding their mutual influences. In Section 2 we outline some concepts in wavelet theory relevant to the stated problem. In Section 3 we propose an algorithm to detect each of the

2

X. Otazu et al.

periodic signals present in a data set by combining wavelet decomposition with classical period determination methods. In Sections 4 and 5 we present some examples of synthetic data we used to test the algorithm and the results we obtained. We summarise our conclusions in Section 6.

f (t) in the set of basis functions ψm,l (t). These basis functions are scaled and translated versions of a general function ψ(t) called Mother Wavelet. To construct the basis functions ψm,l (t), the Mother Wavelet is dilated and translated according to parameters m and l as follows: ψm,l (t) = 2m/2 ψ(2m t − l) .

2

OUTLINE OF THE WAVELET TRANSFORM

Multiresolution analysis based on the wavelet theory introduces the concept of details between successive levels of scale or resolutions (Chui 1992; Daubechies 1992; Meyer 1993; Young 1993; Kaiser 1994; Vetterli & Kovacevic 1995). Wavelet decomposition is increasingly being used for astronomical signal processing (Starck & Murtagh 1994; Starck, Murtagh & Bijaoui 1998) and remotely sensed data (Yocky 1995; Datcu, Luca & Seidel 1996). The method is based on the decomposition of the signal into multiple channels based on their local frequency content. The wavelet transform is an intermediate representation between the Fourier and the temporal one. In the Fourier transform we know the global frequency content of our signal, but we have no information about the temporal localization of these frequencies. However, the wavelet transform gives us an idea of both the local frequency content and the temporal distribution of these frequencies. Since in Fourier space the basis functions are sinusoidal, they extend along all space and do not have temporal concentration. However, wavelets are concentrated around a central point, thus they have a high degree of temporal localization.

Therefore, all the basis functions ψm,l (t) have the same profile, that is, the Mother Wavelet profile. Using (4) we obtain an orthonormal wavelet basis (Daubechies 1992). The inverse discrete wavelet transform is given by the reconstruction formula: f (t) =

Some theory

Next, we will briefly present an outline of the wavelet transform, relying on the multiresolution signal representation concept. Given a signal f (t), we construct a sequence Fm (f (t)) of approximations of f (t). Each Fm (f (t)) is specific for the representation of the signal at a given scale (resolution). Fm (f (t)) represents the projection of the signal f (t) from the signal space S onto subspace Sm . In this representation, Fm (f (t)) is the ‘closest’ approximation of f (t) with resolution 2m . The differences between two consecutive scales m and m+1 are the multiresolution wavelet planes or ‘detail’ signal at resolution 2m : wm (f (t)) = Fm (f (t)) − Fm+1 (f (t)) .

(1)

wm (f (t)) =

Wm,l (f )ψm,l (t) ,

(2)

Wm,l (f )ψm,l (t) .

(5)

l

In summary, the wavelet transform describes at each resolution step the difference between the previous and the current resolution representation. By iterating the process from the highest to the lowest available resolution level we obtain a pyramidal representation of the signal. This usually includes a decimation process, i.e., in each iteration 1 out of 2 points is discarded, which implies that the number of data points at lower frequencies is highly reduced. In the Fourier transform, the noise contribution is spread along all frequencies. In contrast, one of the interesting properties of the wavelet transform (which is also frequency based) is that noise contribution is more important on the higher frequency wavelet planes. The ‘` a trous’ algorithm

In order to obtain a discrete wavelet decomposition for signals, we follow Starck & Murtagh (1994) and we use the discrete wavelet transform algorithm known as ‘` a trous’ (‘with holes’) to decompose the signal into wavelet planes. Given a signal p we construct the sequence of approximations: F1 (p) = p1 , F2 (p1 ) = p2 , F3 (p2 ) = p3 , · · · .

(6)

To construct the sequence, this algorithm performs successive convolutions with a filter obtained from an auxiliary function named scaling function. We use a scaling function which has a B3 cubic spline profile. The use of a B3 cubic spline leads to a convolution with a mask of 5 elements, all elements being scaled up by 16: (1,4,6,4,1). As stated above, the wavelet planes are computed as the differences between two consecutive approximations pi−1 and pi . Letting wi = pi−1 − pi (i = 1, · · · , n), in which p0 = p, we can write the reconstruction formula:

This detail signal can be expressed as:

X

XX m

2.2 2.1

(4)

p=

n X

wi + p r .

(7)

i=1

l

where coefficients Wm,l (f ) are given by the direct wavelet transform of the signal f(t): Wm,l (f ) =

Z

+∞

f (t) ψm,l (t) dt .

(3)

−∞

The coefficients Wm,l (f ) are called wavelet coefficients of f (t). Such coefficients correspond to the fluctuations of the signal f (t) near the point l at resolution level m. Thus, the wavelet transform (3) represents the expansion of signal

In this representation, the signals pi (i = 0, · · · , n) are versions of the original signal p at increasing scales (decreasing resolution levels), wi (i = 1, · · · , n) are the multiresolution wavelet planes and pr is a residual signal (in fact n = r, but we explicitly substitute n by r to clearly express the concept of residual). In our case, we are using a dyadic decomposition scheme. Thus, the original signal p0 has double resolution than p1 , the signal p1 double resolution compared to p2 , and so on. If the resolution of signal p0 is, for example, 10 ∆t (being ∆t the sampling of the data), the resolution of

Detection of superimposed periodic signals p1 would be 20 ∆t, the resolution of p2 would be 40 ∆t, etc. All these pi (i = 0, · · · , n) signals have the same number of data points, in contrast to some more usual wavelet decomposition algorithms which reduce the number of points by a factor 2 when going through increasing scales (process known as decimation step). Later we will see why we are interested in maintaining the number of points.

3

PROPOSED ALGORITHM

We propose to apply the wavelet decomposition using the ‘` a trous’ algorithm to solve our initial stated problem: to isolate each of the periodic signals contained in a set of data and study them separately. Hence, the algorithm we propose is as follows: (i) Choose a value for n and decompose the original signal p into its wavelet planes wi (i = 0, · · · , n). (ii) Detect periods in each of the n obtained wavelet plane wi . Any method can be used to detect the periods present in each wi . In our tests, we used Phase Dispersion Minimization (PDM) (Stellingwerf 1978) and CLEAN (Roberts et al. 1987) methods. PDM belongs to the group of epoch folding methods mentioned in the introduction. Therefore it folds the data points over a set of trial periods and analyses the dispersion of every light curve in the phase space. To do that, the data points are grouped in phase bins and the most probable period is then found by minimizing the sum of the bin variances. On the other hand, CLEAN belongs to the group of Fourier transform decomposition methods. It works in the frequency domain and it is based on the fact that the spectrum obtained from the Fourier transformation of the data set contains artefacts caused by the imcompleteness of the sampling (time gaps) and the finite time span of it. CLEAN uses a non-linear deconvolution algorithm to clean (hence its name) these artefacts from the original (dirty) spectrum. This is done by interatively substracting from the dirty spectrum the expected response of a signal composed by a unique harmonical function (clean component). From the set of clean components, a clean spectrum is obtained which is at the end corrected to preserve the noise level and the frequency resolution. PDM is well suited to identify signals with any profile, not just a sinusoidal one. For instance, if we have a burst-like periodic signal, it is more successfully detected by the PDM algorithm than by CLEAN. On the other hand, if we knew that we are trying to estimate the period of a sinusoidal signal, it would be better to use CLEAN. The algorithm we present here could be seen as an improvement over the usual period detection algorithms, because prior to using them we decompose our signal into its wi wavelet planes, and afterwards we apply these algorithms to each wi . In PDM and CLEAN methods there is only one data set to study and to estimate its period. But in the wavelet based method there are several data sets to study (several wavelet planes), and this improves the detection probabilities. We can search for a certain period in several wavelet planes, which helps to discard some of the usual marginal artefact

3

detections present in PDM and CLEAN and to detect the true periods present in our original data. We can use other wavelet based algorithms (Szatm´ ary, Vink´ o & G´ al 1994), which are based on approximations of a continuous wavelet transform and the study of wavelet space coefficients. However, these algorithms present a non direct inverse wavelet transform: the search for periodicities is based on the fit between the wavelet base function profile and the signal one, and therefore on the values of the wavelet transform coefficients, which highly depend on the wavelet base used. For example, if we took a signal constructed by a Gaussian with 1 day duration which repeats every 30 days, it would be very difficult to find a wavelet base that fits well to this sparse signal. On the contrary, if we take a suitable wavelet base and a decomposition scheme which allows us to perform easily its inverse transform, this base will fit the profile of every Gaussian. Hence, performing the inverse transform of every wavelet scale, we can work in the temporal representation to find directly the periodicities. This is the case of the wavelet decomposition algorithm we use in the method we present here. Another reason to use the ‘` a trous’ algorithm is that the sampling of the wavelet planes is the same as the original data. This allows us to work directly in the temporal space with the frequency content of the corresponding wavelet plane. In other wavelet decomposition algorithms, we are usually forced to work on the decimated wavelet space. Hereafter, and for notational convenience, the wavelet based PDM and CLEAN methods will be called WPDM and WCLEAN respectively.

4

SIMULATED DATA

In order to check the benefit of applying WPDM versus PDM, or WCLEAN versus CLEAN, we generated several sets of simulated data containing two superimposed periodic signals. Each data set is composed of a high amplitude primary sinusoidal function and a secondary low amplitude one. For the latter we used two kinds of functions: sinusoidal and Gaussian. The first are intended to simulate variable sources with pure sinusoidal intensity profile (like precession of accretion discs), and the second burst-like events (like pulses or eclipses). Finally, we added a white-Gaussian noise to this combination of signals. We increased the value of the noise standard deviation, σ, up to the value where detection of periodic signals became statistically insignificant in both the classical and the wavelet based methods. Here we list the characteristics of the signals: (i) Each signal is generated as an evenly spaced data set composed of 1000 points, separated by a unit between them. (ii) The high amplitude sinusoidal function has an amplitude equal to 1 and a period of 108.5 units. The period selected should satisfy two conditions: not be an exact divisor of the number of points in the data set, and allow the presence of more than 5 complete periods in the simulated signal to guarantee the possibility of detection. (iii) The amplitudes of the low amplitude periodic functions (sinusoidal or Gaussian) are 0.1 and 0.5. (iv) The periods used for the secondary functions are 13.13 and 23.11 units. Others could be used, satisfying the

4

X. Otazu et al. w1

SIMULATED DATA

w2

w3

1.2

1.2

1.2

1.0

1.0

1.0

2.00

0.8

0.8

0.8

1.00

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.00 −1.00

0

200

400 600 Time units (∆t)

PDM Normalized Power

Dispersion

1.01 1.00 0.99 0.98 0

10 20 30 40 Trial Period (∆t)

1000

CLEAN

1.02

0.97

800

Dispersion

Flux (count rate)

3.00

50

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

50 100 Trial Period (∆t)

150

5

RESULTS

The simultaneous use of two independent methods, such as CLEAN and PDM, is usually applied to discriminate false period detections from the true ones. A similar procedure can be used with each one of the wavelet planes in the WPDM and WCLEAN methods. Therefore, when comparing the behavior of these different methods, we have to compare the usual PDM-CLEAN method combination for period estimation prior to the new WPDM-WCLEAN combination. The primary period (108.5 units) is always detected by all methods, and does not appear in Table 1 and Table 2. In these tables, the four last columns show the detected low amplitude periods for each data set using PDM, CLEAN, WPDM and WCLEAN methods, respectively. A dash is shown when a period is not detected, and a question mark when the detection is difficult or doubtful. When a period is indicated in the wavelet based methods, we also show in parentheses the wavelet planes where it is detected. 5.1

Sine+sine case

To illustrate the performances of WPDM and WCLEAN methods, we show in Fig. 1 a simulated data set which contains two pure sinusoidal functions without noise, and both their PDM and CLEAN outputs (in PDM we look for minima with

100

150

0

50

w4

100

150

0.2 1.2

1.0

1.0

1.0

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0.0

0.0

0.0

50

50

100

150

100

150

100

150

w6

1.2

0

0

w5

1.2

0 50 100 150 Trial Period (∆t)

w1

Normalized Power

In the first three columns of Table 1 we present the parameters used to generate each simulated data set in the case of sine+sine periodic signals. In Table 2 the first four columns show the parameters used to generate the sine+Gaussian signals.

50

0

50

Figure 2. WPDM output for each wavelet plane of the data at the top of Fig. 1.

Figure 1. Top: simulated data of a sinusoid with 23.11 units period, amplitude=0.1 and σnoise = 0, superimposed on the sinusoid with 108.5 units period and amplitude=1.0. Bottom: PDM and CLEAN outputs of this data set.

condition that they are not harmonics of the high amplitude sinusoidal function period. (v) In the case of the Gaussian signal, we have used two values for the Full Width Half Maximum (FWHM), corresponding to 2 and 6 units, respectively.

0

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

50

w2

100 w4

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

50

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2 150 0

50

w3

100

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2 150 0

50

w5

100

1.2 1.2 1.0 1.0 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0.0 0.0 −0.2 −0.2 150 0 50 100 150 0 Trial Period (∆t)

100

150

100

150

w6

50

Figure 3. WCLEAN output for each wavelet plane of the data at the top of Fig. 1.

subharmonics, while in CLEAN we look for maxima). The 108.5 units period is detected by both methods. In the case of PDM we have zoomed into the range of short trial periods and dispersion to try to detect the 23.11 units minimum. It is clear that even in the case of signals without noise, PDM method is unable to find the secondary signal (period=23.11 units). In contrast, CLEAN is able to find it, although with a low normalized power (less than 0.1, which is the ratio between the amplitudes of both periodic signals). The output of WPDM is shown in Fig. 2 and the one of WCLEAN in Fig. 3. As one would expect, the secondary period is only seen in the lower wavelet planes (w1 to w4 ), while the primary longterm period of 108.5 units is better detected when working with higher wavelet planes (w4 to w6 ). We have worked until w6 , because in this plane the 108.5 units period is clearly detected, and there is no reason to continue with higher values of n. This will be the case in all the analysed data we present here. Some remarks can be made after inspection of Table 1: (i) There is a better performance of CLEAN relative to because the latter has problems when dealing with superimposed low amplitude periodic signals. Only when the PDM,

Detection of superimposed periodic signals

5

Table 1. Periods detected in the sine+sine data sets. The first three columns are the parameters used for the simulated data. A dash is shown when a period is not detected, and a question mark when the detection is difficult or doubtful. When a period is indicated in the wavelet based methods, we also show in parentheses the wavelet planes where it is detected. Amplitude

σnoise

13.13

0.1

0.0 0.1 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 0.0 0.5 1.0 2.0 3.0 4.0 5.0 6.0 8.0 0.0 0.1 0.2 0.4 0.6 0.8 1.0 1.2 0.0 0.5 1.0 2.0 3.0 4.0 5.0

0.5

23.11

0.1

0.5

PDM

CLEAN

WPDM

WCLEAN

− − − − − − − − − − 13.13 13.12 13.12 13.09 13.08? − − − − − − − − − − − − 23.11 23.11 23.10 23.11 − − −

13.13 13.12 13.11 13.09 13.09? − − − − − 13.13 13.12 13.12 13.09 13.08? − − − − 23.11 23.12 23.12 23.12? − − − − 23.11 23.11 23.12 23.11 23.11? − −

13.13 13.13 13.11 13.09 13.09 13.08 13.08 13.07 13.07 − 13.13 13.13 13.13 13.13 13.13 13.13 13.13 13.13 − 23.11 23.11 23.13 23.11 23.17 23.13 23.13 − 23.11 23.11 23.12 23.11 23.17 23.13 −

13.13 13.13 13.11 13.09 13.09 13.08 13.08 13.07 13.07 − 13.13 13.13 13.13 13.13 13.13 13.13 13.13 13.13 − 23.11 23.11 23.13 23.11 23.17 23.13 23.11 − 23.11 23.11 23.13 23.11 23.17 23.13 −

amplitude of the secondary signal is close enough to the primary one (amplitude=0.5) is PDM able to detect it. (ii) In all methods, with high noise-signal ratios the detected periods are slightly different from the simulated ones. (iii) WPDM and WCLEAN perform always better than PDM and CLEAN methods. When CLEAN marginally detects the secondary period, WPDM and WCLEAN have no problems to detect it, and they work properly even with higher noise. In the WPDM case, the results are always much better than with PDM. (iv) As the noise increases, the detection starts to fail in the lower wavelet planes (higher frequencies), and only the higher ones (lower frequencies) are noise-free enough to allow period detection. (v) The maximum ratio between the noise and the signal amplitude to detect a period with WPDM and WCLEAN, is nearly constant for a given period, being around 12 for the 13.13 units period and around 8 for the 23.11 units period. The reason for this difference is that there are 76 complete periods in the case of 13.13 units period, and only 43 complete periods in the case of 23.11 units period.

(1,2,3) (2,3) (2,3) (2,3) (3) (3) (3) (3?) (3?) (1,2,3) (2,3) (2,3) (2,3) (2?,3) (3) (3) (3?) (1,2,3,4) (2,3,4) (3,4) (3,4) (3,4) (3?,4) (4?) (1,2,3,4) (2,3,4) (2?,3,4) (3,4) (3?,4) (4?)

(1,2,3) (2,3) (2,3) (2,3) (2,3) (2,3) (2,3) (2?,3?) (2?,3?) (1,2,3) (2,3) (2,3) (2,3) (2,3) (2,3) (2?,3) (3?) (1,2,3,4) (2,3,4) (3,4) (3,4) (3,4) (3?,4) (4?) (1,2,3,4) (1,2,3,4) (2?,3,4) (3,4) (3?,4) (4?)

w2

w1

w3

0.06

0.04

0.05

0.04

0.03

0.03

0.02

0.02

0.01

0.01 Flux (count rate)

Period

0.00

0.00

−0.02 −0.04 100

−0.01 −0.03

−0.01 150 w4

−0.02 200 100

150 w5

−0.05 200 100

200

300

400

500

700

w6

0.2

0.4

1.0

0.1

0.2

0.6

0.0

0.0

0.2 −0.2 −0.1 −0.2 100

−0.2 300

500

−0.4 700 100

−0.6 300

500

700

−1.0 100

300

Time units (∆t)

Figure 4. Wavelet decomposition of the simulated data which contains a Gaussian with 13.13 units period, a FWHM of 2.0 units, amplitude=0.1 and σnoise = 0, superimposed on the sinusoid with 108.5 units period and amplitude=1.0.

6

X. Otazu et al. Table 2. Periods detected in the sine+Gaussian data sets. The first four columns are the parameters used for the simulated data. A dash is shown when a period is not detected, and a question mark when the detection is difficult or doubtful. When a period is indicated in the wavelet based methods, we also show in parentheses the wavelet planes where it is detected. FWHM

Period

Amplitude

σnoise

2.0

13.13

0.1

0.0 0.1 0.15 0.0 0.5 0.75 1.0 0.0 0.1 0.15 0.0 0.5 0.75 0.0 0.1 0.2 0.3 0.5 0.0 0.5 1.0 1.5 2.0 0.0 0.1 0.15 0.2 0.0 0.5 1.0 1.5 2.0

0.5

23.11

0.1

0.5

6.0

13.13

0.1

0.5

23.11

0.1

0.5

5.2

PDM

CLEAN

WPDM

WCLEAN

− − − − − − − − − − − − − − − − − − 13.13 13.13 13.12? − − − − − − 23.11 23.11 23.12? − −

13.13 − − 13.13 13.14 − − 23.11 − − 23.11 − − 13.13 13.13 13.13? − − 13.13 13.13 13.12? − − 23.11 23.12 23.12? − 23.11 23.11 23.12? − −

13.13 13.13 − 13.13 13.14 13.15 − 23.11 23.11 − 23.11 23.12 − 23.12 23.13 23.14 23.14 − 13.13 13.13 13.13 13.13 − 23.11 23.11 23.12 − 23.11 23.11 23.12 23.12 −

13.13 13.11 − 13.13 13.14 − − 23.11 23.11 − 23.11 − − 13.13 13.13 13.13 13.13 − 13.13 13.13 13.13 − − 13.13 23.11 23.11 − 13.13 23.11 23.11 − −

Sine+Gaussian case

For illustrative purposes, we show in Fig. 4 the wavelet decomposition of a sine+Gaussian signal. In the lower wavelet planes (planes w1 and w2 ) we have isolated the Gaussian contribution. In the w3 we have a mixed contribution from both periods, but in the higher planes we only have the longest period contribution. We must note that the use of two different FWHM for the Gaussian, combined with two different periods (13.13 and 23.11 units), gives 4 different profiles. Hence, the phase duration of the burst-like event ranges from very low to relatively high values in the following order: FWHM=2 and period=23.11, FWHM=2 and period=13.13, FWHM=6 and period=23.11, and finally FWHM=6 and period=13.13. In view of the results displayed in Table 2 we can make the following comments: (i) Again, there is a better performance of CLEAN over for the reason explained above. However, we must note that when the FWHM is only 2 units, PDM never detects the secondary period. Only with FWHM=6 units and an amplitude of 0.5, can PDM detect the low amplitude periodic signals. PDM,

(1,2,3) (2?,3?) (1,2,3) (2,3) (2?,3?) (2?,3) (2?,3?) (2?,3,4) (2?,3?) (1,2,3) (2,3) (2?,3) (3) (1,2,3,4?) (2,3) (2?,3) (3) (1,2,3,4?) (2,3,4?) (2?,3,4?) (1,2,3,4) (1,2,3) (3,4) (3?)

(1,2,3) (3?) (1,2,3) (2?,3)

(3?,4) (4?) (2?,3)

(2,3) (2,3) (2?,3) (2?,3) (1,2,3,4?) (2,3) (2?,3)

(1,2,3,4?) (2,3,4?) (2?,3,4?) (1,2,3,4) (2,3,4) (3,4?)

(ii) In all methods, with high noise-signal ratios the detected periods are slightly different from the simulated ones. (iii) WPDM and WCLEAN perform always better than PDM and CLEAN methods. When CLEAN marginally detects the secondary period, WPDM and WCLEAN have no problems to detect it, and they work properly even with higher noise. In the WPDM case, the results are always much better than with PDM. (iv) As the noise increases, the detection starts to fail in the lower wavelet planes (higher frequencies), and only the higher ones (lower frequencies) are noise-free enough to allow period detection. (v) In all cases with amplitude=0.1, the WPDM performance is very similar to WCLEAN. In the amplitude=0.5 cases, WPDM is always better than WCLEAN because the signal is not sinusoidal and the amplitude is high enough to allow a good detection. (vi) For a given amplitude of the Gaussian signal, the maximum noise-signal ratio achieved with WPDM and WCLEAN increases with the phase duration of the FWHM. Finally, and for illustrative purposes, we show in Fig. 5 and Fig. 8 two simulated data sets generated with the following common parameters: 13.13 units period, FWHM=6.0 units and amplitude=0.5. The only difference is that in

Detection of superimposed periodic signals SIMULATED DATA

4.0

7.0

3.0

5.0

Flux (count rate)

2.0 1.0 0.0 −1.0 −2.0

0

200

400 600 Time units (∆t)

800

PDM

0.98 0.96 10

20 30 40 Trial Period (∆t)

50

w1

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

50 100 Trial Period (∆t)

0.99 0.98 0.97

150

0

10

20 30 40 Trial Period (∆t)

1.2 1.0 0.8 0.6 0

50

100

150

w1

w3

1.05

1.00

1.0

1.000

1.00

1.00

0.99 0.995

0.7

0.85

0.6

0

10 20 30 40 50 w5

0.95

1.2

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.0

0

50

100

150

0.98

0

50

100

150

w6

1.2

0.0

0.990

0.97

0.985

0.96

0

10 20 30 40 50 w4

0

50

100

150

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

50

100

150

10 20 30 40 50 w5

1.2

1.2

1.0

1.0

100

0.85

150

50

100

150

100

150

0.8 0.6

0.90 50

0

w6

1.00

0.6

0

0.85

1.05

0.4

0.4 0

50

100

150

0.2

0.2 0

50

100

150

0.0

0

50

Trial Period (∆t)

Figure 9. WPDM output for each wavelet plane of the data at the top of Fig. 8.

w3

10 20 30 40 50 w5

0

0.8

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

50

w1

100

150

w6 1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

50

100

150

Trial Period (∆t)

Figure 7. WCLEAN output for each wavelet plane of the data at the top of Fig. 5.

Normalized Power

Normalized Power

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

0.90

0.95

w2 1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

150

w3

w2 1.01

w1

10 20 30 40 50 w4

50 100 Trial Period (∆t)

1.005

Figure 6. WPDM output for each wavelet plane of the data at the top of Fig. 5.

0

0

1.1

Trial Period (∆t)

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

50

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

1.05

0.90

1000

Figure 8. Top: simulated data of a Gaussian with 13.13 units period, a FWHM of 6.0 units, amplitude=0.5 and σnoise = 1.5, superimposed on the sinusoid with 108.5 units period and amplitude=1.0. Bottom: PDM and CLEAN outputs of this data set.

Dispersion

Dispersion

10 20 30 40 50 w4

800

CLEAN

1.00

0.8

0

400 600 Time units (∆t)

PDM

0.95

0.4

200

0.9

0.995 0.990

0

1.01

w2

1.000

−5.0

1000

Figure 5. Top: simulated data of a Gaussian with 13.13 units period, a FWHM of 6.0 units, amplitude=0.5 and σnoise = 0.5, superimposed on the sinusoid with 108.5 units period and amplitude=1.0. Bottom: PDM and CLEAN outputs of this data set.

1.005

−3.0

Dispersion

Normalized Power

Dispersion

1.00

0

1.0 −1.0

CLEAN

1.02

0.94

3.0

Normalized Power

Flux (count rate)

SIMULATED DATA

7

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2 1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

0

w2

10 20 30 40 50 w4

50

100

150

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2 1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

0

w3

10 20 30 40 50 w5

50

100

150

1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

50

100

150

100

150

w6 1.2 1.0 0.8 0.6 0.4 0.2 0.0 −0.2

0

50

Trial Period (∆t)

Figure 10. WCLEAN output for each wavelet plane of the data at the top of Fig. 8.

8

X. Otazu et al.

Fig. 5, σnoise = 0.5, while in Fig. 8, σnoise = 1.5. The outputs of PDM, CLEAN, WPDM and WCLEAN are shown from Fig. 5 to Fig. 10. As can be seen in the figures, in the case σnoise = 0.5, the 13.13 units period is clearly detected by all methods. In WPDM and WCLEAN the period is detected in the wavelet planes w2 and w3 . These detections are much better than those of PDM and CLEAN. On the other hand, in the σnoise = 1.5 case, the period is not detected by PDM nor by CLEAN. WCLEAN also fails to detect it: the maximum in the wavelet plane w3 is not indicative of a period detection and, moreover, its position is in the 15.5 units trial period. The only clear detection with a subharmonic is in w3 of WPDM.

6

CONCLUSIONS

When dealing with two superimposed periodic signals, the performance of classical methods, like PDM and CLEAN, depends on the characteristics of these signals. These methods often fail to detect simultaneously both periods, even with high or moderate signal-noise ratios. With the wavelet based methods we present here, WPDM and WCLEAN, it is possible to reach lower values of signal-noise ratios and still detect both periodic signals. This means that the wavelet based methods are less affected by noise than PDM and CLEAN, and they allow us to detect periodic signals with noisier data. Another advantage of the proposed methods is the simultaneous detection of periods in several wavelet planes. If a period is marginally detected in more than one wavelet plane, we are probably detecting a true periodic signal. This allows us to improve the detection confidence. These facts prove the major ability of WPDM and WCLEAN over PDM and CLEAN to deal with noisy signals and to detect superimposed periods even with low signalnoise ratios.

ACKNOWLEDGMENTS We acknowledge partial support from DGI of the Ministerio de Ciencia y Tecnolog´ıa (Spain) under grant AYA20013092. This research has made use of facilities of CESCA and CEPBA, coordinated by C4 (Centre de Computaci´ o i Comunicacions de Catalunya). During this work, M. Rib´ o has been supported by a fellowship from CIRIT (Generalitat de Catalunya, ref. 1999 FI 00199).

REFERENCES Chui C.K., 1992, An Introduction to wavelets, Boston Ac. Press, Boston Corbet R.H.D., Finley J.P., Peele A.G., 1999, ApJ, 511, 876 Datcu M., Luca D., Seidel K., 1996, in Klemm R., Keydel W., Ender J., eds, EUSAR’96. VDE-Verlag, Berlin, Offenbach, p. 375 Daubechies I., 1992, Ten Lectures on Wavelets, SIAM Press, Philadelphia Jurkevich I., 1971, Ap&SS, 13, 154

Kaiser G., 1994, A friendly guide to wavelets, Birkh¨ auser, Berlin Lafler J., Kinman T.D., 1965, ApJS, 11, 216 Levine A., Rappaport S., Putney A., Corbet R., Nagase F., 1991, ApJ, 381, 101 Lomb N.R., 1976, Ap&SS, 39, 447 Meyer Y., 1993, Wavelets. Algorithm and Applications, SIAM Press, Philadelphia Press W.H., Rybicky G.B., 1989, ApJ, 338, 277 Roberts D.H., Leh´ ar J., Dreher J.W., 1987, AJ, 93, 968 Scargle J.D., 1982, ApJ, 263, 835 Scargle J.D., 1989, ApJ, 343, 874 Starck J.L., Murtagh F., 1994, A&A, 288, 342 Starck J.L., Murtagh F., Bijaoui A., 1998, Image and data analysis: the multiscale approach, CUP, Cambridge Stellingwerf R.F., 1978, ApJ, 224, 953 Szatm´ ary K., Vink´ o J., G´ al J., 1994, A&AS, 108, 377 Vetterli M., Kovacevic J., 1995, Wavelets and subband coding, Prentice Hall, Englewood Cliffs, New Jersey Wojdowski P., Clark G.W., Levine A.M., Woo J.W., Nan Zhang S., 1998, ApJ, 502, 253 Yocky D.A., 1995, Journal of the Optical Society of America, series A, 12, 1834 Young R.K., 1993, Wavelet theory and its applications, Kluwer Ac Pub, London This paper has been typeset from a TEX/ LATEX file prepared by the author.