A&A 467, 1353-1371 (2007)
DOI: 10.1051/0004-6361:20066597

SigSpec

I. Frequency- and phase-resolved significance in Fourier space

P. Reegen

Institut für Astronomie, Universität Wien, Türkenschanzstraße 17, 1180 Vienna, Austria

Received 19 October 2006 / Accepted 6 March 2007

Abstract
Context. Identifying frequencies with low signal-to-noise ratios in time series of stellar photometry and spectroscopy, and measuring their amplitude ratios and peak widths accurately, are critical goals for asteroseismology. These are also challenges for time series with gaps or whose data are not sampled at a constant rate, even with modern Discrete Fourier Transform (DFT) software. Also the False-Alarm Probability introduced by Lomb and Scargle is an approximation which becomes less reliable in time series with longer data gaps.
Aims. A rigorous statistical treatment of how to determine the significance of a peak in a DFT, called S IGS PEC, is presented here. S IGS PEC is based on an analytical solution of the probability that a DFT peak of a given amplitude does not arise from white noise in a non-equally spaced data set.
Methods. The underlying Probability Density Function (PDF) of the amplitude spectrum generated by white noise can be derived explicitly if both frequency and phase are incorporated into the solution. In this paper, I define and evaluate an unbiased statistical estimator, the "spectral significance'', which depends on frequency, amplitude, and phase in the DFT, and which takes into account the time-domain sampling.
Results. I also compare this estimator to results from other well established techniques and assess the advantages of S IGS PEC, through comparison of its analytical solutions to the results of extensive numerical calculations. According to those tests, S IGS PEC obtains as accurate frequency values as a least-squares fit of sinusoids to data, and is less susceptible to aliasing than the Lomb-Scargle Periodogram, other DFTs, and Phase Dispersion Minimization (PDM). I demonstrate the effectiveness of S IGS PEC with a few examples of ground- and space-based photometric data, illustratring how S IGS PEC deals with the effects of noise and time-domain sampling in determining significant frequencies.

Key words: methods: data analysis - methods: statistical

  
1 Overview

In this Paper I provide a brief introduction to Fourier methods in astronomical time series analysis (Sect. 2), outline existing statistical approaches (Sect. 3), and address the major weaknesses in available techniques (Sect. 4).

While Sects. 2 to 4 contain previously available (textbook) information, I introduce a new and unbiased reliability criterion (spectral significance) based on theoretical statistics in Sect. 5. This section also addresses the correspondence between the spectral significance and other reliability estimators. Section 6 is devoted to the comparison of the analytically deduced spectral significance to the results of numerical simulations. Finally, the results of comparative tests of spectral significance computation vs. other period detection methods are presented.

An example for the practical application of the new method vs. a widely used standard procedure is provided in Sect. 7. It may be useful for the non-mathematically oriented reader who is mainly interested in how the technique performs in "real life''.

Further topics (Sect. 8) are the application of statistical weights, the fact that a time series consists of individual subsets, and potential problems with colored vs. white noise.

  
2 Introduction

In Fourier Analysis a continuous function of time over a finite time interval is expanded into a series of sine waves. These waves represent a superposition of an oscillation at a fundamental frequency and a discrete, generally infinite set of overtones. The fundamental frequency is determined by the reciprocal time interval width. All other frequencies correspond to integer multiples of this fundamental. The knowledge of the amplitudes and phase angles for all frequency components permits one to entirely recover the given function in the time domain.

Practical applications (such as astronomical observations) generally deal with discrete sets of measurements (time series) rather than continuous functions of time and, on the other hand, consider the Fourier Spectrum as a continuous function of frequency rather than a discrete dataset. This leads to the Discrete Fourier Transform (DFT). It allows one to determine the dominant frequencies of the observed physical process with a higher frequency resolution than is possible with Fourier Analysis.

Motivated by the desire to understand physical oscillations, the scientist is interested in a couple of eigenfrequencies and the exact determination of related amplitudes and phases rather than the complete signal recovery. In practical applications, these are considered to correspond to local maxima (peaks) in the amplitude spectrum. A widely held strategy is to

1.

compute an amplitude spectrum for the given dataset;

2.

identify the maximum amplitude within the frequency range of interest;

3.

decide whether this amplitude is "significant'';

4.

subtract the corresponding sinusoidal signal from the time series; and

5.

use the residuals after subtracting the fit from the time series for the next iteration.

This procedure is to be understood as a loop, terminated if the maximum amplitude is not considered significant any more. The result of these consecutively performed prewhitenings is a list of frequencies, amplitudes, and phase angles, plus a residual time series (hopefully) representing the pure observational noise. In fact, part of this noise is due to measurement errors, but frequently merged with signal components the amplitudes of which are too weak to be detected. As an example, the number of photometrically resolved frequencies in the $\delta$ Sct star 4 CVn increased from 5 to 34 between 1990 and 1999 (Breger et al. 1990, 1999).

In many cases, the results of the prewhitening are subject to a multiperiodic least-squares fit (e.g., Sperl 1998; Lenz & Breger 2005). This represents a fine tuning to adjust frequencies, amplitudes, and phases to a minimum rms residual, but may lead to the exclusion of some terms, or inclusion of new terms.

The eigenfrequencies, amplitudes, and phases of stellar oscillations provide fundamental information on the distribution of mass and temperature, the radiative and convective energy transport, or abundances of elements in the stellar atmosphere and help to determine fundamental parameters such as mass, radius, effective temperature, rotational velocity, and age of a star.

  
3 Statistical aspects

In practical applications, a signal is not only of stellar origin but a superposition of the information received from the star, instrumental (pseudo-)periodicities (e.g., invoked by thermal effects), and transparency variations in the Earth's atmosphere. To eliminate the third component, observations with instruments aboard of spacecraft have been used increasingly by the astronomical community during the past two decades. Unfortunately, stray light scattered from the illuminated surface of the Earth introduces additional quasi-periodic artifacts that are not easy to handle and hence represent the major constraint to the accuracy of space-based data acquisition (Reegen et al. 2006).

An unbiased criterion to decide, whether a peak amplitude is generated by noise or represents an intrinsic variation, is important, because the choice of the most significant peak determines all further iterations in the prewhitening sequence. A falsely identified signal component usually perturbs all results obtained subsequently.

In addition, Scargle (1982) points out that Gaussian noise in the time domain may produce a peak of arbitrary amplitude in the DFT spectrum. Since there is no natural upper limit to amplitudes produced by noise, the danger of misinterpretation is imminent, and the "significance'' of an examined peak needs to be described by a probability distribution.

Of course, the probability that white noise produces a low-amplitude peak is higher than that for a high-amplitude peak. In other words, the False-Alarm Probability that the highest peak in an amplitude spectrum is an artifact due to noise is lower for higher amplitudes. This definition is reasonable because it relies on the highest available peak, since one cannot trust any peak, if even the one with the highest amplitude is unreliable. In this perspective, the False-Alarm Probability appears to be a good criterion of whether to believe in the presence of a signal in a given dataset.

The statistical description of the False-Alarm Probability relies on a Probability Density Function (PDF), which is the continuous version of a histogram, as the bins of which become infinitesimally small. Thus the False-Alarm Probability is an integral over the PDF. In many statistical applications, probability distributions may easily be described by the PDF, but an analytic expression of the integral does not exist, as e.g. for the Gaussian distribution. Hitherto many statistical problems are solved in terms of PDFs rather than of the corresponding cumulative quantities.

The PDF of the amplitude spectrum generated by pure Gaussian noise at equidistant sampling was deduced by Schuster (1898; also Scargle 1982). For non-equidistant sampling, the Lomb-Scargle Periodogram (Lomb 1976) is claimed to provide a better statistical behavior (Scargle 1982). Koen (1990) examined the application of the Fisher test (1929) to the Lomb-Scargle Periodogram to estimate peak significance, and Schwarzenberg-Czerny (1996) combined Lomb's solution with the analysis of variance method (Schwarzenberg-Czerny 1989) to obtain a powerful period search technique for non-equidistant data.

A widespread, but purely empirical approach is the consideration of a peak in the amplitude spectrum to be "real'', if its amplitude signal-to-noise ratio exceeds a given limit. Breger et al. (1993) suggest a value of 4. Numerical simulations for 19 300 synthetic time series by Kuschnig et al. (1997) return an associated False-Alarm Probability of 10^-3 for 1000 data points.

An alternative way of determining periods was introduced by Lafler & Kinman (1965) and statistically examined by Stellingwerf (1978). The Phase Dispersion Minimization method is based on the assumption that the correct period would produce a phase diagram with the lowest intrinsic scatter. This algorithm is a powerful tool especially for non-sinusoidal periodicities, e.g., if the stellar surface structure is observed photometrically as the star rotates. In this situation, the advantage compared to Fourier techniques is the simultaneous treatment of all overtones that determine the shape of the signal.

In the case of multiperiodic variability, PDM loses accuracy, since the phase diagram for one frequency is contaminated by all others, unless they are integer multiples of this frequency.

  
4 Problems

A closer examination of the fundamental properties of the DFT leads to the following issues:

1.

The peak frequency in the amplitude spectrum of a single sinusoidal signal is not recovered correctly (Kovacs 1980). This is a side effect of the restriction to a finite time interval rather than a property of non-equidistant sampling. Also the Lomb-Scargle Periodogram suffers from systematic peak frequency displacements (see 6.2).

2.

The "entire'' spectrum is defined by the Nyquist (Critical) Frequency as the upper limit for equidistantly sampled data (Whittaker 1915; Nyquist 1928; Kotelnikov 1933; Shannon 1949; Press et al. 1992). In the case of non-equidistant sampling, there is no unique definition of such a limit. Horne & Baliunas (1986) suggest the average sampling interval width for the calculation of the Nyquist Frequency; Scargle (1982) and Press et al. (1992) promote the minimum sampling interval; a variant discussed by Eyer & Bartholdi (1999) is to pad the entire dataset with zeros to achieve equidistant, but much denser sampling in the time domain and to take the resulting sampling interval as a Nyquist Frequency estimator. The most promising method was introduced by Sperl (1998): each sampling interval of the non-equidistant series is considered responsible for its own, individual Nyquist Frequency. In a subsequent step, a histogram of all individual Nyquist Frequencies helps to decide where to define a healthy upper frequency limit.

3.

Periodic gaps in the time-domain sampling introduce periodicities in the Fourier Spectrum. In this case, a single sinusoidal signal is represented by a peak that is accompanied by several aliases. In order to overcome the systematic effects of non-equidistant sampling in the frequency domain, Ferraz-Mello (1981) examined the possibility of incorporating the sampling properties into the DFT amplitude spectrum by introducing appropriate statistical weights. His study led to the design of harmonic filters, which may help to improve the selection of peaks. A different approach is the consideration of the entire "comb'' of aliases instead of the peak with the highest amplitude only (Roberts et al. 1987; Foster 1995).

The goals of the present paper are

• to deduce the depedence of the Probability Density Function of non-equidistantly sampled Gaussian noise on frequency and phase (Sect. 5);
• to introduce spectral significance as a measure of False-Alarm Probability (5.5);
• to compare the theoretical solution with the results of numerical simulations (6.1); and
• to quantitatively oppose the accuracy of the peak frequencies using spectral significance (6.2, 6.3) vs. the "simple'' DFT amplitude, the Lomb-Scargle Periodogram, PDM, and the DFT amplitude plus improvement of resulting peak frequencies by least-squares fitting.

A frequent problem in the context of astronomical observations is that all measurements do not necessarily have the same variance. This may be due to different instrumentation for different subsets of the time series or changing environmental conditions, such as thermal noise or transparency variations in the Earth's atmosphere. Statistical weights are introduced into the spectral significance in order to take appropriate account of the variable quality of measurements within a single dataset (8.1).

  
5 Amplitude PDF for non-equidistant sampling

Figure 1: Schematic illustration of a time series with periodic gaps representing Gaussian noise with a standard deviation equal 1 ( top) to visualize the dependency of the corresponding frequency-domain noise on frequency and phase, as well as on the characteristics of the time-domain sampling. In the two examples referred to in the lower panels, a test signal is displayed ( dashed lines) for which the DFT shall be evaluated. If the frequency and phase are combined such that the data points consistently align around the times when the trigonometric function attains a small value, only a small fraction of the time-domain noise will be transformed into the frequency domain ( mid). For a combination of frequency and phase grouping the data points around the maxima and minima of the test signal, the time-domain noise will produce a higher frequency-domain noise, correspondingly ( bottom). The first of these two cases produces a narrower probability distribution in Fourier Space, and consequently, a signal with the same amplitude will be considered more reliable in the first case than in the second.

Open with DEXTER

This section presents the theoretical evaluation of the frequency-domain PDF for a non-equidistantly sampled time series representing Gaussian noise.

A very important fact for the statistical analysis is, that in most practical applications, the mean value of the observable is shifted to a fixed value (frequently zero) before the DFT is evaluated. The statistical description of frequency-domain noise has to take this fact into account.

In statistical terms, the Fourier Coefficients are regarded as weighted sums of random variables. Such sums tend towards a Gaussian distribution as datasets become large enough. Since the Fourier Coefficients for each frequency represent the Cartesian components of the two-dimensional Fourier Space, the resulting Gaussian distribution will be two-dimensional (bivariate) as well. Furthermore, as the Fourier Coefficients are functions of frequency, the considered Fourier-Space probability distribution will depend on frequency, correspondingly.

For a time series with gaps (Fig. 1), the relevant influence of the time-domain sampling on the frequency-domain probability distribution is determined by the phase coverage of the measurements: a combination of frequency and phase for which the intervals containing data are consistently associated to angles where the corresponding trigonometric function ("test signal'') attains low numerical values will result in a lower noise level than a combination allocating the data close to angles where the maxima/minima of the test signal are located. The first case will yield a narrow probability distribution, the second case a broad one. Consequently, a signal with the same amplitude has be considered more significant in the first case. These strong phase dependencies are mitigated for frequencies providing a better phase coverage.

Consequently, the amplitude distribution in Fourier Space will have to be a function of both frequency and phase, which is achieved by expressing the bivariate Gaussian PDF in polar rather than Cartesian coordinates.

The probability of a peak generated by noise to reach a given amplitude level may be evaluated through integration of the PDF over amplitude, which leads to the Cumulative Distribution Function (CDF) and False-Alarm Probability, based on which a more informative quantity, the spectral significance, is defined.

  
5.1 Zero-mean correction

In astronomical applications, magnitudes are usually averaged to a pre-defined constant (zero or non-zero, as obtained by some theoretical concept or calibration) before the amplitude spectrum is evaluated. The following considerations apply to observables adjusted to zero mean. If a non-zero constant is chosen instead, the DFT will change, but the False-Alarm Probabilities will remain the same. In this case, one prefers to evaluate a DFT spectrum for the time series as it is but use a zero-mean corrected version of the dataset for the computation of spectral significances.

Considering a time series $x_k := x\left( t_k\right)$ to be generated by a Gaussian random process with expected value 0 and population variance $\left< x^2\right>$, the time-domain PDF is

$$\phi\left( x_k\right) := \frac{1}{\sqrt{2\pi\left< x^2\right>}}\:{\rm e}^{-\frac{x_k^2}{2\left< x^2\right>}} .$$ (1)

Given a random process that produces an infinite population of Gaussian random variables, the mean of a finite sample of random variables x_k is free to scatter around the population mean. Gauss's law of error propagation returns a variance of the sample mean which is proportional to the inverse number of data points in the sample,

$$\left<\left< x_k\right> ^2\right> = \frac{\left< x^2\right>}{K}\cdot$$ (2)

If the finite sample is artificially adjusted to zero average, the sample mean value is not allowed to scatter any more. Since the standard error of an individual data point implicitly contains the standard error of the mean, too, zero-mean correction will distort the PDF of the random variable. The only exception with an invariant PDF is the Fourier Analysis, i.e., equidistant time-domain sampling and a set of discrete frequencies.

An alternative (and promising) method in this context is the Floating-Mean Periodogram (Cumming et al. 1999), which is based on a least-squares fit of a sinusoid plus a constant to the time series, the latter retaining the free scatter of the sample mean.

  
5.2 Distribution of Fourier coefficients

Incorporating the effect of zero-mean correction into the statistical examination, the zero-mean corrected magnitude values $x_k - \frac{1}{K}\sum_{k=0}^{K-1}x_k$ have to be used for the calculation of Fourier Coefficients according to

  

$$a_{{\rm ZM}}\left(\omega\right) := \frac{1}{K}\sum_{k=0}^{K-1}x_k\cos\omega t_k - \frac{1}{K^2}\sum_{k=0}^{K-1}x_k\sum_{l=0}^{K-1}\cos\omega t_l ,$$ (3)

$$b_{{\rm ZM}}\left(\omega\right) := \frac{1}{K}\sum_{k=0}^{K-1}x_k\sin\omega t_k - \frac{1}{K^2}\sum_{k=0}^{K-1}x_k\sum_{l=0}^{K-1}\sin\omega t_l .$$ (4)

Due to the linearity of the Fourier Transform, the subtraction of a constant in the time domain refers to a subtraction of a spectral window in the frequency domain.

Rearrangement of indices yields

  

$$a_{\rm ZM}\left(\omega\right) = \frac{1}{K}~\sum_{k=0}^{K-1}x_k\left(\cos\omega t_k - \frac{1}{K}~\sum_{l=0}^{K-1}\cos\omega t_l\right) ,$$ (5)

$$b_{\rm ZM}\left(\omega\right) = \frac{1}{K}~\sum_{k=0}^{K-1}x_k\left(\sin\omega t_k - \frac{1}{K}~\sum_{l=0}^{K-1}\sin\omega t_l\right) .$$ (6)

Given pure Gaussian noise in the time domain with a population variance $\left< x^2\right>$, Eqs. (5), (6) allows one to consider both Fourier Coefficients as linear combinations of Gaussian variables with expected values 0 and variances

  

$$\left< a_{\rm ZM}^2\right>\left(\omega\right) = \frac{\left< x^2\... ...1}\left(\cos\omega t_k - \frac{1}{K}~\sum_{l=0}^{K-1}\cos\omega t_l\right) ^2 ,$$ (7)

$$\left< b_{\rm ZM}^2\right>\left(\omega\right) = \frac{\left< x^2\... ...1}\left(\sin\omega t_k - \frac{1}{K}~\sum_{l=0}^{K-1}\sin\omega t_l\right) ^2 .$$ (8)

Thanks to the Central Limit Theorem (de Moivre 1718; Stuart & Ord 1994, p. 310f), the consideration of these coefficients as Gaussian variables holds to a sufficient degree, even if the time-domain noise is not Gaussian, because even short datasets in astronomical applications are long enough compared to the fast convergence of the PDF towards the Gaussian distribution with an increasing number of random variables.

  
5.3 Frequency- and phase-dependent PDF

Since the DFT produces a two-dimensional vector $\left( a,b\right)$, the probability distribution in the frequency domain will also be two-dimensional, so-called bivariate.

The combined probability density of two independent Gaussian variables $\alpha$, $\beta$ with corresponding variances $\left<\alpha ^2\right>$, $\left<\beta ^2\right>$ is given by a bivariate Gaussian PDF,

$$\phi\left(\alpha ,\beta\right) = \frac{1}{2\pi\sqrt{\left<\a... ... ^2\right>}+\frac{\beta ^2}{\left<\beta ^2\right>}\right)}\: ,$$ (9)

if the covariance $\left<\alpha\beta\right>$ vanishes.

The fact that the Fourier Coefficients $a_{\rm ZM}$, $b_{\rm ZM}$ of pure noise are two linear combinations of the same random vector in Fourier Space diminishes the degrees of freedom by 1. Hence they may be considered independent to a sufficient degree, if the sample size is large enough. Consequently, if $\left< a_{\rm ZM}b_{\rm ZM}\right> = 0$, Eq. (9) describes the bivariate distribution of Fourier Coefficients related to noise satisfactorily. According to Appendix A, rotating the Fourier Space coordinates by an angle $\theta _0$ given by

$$\tan 2\theta _0\left(\omega\right) =\frac{K\sum_{k=0}^{K-1}\... ... t_k\right) ^2+\left(\sum_{k=0}^{K-1}\sin\omega t_k\right) ^2}$$ (10)

transforms the Fourier Coefficients $a_{\rm ZM}$, $b_{\rm ZM}$ into coefficients $\alpha$, $\beta$ with zero covariance, as desired.

The DFT of a measured time series x_k contains only an amplitude A but also a phase angle

$$\theta\left(\omega\right) = \frac{\sum_{k=0}^{K-1}x_k\sin{\omega t_k}}{\sum_{k=0}^{K-1}x_k\cos{\omega t_k}}\cdot$$ (11)

This additional information may be taken into account by evaluating the conditional probability density of amplitude for a constant phase angle $\theta$, pre-defined by the DFT of the time series under consideration at the frequency $\omega$.

The transformation of Eq. (9) from Cartesian into polar coordinates $\left( A,\theta\right)$ is performed via ${\rm d}\alpha~{\rm d}\beta = A~{\rm d}A~{\rm d}\theta$ (Appendix B.1) with

$$\alpha = \frac{A}{2}\cos\left(\theta - \theta _0\right) ,$$ (12)

$$\beta = \frac{A}{2}\sin\left(\theta - \theta _0\right) ,$$ (13)

where the fact, that the coordinate system is rotated by a constant angle $\theta _0$, does not contribute to the Jacobian of the transformation. The division by 2 is introduced by collecting the contributions of both $A\left(\omega\right)$ and $A\left( -\omega\right)$ - which are equal for real observables - to the total amplitude. The transformed amplitude PDF becomes

$$\phi\left( A,\theta\left\vert\right.\omega\right) = \frac{A}... ...eft(\theta -\theta _0\right)}{\left<\beta ^2\right>}\right]} .$$ (14)

Of course, $\phi$ does not only depend on amplitude A and phase $\theta$, but also on $\alpha _0$, $\beta _0$, and $\theta _0$, which are determined by the time-domain sampling and are functions of frequency $\omega$. The bar symbol in $\phi\left( A,\theta\left\vert\right.\omega\right)$ is introduced to formally separate random variables to those considered constant.

Changing the normalization condition from $\int\!\!\!\int _{~{\cal R} ^2}{\rm d}A~{\rm d}\theta~\phi\left( A,\theta\left\vert\right.\omega\right) = 1$ into $\int _{~\cal{R}} {\rm d}A~\phi\left( A\left\vert\right.\omega ,\theta\right) = 1$ yields

$$\phi \left( A\left\vert\right.\omega ,\theta\right) = \frac{A}{4R^2}~{\rm e}^{-\frac{A^2}{8R^2}}$$ (15)

with

$$R := \sqrt{\frac{\left<\alpha ^2\right>\left<\beta ^2\right>... ...<\alpha ^2\right>\sin ^2\left(\theta - \theta _0\right)}}\cdot$$ (16)

The difference between Eqs. (14) and (15) is that Eq. (14) is a bivariate PDF of amplitude and phase, whereas in Eq. (15) only the amplitude A is considered as a random variable, and $\theta$ is constant. This relation returns the probability density of amplitude for a fixed frequency and a fixed phase in Fourier Space. Accordingly, the PDF is normalized by the condition that its integral over the entire amplitude range (from 0 to $\infty$) has to be 1. Furthermore, amplitudes being defined $\ge$0 introduce a factor 2 into the argument of the exponential function.

Equation (16) presents $R\left(\theta\right)$ as an ellipse in polar coordinates. The semi-major and semi-minor axes are $\sqrt{\left<\alpha ^2\right>}$ (Eq. (A.3)) and  $\sqrt{\left<\beta ^2\right>}$ (Eq. (A.4)), respectively, and the orientation is determined by $\theta _0$. This ellipse will be called the rms error ellipse. Its orientation and dimensions depend on frequency.

Equation (10) has got a set of solutions for $\theta _0$ assigned to orthogonal directions: if $\theta _0$ is a solution, then the complete set of solutions is $\theta _0 + \frac{z}{2}\pi$ $\forall z\in\cal{Z}$. Whether $\sqrt{\left<\alpha ^2\right>}$ returns the semi-major axis and $\sqrt{\left<\beta ^2\right>}$ the semi-minor axis of the rms error ellipse, or vice versa, depends on the choice of $\theta _0$. This paper consistently uses solutions of $\theta _0$ that assign $\alpha _0$ to semi-minor axes, which yields the maximum spectral significance for all phase angles under consideration.

As shown later (see 5.6), the introduction of the normalized semi-major and semi-minor axes,

  

$$\alpha _0\left(\omega ,\theta _0\right) := \sqrt{2K\frac{\left<\a... ...um_{l=0}^{K-1}\cos\left(\omega t_l - \theta _0\right)\right] ^2\right\rbrace} ,$$ (17)

$$\beta _0\left(\omega ,\theta _0\right) := \sqrt{2K\frac{\left<\be... ...um_{l=0}^{K-1}\sin\left(\omega t_l - \theta _0\right)\right] ^2\right\rbrace} ,$$ (18)

respectively, provides the separation of sampling-dependent quantities from quantities that only depend on the DFT amplitude. In this context, the term "normalized'' means that for arguments $\omega t_k - \theta _0$ uniformly distributed on $\left[ 0,2\pi\right]$, the expected values of both $\alpha _0$ and $\beta _0$ are 1.

Given the orientation of axes, $\theta _0$, and the normalized axes, $\alpha _0$, $\beta _0$, of the ellipse, the standard deviation for an arbitrary phase $\theta$ in Fourier Space is the radius of an ellipse in polar coordinates (R, $\theta$) according to

$$R = \sqrt{\frac{\left< x^2\right>}{2K}\:\frac{\alpha _0^2~\b... ...ght) +\alpha _0^2\sin ^2\left(\theta - \theta _0\right)}}\cdot$$ (19)

The three parameters $\theta _0$, $\alpha _0$, $\beta _0$ describe the ellipticity of the two-dimensional PDF for Gaussian noise in Fourier Space and thus represent the cornerstones of spectral significance evaluation. All subsequent relations will be given in terms of these three quantities.

  
5.4 False-Alarm Probability

The Cumulative Distribution Function (CDF) is obtained by integrating the PDF (Eq. (15)) according to

$$\Phi \left( A\left\vert\right.\omega ,\theta\right) = \int _0... ...ime\phi\left( A^\prime\left\vert\right.\omega ,\theta\right) ,$$ (20)

which yields

$$\Phi \left( A\left\vert\right.\omega ,\theta\right) = 1 - {\rm e}^{-\frac{A^2}{8R^2}} .$$ (21)

Thus the probability for an amplitude to exceed a given limit A is given by

$$\Phi _{\rm FA}\left( A\left\vert\right.\omega ,\theta\right) = {\rm e}^{-\frac{A^2}{8R^2}} ,$$ (22)

which is the False-Alarm Probability of an amplitude level at phase $\theta$ (and frequency $\omega$, since $\theta$, $\theta _0$, $\alpha _0$, and $\beta _0$ are frequency-dependent quantities).

  
5.5 Spectral significance

The frequency- and phase-dependent False-Alarm Probability of an amplitude level was introduced as the probability that random noise in the time domain with the same rms error as the given time series produces a peak in the DFT amplitude spectrum which is at least as high as the corresponding amplitude level for the time series itself: if a peak is assigned a False-Alarm Probability of 0.00 001, its risk of being due to noise is 1:100 000. In this section, the spectral significance as a more informative quantity is introduced (5.5.1). It is the inverse False-Alarm Probability (in this case, 100 000) scaled logarithmically. In the present example, the conversion of a False-Alarm Probability of 0.00 001 into spectral significance returns 5. In this context, the spectral significance is presented as a logarithmic measure for the number of cases in one out of which the considered amplitude would be an artifact.

Plotting the spectral significance vs. frequency yields the significance spectrum, and the identification and consideration of the highest peak in this spectrum may lead to a statement on the significance of the entire spectrum. The argument is similar to many existing significance estimates (e.g., Scargle 1982): if the highest peak in the spectrum is below some limit, the entire spectrum has to be considered insignificant. But instead of using the (statistically biased) signal-to-noise ratio as a threshold, the spectral significance is employed. The application of the spectral significance concept to the highest peak out of a sample is briefly discussed (5.5.2).

The formal correspondence to traditional techniques, namely signal-to-noise ratio and the Lomb-Scargle Periodogram, is of special interest and hence provided subsequently (5.5.3, 5.5.4).

  
5.5.1 Definition

To enhance the compatibility to the popular signal-to-noise ratio criterion (see 5.5.3), the spectral significance of a DFT amplitude is defined as

$${\rm sig}\left( A\left\vert\right.\omega ,\theta\right) := -... ...\Phi _{\rm FA}\left( A\left\vert\right.\omega ,\theta\right) ,$$ (23)

or - using Eqs. (19) and (22) -

$${\rm sig}\left( A\left\vert\right.\omega ,\theta\right) = \f... ...c{\sin ^2\left(\theta - \theta _0\right)}{\beta _0^2}\right] ,$$ (24)

with the normalized axes, $\alpha _0$ and $\beta _0$, as defined by Eqs. (17), (18), and the orientation of the ellipse in Fourier Space according to Eq. (10). Since the angle $\theta _0$ was chosen to refer to $\alpha _0$ as the semi-minor axis of the rms error ellipse, $\alpha _0$ now corresponds to the phase of maximum spectral significance for a given frequency.

The concept of spectral significance computation relies on the analytical comparison of the DFT amplitude generated by a measured time series to noise at the same variance as the time series under consideration. Unless the population variance $\left< x^2\right>$ of the noise used for comparison is given by theory and/or other observations than the ones under consideration, the sample variance $\left< x_k^2\right>$ of the observable may be used as an estimator.

The Cartesian representation of Eq. (24),

 

$${\rm sig}\left( a_{{\rm ZM}},b_{{\rm ZM}}\left\vert\right.\omega\right) \frac{K\log e}{\left< x^2\right>}\left[\left(\frac{a_{{\rm ZM}}\cos\theta _0+b_{{\rm ZM}}\sin\theta _0}{\alpha _0}\right) ^2\right.$$ =

$$+ \left.\left(\frac{a_{{\rm ZM}}\sin\theta _0-b_{{\rm ZM}}\cos\theta _0}{\beta _0}\right) ^2\right] ,$$ (25)

is useful for practical applications and will be employed for the implementation of statistical weights (see 8.1). The subscript "ZM'' indicates zero-mean corrected time series data, in consistency with Eqs. (3), (4).

Thanks to the logarithmic scaling, the spectral significance appears as a product form (as opposed to an exponential function), where one factor (the bracket term) contains all information on the ellipticity of the underlying PDF in Fourier Space and is entirely determined by the time-domain sampling. This term is scaled according to the squared amplitude. The serendipitous consequence of this separation is that the evaluation of the bracket term applies to all datasets with the given sampling. In a prewhitening cascade, the sampling of the time series will not change. Consequently, the bracket term remains valid for the entire sequence and has to be computed only once. In the prewhitening cascade itself, it is sufficient to rescale the bracket term by the squared amplitude, which speeds up the computations considerably.

Thanks to this formal separation, it is possible to pack all the characteristics of the time-domain sampling into an amplitude-independent function of frequency and phase. This will lead to the Sock Diagram (5.6).

A further practical advantage of this separation is the occurrence of the population variance $\left< x^2\right>$ independently of frequency and phase. For small samples, the higher uncertainty of the estimated population variance may be overcome by using the sample variance $\left< x_k^2\right>$ instead of the population variance $\left< x^2\right>$ and increasing the spectral significance limit for peak acceptance accordingly.

  
5.5.2 Spectral significance for a statistically independent sample

One may desire to evaluate the spectral significance for the highest out of a sample of peaks in the significance spectrum - in analogy to the procedure presented by Scargle (1982). His arguments may be applied to the spectral significance directly.

For a given spectral significance level ${\rm sig}$, the probability of an amplitude level generated by a noise process to exceed the spectral significance limit ${\rm sig}$ is the False-Alarm Probability  $\Phi _{\rm FA}$. It is linked to the spectral significance via Eq. (23). The complementary probability that such an amplitude level is below ${\rm sig}$ is $1 - \Phi _{\rm FA}$. Given a sample of N such amplitude levels, which are statistically independent, the probability of none exceeding ${\rm sig}$ is $\left( 1 - \Phi _{\rm FA}\right) ^N$. Again, the complement, $1 - \left( 1 - \Phi _{\rm FA}\right) ^N$, returns the probability for at least one peak out of the sample to exceed ${\rm sig}$. In other words (and using Eq. (23) to substitute for  $\Phi _{\rm FA}$), the False-Alarm Probability for the maximum of a statistically independent sample of N peaks is

$$\widehat{\Phi _{\rm FA}} = 1 - \left( 1 - 10^{-{\rm sig}}\right) ^N .$$ (26)

This False-Alarm Probability may be transformed into a spectral significance (using Eq. (23) again), which yields

$$\widehat{{\rm sig}} = -\log\left[ 1 - \left( 1 - 10^{-{\rm sig}}\right) ^N\right]$$ (27)

for the spectral significance of the maximum out of a sample of N statistically independent peaks in the significance spectrum.

Solving Eq. (27) for ${\rm sig}$ yields

$${\rm sig} = -\log\left( 1 - \sqrt[N]{1 - 10^{-\widehat{{\rm sig}}}}\right) ,$$ (28)

which allows one to immediately convert a chosen threshold for maximum spectral significance into "individual'' spectral significance (as given by Eq. (24)). In most practical applications, the approximation

$${\rm sig}\approx \widehat{{\rm sig}} +\log N$$ (29)

is sufficiently accurate.

For an equidistant time series consisting of K data points, the number of statistically independent DFT amplitudes is $\frac{K}{2}$, if K is an even number. One may set $N := \frac{K}{2}$ as a rough estimate also for the non-equidistant case, which performs quite reliably in practical applications.

For example, if a maximum spectral significance threshold of 5.46 shall be applied to a time series consisting of 1000 data points, the corresponding "individual'' spectral significance threshold would be 2.76, which is in good agreement with the numerical results by Kuschnig et al. (1997), who obtained a significance of $\approx$3 in this case, examining 19 300 synthetic time series.

  
5.5.3 Connection between spectral significance and signal-to-noise ratio

The correspondence between spectral significance and signal-to-noise ratio is obtained through substitution of Eq. (24) by the amplitude signal-to-noise ratio $\frac{A}{\left< A\right>}$ according to Eq. (B.11). This yields

$${\rm sig}\left( A\left\vert\right.\omega ,\theta\right) =\fr... ...sin ^2\left(\theta - \theta _0\right)}{\beta _0^2}\right]\cdot$$ (30)

Given uniformly distributed arguments of the trigonometric functions, the expected values of both $\alpha _0$ and $\beta _0$ evaluate to 1. Thus an approximation for the correspondence between amplitude signal-to-noise ratio and spectral significance is obtained by

$${\rm sig}\left( A\right)\approx\frac{\pi\log e}{4}~\left(\frac{A}{\left< A\right>}\right) ^2\cdot$$ (31)

For example, an amplitude signal-to-noise ratio of 4 - as a suggested significance estimator by Breger et al. (1993) - roughly corresponds to a spectral significance of

$${\rm sig}\left( 4\left< A\right>\right)\approx 4\pi\log e\approx 5.4575 .$$ (32)

A numerical simulation for 42 597 time series, each consisting of 14 400 equidistantly sampled data points representing a single sinusoidal signal with randomly chosen amplitude, frequency, and phase, plus Gaussian noise with randomly chosen rms error was performed. The agreement with Eq. (31) for spectral significances <200 and a corresponding signal-to-noise ratio of $\approx$25 is excellent. At higher spectral significances, the following effect has to be taken into account:

Figure 2: Spectral significance (sig) associated with an amplitude signal-to-noise ratio of 4, without ( dashed line) and with ( solid line) prewhitening of the considered peak, depending on the number of time series data points. If the amplitude noise is calculated after prewhitening, the spectral significance associated with a signal-to-noise ratio of 4 decreases below the equivalent limit of 5.46 for short datasets.

Open with DEXTER

The spectral significance is related to a peak generated by noise with the same variance as the observable under consideration, i.e. without prewhitening. Frequently, a signal-to-noise ratio calculation relies on the rms residual, i.e. with the peak under consideration prewhitened. Since (on average) a sinusoidal signal with amplitude A contributes $\frac{A^2}{2}$ to the variance in the time domain, the population variance of residuals in the time domain evaluates to

$$\left< r^2\right> \approx \left< x^2\right> - \frac{A^2}{2}\cdot$$ (33)

Using this relation, Eq. (B.11) permits one to calculate an amplitude noise level with the considered peak prewhitened, according to

$$N\left( A\right) \approx \sqrt{\frac{\pi}{K}\left(\left< x^2\right> - \frac{A^2}{2}\right)}\: ,$$ (34)

and writing Eq. (31) in terms of $N\left( A\right)$ leads to

$${\rm sig}\left( A\right)\approx\frac{\log e}{2}\frac{K\pi}{\pi +2K\left[\frac{A}{N\left( A\right)}\right] ^{-2}}\cdot$$ (35)

Figure 2 shows the difference in spectral significance for amplitude noise with and without prewhitening as a function of time series length. This difference increases dramatically as datasets become shorter. When relying on signal-to-noise ratio estimation, the issue of prewhitening is crucial and in some way philosophical. If the considered peak is believed to be "true'' a priori, then it will have to be prewhitened for the noise calculation. If it is assumed to be an artifact, noise will have to be calculated without prewhitening. However, the answer in terms of spectral significance is unique and clear: since the spectral significance analytically compares the considered time series with a randomized one, i.e. keeping the time-domain rms deviation the same, the correspondence to a non-prewhitened signal-to-noise ratio would be correct - if desired at all. The computation of the spectral significance allows one to completely omit the potentially ambiguous computation of noise levels by averaging amplitude over a frequency interval about a considered peak. In general, the latter are not objective, since the resulting noise level depends on the choice of the interval width and whether unresolved peaks are encountered.

  
5.5.4 Spectral significance and Lomb-Scargle Periodogram

If the influence of the zero-mean correction on the statistical characteristics of the Fourier Coefficients is neglected, Eqs. (17), (18) simplify to

  

$$\alpha _0\left(\omega ,\theta _0\right) := \sqrt{\frac{2}{K^2}\left[ K\sum_{k=0}^{K-1}\cos ^2\left(\omega t_k - \theta _0\right)\right]} ,$$ (36)

$$\beta _0\left(\omega ,\theta _0\right) := \sqrt{\frac{2}{K^2}\left[ K\sum_{k=0}^{K-1}\sin ^2\left(\omega t_k - \theta _0\right)\right]} ,$$ (37)

respectively. In this case, the evaluation of $\theta _0$ - as performed in Appendix A - transforms Eq. (10) into

$$\tan 2\theta _0\left(\omega\right) = \frac{\sum_{k=0}^{K-1} \sin 2\omega t_k}{\sum_{k=0}^{K-1} \cos 2\omega t_k} ,$$ (38)

and Eq. (24) becomes fully compatible to the definition of the Lomb-Scargle Periodogram. In this context, the improvement of the spectral significance compared to the Lomb-Scargle Periodogram is the correct statistical implementation of the artificially fixed time series mean. Remembering that the initial condition that led to the Lomb-Scargle Periodogram was a least-squares solution (unfortunately without handling the zero-mean correction appropriately), the significance spectrum will satisfy the corresponding least-squares condition for the zero-mean corrected data as well. This is in perfect agreement with the results of simulations, as presented in 6.2.

Since the Fourier-Space effect of the zero-mean correction tends to vanish for infinitely long time series, one would expect the difference between the Lomb-Scargle Periodogram and the spectral significance to become small (or even negligible) for long datasets. The performance of DFT, Lomb-Scargle Periodogram, and spectral significance is compared for time series of different length in 6.4.

On the other hand, the spectral significance for the Floating-Mean Periodogram (Cumming et al. 1999), the statistic of which appears not to suffer from zero-mean correction problems, is directly obtained by using Eq. (24) with $\alpha _0$, $\beta _0$ and $\theta _0$ as given above (Eqs. (36), (37), (38)).

  
5.6 The Sock Diagram

In 5.5.1, the spectral significance was introduced as a representation of the statistical properties of time-domain sampling in Fourier Space, applied to the amplitude spectrum of a given observable (Eq. (24)). Selecting a single frequency and phase angle, one finds the spectral significance to be proportional to the squared amplitude. This property makes it easy to determine an analogy to the spectral window.

For classical DFT-based methods, the spectral window is frequently used to determine the effects of time series sampling in the frequency domain. It is defined as the DFT amplitude spectrum of a constant in the time domain, normalized to an amplitude of 1 at zero frequency. In the case of non-equidistant sampling, peaks in the amplitude window indicate periodicities in the sampling of the time series corresponding to the frequencies where these peaks occur. A frequently returned signature in the spectral window of astronomical single-site measurements is a set of peaks at integer multiples of 1 d^-1. This is the Fourier-Space representation of periodic data gaps due to daylight and termed "1 d^-1 aliasing''.

Since the spectral significance is a more subtle and sensitive quantity taking into account more information on the time-domain sampling, it is possible to introduce a more sensitive analogy to the spectral window, the Sock Diagram. As with all formalism in terms of spectral significance, it is frequency- and phase-resolved.

In analogy to the spectral window, the Sock Diagram represents the spectral significance variations with frequency and phase for a constant amplitude, or amplitude signal-to-noise ratio. It provides quantitative information on the influence by the time-domain sampling on the spectral significance. Furthermore, it displays the quality of signal-to-noise ratio-based estimation for all possible frequencies and phases at a glance. Providing both information on gaps in the sampling and frequency regions with poor accuracy of the DFT amplitude, the Sock Diagram shows where DFT-based signal-to-noise ratio estimation fails.

The normalization of the Sock Function,

$${\rm sock}\left(\omega ,\theta\right) := \left[\frac{\cos ^2... ...c{\sin ^2\left(\theta - \theta _0\right)}{\beta _0^2}\right] ,$$ (39)

provides an expected value 1 on the assumption of uniformly distributed arguments of the trigonometric functions. Thus

$${\rm sig}\left( A\left\vert\right.\omega ,\theta\right) = \fr... ...}{4\left< x^2\right>}\:{\rm sock}\left(\omega ,\theta\right) ,$$ (40)

as obtained from Eq. (24), permits to compute the spectral significance associated with an amplitude A based on the Sock Diagram. In terms of signal-to-noise ratio, Eq. (30) evaluates to

$${\rm sig}\left( A\left\vert\right.\omega ,\theta\right) = \fr... ...t< A\right>}\right) ^2~{\rm sock}\left(\omega ,\theta\right) .$$ (41)

Since the orientation $\theta _0$ of the elliptical PDF in Fourier Space appears only in the form $\theta - \theta _0$ in all equations related to spectral significance, phase angles in Figs. 3 and 4 consistently refer to the position of the semi-minor axis of the rms error ellipse to achieve better visibility: using $\theta - \theta _0$ instead of $\theta$ provides a constant alignment of the spectral significance maxima in phase for all frequencies.

Figure 3: Sock Diagram (in cylindrical coordinates) for the V measurements of IC 4996 # 89, displaying the relative variations of the spectral significance (radial coordinate) with frequency (height coordinate) and phase (azimuthal coordinate) for a constant signal-to-noise ratio, and hence providing an overview of the effect of time-domain sampling properties in Fourier Space. Better visibility is achieved by additional color coding, referring to the color bar in the lower panel. The susceptibility of DFT spectra to 1 d-1 aliasing shows up in spectral significance variations of up to a factor $\approx$35.

Open with DEXTER

Figure 4: Sock Diagram for MOST ("Microvariability and Oscillations of STars''; Walker et al. 2003) measurements of $\zeta$ Oph. Due to practically equidistant sampling (duty cycle $\approx$$99.9\%$, the relative variations of the spectral significance with frequency and phase are $\approx$10-3, even close to the orbital period of the satellite ($\equiv$14.2 d-1). As illustrated by the color bar in the lower panel, the color coding scale differs from Fig. 3 considerably.

Open with DEXTER

Figure 3 displays the Sock Diagram for typical non-equidistant sampling representing 10 nights yielding 381 data points of single-site V photometry of star # 89 in the young open cluster IC 4996 (Zwintz et al. 2004; Zwintz & Weiss 2006). Regardless of the astrophysical dimension, these observations have been chosen as the primary test dataset, because they impressively show all the characteristics typical for single-site measurements that make multifrequency analysis a puzzling task.

The Sock Diagram uses three-dimensional polar coordinates: for each frequency, the angular coordinate refers to phase, and the radial component is associated with the spectral significance normalized to an expected value of 1 (according to Eq. (39)). To enhance the visibility, the radial information is color-coded additionally.

The spectral significance variations close to 0.5, 1, 1.5, and 2 d^-1 clearly indicate frequencies where signal-to-noise ratio estimation of False-Alarm Probability is potentially misleading.

An example for excellent sampling is shown in Fig. 4, representing 72 055 MOST data points of $\zeta$ Oph Fabry Imaging photometry (Walker et al. 2004, 2005) with a duty cycle of $\approx$$99.9\%$, obtained between May 17 and June 14, 2004. The data reduction is performed according to the procedure described by Reegen et al. (2006). Even for frequencies close to 14.2 d^-1 - corresponding to the orbital period of the satellite (101.4 min; Walker et al. 2003) - the spectral significance varies with frequency and phase by only $\approx$$0.1\%$.

  
5.7 The marginal distribution of phase angles

An alternative method to the computation of spectral significance as a function of amplitude and phase simultaneously would be to use the phase information separately in addition to, e.g., some traditional signal-to-noise ratio-based reliability estimate. The idea is that for white noise in the time domain, the Fourier phases are not uniformly distributed and that the additional incorporation of the phase information may provide a more detailed criterion on the reliability of a peak.

The overall distribution of phases at a given frequency is immediately obtained by integrating the bivariate PDF (Eq. (14)) over amplitude, which yields

$$\phi\left(\theta\right) = \frac{K\alpha _0~\beta _0}{\beta _... ...0\right) - \alpha _0^2\sin ^2\left(\theta -\theta _0\right)} ,$$ (42)

normalized for phases $\theta$ on an interval of width $\pi$. This avoids the ambiguity of solutions for $\theta _0$ modulo $\pi$, as discussed in 5.3.

The expected value of this probability distribution is $\theta _0 - \frac{\pi}{2}$: phases associated to a low spectral significance (for given amplitude) generally occur more frequently than phases for which the spectral significance is high. The special case $\alpha _0 = \beta _0$ yields an upper limit of $\frac{\pi}{2\sqrt{3}} \approx 52^{\circ}$ for the standard deviation obtained by Eq. (42).

In statistical terms, the phase distribution provided above is a marginal distribution of the bivariate PDF given by Eq. (14). Instead of statistically examining a bivariate distribution, one may use the two marginal distributions (in the present case, the amplitude and phase distributions) instead, but encountering a loss of accuracy, since correlations between the two random variables remain unresolved. Therefore, from the theoretical point of view, it is advisable to use the bivariate form rather than the marginal distributions. An additional problem is that there is no analytical solution for the integral over phase, which would transform the bivariate PDF into the marginal distribution of amplitudes. One would have to employ classical techniques relying on the amplitude signal-to-noise ratio instead.

In addition to these theoretical objections to the examination of marginal distributions, there is a major practical constraint: Fig. 5 displays the expected values and standard deviations of Fourier phase for frequencies between 1.3 and 2.7 d^-1, given the sampling of IC 4996 # 89. In the considered frequency range, the rms scatter of phases is greater than  $44^{\circ}$ for the entire frequency range under consideration. Compared to the upper limit of  $52^{\circ}$, this scatter will probably be much too high to reveal additional information on the reliability of a peak, if the marginal distribution of phases is considered.

The PDF of phases provided by Eq. (42) are in perfect agreement to the results of numerical simulations.

	Figure 5: Expected values ( solid lines) and rms errors ( gray-shaded areas) of the probability distribution of phase angles in Fourier Space for white noise on the IC 4996 # 89 sampling.
Open with DEXTER

  
5.8 Spectral significance-based signal recovery

Practical applications are frequently based on a cascade of consecutive prewhitenings. In this case, the frequency at maximum spectral significance, $\hat{\omega}$, is considered as that of the strongest signal component, which is understood to be the next candidate for prewhitening. The next step is to determine the best fit for amplitude and phase at $\hat{\omega}$. The least-squares condition for the best sinusoidal fit at constant frequency $\hat{\omega}$ is obtained by

$$\frac{\partial}{\partial\hat{A}}\sum_{k=0}^{K-1}\left\lbrace x_k ... ...left(\hat{\omega} t_l-\hat{\theta}\right)\right>\right]\right\rbrace ^2 = 0\: ,$$ (43)

$$\frac{\partial}{\partial~\hat{\theta}}\sum_{k=0}^{K-1}\left\lbrac... ...left(\hat{\omega} t_l-\hat{\theta}\right)\right>\right]\right\rbrace ^2 = 0\: ,$$ (44)

where the term $\left<\cos\left(\hat{\omega} t_l-\hat{\theta}\right)\right> = \frac{1}{K}\sum_{l=0}^{K-1}\cos\left(\hat{\omega} t_l-\hat{\theta}\right)$ takes into account that the discrete fit has to be zero-mean corrected, if applied to zero-mean corrected data x_k. The derivatives lead to

 

$$\sum_{k=0}^{K-1} x_k\cos\left(\hat{\omega} t_k -\hat{\theta}\righ... ...}\Bigg\lbrace\sum_{k=0}^{K-1}\cos ^2\left(\hat{\omega} t_k -\hat{\theta}\right)$$

$$- \frac{1}{K}\left[\cos\left(\hat{\omega} t_k -\hat{\theta}\right)\right] ^2\Bigg\rbrace ,$$ (45)

$$\sum_{k=0}^{K-1} x_k\sin\left(\hat{\omega} t_k -\hat{\theta}\righ... ...\omega} t_k -\hat{\theta}\right)\sin\left(\hat{\omega} t_k -\hat{\theta}\right)$$

$$- \frac{1}{K}\left[\cos\left(\hat{\omega} t_k -\hat{\theta}\right)\right]\left[\sin\left(\hat{\omega} t_k -\hat{\theta}\right)\right]\Bigg\rbrace .$$ (46)

Solving for $\hat{\theta}$ yields

$$\sum_{k=0}^{K-1}\sum_{l=0}^{K-1}\sum_{m=0}^{K-1}\left[\sin\ha... ...right] x_k\cos\left(\hat{\omega} t_l-\hat{\theta}\right) = 0 ,$$ (47)

which finally reduces to

$$\tan\hat{\theta} = \frac{P_1\sum _{k=0}^{K-1}\cos\hat{\omega... ...in\hat{\omega}t_k - P_3\sum _{k=0}^{K-1}\cos\hat{\omega}t_k} ,$$ (48)

using

$$K\sum _{k=0}^{K-1}\cos\hat{\omega} t_k\sin\hat{\omega} t_k -\left... ...\cos\hat{\omega} t_k\right)\left(\sum _{k=0}^{K-1}\sin\hat{\omega} t_k\right) ,$$                                      P₁ := (49)

$$K\sum _{k=0}^{K-1}\cos ^2\hat{\omega} t_k - \left(\sum _{k=0}^{K-1}\cos\hat{\omega} t_k\right) ^2 ,$$ P₂ := (50)

$$K\sum _{k=0}^{K-1}\sin ^2\hat{\omega} t_k - \left(\sum _{k=0}^{K-1}\sin\hat{\omega} t_k\right) ^2 .$$ P₃ := (51)

Equation (48) provides two solutions for $\hat{\theta}$. To pick the least-squares related solution is an easy task for a program.

Once frequency $\hat{\omega}$ and phase $\hat{\theta}$ of the signal are evaluated, Eq. (45) immediately yields the best-fitting amplitude,

$$\hat{A} = \frac{\sum_{k=0}^{K-1} x_k\cos\left(\hat{\omega} t_... ...\cos\left(\hat{\omega} t_k-\hat{\theta}\right)\right] ^2}\cdot$$ (52)

  
6 Numerical tests

Two sets of numerical simulations have been performed, the first one to confirm the agreement between the theoretically evaluated spectral significance and a straight-forward histogram analysis, and the second one to quantitatively compare the accuracy of frequencies returned by the various methods discussed in this paper.

   
6.1 Comparison of analytical and numerical solutions

Extensive numerical simulations have been performed in order to examine the validity of the above theoretical considerations. The simulations are set up in a way that compares closely to real life. The algorithm to generate Gaussian noise is based on the Central Limit Theorem (Stuart & Ord 1994, p. 310f), providing fast convergence of a mean value of uniformly distributed random variables towards a Gaussian distribution with increasing number of variables. All of the following numerical applications rely on Gaussian noise produced by summing up 10 values from the random number generator. A comprehensive compilation on alternative methods to generate Gaussian noise is provided by Firneis (1970).

An all-in-one simulation for many frequencies, many phases, and many amplitudes would result in tremendously time-consuming computations. In order to keep the effort reasonable, only a single frequency is picked to examine the phase dependency of the spectral significance for different amplitudes.

Figure 6 displays the spectral significance associated with an amplitude signal-to-noise ratio of 4 - according to Eq. (30) - for the V data of IC 4996 # 89. The blue and red graphs represent two orthogonal phases in Fourier Space. The comparison of numerical and theoretical results is performed for a frequency of 1.956 d^-1, where the deviation of spectral significances from the expected value 5.46 is $\approx$5 for selected phase angles. Since it is desired to examine the presence of spectral significance variations with phase in numerical results where predicted by theory, this frequency is a reasonable choice for a test.

Figure 6: Spectral significance associated with an amplitude signal-to-noise ratio of 4 for the V measurements of IC 4996 # 89. The blue and red graphs refer to orthogonal phases in Fourier Space. The expected spectral significance of $\approx$5.46 is displayed by the dashed black line. The vertical dashed-dotted green line indicates the frequency of 1.956 d-1, which was selected for a numerical simulation to check the validity of the theoretical approach.

Open with DEXTER

The procedure consists of five steps.

1.

Zero-mean Gaussian noise with a standard deviation of 1 is imposed upon the IC 4996 # 89 sampling.

2.

The zero-mean correction is performed to avoid scatter in the mean of the finite time series about the population mean (see 5.1).

3.

The Fourier-Space phase angle for the synthetic time series at 1.956 d^-1 is evaluated.

4.

The Fourier-Space noise is evaluated upon the variance of the time series according to Eq. (B.11).

5.

If the Fourier Amplitude exceeds the signal-to-noise ratio under consideration, the resulting histogram is updated correspondingly.

The number of amplitudes exceeding the preselected limit relative to the total number of synthetic datasets provides an estimator for the False-Alarm Probability (and consequently spectral significance) associated with the chosen signal-to-noise ratio.

Figures 7 to 9 display the agreement between theory and simulations for amplitude signal-to-noise ratios of 1, 2, and 3. The corresponding numbers of synthetic datasets computed are 1.5 million, 22 million, and 250 million, correspondingly. In all three cases, phase bins of width $1^{\circ }$ were used for the histograms. For a signal-to-noise ratio of 4, the number of synthetic datasets required for an acceptable numerical accuracy would exceed the capabilities of computational performance, especially those of system random number sequences, by far.

Figure 7: Phase-dependent spectral significance (radial coordinate, referred to by the vertical axis) for the V measurements of IC 4996 # 89 at a constant frequency of 1.956 d-1 as a function of phase, referring to an amplitude signal-to-noise ratio of 1. The solid line represents the theoretical result. A numerical simulation for 1.5 million synthetic datasets ( dots) - counted in phase bins of width $1^{\circ }$ - illustrates the excellent agreement. The systematically higher deviation of the numerical results from the theoretical solutions at higher spectral significances is due to sample-size effects: phase bins associated with a high theoretical spectral significance are hit less often than others. Hence the total number of Fourier Amplitudes to be examined differs from phase bin to phase bin. The dashed line represents the expected spectral significance of $\approx$0.34 associated with a signal-to-noise ratio of 1 for uniformly distributed arguments of all trigonometric functions.

Open with DEXTER

	Figure 8: Same as Fig. 7, but for an amplitude signal-to-noise ratio of 2. The dots refer to a numerical simulation for 22 million synthetic datasets. The expected spectral significance level is $\approx$1.36 ( dashed line).
Open with DEXTER

	Figure 9: Same as Fig. 7, but for an amplitude signal-to-noise ratio of 3. The dots refer to a numerical simulation for 250 million synthetic datasets. The expected spectral significance level is $\approx$3.07 ( dashed line).
Open with DEXTER

Phases associated with a high spectral significance occur less frequently than others (see 5.7), whence the scatter of numerical results is systematically higher at higher spectral significance levels due to sample-size effects. Taking this into account, the overall quality of fits is good. In neither of the three plots, a systematic deviation of the numerical results from the analytical functions is visible. However, the theoretical prediction for the overall shape is recovered by the numerical results, indicating that the ellipticity given by the normalized axes $\alpha _0$ and $\beta _0$ matches the "reality'' of simulation. Also the orientation of the ellipse ($\theta _0$) appears consistent with the synthetic data. Finally, the consecutive comparison of Figs. 7 to 9 offer convincing evidence of the spectral significance at a constant frequency and phase angle to be proportional to the squared amplitude, as predicted by the theoretical solution.

  
6.2 Accuracy of peak frequencies

One of the major issues in period searches is the accuracy of the resulting signal frequencies. Classical DFT with or without improvement of peak frequencies by a subsequent least-squares fit, Lomb-Scargle Periodogram, and spectral significance analysis are based on the identification of the highest peak in a chosen frequency interval.

A quantitative description of the quality of resulting frequencies with and without noise was obtained by means of simulations. For the sampling of the IC 4996 # 89 dataset (V), the following procedure is performed:

1.

a synthetic signal of given frequency f₀ and amplitude with random phase (uniformly distributed on $\left[ -\pi, \pi\right]$) is generated;

2.

Gaussian noise with given standard deviation is added (optional);

3.

zero-mean correction of the resulting dataset is performed;

4.

a DFT amplitude spectrum, a Lomb-Scargle Periodogram, a spectrum of phase dispersion, and a significance spectrum are computed for a predefined frequency range;

5.

an interpolation routine is performed to find the maximum (or minimum phase dispersion, respectively) in each of the four spectra; and

6.

the deviation of the resulting frequency from the corresponding input frequency, $\Delta f$, is evaluated.

For each frequency, deviations $\Delta f$ were collected to obtain a frequency-dependent rms error of recovered frequencies, where only attempts with $\left\vert\Delta f\right\vert \le 0.5$ d^-1 were taken into account. Attempts resulting in $\left\vert\Delta f\right\vert > 0.5$ d^-1 were considered alias.

The example illustrated by Figs. 10, 11 uses the IC 4996 # 89 data.

Figure 10: Frequency accuracy (rms scatter of resulting frequencies about the initial frequency for a single sinusoidal signal with uniformly distributed phase) vs. signal frequency of five methods: DFT ( solid orange line), DFT plus least-squares fitting ( solid blue line), Lomb-Scargle Periodogram ( solid green line), PDM ( dotted black line), and spectral significance ( solid red line). The time-domain sampling represents the V measurements of IC 4996 # 89. The top panel is the frequency accuracy for a pure signal. Towards the bottom panel, Gaussian noise with increasing standard deviation is added to the signal in two steps, corresponding to amplitude signal-to-noise ratios of 40 and 4, respectively. Only those attempts where the distance between resulting frequency and input frequency does not exceed 0.5 d-1 were taken into account, the rest was considered alias (see 6.3). The results are based on a numerical simulation investigating 1000 datasets for every frequency.

Open with DEXTER

Figure 11: Relative number of aliases vs. signal frequency of four methods: DFT plus least-squares fitting ( dashed blue line), Lomb-Scargle Periodogram ( solid green line), PDM ( dotted black line), and spectral significance ( solid red line). In this context, a result is considered alias, if the absolute difference between resulting frequency and signal frequency exceeds 0.5 d-1. The DFT without least-squares fitting is not displayed here. It produces essentially the same fraction of alias, because the least-squares fit allows a fine tuning of resulting peak frequencies only, which is negligible considering frequency errors of 0.5 d-1 and more. The time-domain sampling represents the V measurements of IC 4996 # 89. The top panel is the fraction of alias peaks for a pure signal. Towards the bottom panel, Gaussian noise with increasing standard deviation is added to the signal in two steps, corresponding to amplitude signal-to-noise ratios of 40 and 4, respectively. In the upper two panels, the DFT graph is offset vertically by -0.01 and the Lomb-Scargle graph by +0.01 to provide better visibility at and close to zero. The results are based on a numerical simulation investigating 1000 datasets for every frequency.

Open with DEXTER

Figure 10 shows the comparison of the five methods. For each frequency, 1000 datasets were examined. The figure displays the rms deviation (denoted $\left\vert\Delta f\right\vert$) of resulting frequencies as a function of the signal frequency f₀. The top panel refers to a clean sinusoidal signal without noise. Towards the bottom panel, Gaussian noise with increasing standard deviation is added (corresponding to signal-to-noise ratios of 40 and 4, respectively, using Eq. (35)).

For a signal without noise (top panel), the frequency accuracies of the Lomb-Scargle Periodogram exceed that of the DFT by a factor $\approx$10, and both profiles show accuracy variations with frequency. The Lomb-Scargle Periodogram does not take into account zero-mean correction, as shown in 5.5.4. This leads to systematic effects not only for short datasets, but also in frequency regions where the phase coverage becomes poor. In the absence of noise, frequencies at maximum spectral significance are $\approx$100 000 times more accurate than the DFT results, which is in the accuracy domain of the least-squares solutions. Apparently the accuracy of the spectral significance peaks is only limited by the internal accuracy of the computer for double-precision floating-point numbers. In addition, the accuracy of spectral significance is, in this respect, practically independent of frequency.

Once noise is added to the signal, the peak frequency accuracy of spectral significance becomes much poorer (in consistency with O'Donoghue & Montgomery 1999), but as illustrated by Fig. 10, the method does not get worse than either alternative procedure. The improvement of the frequency accuracy by the spectral significance solution for high signal-to-noise ratio is valuable, since the exact knowledge of frequency provides exact prewhitening. Only exact prewhitening guarantees that weaker frequency components are not contaminated by spurious residuals of the higher peaks.

  
6.3 Aliasing

One potential weakness of the DFT, which usually becomes puzzling in the analysis of non-equidistantly sampled time series, is the susceptibility to aliasing. Whereas measurements over infinite time would lead to perfectly delta-shaped signal components, the spectral window of "real'' data is convolved into this "ideal'' spectrum. For typical single-site observations, 1 d^-1 sidelobes are visible around every signal peak, and in many cases it is hard to decide which peak out of the "comb'' of aliases is the right one. For multiperiodic signal, aliases of individual components may interfere. If this interference leads to an amplification, an erroneous component identification is inevitable, when only relying on the highest amplitude. This applies to the spectral significance as well and reflects a major weakness of the step-by-step prewhitening technique rather than the calculation of the spectral significance. Strategies to overcome the aliasing problem have been examined for many years (e.g. Ferraz-Mello 1981; Roberts et al. 1987; Foster 1995).

Figure 11 displays a comparison of DFT, Lomb-Scargle Periodogram, PDM, and spectral significance in terms of fraction of aliases among the 1000 simulated test datasets used in 6.2. As in Fig. 10, the top panel refers to signal without noise, and Gaussian noise of increasing standard deviation is added towards the bottom panel. Not surprisingly, the susceptibility of the DFT to aliasing does not improve, if a least-squares algorithm is appended. The portion of aliases obtained using the Lomb-Scargle Periodogram is fairly below that of the DFT. Finally, the figure clearly shows that the spectral significance analysis is more stable against aliasing than either compared strategy. Without noise, not a single alias peak occurs among altogether 300 000 simulated datasets. For a signal-to-noise ratio of 40, the percentage of alias peaks is less than a third of the corresponding percentage of the alternative methods, even at 1 d^-1. Only for the highest noise level, to which the bottom panel refers, all results are quite comparable.

Since the spectral significance is not initially intended to correct for aliasing, its capabilities to avoid potential misidentification of peaks are limited to the extent of spuriously amplified peaks by systematic inhomogeneities of the frequency-domain noise. Simulations and practical experience show that the systematic errors of spectral noise are an inferior error source compared to the interference of aliases of different signal components. Presently investigations are performed to apply the spectral significance technique to simultaneous multi-frequency solving algorithms for an improved treatment of aliases. The results are planned to be presented in Paper II.

  
6.4 Sample size effects

A further example is provided to compare the performance of DFT, Lomb-Scargle Periodogram, and spectral significance for datasets of equal (or at least comparable) characteristics, but different length.

The "long'' dataset represents 34 nights (not consecutive, but covering 81 d) of single-site photometry of the Delta Scuti star 44 Tau (Strømgren y, 2981 data points; Antoci et al. 2006) acquired by the Vienna University Automatic Photoelectric Telescope (APT; Strassmeier et al. 1997; Granzer et al. 2001). More details on the data are provided in Sect. 7. The "short'' dataset is a subset of 7 consecutive nights (619 points). For the subsequent investigations, only the sampling of these two datasets is used. Peak frequency accuracies are computed according to the procedure described in 6.2.

Figure 12 shows the frequency accuracies of the compared methods for the short dataset in the top panel. The display is exactly according to Fig. 10. The bottom panel contains the identical analysis of frequency accuracies for the long dataset and indicates an improved overall precision of the DFT and Lomb-Scargle Periodogram, if the number of data increases.

Figure 12: Frequency accuracy (rms scatter of resulting frequencies about initial frequency for a single sinusoidal signal with uniformly distributed phase) vs. signal frequency of DFT ( orange) and Lomb-Scargle Periodogram ( green). The simulated time series represent pure signal without noise at uniformly distributed phase angles. The top panel represents the frequency accuracy for seven consecutive nights of single-site photometry of 44 Tau (Strømgren y, 619 data points). These data are a subset of the 81 days long time series presented in the bottom panel (34 nights, 2981 data points). The plots illustrate the improvement of the overall frequency accuracy of both methods for the longer dataset. The corresponding deviations obtained using the spectral significance are practically zero, similar to the top panel in Fig. 10, and thus not presented here. Only those attempts where the distance between resulting frequency and input frequency does not exceed 0.5 d-1 were taken into account, the rest was considered alias. The results are based on a numerical simulation investigating 1000 datasets for every frequency.

Open with DEXTER

However, 1 d^-1 aliasing persists also for the long time series (Fig. 13), which is compatible with the persisting alias peaks in the spectral window at integer multiples of 1 d^-1. Since the zero-mean correction is represented by the subtraction of a constant from all observables in the entire time series (referring to Eqs. (3), (4)), it is not surprising that the techniques which ignore the statistical effects of zero-mean correction retain these spectral window-related weaknesses also in case of a large number of data points.

Figure 13: Same as Fig. 12, but displaying the relative number of aliases vs. the signal frequency of DFT ( solid orange) and Lomb-Scargle Periodogram ( dashed green). In this context, a result is considered alias, if the absolute difference between resulting frequency and signal frequency exceeds 0.5 d-1. The spectral significance does not produce a single alias maximum in the simulations, similar to the top panel in Fig. 11. The corresponding graphs are thus not presented here. The DFT graph is offset vertically by -0.01 and the Lomb-Scargle graph by +0.01 to provide better visibility at and close to zero.

Open with DEXTER

The result of spectral significance calculations is not provided in these graphs: the deviation of frequencies is practically equal to zero independently of frequency - as for the IC 4996 # 89 data (top panel in Fig. 10).

  
7 SigSpec: practical application

The formal concept of spectral significance does not only provide reliable information on the sampling characteristics in the time domain, but also allows to include the computations into a prewhitening sequence. S IGS PEC realizes the entire procedure as a high-performance algorithm. Implementations for various operating systems and CPUs for free download at http://www.astro.univie.ac.at/SigSpec. The ANSI-C source code is available on request (reegen@astro.univie.ac.at).

The S IGS PEC technique was applied to numerous datasets and has proven its value also in other scientific problems. A practical example for a situation where S IGS PEC performs superior to classical DFT-based signal-to-noise ratio estimation is given in Table 1. Columns 1 and 2 represent 13 identified eigenfrequencies plus 16 excited linear combinations of the $\delta$ Sct star 44 Tau (Antoci et al. 2006). The corresponding Strømgren v and y amplitudes as obtained from a multisite campaign in 2001/02 are displayed in Cols. 3 and 4. The underlying light curves consist of 3890 points in v and 3582 points in y. For the analysis of the campaign data, the signal-to-noise ratio threshold for peak acceptance was chosen to be 3.5.

Table 1: Frequencies detected in multisite data of 44 Tau, according to Antoci et al. (2006). A_v and A_y are the published amplitudes (mmag) in Strømgren v and y for the 2001/02 multisite data (3890 points in v, 3582 points in y). The corresponding columns indicated by "S/N'' contain a cross-identification for the data obtained by the Vienna University Automatic Photoelectric Telescope (Strassmeier et al. 1997; Granzer et al. 2001) in the course of this campaign (3280 points in v, 2981 points in y). An amplitude signal-to-noise ratio >4 was used. The columns indicated by "sig'' represent the cross-identification of the Vienna APT data with the complete campaign dataset using the result of the S IGS PEC analysis, based on a spectral significance limit of 5.46. Amplitudes in italic print refer to 1 d^-1 alias peaks. Since the peaks for the APT subset were assigned to corresponding peaks in the total dataset, individual frequencies are not displayed.

The major part of the data (3280 points in v, 2981 points in y) was acquired by the Vienna University Automatic Photoelectric Telescope (APT; Strassmeier et al. 1997; Granzer et al. 2001). The data represent 34 nights covering a time interval of 81 d. As a comparative test of the practical performance of S IGS PEC, the APT data alone were analysed using a DFT-based prewhitening sequence relying on a signal-to-noise ratio limit of 4. Except for the higher threshold, this is exactly the same technique as applied to the campaign data. In any case, the frequencies found in the APT subset were cross-identified with the frequencies published for the full dataset (Antoci et al. 2006). The result for the APT data is displayed in Cols. 5 and 6. In v, 20 of the 29 signal components are found, in y only 15. With each filter, 1 d^-1 aliasing is found for 2 frequencies, which are indicated by italc print.

Columns 7 and 8 refer to the corresponding S IGS PEC analysis using a spectral significance limit of 5.46 and reproducing all frequencies except for one in v and 5 in y. The number of aliases is 8 in v and 4 in y, but the majority of aliases in the S IGS PEC result occurs for frequencies not resolved by the alternative technique at all.

For additional practical examples, the reader is referred to Reegen (2005).

  
8 Further topics

  
8.1 Spectral significance for statistically weighted time series

Astronomical measurements are generally influenced by instrumental and environmental conditions changing with time. Thus different accuracies for different data points are involved, which are frequently desired to be taken into account by applying weights to the observables. If the accuracy is poor, the corresponding weight is low, and high-accuracy data points are assigned a high weight, respectively. Furthermore, multisite campaigns employ different telescopes with different instrumental parameters, which may also require an appropriate weighting.

This section refers to a set of statistical weights, $\gamma _k$, k = 0,  ...  K-1, normalized according to

$$\sum _{k=0}^{K-1}\gamma _k =: K .$$ (53)

Weighting may be considered as counting an item $\xi _k$ of, e.g., an arithmetic mean repeatedly instead of only once. If the number of counts for an arbitrary $\xi _k$ is n_k, then the arithmetic mean over all $\xi _k$, k=0,1,...,K-1 is expressed as

$$\left<\xi _k\right> = \frac{\sum _{k=0}^{K}n_k\xi _k}{\sum _{k=0}^{K}n_k}\cdot$$ (54)

In this case, Eq. (53) demands the normalized weights to be introduced according to

$$\gamma _k := \frac{K\xi _k}{\sum _{k=0}^{K}n_k}\cdot$$ (55)

These arguments may consistently be followed to obtain the subsequent relations (Eqs. (56) to (62)).

The implementation of statistical weights has to be performed for all elements of the spectral significance relation (Eq. (24), or (25), respectively).

1.

The weighted mean value to be subtracted from every data point is calculated according to

$$\left< x_k\right> = \frac{1}{K}\sum _{k=0}^{K-1}\gamma _k x_k .$$ (56)

Further considerations assume x_k to be already zero-mean corrected.

2.

Also the weighted mean variance of the data,

$$\left< x_k^2\right> = \frac{1}{K}\sum _{k=0}^{K-1}\gamma _k x_k^2 ,$$ (57)

can be used as an estimator for the variance of the corresponding population.

3.

The DFT is determined by the weighted Fourier Coefficients

$$a_{{\rm ZM}}\left(\omega _n\right) := \frac{1}{K}\sum_{k=0}^{K-1}\gamma _k x_k\cos\omega _nt_k ,$$ (58)

$$b_{{\rm ZM}}\left(\omega _n\right) := \frac{1}{K}\sum_{k=0}^{K-1}\gamma _k x_k\sin\omega _nt_k .$$ (59)

4.

The parameters characterizing the sampling properties may be generalized according to

   

$$\alpha _0^2\left(\omega ,\theta _0\right) = \frac{2}{K^2}$$

$$\left\lbrace K\sum_{k=0}^{K-1}\gamma _k\cos ^2\left(\omega t_k - ... ...^{K-1}\gamma _k\cos\left(\omega t_l - \theta _0\right)\right] ^2\right\rbrace ,$$ (60)

$$\beta _0^2\left(\omega ,\theta _0\right) = \frac{2}{K^2}$$

$$\left\lbrace K\sum_{k=0}^{K-1}\gamma _k\sin ^2\left(\omega t_k - ... ...^{K-1}\gamma _k\sin\left(\omega t_l - \theta _0\right)\right] ^2\right\rbrace ,$$ (61)

$$\tan 2\theta _0\left(\omega\right) =$$

$$\left[ K\sum_{k=0}^{K-1}\gamma _k\sin 2\omega t_k - 2\sum_{k=0}^{K-1}\gamma _k\cos\omega t_k\sum_{k=0}^{K-1}\gamma _k\sin\omega t_k\right]\:$$

$$\times\left[ K\sum_{k=0}^{K-1}\gamma _k\cos 2\omega t_k - \left(\sum_{k=0}^{K-1}\gamma _k\cos\omega t_k\right) ^2\right.$$

$$+ \left.\left(\sum_{k=0}^{K-1}\gamma _k\sin\omega t_k\right) ^2\right] ^{-1} .$$ (62)

These parameters provide Eq. (25) to return the weighted significance spectrum of the dataset.

In the case of $\gamma _k =: 1$ $\forall k$, these relations consistently reduce to Eqs. (10), (17), (18) for unweighted analysis, respectively. Given the above re-definitions, Eq. (24) and the corresponding considerations - conversion between signal-to-noise ratio and signficance (5.5.3) and Sock Diagram (5.6) - readily apply to weighted time series as well.

  
8.2 Spectral significance for zero-mean corrected subsets

In some cases one may demand to correct the mean magnitude for individual subsets instead of the entire time series, i.e. if the data were obtained by different instruments, or if the subtraction of nightly averages is performed. In this case the influence on the statistical description of the noise to which an amplitude level is compared has to reflect the subdivision correspondingly. Given a subset index, s = 1,2,...,S, where S denotes the total number of subsets, Eqs. (3), (4) have to be re-written according to

$$a_{{\rm ZM}}\left(\omega\right) = \frac{1}{K}\left[\sum_{s=1}^{S}... ...K_s}\sum_{k=0}^{K_s-1}x_{sk}\sum_{l=0}^{K_s-1}\cos\omega t_{sl}\right)\right] ,$$ (63)

$$b_{{\rm ZM}}\left(\omega\right) = \frac{1}{K}\left[\sum_{s=1}^{S}... ...K_s}\sum_{k=0}^{K_s-1}x_{sk}\sum_{l=0}^{K_s-1}\sin\omega t_{sl}\right)\right] .$$ (64)

In this context, each registration time and observable is assigned two indices, writing x_sk. The first index refers to the subset, and the second index is the identifier of the data point within the subset. The number of data points within the subset s is denoted by K_s, and finally $K := \sum_{s=1}^{S}K_s$ for the total number of data points in the entire time series.

Further calculations follow the procedure described in Sect. 5 exactly and yield

$$\tan 2\theta _0\left(\omega\right) = \left[\sum_{s=1}^{S}K_s\sum_... ...-\sum _{k=0}^{K_s-1}\cos\omega t_{sk}\sum_{k=0}^{K_s-1}\sin\omega t_{sk}\right]$$

$$\times\left[\sum_{s=1}^{S}K_s\sum_{k=0}^{K_s-1}\cos 2\omega t_{sk... ...{sk}\right) ^2+\left(\sum_{k=0}^{K_s-1}\sin\omega t_{sk}\right) ^2\right] ^{-1}$$ (65)

for the orientation of the rms error ellipse, in analogy to Eq. (10). Correspondingly, the normalized axes of the ellipse (Eqs. (17), (18)) evaluate to

$$\alpha _0\left(\omega ,\theta _0\right) \left(\frac{2}{K}\sum_{s=1}^{S}\left\lbrace\sum_{k=0}^{K_s-1}\cos ^2\left(\omega t_{sk} - \theta _0\right)\right.\right.$$ =

$$- \left.\left.\frac{1}{K_s}\left[\sum_{l=0}^{K_s-1}\cos\left(\omega t_{sl} - \theta _0\right)\right] ^2\right\rbrace\right) ^{\frac{1}{2}} ,$$ (66)

$$\beta _0\left(\omega ,\theta _0\right) \left(\frac{2}{K}\sum_{s=1}^{S}\left\lbrace\sum_{k=0}^{K_s-1}\sin ^2\left(\omega t_{sk} - \theta _0\right)\right.\right.$$ =

$$- \left.\left.\frac{1}{K_s}\left[\sum_{l=0}^{K_s-1}\sin\left(\omega t_{sl} - \theta _0\right)\right] ^2\right\rbrace\right) ^{\frac{1}{2}} .$$ (67)

  
8.3 Colored noise

S IGS PEC does not take into account colored noise. Since there are both instrumental effects (e.g. CCD readout, stability of the spacecraft position) and stellar variations (e.g. low-amplitude modes, granulation noise) to be considered and neither of these two sources may be determined unambiguously, it is presently impossible to deduce a reliable amplitude noise profile for measurements in many cases.

The heuristic approach to generate a noise spectrum by means of (weighted) moving averages suffers from the presence of unresolved peaks, which increases the risk to miss intrinsic signal. Hence it is advisable to use S IGS PEC for the detection of significant sinusoidal signals and to take effects due to colored noise (however the latter may have been determined) into account by choosing different spectral significance thresholds for different frequency regions.

However, there is evident demand for further investigation in the presence of colored noise. The statistical description of colored noise and its implementation into the evaluation of the spectral significance is the subject of present and ongoing investigation and will be discussed in a dedicated publication (Paper III).

  
9 Conclusions

S IGS PEC exceeds the diagnostic capabilities of Discrete Fourier Transform, Lomb-Scargle Periodogram, Phase Dispersion Minimization, and least-squares fits in various respects.

1.

S IGS PEC does not compare one peak amplitude to another. Thus it avoids implications to the statistical independence of peaks in the amplitude spectrum at a single prewhitening stage and the ambiguity in defining a Nyquist Frequency, which all compared methods suffer from. The implicit comparison of peaks by selecting the dominant one for prewhitening may be performed free of the frequency and phase dependencies of Fourier amplitudes, if the spectral significance is used instead (6.1). One remaining issue in this context is the effect of multiple hypothesis testing on step-by-step identification and prewhitening. The propagation of uncertainties from one prewhitening step to another will be a topic in Paper II examining the DFT on multi-periodic signals.

2.

S IGS PEC is based on an analytically straight determination of the amplitude Probability Density Function (5.3). Provided the noise is white (but not necessarily Gaussian), S IGS PEC returns exact results instead of signal-to-noise ratio estimates (6.1). A more realistic approach is the introduction of colored noise in the sense of serially correlated measurements. Paper III is planned to discuss the appropriate incorporation of the serial correlation into the spectral significance computation.

3.

S IGS PEC is the first technique in astronomical time series analysis to use both frequency and phase angle in the computation of the False-Alarm Probability (Sect. 5), appropriately encountering the fact that the amplitude noise level is systematically different for different frequencies and phases in Fourier Space (6.2). The statistics employed by S IGS PEC takes into account all available information provided by the Discrete Fourier Transform.

4.

The performance of S IGS PEC for a single signal plus noise is discussed in this paper (6.2, 6.3). For strong signal components, the peak frequencies returned by S IGS PEC are considerably more accurate than the solutions of all compared methods. For weak signal components, S IGS PEC still provides a slightly higher accuracy. Furthermore, there is indication that the S IGS PEC frequencies are at least as accurate as the least-squares fits. The systematic examination of multiperiodic signals and a comparison to non-linear least-squares fits will be provided in Paper II.

5.

Signal-to-noise ratio estimation provides a valid significance estimate, but becomes very poor for low frequencies and in case of periodicities in the time-domain sampling (6.4).

6.

Signal-to-noise ratio-based period detection is a matter of experience and - to some extent - personal taste. S IGS PEC provides a solution that is completely free of subjective influence. Interpretation is no longer a part of the analysis, and the results remain exactly the same, if two different persons apply the calculations to the same data.

Tests on synthetic and real data confirm that the S IGS PEC analysis is superior to signal-to-noise ratio estimates (Sect. 6).

The Sock Diagram (5.6), as a generalized analogy to the spectral window, provides increased information on the properties of the time-domain sampling in Fourier Space, since incorporating phase-resolution as well. Periodicities in the sampling of a time series are recovered by the spectral window and provide a measure of the susceptibility of DFT to aliasing. In addition, the Sock Diagram indicates frequency and phase regions with high systematic errors of DFT amplitudes. Both error sources may cause the classical signal-to-noise ratio estimation to fail.

The incorporation of statistical weights into the spectral significance calculation is provided (8.1), as well as the modification of the spectral significance computation, if a time series is assumed to consist of several zero-mean corrected subsets (8.2). Both features also implemented in the present software version.

The S IGS PEC software has proved its validity and excellent performance in numerous practical applications, mainly in connection with the data from the MOST satellite (Walker et al. 2005; Guenther et al. 2005), but also for the photometric search for new pre-main sequence Delta Scuti stars in young open clusters (Zwintz et al. 2004; Zwintz & Weiss 2006). In addition, S IGS PEC was successively employed for the photometric and spectroscopic determination of the rotation periods of selected roAp stars (Sachkov et al. 2004; Ryabchikova et al. 2005).

However, although the statistically unbiased approach used by S IGS PEC improves the reliability of results considerably, the software does not provide an invitation to blind belief. It successfully separates the detection of formally significant sinusoidal signal components in a time series from their interpretation in a physical context, nothing more. The practical application of the program still requires a careful examination of results, namely the decision whether a formally significant frequency is related to the observed physical process or an (e.g. instrumental) artifact.

Acknowledgements

P.R. received financial support from the Fonds zur Förderung der wissenschaftlichen Forschung (FWF, projects P14546-PHY, P14984-PHY) Furthermore, it is a pleasure to thank T. Appourchaux (IAS, Orsay), A. Baglin (Obs. de Paris, Meudon), T. Boehm (Obs. M.-P., Toulouse), M. Breger, R. Dvorak, M. G. Firneis, D. Frast (Univ. of Vienna), R. Garrido (Inst. Astrof. Andalucia, Granada), M. Gruberbauer (Univ. of Vienna), D. B. Guenther (St. Mary's Univ., Halifax), M. Hareter, D. Huber, T. Kallinger (Univ. of Vienna), R. Kuschnig (UBC, Vancouver), S. Marchenko (Western Kentucky Univ., Bowling Green, KY), M. Masser (Univ. of Vienna), J. M. Matthews (UBC, Vancouver), E. Michel (Obs. de Paris, Meudon), A. F. J. Moffat (Univ. de Montreal), E. Paunzen, D. Punz (Univ. of Vienna), V. Ripepi (INAF, Naples), S. M. Rucinski (D. Dunlap Obs., Toronto), T. A. Ryabchikova (Inst. Astr. RAS, Moscow), D. Sasselov (Harvard-Smithsonian Center, Cambridge, MA), S. Schraml (Univ. of Technology, Vienna), G. A. Wade (Royal Military College, Kingston), G. A. H. Walker (UBC, Vancouver), W. W. Weiss, and K. Zwintz (Univ. of Vienna) for valuable discussion and support with extensive software tests. Moreover, M. G. Firneis, D. B. Guenther, J. M. Matthews, G. A. H. Walker, and W. W. Weiss were a great help in clarifying the exposition of this paper.

Finally, I address my very special thanks to the referee of this paper, J. D. Scargle, for his contribution that substantially improved the overall quality of the publication.

  
Appendix A: Orientation of the rms error ellipse of Fourier coefficients

Rotating the coordinate system in Fourier Space by a phase angle $\theta _0$ transforms $a_{\rm ZM}$, $b_{\rm ZM}$ from Eqs. (5), (6) into

$$\alpha\left(\omega ,\theta _0\right) = \frac{1}{K}\sum_{k=0}^{K-1... ...-\! \frac{1}{K}\!\sum_{l=0}^{K-1}\cos\left(\omega t_l-\theta _0\right)\right] ,$$ (A.1)

$$\beta\left(\omega ,\theta _0\right) = \frac{1}{K}\sum_{k=0}^{K-1}... ...! -\! \frac{1}{K}\sum_{l=0}^{K-1}\sin\left(\omega t_l-\theta _0\right)\right] .$$ (A.2)

In terms of the population variance $\left< x^2\right>$ of the time-domain data x_k, the population variances of $\alpha$ and $\beta$ are

  

$$\left<\alpha ^2\right>\left(\omega ,\theta _0\right) = \frac{\lef... ...\sum_{l=0}^{K-1}\cos\left(\omega t_l - \theta _0\right)\right] ^2\Bigg\rbrace ,$$ (A.3)

$$\left<\beta ^2\right>\left(\omega ,\theta _0\right) = \frac{\left... ...m_{l=0}^{K-1}\sin\left(\omega t_l - \theta _0\right)\right] ^2\Bigg\rbrace\cdot$$ (A.4)

The covariance

$$\left<\alpha\beta\right> \frac{1}{K}\sum_{k=0}^{K-1}x_k^2\left[\cos\left(\omega t_k-\theta... ...ight) - \frac{1}{K}\sum_{l=0}^{K-1}\cos\left(\omega t_l-\theta _0\right)\right]$$ =

$$\times\left[\sin\left(\omega t_k-\theta _0\right) - \frac{1}{K}\sum_{l=0}^{K-1}\sin\left(\omega t_l-\theta _0\right)\right]$$ (A.5)

vanishes for a random population, if

 

$$\sum_{k=0}^{K-1}\left[\cos\left(\omega t_k-\theta _0\right) - \fr... ...- \frac{1}{K}\sum_{l=0}^{K-1}\sin\left(\omega t_l-\theta _0\right)\right] = 0 .$$ (A.6)

The constant $\theta _0$ may be separated from the sums, which yields

 

$$\left(\cos ^2\theta _0\!-\! \sin ^2\theta _0\right)\left(K\!\sum_... ...k \!-\! \!\sum_{k=0}^{K}\!\cos\omega t_k\!\sum_{k=0}^{K}\sin\omega t_k\right) =$$

$$\cos\theta _0\sin\theta _0\left[K\sum_{k=0}^{K}\left(\cos ^2\omega t_k-\sin ^2\omega t_k\right)\right.$$

$$-\left.\left(\sum_{k=0}^{K}a\cos\omega t_k\right) ^2+\left(\sum_{k=0}^{K}\sin\omega t_k\right) ^2\right] ,$$ (A.7)

where it is not necessary to distinguish between indices k and l any more. Equation (A.7) may be expressed in terms of $2~\theta _0$ and  $2~\omega t_k$, respectively, according to

$$\cos 2~\theta _0\left[K\sum_{k=0}^{K}\sin 2\omega t_k-2\sum_{k=0}^{K}\cos\omega t_k\sum_{k=0}^{K}\sin\omega t_k\right] =$$

$$\sin 2~\theta _0\left[K\sum_{k=0}^{K}\cos 2\omega t_k-\left(\sum_... ...}\cos\omega t_k\right) ^2+\left(\sum_{k=0}^{K}\sin\omega t_k\right) ^2\right] .$$ (A.8)

This immediately leads to Eq. (10).

  
Appendix B: PDF transformation

  
B.1 General PDF transformation

According to, e.g., Stuart & Ord (1994), p. 20ff, the PDF conversion between a random variable x and a random variable y connected via $y := y\left( x\right)$ is given by

$$\phi_x\left( x\right) = \phi_y\left( y\right)~\left\vert\frac{{\rm d}y}{{\rm d}x}\right\vert\cdot$$ (B.1)

For multivariate distributions, the generalization of differentials yields

$$\phi_x\left( x\right) = \phi_y\left( y\right)~\left\vert\frac{{\rm d}\vec{y}}{{\rm d}\vec{x}}\right\vert ,$$ (B.2)

with the Jacobian

$$\left\vert\frac{{\rm d}\vec{y}}{{\rm d}\vec{x}}\right\vert :... ...\\ \vdots & \vdots & \vdots & \ddots \end{array}\right\vert .$$ (B.3)

   
B.2 Example: transformations of PDF and expected value for the exponential distribution

If y is an exponentially distributed random variable with the PDF

$$\phi_y\left( y\right) =: \kappa~{\rm e}^{-\kappa y} ,$$ (B.4)

Eq. (B.1) for $x:=\sqrt{y}$ evaluates to

$$\phi_x\left( x\right) = 2\kappa x~{\rm e}^{-\kappa x^2} ,$$ (B.5)

where both x, y are considered positive semi-definite.

The expected value of y is

$$\left< y\right> =: \int _{0}^{\infty}{\rm d}y~ y\phi_y\left( y\right) = \frac{1}{\kappa} ,$$ (B.6)

as obtained from Eq. (B.4). Correspondingly, Eq. (B.5) yields

$$\left< x\right> =: \int _{0}^{\infty}dx~ x\phi_x\left( x\right) = \sqrt{\frac{\pi}{2~\kappa}} ,$$ (B.7)

so that

$$\left< y\right> = \frac{4}{\pi}~\left< x\right> ^2$$ (B.8)

describes the transformation of expected values corresponding to y = x².

For equidistantly sampled noise x_k in the time domain, the power spectrum associated with the Fourier Series is exponentially distributed (Schuster 1898; Scargle 1982) according to

$$\phi\left( A^2\right) = \frac{K}{2\left< x^2\right>}~{\rm e}^{-\frac{KA^2}{2\left< x^2\right>}} ,$$ (B.9)

where $\left< x^2\right>$ denotes the population variance of the noise. The previous considerations apply through the substitution $\kappa =: \frac{K}{2\left< x^2\right>}$, normalized to single-sided power spectral density (which introduces the factor 2). According to Eq. (B.6), the noise level in the spectrum of squared amplitudes evaluates to

$$\left< A^2\right> = \frac{2\left< x^2\right>}{K}\cdot$$ (B.10)

For the expected value of the amplitude (i.e. the amplitude noise level), Eq. (B.7) yields

$$\left< A\right> = \sqrt{\frac{\pi}{K}~\left< x^2\right>} .$$ (B.11)

Equation B.8 leads to the conversion between the expected values of amplitude, $\left< A\right>$ and squared amplitude, $\left< A^2\right>$, according to

$$\left< A^2\right> = \frac{4}{\pi}~\left< A\right> ^2.$$ (B.12)

This is in agreement with the observational results by Horne & Baliunas (1986).

An amplitude signal-to-noise ratio of 4 - introduced as an approximate significance criterion by Breger et al. (1993) - corresponds to a power signal-to-noise ratio of $4~\pi\approx 12.57$.

References

• Antoci, V., Breger, M., Bischof, K., & Garrido, R. 2006, ASPC, 349, 181 [NASA ADS] (In the text)
• Breger, M., McNamara, B. J., Kerschbaum, F., et al. 1990, A&A, 231, 56 [NASA ADS] (In the text)
• Breger, M., Stich, J., Garrido, R., et al. 1993, A&A, 271, 482 [NASA ADS] (In the text)
• Breger, M., Handler, G., Garrido, R., et al. 1999, A&A, 349, 225 [NASA ADS] (In the text)
• Cumming, A., Marcy, G. W., & Butler, R. P. 1999, ApJ, 256, 890 [NASA ADS] [CrossRef] (In the text)
• Eyer, L., & Bartholdi, P. 1999, A&AS, 135, 1 [NASA ADS] (In the text)
• Firneis, M. G. 1970, Zufallszahlengeber und stellarstatistische Phänomene (Ph.D. Thesis, Univ. of Vienna) (In the text)
• Foster, G. 1995, AJ, 109, 1889 [NASA ADS] [CrossRef] (In the text)
• Granzer, T., Strassmeier, K. G., & Reegen, P. 2001, AN, 322, 325 [NASA ADS] (In the text)
• Guenther, D. B., Kallinger, T., Reegen, P., et al. 2005, ApJ, 635, 547 [NASA ADS] [CrossRef] (In the text)
• Horne, J. H., & Baliunas, S. L. 1986, ApJ, 302, 757 [NASA ADS] [CrossRef] (In the text)
• Kotelnikov, V. A. 1933, On the carrying capacity of the ether and wire in telecommunications, Material for the First All-Union Conference on Questions of Communication (Moscow: Izd. Red. Upr. Svyazi RKKA) (In the text)
• Koen, C. 1990, ApJ, 348, 700 [NASA ADS] [CrossRef] (In the text)
• Kovacs, G. 1980, Ap&SS, 78, 175 [NASA ADS] (In the text)
• Kuschnig, R., Weiss, W. W., Gruber, R., Bely, P. Y., & Jenkner, H. 1997, A&A, 328, 544 [NASA ADS] (In the text)
• Lafler, J., & Kinman, T. D. 1965, ApJS, 11, 216 [NASA ADS] [CrossRef] (In the text)
• Lenz, P., & Breger, M. 2005, CoAst, 146, 53 [NASA ADS] (In the text)
• Lomb, N. R. 1976, ApSS, 39, 447 [NASA ADS] (In the text)
• Marconi, M., & Palla, F. 2004, in The A-Star Puzzle, Proceedings of IAUS 224, ed. J. Zverko, J. Ziznovsky, S. J. Adelman, W. W. Weiss (Cambridge: Cambridge Univ. Press), 69
• de Moivre, A. 1718, The Doctrine of Chances (London: Millar) (In the text)
• Nyquist, H. 1928, Trans. AIEE, 47, 617 (In the text)
• Montgomery, M. H., & O'Donoghue, D. 1999, Delta Scuti Star Newsletter, 13, 28 (University of Vienna, Austria)
• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in C, 2nd ed. (Cambridge: Cambridge Univ. Press) (In the text)
• Reegen, P. 2005, in The A-Star Puzzle, Proceedings of IAUS 224, ed. J. Zverko, J. Ziznovsky, S. J. Adelman, & W. W. Weiss (Cambridge: Cambridge Univ. Press), 791 (In the text)
• Reegen, P., Kallinger, T., Frast, D., et al. 2006, MNRAS, 367, 1414 [NASA ADS] (In the text)
• Roberts, D. H., Lehar, J., & Dreher, J. W. 1987, AJ, 93, 968 [NASA ADS] [CrossRef] (In the text)
• Ryabchikova, T., Wade, G. A., Auriere, M., et al. 2005, A&A, 429, 55 [NASA ADS] [CrossRef] (In the text)
• Sachkov, M., Ryabchikova, T., Kochukhov, O., et al. 2004, ASPC, 310, 208 [NASA ADS] (In the text)
• Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] (In the text)
• Schuster, A. 1898, Terr. Mag., 3, 13 (In the text)
• Schwarzenberg-Czerny, A. 1989, MNRAS, 241, 153 [NASA ADS] (In the text)
• Schwarzenberg-Czerny, A. 1996, ApJ, 460, L107 [NASA ADS] [CrossRef] (In the text)
• Shannon, C. E. 1949, Proc. Institute of Radio Engineers, 37, 10 (In the text)
• Sperl, M. 1998, CoAst, 112, 1 [NASA ADS] (In the text)
• Stellingwerf, R. F. 1978, ApJ, 224, 953 [NASA ADS] [CrossRef] (In the text)
• Strassmeier, K. G., Boyd, L. J., Epand, D. H., & Granzer, T. 1997, PASP, 109, 697 [NASA ADS] [CrossRef] (In the text)
• Stuart, A., & Ord, J. K. 1994, Kendall's Advanced Theory of Statistics, vol. 1 Distribution Theory, 6th ed. (London: Arnold) (In the text)
• Walker, G., Matthews, J., Kuschnig, R., et al. 2003, PASP, 115, 1023 [NASA ADS] [CrossRef] (In the text)
• Walker G. A. H., Matthews, J. M., Kuschnig, R., et al. 2004, BAAS, 36, 1361 [NASA ADS] (In the text)
• Walker G. A. H., Kuschnig, R., Matthews, J. M., et al. 2005, ApJ, 623, L145 [NASA ADS] [CrossRef] (In the text)
• Whittaker, E. T. 1959, Proc. Royal Soc. Edinburgh, A35, 181
• Zwintz, K., & Weiss, W. W. 2006, A&A, 457, 237 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Zwintz, K., Marconi, M., Kallinger, T., & Weiss, W. W. 2004, in The A-Star Puzzle, Proceedings of IAUS 224, ed. J. Zverko, J. Ziznovsky, S. J. Adelman, & W. W. Weiss (Cambridge: Cambridge Univ. Press), p. 353 (In the text)