3 Simple tests

We tested MPA with two groups of artificially generated data. For both data groups we used a light curve model of two harmonics

 

$$A_0^c + A_1^c\cos(\frac{2 \pi t_i}{P}) + B_1^c\sin(\frac{2 \pi t_i}{P})$$ M^c(t_i) =

$$+ A_2^c\cos(\frac{4 \pi t_i}{P}) + B_2^c\sin(\frac{4 \pi t_i}{P})+\epsilon_i^c.$$ (25)

Specific parameters for different models were chosen as follows. Time points $t_i, i = 1, \dots, N$ (where the number of observations N = 211) were obtained from the yet unpublished photometric data set of the RS CVn binary ER Vul obtained at Mount Maidanak Observatory in 1990-1996 with the 60 cm telescope. The typical errors for this UBVR photometry are about $0\hbox{$.\!\!^{\rm m}$ }030$ for the U channel and $0\hbox{$.\!\!^{\rm m}$ }015$ for the others (for details see Shevchenko 1980). In this way, the distribution of time points for our simulations contains complexities which are characteristic for real observational sequences. For all the experiments described here, we used the fixed period P=0.6980993, which comes from our preliminary analysis of the ER Vul data. In the first group of the experiments, we chose amplitudes A_r^c and B_r^c arbitrarily (see Table 1 and Fig. 1) to get a rather complex waveform as the light curve.

In the second group of experiments we used the parameters of a two-harmonic fit to the real light curve of ER Vul, which were estimated previously (You 1999). The corresponding amplitudes and error levels are given in Table 2 and the light curves are shown in Fig. 2. The most important difference between the two groups of simulations is the relatively low amplitude of the first harmonic in the second group. Finally, to the computed trigonometric polynomials we added normally distributed uncorrelated noise N_i(0,1) which we computed using a random number generator:

$$\epsilon_i^c = a \cdot E^c \cdot N_i(0,1),$$ (26)

where a is a noise amplification factor which we choose to be in the range 0.01 - 10 and E^c are the error levels in different channels ( c = 1, 2, 3, 4).

To test the new method, both the multichannel PDM and the single channel PDM were applied for all channels of the two groups of simulations. In the first place, we wanted to compare the detection capabilities of the multichannel method to that of the single channel analysis. This is why we chose a relatively narrow period search range $P_{\rm min} = 0.65$ and $P_{\rm max} = 0.75$. In this way the real period buried in the noise will basically "compete'' only with the random fluctuations. A more comprehensive search where spurious periods are taken into account seriously is described afterwards.

For the smoothing criterion L(t_ij), we adopted $D_{\rm min} = 0.9 \cdot P_{\rm min}, D_{\rm max} = 100 \cdot P_{\rm max}$ and for the selection criterion g(t_ij,P) we set $\tau = 0.125$ according to the above adopted models. These are typical values and do not depend strongly on the data to be analysed.

 

 

Table 3: Results of the three stage weighted MPA for the first group of simulations in comparison with the breakdown noise levels of the single channel search (see Sect. 3 for explanations). The true period is 0.6980993

Set	Noise	$S/N^{\rm U}$	$S/N^{\rm B}$	$S/N^{\rm V}$	$S/N^{\rm R}$	$P_{\rm PDM}$	$P_{\rm LM}$	$P_{\rm NLM}$
1	0.01Ec	1.43	1.53	1.63	1.78	0.698	0.69810094	0.69809942 $\pm$ 0.00000014
2	0.1Ec	0.41	0.49	0.60	0.77	0.698	0.69810094	0.6981014 $\pm$ 0.0000014
3	0.25Ec	0.03	0.08	0.22	0.35	0.698	0.69809550	0.6980980 $\pm$ 0.0000036
4	0.5Ec	-0.25	-0.19	-0.05	0.06	0.698	0.69810094	0.6980991 $\pm$ 0.0000075
5	0.75Ec	-0.45	-0.34	-0.22	-0.10	0.698	0.69810094	0.698099 $\pm$ 0.000010
6	1.0Ec	-0.53	-0.49	-0.39	-0.23	0.698	0.69807374	0.698075 $\pm$ 0.000013
7	1.25Ec	-0.68	-0.57	-0.47	-0.36	0.698	0.69809006	0.698088 $\pm$ 0.000019
8	1.5Ec	-0.75	-0.65	-0.55	-0.40	0.697	0.69811130	0.698109 $\pm$ 0.000020
9	1.75Ec	-0.79	-0.71	-0.57	-0.48	0.697	0.69810586	0.698108 $\pm$ 0.000032
10	2.0Ec	-0.84	-0.82	-0.65	-0.53	0.698	0.69816081	0.698162 $\pm$ 0.000029
11	2.25Ec	-0.88	-0.83	-0.72	-0.58	0.697	0.69816570	0.698166 $\pm$ 0.000037
12	2.5Ec	-0.94	-0.88	-0.76	-0.66	--	--	--

 

 

Table 4: Results of the three stage weighted MPA for the second group of simulations (see Sect. 3 for explanations)

Set	Noise	$S/N^{\rm U}$	$S/N^{\rm B}$	$S/N^{\rm V}$	$S/N^{\rm R}$	$P_{\rm PDM}$	$P_{\rm LM}$	$P_{\rm NLM}$
1	0.01Ec	1.75	1.98	1.92	1.91	0.698	0.69809688	0.698099296 $\pm$ 0.000000023
2	0.1Ec	0.96	1.33	1.29	1.27	0.698	0.69809688	0.69809935 $\pm$ 0.00000023
3	0.25Ec	0.61	0.97	0.91	0.90	0.698	0.69810232	0.69810025 $\pm$ 0.00000057
4	0.5Ec	0.32	0.66	0.60	0.61	0.698	0.69809688	0.6980982 $\pm$ 0.0000011
5	1.0Ec	-0.01	0.32	0.28	0.31	0.698	0.69809688	0.6980982 $\pm$ 0.0000022
6	2.0Ec	-0.31	0.02	-0.02	0.02	0.698	0.69809688	0.6980978 $\pm$ 0.0000042
7	3.0Ec	-0.45	-0.14	-0.12	-0.19	0.698	0.69810232	0.6981013 $\pm$ 0.0000068
8	4.0Ec	-0.66	-0.24	-0.27	-0.31	0.698	0.69810776	0.698108 $\pm$ 0.000010
9	5.0Ec	-0.71	-0.35	-0.39	-0.39	0.698	0.69811321	0.698111 $\pm$ 0.000013
10	6.0Ec	-0.75	-0.44	-0.44	-0.50	0.698	0.69808464	0.698083 $\pm$ 0.000013
11	7.0Ec	-0.83	-0.53	-0.53	-0.53	0.698	0.69810097	0.698100 $\pm$ 0.000019
12	8.0Ec	-0.92	-0.59	-0.60	-0.55	0.698	0.69811865	0.698116 $\pm$ 0.000018
13	9.0Ec	-0.95	-0.64	-0.63	-0.67	0.698	0.69810232	0.698104 $\pm$ 0.000018
14	10.0Ec	-0.97	-0.67	-0.70	-0.72	--	--	--

Tables 3 and 4 show the recovered periods $P_{\rm PDM}$ from the multichannel PDM for the first and the second groups of simulations, respectively. As a noise level indicator, we define a logarithmic signal to noise ratio as

$${S/N}^c = \log\left({\sigma_m^c \over \sigma_n^c}\right), \ \ \ c=1,2,3,4,$$ (27)

where $\sigma_m^c$ is the standard deviation of the model data points around the mean and $\sigma_n^c$ is the standard deviation of the simulated data points around the model curve. $\sigma_m^c$ is chosen for two reasons: (1) The true amplitude of the model curve would be a natural choice, however, for complex shapes with many harmonics the definition of the amplitude is not straightforward; (2) We are dealing with only the discrete time points and there is no guarantee that the true maximum/minimum of the model curve appears among these points. When the noise level increases, the S/N^c decreases. Starting from a certain noise level, the correct peak in the spectrum will be less deep than some of the nearby random fluctuations. We marked in the Tables 3 and 4 with boldface those noise levels where the correct peak occured for the last time as the strongest minimum. It is well seen that the multichannel version of the spectrum reveals the true period for higher noise levels compared to the single channel spectra. Due to the random nature of the simulations and different noise levels in the different channels, the breakdown noise levels cannot be predicted precisely using analytical expressions.

Figure 3: Detection levels in the PDM search for the first group of simulations. Vertical bars denote the S/N range over four channels. The last successful detections are marked with U, B, V and R for computations of a single channel and with UBVR for the multichannel method

Figure 4: Detection levels in the PDM search for the second group of simulations. For details, see Fig. 3

Figure 5: Set of PDM spectra computed for the single channels and for the combined multichannel of set 7 (noise level 3.0 Ec) in the second group of simulations. The numbers label the depth ranks of the valleys in the spectra

Figure 6: Set of PDM spectra computed for the single channels and for the combined multichannel of set 12 (noise level 8.0 Ec) in the second group of simulations

In Figs. 3 and 4, these results are depicted graphically. We can see the clear difference between detection levels for the PDM search. The main peak in the spectra for the U, B, V and R channels "sinks'' into the noise significantly earlier than the main peak in the combined spectrum. As an illustration, we depict some actual spectra in Figs. 5 and 6. Here, it needs to be remarked that the deepest minimum at 0.698 in the (V) spectrum of Fig. 6 was produced accidentally by random fluctuations in the data set 12. In the second group of simulations, the V channel has already stopped working from the set 7 (Fig. 4), for which Fig. 5 shows its last correctly working spectrum.

Our simple tests with the restricted period search range well demonstrate that the expected improvement in precision and in period detection capability indeed takes place.

For a wider search range this nice picture is spoilt by spurious periods. We applied PDM for the range of periods $P_{\rm min} = 0.3$ and $P_{\rm max} = 1.9$ to obtain the single channel and the combined spectra which are depicted in Fig. 7. All five spectra show the four most significant minima at periods 0.349, 0.537, 1.161 and 0.698. The correct period P=0.698 is always detected but it is far from the best one! If we look at the corresponding data window which is computed as in Deeming (1975), we see that our data spacing is nearly periodic with $\delta = 0.99728738$ and consequently every peak in the spectrum is echoed multiple times. Using Eq. (19), we can identify spurious periods which are actually echoes of the real period. We summarize the results of the identification of all the nine spurious periods shown in Fig. 7 in Table 5.

 

 

Table 5: Identification of spurious periods with the Deeming window of set 5 in the second group of simulations. No. is the rank of the period shown in the combined spectrum of (UBVR) in Fig. 7, $P_{\rm PDM}$ is the value found with the multichannel PDM and P_r,l,s are computed from Eq. (19) with the corresponding (r,l,s) and $\delta = 0.99728738$

No.	$P_{\rm PDM}$	Pr,l,s	(r,l,s)
1	0.3490423	0.3490497	(1, 2, 0)
2	0.5365817	0.5369987	(1, 2, -1)
3	1.1614638	1.1634916	(1, 2, -2)
5	1.0733655	1.0739974	(2, 2, -1)
6	0.5169787	0.5171110	(2, 2, 1)
7	0.8743941	0.8726323	(1, -2, 4)
8	0.4107931	0.4106471	(1, 1, 1)
9	1.7493248	1.7452646	(1, -1, 2)
10	0.4662043	0.4654024	(1, -2, 5)

Figure 7: PDM spectra of set 5 (noise level 1.0Ec) in the second group of simulations. In all five spectra the real period is fourth ranked according to the depth of a peak. All others are spurious and can be identified according to Eq. (19)

The spurious period of 0.349 with the deepest minimum occurs at the location where P_{1, 2, 0} = 0.34904965; the spurious period 0.537 with the second deepest minimum occurs at the location where P_{1, 2, -1} = 0.53699869 and finally the spurious period of 1.161 with the third deepest minimum occurs at P_{1, 2, -2} = 1.1634916. Because the amplitudes of the second sinusoidal function is much larger than those of the first one in the simulated model (Table 2), the model curve resembles that of only one sinusoidal function with half the correct period. Thus, this results in the large significance of the first three spurious periods and also for P_2,2,-1 and P_2,2,1which occur due to fitting with only the (r = 2) sinusoidal function (see explanation in Sect. 2.4).

This is a very important observation. The interplay of the real periodicity with a periodicity in the time point distribution can totally change the general appearance of the spectra. Even when we use all available information (as in the case of the multichannel search), the spurious periods can still complicate our analysis. Paradoxically enough, the improved detection capability of multichannel methods increases also the ability to detect spurious periods. So, when working with real data we must always carefully inspect the spectra using the information from the corresponding data windows to reveal the real periods.

In the last two columns of Tables 3 and 4 the refined periods obtained with LM and NLM are given. Figures 8 and 9 show the final refined periods computed with MPA for the two groups of simulations. For comparison, the results from the search using the best single channel (with highest S/N) are also given (if available). We see that MPA gives as a rule more precise estimates for the main period than those from only the best single channel.

Because the differences between standard least squares methods for the single channels and for the multichannel are minimal, we are not going here into the details of computation. There are only two important additional considerations to be taken into account. The first one is the correct choice of the frequency step for the LM and the correct frequency bracketing for the final NLM refinement. Because of the improved detection capability of the multichannel method, we can in principle choose somewhat larger step sizes for LM and somewhat narrower brackets for a final period improvement. This is certainly true for exact periodicities. But in real life situations, the periods can be slightly varying and, as a result of this, the minima around the periods can have less regular profiles.

Figure 8: Periods (filled circles) and their errors ($3\sigma$ bars) recovered with MPA for the first group of simulations. The x-axis is the data set number and the y-axis is the period in 10-3. The horizontal line is at the correct period 0.6980993. Additionally (if available), the best estimates for periods from the single channel search are depicted as open squares with $3\sigma$ error bars

Figure 9: Multi and single channel search results for the second group of simulations. See detailed explanations in Fig. 8

The second consideration is the error estimation for the final computed periods. In the ideal case of exact periods with reasonably even sampling and full time coverage, the error estimates based on the curvature of WRSS(P) can give good estimates for period errors. But for more complex cases, bootstrap type procedures must be involved (see Jetsu & Pelt 1999). The bootstrap for multichannel methods must take into account the special structure of the input data and will only work if the residuals are uncorrelated. It is reasonable to reshuffle residuals for each channel separately with correct reweighting to ensure that the statistical structure of each bootstrap run will be analogous to the original data.