1

A study of periodograms standardized using training
data sets and application to exoplanet detection

arXiv:1702.02049v1 [] 7 Feb 2017

Sophia Sulis, David Mary and Lionel Bigot

Abstract—When the noise affecting time series is colored
with unknown statistics, a difficulty for sinusoid detection is to
control the true significance level of the test outcome. This paper
investigates the possibility of using training data sets of the noise
to improve this control. Specifically, we analyze the performances
of various detectors applied to periodograms standardized using
training data sets. Emphasis is put on sparse detection in the
Fourier domain and on the limitation posed by the necessarily
finite size of the training sets available in practice. We study the
resulting false alarm and detection rates and show that standardization leads in some cases to powerful constant false alarm
rate tests. The study is both analytical and numerical. Although
analytical results are derived in an asymptotic regime, numerical
results show that theory accurately describes the tests’ behaviour
for moderately large sample sizes. Throughout the paper, an
application of the considered periodogram standardization is
presented for exoplanet detection in radial velocity data.
Index Terms—Multiple sinusoids’ detection, colored noise,
periodogram standardization, sparse detection, asymptotic.

I. I NTRODUCTION
A. Considered problem and application

D

ETECTING sinusoids in noise is one of the most studied
problems in signal processing. Sinusoid detection is
classically based on Fourier analysis, which benefits from a
considerable set of statistical results. Often, the assumptions
required for these results to hold are :
(i) Under the null hypothesis (H0 : noise only), the random
process X is stationary Gaussian with known statistics;
(ii) The time series {X(tj )}j=1,··· ,N is regulary sampled,
in which case Fourier analysis is performed through the
periodogram [1], [2]
P (ν) :=

N
2
1 X
X(tj )e−i2πνj∆t ,
N ∆t j=1

(1)

where ν is the considered frequency, N the number of samples
and ∆t the time sampling step;
(iii) The number of samples N is large (asymptotic regime);
(iv) The periodogram is evaluated at Fourier frequencies.
Assumptions (i)−(iv) allow to characterize the statistical
properties of the periodogram as an estimate of the power
spectral density (PSD) [2]–[9]. In particular, any finite set
of periodogram ordinates near some fixed frequency are
This work was supported by Thales Alenia Space, PACA region, CNRS
project DETECTION/Imag’In and Programme National de Physique Stellaire
(PNPS) of CNRS/INSU, France.
The authors are with the Université Côte d’Azur, OCA, CNRS, Laboratoire Lagrange, France (e-mail: Sulis.Sophia@oca.eu; David.Mary@unice.fr;
Lionel.Bigot@oca.eu).

asymptotically independent and exponentially distributed with
parameters dependent on the noise PSD at the relevant frequencies (this assumption actually extends to linear processes
other than Gaussian, see Th. 4 of [9]). In detection, such
properties can be exploited to obtain the probabilities of false
alarm and detection of tests based on periodogram ordinates.
Despite a large literature on the subject, the detection of
periodic signals remains an active field of research because
in practical situations some (or all) assumptions (i)-(iv) above
may not be met. Before reviewing some such situations, it is
worth introducing the application that motivated the present
study (the proposed detection approach may cover other
applications, however).
With almost 3 000 confirmed exoplanets and nearly 2 500
candidates (in September 2016, www.exoplanets.org), the
detection of extrasolar planets is an extremely active research
field in Astrophysics since two decades. This field benefits
from constant technological improvements allowing extremely
low noise detectors [10] and long (months to years) spaceborne
observations with high (a few tens of seconds) sampling rates
[11]–[13]. The Radial Velocity (RV) technique is one method
for exoplanet detection. When planets orbit a host star, the
resulting gravitational force creates a periodic displacement
of the star. This induces a modulation of the RV of the star
with respect to (w.r.t.) Earth, which translates into a periodic
Doppler shift of the stellar light. The RV technique consists in
detecting such variations in stellar RV time series [14], [15].
Modern instrumental performances have reached signal to
noise ratios allowing, in principle, to detect exoplanets comparable to the Earth by the RV method. However, at such low
levels of instrumental noise, a new critical and limiting issue
appears. The stellar surface can be seen as a boiling fluid, with
millions of convection cells generating upward and downward
plasma flows visible under the form of granules having typical
lifetime of a few minutes. These motions generate random
fluctuations in the measured RV of the star that can mimic,
or hide, exoplanetary signatures. The case of α Centauri Bb
planet, a detection claimed in 2012 [16] (with an evaluated
P -value of 0.02%) and since then subject to controversy [17],
[18], is one example showing how our incomplete knowledge
of the noise affects the reliability assigned to exoplanet detection claims. This specific issue motivated the present study.
A key point is that in parallel to technological advances,
astrophysicists have continuously improved stellar models and
elaborated numerical simulation codes able to account for
the complex interplay of various astrophysical processes in
the star’s interior and surface. Recent works demonstrate that
granulation noise can be reproduced by large scale numerical

©2017 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or
reuse of any copyrighted component of this work in other works. DOI: 10.1109/TSP.2017.2652391

2

simulations in a reliable way [19]. This suggests the possibility
of using such simulations to calibrate the detection process, as
RV time series are strongly affected by this stellar noise. The
present study shows that one such calibration indeed leads to
improved control of the statistical significance and to detection
tests with increased power.
B. Unknown noise statistics: related works
Turning back to the deviations encountered in practice
w.r.t. assumptions (i)−(iv) above, this work will be primarily
interested in the case of an incomplete knowledge of the noise
statistics under the null hypothesis (i.e., relaxing condition (i);
we mention in the last section perspectives to relax the other
conditions). In this situation, the distribution of P under the
null is not known and consequently the significance level (the
size) of the test is not known either. Constant false alarm
rate (CFAR) detectors have been devised when the noise is
white [20]–[22]. When the noise is colored, the detection
problem is more complicated. In practice, two approaches can
be followed. The first approach is simply to ignore possible
noise correlations and to apply tests designed for white noise.
However, as will be illustrated in this study, the statistical
behavior of the resulting testing procedure may be hazardous,
with unpredictable significance level and poor power (see also
[23] on this point).
A more sophisticated approach consists of estimating the
noise PSD (called SE below) so that condition (i) above is
considered to hold approximately. This estimate can then be
used to calibrate the periodogram of the data P (ν), leading to
a frequency-wise standardized periodogram of the form
P (ν)
Pe (ν | SbE ) :=
.
b
SE (ν)

(2)

Note that the classical Fisher’s test [20] standardizes the
periodogram ordinates by the estimated PSD of a white noise.
Standardization (2) can be seen as a generalization of this
approach (see [24] for a recent review).
The estimate SbE (ν) can be parametric or non-parametric. Non
parametric approches originate from seminal works of Whittle
[25] and Bartlett [26]. Parametric methods often proceed by
fitting AutoRegressive (AR) or ARMA (AR Moving Average)
processes to the time series.
A further complication arises when multiple sinusoids are
present under the alternative, as they perturb the estimation
of the noise PSD [5], [27]. Standardized tests for this case
can be found in [21], [22], [24], [27]–[29]. These tests are
however non adaptive in the number of sinusoids (which must
be set a priori) and designed for white noise. For adaptive
procedures for colored noise see Chap. 8 of [5], and [30]–
[38]. Techniques reducing the influence of signal peaks under
the alternative are proposed in [39], [40]. Different approaches,
related to standardization (2), can be found in [41]–[44].
When following Generalized Likelihood Ratio (GLR) approaches for detecting multiple sinusoids in unknown number,
the GLR must be combined with model selection procedures.
While sharp model selection criteria and CFAR detectors
exist under white noise, the correlated case remains an open

problem [45]. In the particular field of exoplanet detection
using RV, we find similar families of techniques [17], [46]–
[48].
In conclusion, regarding the problem of assessing tests’
significance levels for multiple sinusoids detection in noise,
an inspection of the literature shows that:
• For white noise of unknown variance, several studies provide
accurate results for standardized test statistics of the form (2),
e.g., [21], [22], [27]–[29].
• For colored noise with unknown PSD, we are not aware
of works studying the false alarm rate when AR/ARMA or
other models are used for test standardization as in (2). The
difficulty in this case is the dependence of the distribution of
SbE on estimated parameters, which complicates the analytical
characterization of the distribution of Pe.
The procedures described in [5], [32], [37] provide asymptotic
control of the false alarm rate. These procedures do not operate
explicitly on test statistics of the form (2) but rather on
windowed periodograms that depend on several parameters.
We found that these tests are in practice sensitive to parameter
setting and that estimating these parameters impacts the
significance level at which the tests are conducted. This
level can of course be approximated by simply neglecting
the influence of such a ‘preprocessing stage’ (dealing, e.g.,
with model order selection, filtering, adaptive window design,
or standardization). For instance, we might pretend that
SbE = SE in (2). As will be highlighted in Sec. VI and VII, the
actual significance level obtained when doing so may however
be far from the assumed one. This leaves open the question of
designing both powerful and CFAR tests for unknown colored
noise and we propose such tests in this paper.
Before closing this literature survey, we mention a few
tests designed for a particularly interesting configuration of
the detection problem, which is the so-called rare and weak
setting. In this setting, the sinusoids are both of small amplitudes w.r.t. the noise level and in small (and unknown)
number w.r.t. the number of samples N . When viewed in the
Fourier domain, sinusoid detection can be casted as a sparse
heterogeneous mixtures problem, which has attracted much
attention in the last decade [49]–[53]. We will see that while
fixed and adaptive (in the number of sinusoids) procedures
lead to inconsistencies when noise correlations are ignored
(because then the statistics under the null hypothesis are
wrongly specified), such tests keep their nominal properties,
with the CFAR property added, when applied to periodograms
suitably standardized with training data sets. In the particular
case of the Higher Criticism [49], standardization (2) is an
alternative approach to that of [54].
C. This study
The present study (an extended version of [55]) focuses on
the statistical characterization of test statistics when both the
PSD of the colored noise and the parameters of the sinusoids
are unknown. We propose a detailed analysis of the effects
of periodogram standardization by means of, say, L training
time series, which are independent realizations of the noise
process alone, and of the gain that can be expected by using

3

such training signals in a detection framework.
We, of course, make the important assumption that such a
training data set is available. Beyond the case of exoplanet
detection considered here, one may imagine various situations
where training signals can be obtained. In astronomical instruments for instances, secondary optical paths are often devoted
to monitor ‘empty’ regions of the sky or calibration stars [56].
Note that training noise vectors are routinely used for detection
in radar systems, with however, an important difference.
Adaptive test statistics in radar typically use estimates of the
covariance matrix of the training vectors and therefore require
L > N for this matrix to be nonsingular. This is a very
different regime from that considered here, where L ≪ N .
In the present study, we assume that the training data set is
unbiased, in the sense that an averaged periodogram obtained
from an infinitely large batch would converge uniformly to the
true noise PSD. In practice, finite (possibly small) batch sizes
can be encountered. For this reason, we say below that the
noise is partially unknown and we address the effect of small
values of L on the detection performances.
Because one important objective of this study is to obtain
analytical characterization of the test performances, we consider here a regular sampling. Comments on how to relax this
assumption are discussed at the end of the paper. Also, our
results are asymptotic in the number of samples N , which
is characteristic of time series analysis. However, we will
also pay attention to whether asymptotic theory accurately
describes reality for finite sample sizes through simulations.
In the considered application framework of exoplanet detection
in RV data the working hypotheses are justified because
accurate simulations of stellar noise can be produced to form
training data sets. These simulations are however computationally demanding. Obtaining a simulation of 100 days,
for a star similar to the one shown in [19], takes about 3
months of computing time on 120 cores on modern clusters.
Consequently, realistic values of L are in the range of one to,
say, a hundred at most. This motivates the study of the impact
of estimation noise in the proposed standardization approach.
We proceed as follows. Sec. II presents the model and the
detection approach. Sec. III recalls classical results regarding
periodogram’s distribution. Sec. IV and V derive false alarm
and detection rates for several tests. Sec. VII is a numerical
study. Table I summarizes the main notations used in the paper.

II. S TATISTICAL MODEL AND DETECTION APPROACH
We consider the two hypotheses:

H : X(tj ) = E(tj )


 0
Ns
X
(3)

αq sin(2πfq tj + ϕq ) + E(tj )
H
:
X(t
)
=
1
j


q=1

where X(tj = j∆t), j = 1, . . . , N , is an evenly sampled
data time series, E(tj ) is a zero-mean second-order
stationary
P
Gaussian noise, with inf(SE (ν)) > 0 and u∈R |rE (u)| <
∞. This is the noise of which we assume a set of training time
series is available. To simplify the presentation, we consider
for the rest of the paper a unit sampling step ∆t = 1 in (1).

TABLE I
TABLE OF NOTATIONS
N
X(tj )

Number of data points
Regularly sampled time series

E(tj )

Zero-mean stationary Gaussian colored noise

ν, νk

Continuous frequency, Fourier frequency

SE (ν)

Noise PSD

rE

Noise autocorrelation function

Ns , αq , fq , ϕq

Parameters of model (3): Number of sines,
sines’ amplitude, frequency and phase

Np

Number of exoplanets orbiting target star

T p , K p , Mp
ep , ωp , t0 , γ0
NC

Planet period and its six Keplerian parameters

P (ν)

Classical periodogram

Proxy for Ns

L

Number of available training data sets

P L (ν)
Pe(ν|P L )

Periodogram averaged with L training data sets

Ω

Indices set of considered Fourier frequencies

λk := λ(νk )

Non centrality parameter

Fλk (d1 , d2 )

Non central Fisher-Snedecor distribution
with d1 and d2 degrees of freedom

Z, z

Scalar random variable, one realization of Z

Periodogram standardized by P L

vZ (z)

Pr(Z > z) (observed p-value)

VZ
Z = [Z1 , . . . , ZN ]⊤

P-value as a random variable (VZ ∼ U[0 1] )
Vector of random variables

VZ,k

P-value associated to component Zk of Z

VZ,(k)

k-th ordered p-value of Z

Under the alternative, the Ns amplitudes αq ∈ R∗+ ,
frequencies fq ∈ R∗+ and phases ϕq ∈ [0, 2π[ of the
deterministic part are unknown. In RV exoplanet detection,
this deterministic part represents the planetary signature(s).
Note that model (3) can actually be useful for the detection
of periodic signals more general than only pure sinusoids. In
such cases, the Fourier spectrum may contain many harmonics.
Because the fundamental frequency has zero probability to
coincide with a Fourier frequency, the corresponding number
of nonzero Fourier coefficients (i.e., of deviations under the
alternative) always equals N . However, most of the energy
of periodic signals is captured by a small fraction of Fourier
coefficients, so that model (3) would often be accurate for such
signals with some Ns ≪ N .
In the case of RV signals for instance, a study of the influence
of the Keplerian parameters (planets’ orbital parameters)
shows that this is indeed the case [57]. In all cases except
perhaps very rare and exotic systems, the RV spectral signatures exhibit only a small fraction of significant harmonics
because the planets tend to have low eccentricities and are in
small number (Np ). In short, the spectrum is sparse (though
not strictly sparse) and RV signals can be modeled by a sum
of a small number of pure sinusoids. We call this number Ns ,
and we say that Ns ≪ N . When there is one planet with
frequency close to the Fourier grid, Ns will be essentially 1.
In our simulations for Sec. VII, we found that multiplanetary
systems with 5 eccentric planets behave like model (3) with
Ns not exceeding, say, 20 at most.

4

III. P ERIODOGRAMS ’

STATISTICS : ASYMPTOTICS

A. Classical (Schuster’s) periodogram
The frequencies considered in (1) will be Fourier frequenk
cies {νk := N
}k=0,...,N −1 and N is considered even. For
simplicity but without loss of generality we will often consider
the subset of ( N2 − 1) Fourier frequencies corresponding to
k ∈ Ω := {1, . . . , N2 − 1}. Asymptotically, the periodogram
P in (1) is an unbiased but inconsistent estimate of the
PSD [2]. Under the above assumptions on E, the periodogram
ordinates at different frequencies νk and νk′ are asymptotically
independent [24].
Under H0 , the asymptotic distribution of P is (Th. 5.2.6, [2]):

S (ν )

 E k χ22 ,
∀k ∈ Ω,
2
P (νk |H0 ) ∼
(4)

 SE (νk )χ2 , for k = 0, N .
1
2
Under H1 , the distribution of P is known when cisoids are
present in (3) ( [24], Corollary 6.2(b)). The real case of model
(3) can be treated similarly (see Appendix A). This leads to
the asymptotic distribution:

S (ν )

 E k χ22,λk , ∀k ∈ Ω,
2
(5)
P (νk |H1 ) ∼

 SE (νk )χ2 , for k = 0, N .
1,λk
2
The {λk := λ(νk )} are non centrality parameters given for
k ∈ Ω by
λk =

Ns h
Ns
i
X
X
N
α2q κ2q +2αq κq
αℓ κℓ cos(θq −θℓ ) , (6)
2SE (νk ) q=1
ℓ=q+1

and for k = 0, N2 this expression is halved. The terms κq and
θq , given by (39) and (40) in Appendix A, arise from spectral
leakage through the spectral window KN (36). Owing to the
fast decay of KN (ν), the proportion of parameters λk that
significantly differ from 0 is small if Ns ≪ N .
B. Averaged periodogram
We assume that a training data set T of independent
realisations of the colored noise is available. This set is
obtained by L independent simulations corresponding to L
time series
Xℓ sampled on the same grid as the observations:

T = {Xℓ (tj )}j=1,...,N ℓ=1,...,L . A straightforward estimate
of the noise PSD is the averaged periodogram [3]:
P L (νk |H0 ) :=

N
L
2
1X 1 X
Xℓ (tj )e−i2πνk j .
L
N j=1

The effect of stochastic estimation noise caused by the
finiteness of T is encapsulated in L. This clearly impacts the
distribution of P L in (7) and in turn the efficiency of the
subsequent standardization by P L .
C. Periodograms standardized with P L
We now turn to statistical properties of standardized periodograms of the form (2). When the averaged periodogram
P L is used, this yields:
P (νk )
Pe(νk | P L ) :=
.
P L (νk )

As the numerator and denominator are independent variables with known asymptotic distributions, assessing the distribution of their ratio is straightforward. The ratio of two
independent random variables (r.v.) V1 ∼ χ2d1 and V2 ∼ χ2d2
follows a Fisher-Snedecor law noted F (d1 , d2 ) with (d1 , d2 )
1
degrees of freedom: VV12 /d
/d2 ∼ F (d1 , d2 ) [59].
Consequently, from (4) and (7), the asymptotic distribution of
this standardized periodogram under H0 is:

SE (νk )χ22 /2


 S (ν )χ2 /2L ∼ F (2, 2L), ∀k ∈ Ω,
E k
2L
Pe(νk |P L , H0 ) ∼

N
SE (νk )χ21


∼ F (1, L), for k = 0, .
2
SE (νk )χL /L
2
(9)
Similarly, from (5) and (7), we have under H1 :
 2

 χ2,λk /2 ∼ F (2, 2L), ∀k ∈ Ω,

λk
 χ2 /2L
2L
(10)
Pe (νk |P L , H1 ) ∼

χ21,λk
N


∼ Fλk (1, L), for k = 0, ,
 2
χL /L
2

where Fλk denotes a non-central F distribution with non
centrality parameter λk given by (6).
Under the null hypothesis, (9) shows that the distribution of
the standardized periodogram is independent of the nuisance
signal, i.e., of the partially unknown noise PSD. This property
is important as it will be inherited by some of the test statistics
investigated in Sec. VI, leading to CFAR tests.
IV. C ONSIDERED

TESTS

A. Preliminary notations
These tests are better presented using P -values1 and order
statistics. When necessary the notation will distinguish between a r.v. Z and its realization z. We recall that the observed
P -value vZ is defined as:
vZ (z) := Pr (Z > z)

ℓ=1

Using (4), the asymptotic distribution of P L is:

S (ν )

 E k χ22L ,
∀k ∈ Ω,
2L
P L (νk |H0 ) ∼

N
 SE (νk ) χ2 ,
for k = 0, .
L
L
2

(8)

(7)

P L is a consistent and unbiased estimate of SE (ν) when both
L → ∞ and N → ∞ ( [58], Chap.14).

and vZ is one realization of the r.v. VZ , which is uniformly distributed. Similarly, for a vector of r.v. Z = [Z1 , Z2 , . . . , ZN ]⊤
of which z = [z1 , z2 , . . . , zN ]⊤ is one realization we will
denote by
min zk := z(1) < z(2) < . . . < z(N ) := max zk
k

k

1 We will denote the P -values by V because P denotes the periodogram
in this work.

5

the ordered values of z and by Z(1) , . . . , Z(N ) the order
statistics of Z. The observed P -values corresponding to z will
be denoted by vZ,k (with vZ,k := Pr (Zk > zk )) and the
observed ordered P -values by vZ,(k) . The corresponding r.v.
will be denoted by VZ,k and VZ,(k) . The ordered P -values
VZ,(k) are not uniform, because they are order statistics from
a uniform distribution, and obviously dependent [60].
We now present some tests discussed in the Introduction and
selected for reference as they cover different classical models.
Under assumptions specified below on the distributions of the
variates {Zi }, the properties of these tests are known and we
shall summarize them.
B.

Test statistics

1
All tests below are of the form T(z) ≷H
H0 γ, with z the
data (which below will be a set of ordinates of one of the
periodograms discussed above), T(·) the test statistic and
γ ∈ R+ a threshold that determines the false alarm rate.

1) Test of the maximum:
TM (Z) := Z(N ) .

(11)

For independent variates Zi of Cumulative Distribution
Function (CDF) ΦQ
Zi , the false alarm of the test is
N
PFA (γ) = 1 − i=1 ΦZi (γ). For a model involving
under H1 a single sinusoid with unknown frequency
(but on the Fourier grid) and under H0 a white Gaussian
noise (WGN) of known variance σ 2 , the ordinates of the
periodogram P evaluated at successive Fourier frequencies
are under H0 independent and identically distributed (i.i.d.)
with distributions given by (4), where SE = σ 2 . For this
model TM (P) corresponds to the GLR test [61].
2) Fisher’s test:
Z(N )
TF (Z) := X
.
Z(k)

(12)

k

When applied to periodogram of WGN of unknown
variance, Fisher’s test is CFAR (the distribution of the
test statistics is established in [20]). This test is actually
the GLR test under the model of a single sinusoid on the
Fourier grid and WGN of unknown variance (see [62], who
also covers the case of more than one sinusoid). Examples
of works using this test in Astronomy are [46], [47], [63]–[65].
3) A test inspired by the tests of Chiu and Shimshoni:
TC (Z) := Z(N −NC +1) ,

(13)

with NC ≥ 1 a parameter related to the number of sinusoids.
This test statistic is justified by the observation made in [21],
[22] that for multiple sinusoids, order statistics different from
the maximum may be more discriminative than Z(N ) against
the null. These tests are designed for periodograms of white
noise of unknown variance and their asymptotic false alarm
rates are given in [21], [22]. As Fisher’s test, they involve
denominators whose purpose is normalization by consistent

estimates of the noise variance. TC is a simplification of
these tests: the normalization is ignored, because it will
not be necessary to yield a CFAR detector once applied to
periodograms standardized by P L . The false alarm rate of
test TC for white noise can be deduced from the expression
obtained in Sec. V-C.
As for the tests [21], [22], we expect TC to have decreasing
power in case of strong mismatch between the value of
parameter NC and the number of detectable deviations under
H1 (roughly speaking, Ns ). Not fixing NC in advance but
estimating this parameter from the data may lead to more
powerful tests, but at the cost of a more difficult control
of the FA rate (as NC becomes random). This suggests to
consider other approaches that are adaptive in the number of
sinusoids, which is the case of the last two tests (14) and (16).
4) Higher Criticism: This test statistic is defined by [49]:
√
N (k/N − vZ,(k) )
⋆
HC (Z) := max p
,
(14)
1≤k≤α0 N
vZ,(k) (1 − vZ,(k) )

where vZ,(k) are ordered P -values and parameter α0 ∈ [ N1 , 1].
HC is designed under the assumption that under the null
hypothesis the ordinates {Zk } are i.i.d. with known marginal
distribution. When Z = σ2P
2 , with P the periodogram of
E
2
a white noise of known variance σE
under the null, this
2
distribution is given by (4) with SE = σE
. The ordered P values involved in (14) are thus ordered values of

 1 − Φχ22 (P (νk )), ∀k ∈ Ω,
(15)
v2P/σE2 , k
 1 − Φχ2 (P (νk )), for k = 0, N .
1
2
Under the alternative, a fraction of the ordinates contain a
2
deterministic part, and hence follow (5), with SE = σE
(see
Sec. 1.7 of [49]). The frequencies at which these deviations
occur are unknown. Their magnitudes ( related to λk ) are also
unknown and weak, in the sense that they are comparable to
the expected magnitude of the periodogram maximum under
the null hypothesis.
Optimality2 results of HC are asymptotic and established
for a specific sparsity vs amplitude model. When the deviation
occurs at a single frequency or for extremely sparse signatures
in the Fourier domain, [49] showed that the test based on
the maximum (11) is asymptotically optimal in the sense that
whenever the Neyman-Pearson test has full power, test (11)
has full power as well. In such cases, the largest periodogram
ordinate (or, equivalently, the smallest P -value) is the most
powerful test statistic to discriminate against the null.
As discussed about the tests [21], [22] and TC above, better
regions than the maximum ordinate/smallest P -value may be
useful for sparse but not extremely sparse signatures. [49]
showed that there exists sparse signals of weak amplitude
that can be optimally detected (in the sense of full power) by
HC but not by TM . This makes HC particularly interesting in
the present study.
2 A test is said optimal in [49] if the sum of the probability of error and
the probability of missed detection tends to zero.

6

5) Berk-Jones : A test statistic related to HC is that of
Berk-Jones [52], [66]–[69] defined by:
BJ(Z) := max I1−vZ,(k) (N − k + 1, k),
1≤k≤α0 N

(16)

where I denotes the regularized incomplete beta function [59].
The deviations in form of Z−scores in (14) for HC are
established using the asymptotic convergence of a binomial
distribution to a Gaussian distribution. In the tails, however,
this convergence is very slow (see [52], [67], [70] for illustrations). For this reason the test statistics (16), which compares
favorably to HC and other goodness-of-fit (GOF) tests in some
cases, was recently (and almost simultaneously) proposed by
[52], [66], [68], [69]. As noted in [52], [53], this test was
initially proposed by Berk and Jones (and called Mn+ ) in [71].
This test is based on the exact significance reflected by
the P -values, that is, on the P -values of the ordered P values. Since the ordered P -values are Beta distributed, with
V2P/σE2 ,(k) ∼ Beta(k, N − k + 1), their P -values involve the
CDF of Beta variables, which is an incomplete Beta function:
Pr(Beta(k, N − k + 1) ≤ x) = Ix (k, N − k + 1) and leads
to test (16) with the P -values computed as in (15).
BJ presents the same adaptive optimality as HC for sparse
mixture detection, and the asymptotic distribution of the BJ
test statistic can be found in Th. 4.1 of [52]. As for HC, convergence to the asymptotic distribution may be slow. Efficient
algorithms for computing significance levels (and hence the
function γ 7→ PFA (γ)) of HC and BJ for finite (but possibly
large) values of N can be found in [52], [72]. BJ and HC are
studied from the viewpoint of local levels in [53].
V. T ESTS

APPLIED TO

e | PL
Z=P

To simplify the presentation of the results we restrict in the
following to the frequency set Ω, i.e., to standardized vectors
#⊤
"
P (ν N −1 )
P
(ν
)
1
2
e | P :=
.
P
,...,
L
P L (ν1 )
P L (ν N −1 )
2

Extension of the results to ν0 , ν N can be obtained using the
2
distributions (9) and (10).
In the following we evaluate false alarm and detection
rates by postulating independence of the considered ordinates. While the asymptotic independence of the periodograms
ordinates at positive Fourier frequencies is well established
(e.g. [73], Th. 2.14), theoretical results regarding the joint
distribution of periodogram ordinates under departures from
whiteness (and Gaussian) assumptions are lacking, however.
For some results on the largest order statistics, see [74], who
consider MA Gaussian processes, and [75], who consider nonGaussian sequences. For finite values of N , the approaches
of [76], [77] might be followed to better characterize the
performances of some tests considered below.
In practical situations, the marginal distributions of the considered ordinates are only approximated by their asymptotic
distribution (this is visible in Eq. (29) for instance). Besides,
as N grows the correlation between the periodogram ordinates
at the signal frequencies approaches zero at a slower rate than
the correlation between other ordinates ([24], Theorem 6.5),

indicating that departure to the independence assumption may
be more pronounced under H1 than H0 . Consequently, we
are using the independence as an operational assumption to
quantify the tests’ performances, and the validity of the resulting expressions below should be checked against numerical
simulations. Sec. VII provides several examples.
A. TM .
e | PL )
Under H0 and in the asymptotic regime, TM (P
is the maximum of independent variables given by (9). For
k ∈ Ω these N2 − 1 variables follow an F (2, 2L) distribution
with generalR density given in [59]. Using the beta function
1
B(1, L) := 0 (1 − t)L−1 dt = L1 , this density ϕF (γ, 2, 2L)
can be expressed as

γ −L−1
1
1 
γ −L−1 
= 1+
.
· · 1+
B(1, L) L
L
L
R∞
It can be checked that 0 ϕF (γ, 2, 2L)dγ = 1. The corresponding CDF ΦF (γ, 2, 2L) is obtained by integration of ϕF :
ϕF (γ, 2, 2L) =

ΦF (γ, 2, 2L) =

Z

0

γ

L
ϕF (γ, 2, 2L)dγ = 1−
L+γ

!L

. (17)

The probability of false alarm (PFA ) can be computed thanks
e | PL :
to the asymptotic independence of the ordinates of P

e | PL ), γ) := Pr (TM (P
e | PL ) > γ|H0)
PFA (TM (P


 N2 −1

Y
=1−
Pr Pe (νk ) ≤ γ|H0 , P L = 1 − ΦF (γ, 2, 2L)
k∈Ω



=1− 1−



L L  N2 −1
.
γ+L

(18)
The probability of false alarm is (asymptotically) indepene | PL ) a CFAR
dent of the noise PSD, which makes TM (P
detector.
Under H1 , using (10) and the approximate independence of
periodogram ordinates ( [24], Theorem 6.5), the probability of
e | PL ) can be approximated as:
detection (PDET ) of TM (P


e | PL ), γ) := Pr TM (P
e | PL ) > γ|H1
PDET (TM (P
Y
(19)
≈1−
ΦFλk (γ, 2, 2L).
k∈Ω

With (18), the relationship γ(PFA ) for TM can be derived as:
e | PL ), PFA ) = L
γ(TM (P

h

1

1 − (1 − PFA ) η

− L1

i
− 1 , (20)

where η := N2 −1. With (19) and (20), we deduce PDET (PFA )
which can be used to compute ROC (Receiver Operating
Characteristics) curves:

e | PL ), PFA )
PDET (TM (P
Y
1
1
≈1−
ΦFλk (L[(1 − (1 − PFA ) η )− L − 1], 2, 2L) (21)
k∈Ω

7

B. TF .

D. HC⋆ and BJ.

e | PL looks for
Under H0 , the Fisher’s test TF applied to P
the maximum of identically distributed variables, which from
(9) and (12) correspond to the ratio of a F variable over a
sum of F variables. To our knowledge, there is no analytical
characterization of the resulting distribution for finite values
of L. Hence, although this test is CFAR, computing the
false alarm rate is problematic. Resorting to Monte Carlo
simulations to evaluate the function γ 7→ PFA is not possible,
owing to the limited number of available noise realizations.
Similar remarks can be made about standardized versions of
other tests like [21], [22], [27].

e | PL , the P -values involved in the HC
When applied to P
test (14) and in the BJ test (16) should be computed according
to the distribution under the null given by (9). Hence, (15) is


replaced by: 
P (νk )

 1 − ΦF
, 2, 2L , ∀k ∈ Ω,

P L (νk )


:=
vP|P
e L, k

P (νk )
N

 1 − ΦF
, 1, L , for k = 0, .
2
P L (νk )
The properties of the two tests are otherwise left unchanged,
with the CFAR property added: thanks to the standardization
of P by PL , the P -values are independent of the noise PSD.

C. TC .
e | PL , NC ), where parameter NC allows
We turn to TC (P
to focus on the regions of order statistics where deviations
under the alternative are likely to be significant. The PFA can
be obtained by observing that, owing to (9), the number K
of standardized ordinates larger than γ under H0 follows a
binomial distribution: K ∼ Bin(N/2 − 1, 1 − ΦF (γ, 2, 2L)),
whose CDF is: Pr (K ≤ k) = IΦF (γ,2,2L) (N/2 − NC , NC ),
see [59]. Using (13), (17) and noting
u := 1 − ΦF (γ, 2, 2L) =



L L
,
γ+L

e | PL can be expressed as:
the PFA of this test applied to P
e L , NC )> γ|H0)
e L , NC ), γ) := Pr (TC (P|P
PFA (TC (P|P
=1−

NX
C −1
k=0

= Iu (NC ,

Pr (K = k | H0 ) = 1 − I1−u (

N
− NC , NC )
2

N
− NC ),
2

(22)
where the last equation uses Ix (a, b) = 1−I1−x (b, a) (cf Prop.
e | PL ), this test is CFAR.
6.6.3 in [59]). As TM (P
As a side remark, note that the false alarm rate of the test
TC applied to the periodogram of a white noise of known
variance is obtained as above, by replacing ΦF with the CDF
of χ2 variables according to (4).
e | PL , NC ) can be
The function γ 7→ PDET (γ) of TC (P
deduced similarly to (22), with the difference that K is no
longer binomial owing to the {λk }. Appendix B shows that
e | PL , NC ),γ) := Pr (TC (P
e | PL , NC ) > γ|H1 )
PDET (TC (P
NX
C −1

≈1−

i 
X Y

i=0 Ω(i) ∈Ωi k=1

N

−1−i
2 Y
ϕFλ
1 − ϕFλ (γ,2,2L)
Ω

(i)
k

k′ =1

Ω(i)
k′

(γ,2,2L) ,

(23)
which can be used with (22) to compute ROC curves. The
non centrality parameters {λΩ(i) := λ(νΩ(i) )} and {λΩ(i) :=
k
k
k′
λ(νΩ(i) )} are given by (6) with the notation defined in (43).
k′
Note that the tests TM and TC are both working on order
statistics of P -values, which are Beta random variables, so in
both cases the PFA is given by the tail of a Beta distribution.

e |S
bE
VI. T ESTS APPLIED TO Z = P
Estimates of SbE different from the averaged periodogram
P L can be used for standardization in (8). Sec. I-B has
reviewed some methods to obtain such estimates. For the
eS
b E with P|P
e L , we opt for parametric
purpose of comparing P|
estimates allowing for an automatic parameter setting, as this
approach is commonly used in practice. For instance, when
using an estimate based on an AR process of order o, this
order can be estimated using many criteria, e.g. [78]–[82].
In our studies, we found that the selected order is often
different from the true order (as expected, see e.g. p. 211
of [78] and [83] for similar conclusions) but, as far as
detection results are concerned, these criteria have very similar
behaviour for sufficiently large N .
In any such method, let us denote respectively by obAR ,
{b
cj }j=,1,··· ,boAR and σ
b2 (b
oAR ) the selected order, AR coefficients and corresponding estimated prediction error variance.
The PSD estimate and resulting standardized periodogram are
SbE,AR (ν) :=

P (νk )
σ
b2 (b
oAR )
, Pe(νk |SbE,AR ) :=
.
o
b
AR
b
X
2
SE,AR (νk )
−2πijν
cj e
b
1+
j=1

(24)

Even if such approaches are relatively straightforward to
implement, characterizing the distribution of Pe (νk |SbE,AR ) is
more difficult than in the case of Pe (νk |P L ), owing essentially
to the stochastic nature of obAR in (24). In practice, the
‘whitening’ effect of SbE,AR is efficient because the selection
procedures have good fitting properties (they are approxi≈
mately consistent, i.e., ŜE,AR −→ SE as N −→ +∞).
This leads to consider as reasonable the assumptions that,
in effect, Sb P (νk(ν) ) ≈ SPE(ν(νkk)) and, with (4), that SPE(ν(νkk)) is
E,AR

k

approximately a χ22 /2 r.v. for k ∈ Ω and a χ21 r.v. for k = 0, N2 .
An approximate false alarm rate can then be evaluated from
these assumptions. For example, for the TC test applied to
P
, following the lines of (22) leads to:
b
S
E,AR

eS
b E,AR , NC ), γ) ≈ PFA (TC (P|S
e E , NC ), γ)
PFA (TC (P|
N
− NC ).
= Iu (NC ,
2
(25)
with u := Φχ22 (2γ) = e−γ . We will evaluate the reliability
of this approximation by numerical simulations in Sec.VII-D.

8

VII. N UMERICAL

as the estimation noise decreases with the increasing size of
the training data set.

STUDY

A. Simulation setting
We consider under H0 two PSD models for the noise E.
The first model comes from real RV data of the Sun obtained
from the GOLF spectrophotometer on board SoHO satellite
[84]. This instrument has been observing the Sun for 18 years
with a sampling rate of 20 s. As several gaps are present in
the resulting time series, we selected some (158) regularly
sampled data blocks of T ≈ 23 days, of which we averaged
and smoothed the periodograms (see Fig. 1, left panel). (Note
that the data we used are filtered at low frequencies so that
the resulting PSD estimate does probably not accurately reflect
the solar PSD at low frequencies).
The second model (Fig. 1, right panel) corresponds to a zeromean second-order stationary Gaussian AR(6) process. The
coefficients were chosen to yield a correlated process exhibiting higher energy at low frequencies and local variations, as
in some stars, but the choice of this PSD is not intended to
reflect the reality of a particular star.
AR(6)

log P(ν) [m2 .s−2 .H−1
z ]

GOLF
10

1

10

0

10

-1

10

-2

10

-3

10

-4

10

2

10

0

10

10 -6

10 -5

ν [Hz ]

10 -4

e | PL ) and TC (P
e | PL )
Fig. 2. Theoretical vs empirical results for TM (P
by MC simulations. Parameters: ∆t = 30 min, Ns = 3, αq = 0.2 m.s−1
for the three RV signatures of respective periods 11 h, 2.45 d and 6.61 d. The
curves in dashed lines show the inconsistency of ignoring noise correlations
(Sec. VII-C).

-2

10 -4

10 -3

ν [Hz ]

Fig. 1.
Left: Estimated PSD of the solar noise with part of GOLF data
(blue) and periodogram of one of the 158 data blocks (grey). Parameters:
N = 1110, ∆t = 30 min. This PSD is used to generate the noise for Figs.
2 and 4. Right: Theoretical PSD of the AR(6) noise (blue) and one noise
periodogram (grey). This PSD is used for Fig. 3. Parameters:N = 1024,
∆t = 1 min.

Under H1 , several cases of exoplanetary RV signatures will
be considered. The tests’ performances will be illustrated by
ROC curves representing PDET as a function of PFA . For
Monte Carlo (MC) simulations, 104 realizations have been
used.
e | PL
B. Analytical expressions for tests based on P
e | PL ).
e | PL ) and TC (P
We first consider the tests TM (P
The first panel of Fig.2 compares to empirical results obtained
from MC simulations to the expressions obtained for the
PFA (γ) for both tests (see (18) and (22)). The second panel
regards the corresponding expressions for PDET (γ) ((19) and
(23)) and the last panel the expression (21) for PDET (PFA ).
All the theoretical expressions are shown by color dots and
the empirical results from MC simulations are plotted in full
lines. Different values of L are considered to illustrate the
improvement brought by larger training data sets. The figure
shows a fair agreement between theoretical and empirical
results, even for the not so large value of N considered here
(N = 1110). The test performances logically increase with L,

C. Effects of standardization by PL
We compare the detection performances of the tests TM
e L ) with
and TC applied to the standardized periodogram (P|P
their unstandardized versions TM (P) and TC (P), as described
in Sec. IV-B. We evaluate first how accurate would be the false
alarm rate obtained by neglecting noise correlation for tests TC
and TM . For this we assume the detectors consider the noise
2
is white and have knowledge of σE
, the exact variance of E.
2
When Z = 2P/σE , it is easy to show that the false alarm
rates assumed by these two tests are given by

N
 PFA TM (2P/σ 2 ), γ  = 1 − Φ 22−1 (γ),

E
χ2

N

2
 PFA TC (2P/σE
, NC ), γ = IΦχ2 (γ) (NC ,
− NC ).
2
2
(26)
These expressions are compared to the true false alarm rates in
the top left panel of Fig. 2. The cyan dashed and blue dotted
curves show respectively approximation (26) and empirical
PFA for TM , while the yellow dashed and red dotted curves
show respectively approximation (26) and empirical PFA for
TC . Clearly, the correspondence between the thresholds values
and the target false alarm rates is destroyed because of noise
correlation in absence of standardization.
b E,AR
D. Effects of standardization by S

As discussed in Sec.VI, a way to deal with the frequency
dependence of the noise is to estimate its PSD by parametric

9

models. A difficulty with such methods is the injection of
estimation noise in the detection process. A standard approach
discussed in Sec.VI is to consider that the estimates are
sufficiently accurate for their intrinsic error to be negligible.
We study the performances of this approach here. For this we
consider the case of five sinusoidal signals with frequencies
falling into a ‘valley’ of the noise PSD (Fig.1, right panel). We
assume the noise PSD follows an AR model (which is indeed
the case here) and we estimate SbE,AR as described in Sec.VI.
e |S
b E,AR , and in
The question is the reliability of tests using P
particular how accurate is expression (25) with this approach.
Fig.3 compares, as a function of the test threshold, the PFA
assumed by approximation (25) (blue curve) with the actual
false alarm rates obtained for 1000 MC simulations.

e L , NC ) and TM (P|P
e L ) vs adaptative approaches
E. TC (P|P

We do not attempt to apply HC⋆ and BJ to P in the
correlated case as this leads to the same inconsistencies as
those illustrated in Sec. VII-C and Fig. 2, top left panel, dotted
curves. We compare here the adaptative approaches HC⋆ and
BJ standardized by PL as in Sec.V-D with the corresponding
tests TM and TC for L = 1 and L = 50. (Note that when
NC = 1, TC reduces to TM , cf (13)).
For this comparison, we consider two different cases of
Keplerian signals under H1 , which are typical of ‘super-Earth’
planets (Fig. 4): Np = 1 planet with null eccentricity in Case
1 (left column), and Np = 5 eccentric planets in Case 2 (right
column).
Case 1

1

PDET

0.8

0.8

0.6

0.6
e 1)
TM (P|P
e 1 , NC = 2Np )
TC (P|P
⋆ e
HC (P|P1 )
e 1)
BJ(P|P
e 50 )
TM (P|P
e 50 , NC = 2Np )
TC (P|P
⋆ e
HC (P|P50 )
e 50 )
BJ(P|P

0.4

0.2

0

b E,AR (L) was obIn each such simulation, an estimate S
tained with Akaike’s Final Prediction Error (FPE) [78] from
L noise time series and used to calibrate the periodogram. For
e |S
b E,AR , NC =
each such estimate, the true PFA of test TC (P
Ns ) was evaluated using 100 MC simulations. The figure
plots, respectively for L = 1, 20 and 100 the average PFA ,
respectively in black, green and red solid lines with dots.
We see that (25) is accurate, on average, only when L is
large. This figure also indicates the variability of the true
false alarm rate. For each value of L, the figure shows the 3σ
dispersion of the true PFA w.r.t. its empirical average (shaded
regions in grey for L = 1, green for L = 20 and red for
L = 100). Even when N is large, the true significance levels
at which such tests are conducted can undergo wild (and in
practice unknown) variations. The right panel is a zoom on
the 3σ region for L = 100. For a threshold γ = 5.9 for
instance, the PFA ≈ 0.013 from (25). In the right panel,
we see that the true false alarm rate for this threshold varies
in reality in the range [0.001 0.3]. For smaller values of L,
the excursions of the true false alarm rates are so large that
(25) is simply useless. Conclusions drawn from tests based
on parametric estimation may thus be very hazardous, even
when the parametric model is true and when large data sets
are available for PSD estimation.

P(ν) [m2 /s2 /Hz]

0

Fig. 3.
Comparison of the approximated PFA (25) (blue curve) with
true PFA for L = 1, 20 and 100. The solid lines with dots represent the
average actual PFA . The shaded regions are the corresponding 3σ enveloppes.
The right panel is a zoom in log-scale on the violet square in left panel.
Parameters: Ns = 5, αq = 0.07 m.s−1 , fq = [5.0, 5.5, 5.75, 6, 6.5] mHz.

Case 2

1

0.2

0.4

0.6

0.8

1

PFA

0.4

0.2

0
0

0.2

0.4

0.6

0.8

1

PFA
Signal for L = 1
Signal for L = 50

10

100

0

1.5

2

2.5

ν [µHz]

3

3.5

5

10

15

20

25

30

ν [µHz]

Fig. 4.
Top panels: Empirical ROC curves comparing classical and
adaptative tests for different signals. Case 1: Np = 1, Mp = 3.5M⊕
(for L = 1), Tp = 5.7813 d (signal frequency 1/Tp on-grid), ep = 0,
ωp = 0, T0 = 0, γ0 = 0. For L = 50, Mp = 0.8M⊕ . Case
2: Np = 5, Mp = [0.15, 0.15, 0.25, 0.25, 0.25]M⊕ (for L = 1),
Tp = [11.21 h, 1.33 d, 2.45 d, 6.91 d, 9.25 d] (signal frequencies off-grid),
ep = 0.9, ωp = π, T0 = 0, γ0 = 0. For L = 50, Mp (Tp = 11.21 h) =
0.07M⊕ . In the adaptative tests, α0 = 1/2. Bottom panels: periodograms
(logscale) of the signals under H1 .

The lower panels illustrate the periodogram of the (noiseless) Keplerian signatures for the two cases. We see the apparition of significant harmonics in the case of off-grid signal
frequencies and highly eccentric orbits. For each case, the
planet masses have been adapted depending on the considered
value of L for a better display of the ROC curves.
From L = 1 (dashed lines) to L = 50 (solid lines), the
performances of all tests increase in both cases. When the
signal is extremely sparse in the Fourier domain (Case 1),
TM is more powerful than the considered adaptive approaches.
This situation changes when the spectrum is less sparse (Case
2, compare the bottom panels), which is expected.
Results of [49] show that for a proportion of deviations in
1
1
the range [( N2 − 1) 4 ( N2 − 1) 2 ] an adaptive procedure such
as HC may have better asymptotic power than TM . For the
case N = 1110 considered here, this corresponds to the range
[5 23] (considering Ω). The situation should be opposite for
very sparse signals (in the range [1 4] here). It turns out that

10

the superiority of adaptive procedures is confirmed in Fig.
4, where the spectrum is 1-sparse in Case 1 and about 10sparse in Case 2. Note, however, that while the theory of [49]
may be used as a guideline for guessing the sparsity range in
which each test should work better, we should not expect a
too tight agreement with this theory. Indeed, in the framework
of [49], all deviations under the alternative have the same
amplitude, while RV signals lead to different amplitudes in
general. Moreover, these theoretical results are asymptotical
(while N is not so large here).
An interesting point is the comparison of HC⋆ and BJ with
e L , NC ), for which NC is a proxy for the number
TC (P|P
of significant deviations in the Fourier spectrum (here we
considered NC = 1 in Case 1 and NC = 10 in Case 2).
TC represents a kind of Oracle, which knows in which region
of the P -values to look at in order to ‘make the case’ against
H0 . The right panel shows that adaptive procedures have better
power than TC , yet without prior knowledge.

probability of detection of 0.9 while ensuring a false alarm rate
of 0.01. This situation is indicated by the black square. With
only L = 5 training time series, the probability of detection
would fall to about 0.1, all other parameters equal (red square).

F. A detectability study

Fig. 5. Example of a detectability study for a planet (relevant Keplerian
parameters: K = 0.54 m.s−1 , Tp = 3.23 d and e = 0; Mp = 1.241M⊕ ,
orbit inclination of 90°, semi-major axis a = 0.0425 Astronomical Unit)
orbiting a G2 type star. The plot shows the achievable PDET for different
configurations of PFA budgets and numbers of available training light curves.

A direct application of the previous results is detectability
studies, which can be used for the design of observational
strategies for instance. To illustrate this, we consider under
H1 a planet with the parameters of α Centauri B’s exoplanet
candidate as estimated in [16] (see the legend of Fig.5). As the
eccentricity is supposed null, the signal can be considered very
sparse in the Fourier domain. We consider the TM test and a
time sampling of ∆t = 4 hours. It was allowed to slightly vary
from one value of N to another in order to guarantee that the
planet’s period yields a frequency exactly on the Fourier grid,
in which case the spectrum is 1-sparse on Ω and TM is the
test that yields the best performances.
As for the noise PSD, we used a model based on HD
simulations of a star with similar spectral type as that of
αCenB. There is some mismatch between the simulated RV
time series of the spectral type (G2) and the true spectral
type (K1) of αCenB. Because spectral type affects the noise
properties [85], these results should not be considered to reflect
perfectly the case of the candidate planet orbiting αCenB.
In Fig.5, we illustrate the feasible performance compromises
(PDET , PFA ) as a function of N , for three target PFA (0.5, 0.1
and 0.01, indicated respectively by the dotted, solid and dashed
lines) and for different sizes of available training data sets L
(∞, 100, 20, 5, shown respectively in black, green, blue and
red). These curves were built using the expressions (20) and
(21) for the TM test and we checked that they are accurate
in separate MC simulations (not shown). In the case L →
+∞, we use the fact that F (2, 2L) −→ χ22 to calculate the
L→+∞

theoretical PDET (PFA ).
The study presented in Fig.5 allows to quantify interesting
facts. First, of course, PDET is larger if the allowed PFA is
larger. Second, for a fixed PFA , PDET is larger for a larger
value of L. Going to specific cases, we see that if a planet
similar to the considered candidate was orbiting a star of the
considered spectral type (G2), of which 100 training time
series are available, it would require 250 days (1500 × 4h)
of observations with 1 point every 4 hours to guarantee a

1

PFA = 50%, L → ∞
L = 100
L = 20
L=5
PFA = 10%, L → ∞
L = 100
L = 20
L=5
PFA = 1%, L → ∞
L = 100
L = 20
L=5

0.9
0.8

PDET

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
2
10

10

3

1500

10

4

N

VIII. S UMMARY

AND PERSPECTIVES

This paper has investigated the possibility of using training
data sets to standardize periodograms in order to improve
the control of the resulting false alarm rate. The paper first
provided an extended (though unavoidably selective) overview
of classical and recent techniques in sinusoid detection, with
emphasis on the problem of designing CFAR detectors for the
composite hypotheses of multiple sinusoids and colored noise.
We proposed an asymptotic analysis of the periodograms
statistics after standardization for a model involving an unknown number of sinusoids with unknown parameters in
partially unknown colored noise. This analysis allowed to
characterize the performances of some standardized tests in
terms of false alarm and detection rates. We showed that
when standardization is performed with a simple averaged
periodogram as a noise PSD estimate these tests are CFAR
for all sizes of training data sets. In contrast, we pointed out
that standardization based on parametric estimates of the noise
PSD may present actual false alarm rates that may be very far
from the assumed ones, even for large data sets.
The tests we considered include classical approaches and also
more recent adaptive tests designed for the rare and weak
setting. For the latter tests, the standardization (by PL ) offers
the same benefits as if the statistics of the noise were known
a priori, with the CFAR property added. We also showed that
such tests can present better power comparing to procedures
for which the number of sinusoids would be known in advance.
In practical situations, some of the assumptions (i-iv) in Sec.
I-A may not be met. This can be the case in exoplanet detection using RV time series, owing for instance to astrophysical
effects linked to magnetic activity (like spots, affecting (i)) or
instrumental defects / observational constraints (affecting (ii)

11

and (iii)). Comparing theory to pratice, the present study is
useful to make feasibility studies (as in Sec. VII-F), which then
describe best achievable performances (i.e., around “quiet”
stars and in absence of other unmodelled perturbations) for
a regular sampling.
An important extension to the considered framework regards
the case when, for periodicity analysis, the time series is not
correlated to orthogonal exponentials. This case encompasses
situations where (a) the sampling is irregular, (b) P is not evaluated on the Fourier grid (as in oversampled periodograms)
(c) P is modified in the form of “generalized periodograms”,
which correlate the time series with highly redundant dictionaries of specific features [28], [86]–[93]. In such cases, the
considered ordinates exhibit strong dependencies. With the
additional complication of partially unknown colored noise,
analytical evaluation of the false alarm rate for the considered
tests seems out of reach. We conducted however preliminary
studies suggesting that it might be possible to obtain accurate
estimates of the false alarm rate when the noise is colored, the
sampling irregular, the considered frequencies not restricted to
Fourier grid and N small, by combining the standardization
proposed in the present paper with bootstrap procedures [94]
and maximum likelihood estimation of Generalized Extreme
Values distributions’ parameters [87], [95]. This numerical
approach, still under study, would allow to address important
questions, like that of the impact of the sampling distribution
on the detection performances.

OF EXPRESSIONS

(5) AND (6)

We prove here that for model (3) the periodogram is asymptotically distributed as in (5) with non centrality parameters as
in (6). The proof is adapted from Theorem 6.2 of [24], which
considers the complex case. We first prove (5) and then turn
to (6). The time series of model (3) can also be written as
X(j) =

Ns
X

αq sin(2πfq j + ϕq ) + E(j) = R(j) + E(j),

q=1

PNs
with R(j) := q=1
αq sin(2πfq j + ϕq ) a deterministic part,
which using Euler formulae can be written as
R(j) =

Ns
X
αq
q=1

2

e

i(ϕq − π
2 ) 2πifq j

e

αq −i(ϕq + π ) −2πifq j
2 e
e
−
.
2


i2πν
, . . . , ei2πN ν ]⊤ ,

 f (ν) := [e
π
By introducing:
f + (fq ) := ei(ϕq − 2 ) f (fq ),

π
 −
f (fq ) := e−i(ϕq + 2 ) f (fq ),
the time series writes in vector form :
Ns

X
αq  +
X=
f (fq ) − f − (fq ) + E = R + E
2
q=1

and ǫk :=

1 H
f (νk )E.
N

(28)

Because E is a zero mean Gaussian process, the random
variable yk = µk + ǫk is Gaussian with mean µk and variance
σk2 . This is a complex variable for all Fourier frequencies
except for ν0 and ν N because f (ν0 ) and f (ν N ) are real.
2
2
The distribution of the periodogram requires to investigate
the variance of ǫk which, with (28), writes:
N var ǫk =


1
E f H (νk )EE⊤ f (νk )
N

N
1 X
rE (t − s)e−i2πνk (t−s)
N t,s=1
1 X
(N − |u|)rE (u)e−i2πνk u
=
N

=

|u|<N

= SE (νk ) −

X |u|
X
rE (u)e−i2πνk u
rE (u)e−i2πνk u −
N

|u|<N

|u|≥N

= SE (νk ) + O(rN ),

since,
using the absolutely summable autocorrelation function
P
u |rE (u)| < ∞ and the dominated convergence theorem we
have
X
|u|
)|rE (u)| → 0 as N → ∞.
(29)
rN :=
min(1,
N
u

1 H
1
1
f (νk )X = f H (νk )R + f H (νk )E.
N
N
N

σk2 := var ǫk = N −1 SE (νk ) + O(N −1 rN ).

(30)

By Lemma 12.2.1(b) of [24], we obtain

1 2


∀k ∈ Ω,
 2 χ2,2 |µk |2 ,

σ2
2
2
k
|yk | /σk ∼
N

2


 χ1, |µk |2 , for k = 0, 2 .
2
σ
k

Hence, for the periodogram this implies

P (νk |H1 ) = N |yk |2 = N σk2 (|yk |2 /σk2 )

1
2

 ρ−1
k SE (νk )χ2,2ρk γk , ∀k ∈ Ω,
2
∼
N

 ρ−1 SE (νk )χ2
,
1,ρk γk , for k = 0,
k
2
where
ρk := SE (νk )/(N σk2 )

(31)

and γk := N |µk |2 /SE (νk ). (32)

With (30), we see that
ρk =
(27)

and its discrete Fourier transform (DFT) yk at frequency νk
is
yk =

1 H
f (νk )R
N

µk :=

Hence, for all Fourier frequencies,

A PPENDIX A
D ERIVATION

The DFT is composed of a deterministic part, µk , and a
stochastic part, ǫk , defined as

SE (νk )
= 1 + O(rN )
SE (νk ) + O(rN )

(33)

for all Fourier frequencies. Owing to (29), an approximated
distribution of (31) can be obtained by neglecting the O(rN )
in (33). The distribution (5) follows by noting
λk := 2γk for k ∈ Ω

and λk := γk for k = 0,

N
. (34)
2

12

We now turn to the computation of the non centrality parameters. We have from (27), (28) and (32)
γk =
=

N
1
| f H (νk )R|2
SE (νk ) N
Ns
N X
X
1
αq 
N SE (νk )

j=1 q=1

2

γk =

Ns
Ns h
i
X
N X
α2q κ2q+2αq κq
αℓ κℓ cos(θq − θℓ ) (41)
4SE (νk )q=1
ℓ=q+1

i(ϕq − π
) 2πi(fq −νk )j
2

e

(35)

e

π

. . . −e−i(ϕq + 2 ) e−2πi(fq +νk )j



2

.

Introducing the Dirichlet Kernel (cf Lemma 12.1.3 of [24]):
DN (ν) :=

Consequently, the expression of the {γk } in (37) becomes

N
1 X i2πνj
sin(N πν) i(N +1)πν
e
=
e
,
N j=1
N sin(πν)

and the corresponding Fejér kernel (or spectral window)
!2
sin(N
πν)
,
(36)
KN (ν) := |DN (ν)|2 =
N sin(πν)

and the non centrality parameters {λk } of (6) follow from
(34),with κq and θq given by (39) and (40).

Note that if all signal frequencies {fp } fall on the Fourier
frequency grid, the crossed term in (41) vanish owing to the
orthogonality of the Fejér kernels centered at different signal
frequencies. In this case, expression (41) precisely reduces to
expression given in Remark 6.6 of [24].
We finally wish to mention that the expression of the non
centrality parameters is erroneously reported in exp. (5) of [55]
(sign error and crossed terms missing).

A PPENDIX B
D ERIVATION OF EXPRESSION (23)

expression (35) becomes
γk =

Ns

X
π
N
αq DN (fq − νk )ei(ϕq − 2 )
4SE (νk ) q=1

2
π
. . . − DN (fq + νk )e−i(ϕq + 2 ) .

This equation can also be written as :
γk =

Ns
X
2
N
αq zq (νk ) ,
4SE (νk ) q=1

(37)

e | PL larger than
Let K denote the number of ordinates of P
γ under H1 , and pi := Pr (K = i | H1 ). From the definition
(13), we have:
e L ), γ, NC ) := Pr (TC (P|P
e L , NC )> γ|H1)
PDET (TC (P|P
= Pr (K ≥ NC | H1 )
=1−

with

NX
C −1
i=0

Pr (K = i | H1 ) = 1 −

NX
C −1

pi .

i=0

(42)
e L )i := Pe(νi ) has probaOwing
to
(10)
each
ordinate
(
P|P
P L (νi )
= x+ eiθ+ − x− eiθ−
bility 1 − ΦFλi (γ) to be larger than γ. These variates can be
= (x+ cos θ+ − x− cos θ− ) + i(x+ sin θ+ − x− sin θ− ), considered approximately independent but not i.i.d. Hence, the
(38) variable K is not binomially distributed (as it is under H ) and
0
where
the
probabilities
{p
}
require
further
investigation.
We
proceed
i

sin(N π(fq − νk ))
by induction. In the following, all probabilities are under H1 .


,
x+ = x+ (νk , q) :=


N
sin(π(f
−
ν
))
The
first probability p0 can simply be approximated as

q
k




sin(N
π(f
+
ν
))
N
 N
q
k
 x = x (ν , q) :=
,
−1
−1
−
−
k
2
 2\

Y
N sin(π(fq + νk ))
e

P
P
p
=
Pr
(
ΦFλk .
|
)
≤
γ
≈
0
L k
π





θ+ = θ+ (νk , q) := +[(N + 1)π(fq − νk ) + (ϕq − )],

k=1
k=1

2



 θ = θ (ν , q) := −[(N + 1)π(f + ν ) + (ϕ + π )].
−
−
k
q
k
q
The probability p1 = Pr(K = 1) is similarly
2

The modulus κq of zq may be written as

N
−1 
2[

\
1

e | PL )k > γ
e | PL )j ≤ γ
(39)
(P
(P
p1 = Pr
κq := |zq | = x2+ + x2− − 2x+ x− cos (θ+ − θ− ) 2 ,


k=1
j6=k

2
x
=
K
(f
−
ν
),

N
N q
k
N
+

−1
2 −1 h
2
i
X
Y
with
x2− = KN (fq + νk ),
.
Φ
(1
−
Φ
)
≈
F
F
λ
λ

j
k

j=1,
θ+ − θ− = 2π(N + 1)fq + 2ϕq ,
k=1
π

π

zq (νk ) := DN (fq − νk )ei(ϕq − 2 ) −DN (fq + νk )e−i(ϕq + 2 )

j6=k

θq := ∠ zq , θq ∈ ] − π, π],

and

(40)

the phase of zp obtained from the real and imaginary parts of
(38) [96]. With these notations, it is easy to show that
Ns
X
q=1

zq

2

Ns
Ns h
i
X
X
2
κq +2κq
κℓ cos(θq −θℓ ) for Ns > 1.
=
q=1

ℓ=q+1

To generalize further, denote by Ω(i) one particular com(i)
bination of i indices taken in Ω and Ω
:= Ω\Ω(i)
(i)
(i)
the set of remaining indices. Let {Ω1 , . . . , Ωi } (resp.
(i)
(i)
{Ω1 , . . . , Ω N −1−i }) denote the indices in two such combi2
nations, and let Ωi be the set of all the {Ω(i) }. With these

13

notations we obtain for i > 1 :

i
[ \
e | PL ) (i) > γ
(P
pi = Pr
Ωk

(i)
i

[28] E. Bölviken, “New tests of significance in periodogram analysis,”
Scandinavian J. Stat., vol. 10, no. 1, pp. 1–9, 1983.
E. Bölviken, “The distribution of certain rational functions of order
\
e | PL ) (i) ≤ γ
statistics from exponential distributions,” Scandinavian J. Stat., vol. 10,
(P
Ωk′

no. 2, pp. 117–123, 1983.
k=1
k′ =1
Ω ∈Ω
[30] R. Von Sachs, “Estimating the spectrum of a stochastic process in the
N
−1−i
presence of a contaminating signal,” IEEE Trans. Signal Process., vol.
i 
 2Y
XY
41, no. 1, pp. 323, 1993.
ΦFλ (i)(γ,2,2L).
≈
1− ΦFλ (γ,2,2L)
[31] R. Von Sachs, “Peak-insensitive non-parametric spectrum estimation,”
(i)
Ω ′
Ω
k
k
k′ =1
Ω(i) k=1
J. Time Series Anal., vol. 15, no. 4, pp. 429–452, 1994.
(43) [32] R.J. Bhansali, “A mixed spectrum analysis of the lynx data,” J. R. Stat.
Soc. Ser. A, , no. 142, pp. 199–209, 1979.
Expression (23) follows by combining (42) and (43).

[33] B. Truong-Van, “A new approach to frequency analysis with amplified
harmonics,” J. R. Stat. Soc. Ser. B, vol. 52, pp. 203–221, 1990.
ACKNOWLEDGEMENT
[34] B.G. Quinn and J.M. Fernandes, “A fast efficient technique for the
estimation of frequency,” Biometrika, vol. 78, no. 3, pp. 489–497, 1991.
The GOLF instrument onboard SoHO is a cooperative effort [35] B.G. Quinn, “A fast efficient technique for the estimation of frequency:
of scientists, engineers, and technicians, to whom we are
Interpretation and generalisation,” Biometrika, vol. 86, pp. 213–220,
1999.
indebted. SoHO is a project of international collaboration
[36] L. Kavalieris and Hannan E.J., “Determining the number of terms in a
between ESA and NASA.
trigonometric regression.,” J. Time Series Analysis, vol. 15, pp. 613–
625, 1994.
[37] E.J. Hannan, “Testing for a jump in the spectral function,” J. R. Stat.
R EFERENCES
Soc. Ser. B, vol. 23, no. 2, pp. 394–404, 1961.
[38] D.F. Nicholls, “Estimation of the spectral density function when testing
[1] A. Schuster, “On the investigation of hidden periodicities,” J. Geophys.
Res., vol. 3, pp. 13, 1898.
for a jump in the spectrum,” Austr. J. Stat., vol. 9, pp. 103–108, 1967.
[2] D.R. Brillinger, Time Series : Data Analysis and Theory, Holden Day,
[39] S.T. Chiu, “Peak-insensitive parametric spectrum estimation,” Stochastic
San Francisco, 1981.
Processes and their Applications , vol. 35, pp. 121–140, 1990.
[3] M.S. Bartlett,
“Periodogram analysis and continuous spectra,”
[40] J.K. Gryca, “Detection of multiple sinusoids buried in noise via balanced
Biometrika, vol. 37, pp. 1–16, 1950.
model truncation,” IEEE Instrum. Meas. Conf., pp. 1353–1358, 1998.
[4] U. Grenander and M. Rosenblatt, Statistical analysis of stationary time
[41] L.B. White, “Detection of sinusoids in unknown coloured noise using
series, John Wiley and Sons, 1957.
ratios of ar spectrum estimates,” Proc. Inform., Decision and Control,
[5] M.B. Priestley, Spectral Analysis and Time Series, Academic Press,
pp. 257–262, 1999.
[42] N. Lu and D. Zimmerman, “Testing for directional symmetry in spatial
San Diego, 1981.
dependence using the periodogram,” J. Stat. Planning and Inference,
[6] P.J. Brockwell and R.A. Davis, Time series : theory and methods,
Springer, 1991.
vol. 129, pp. 369–385, 2005.
[7] P. Bloomfield, Fourier Analysis of Time Series , Wiley-Intersci., 2000.
[43] A.P. Liavas et al., “A periodogram-based method for the detection of
[8] P. Stoica and R. Moses, Spectral analysis of signals, Prentice Hall,
steady-state visually evoked potentials.,” IEEE Trans. Biom. Eng., vol.
2005.
45, no. 2, pp. 242–248, 1998.
[9] B. G. Quinn and E.J. Hannan, The Estimation and Tracking of [44] C. Zheng, “Detection of multiple sinusoids in unknown colored noise
Frequency, Cambridge Univ., 2001.
using truncated cepstrum thresholding and local signal-to-noise-ratio,”
[10] F. Pepe et al., “Instrumentation for the detection and characterization of
Applied Acoust., pp. 809–816, 2012.
exoplanets,” Nature, vol. 513, pp. 358–366, 2014.
[45] B. Nadler and A. Kontorovich, “Model selection for sinusoids in noise:
[11] N.M. Batalha, “Exploring exoplanet populations with NASA’s Kepler
Statistical analysis and a new penalty term,” IEEE Trans. Signal Process.,
vol. 59, no. 4, pp. 1333–1345, 2011.
Mission,” Proc. Nat. Acad. Sci., vol. 111, pp. 12647–12654, Sept. 2014.
[12] M. Auvergne et al., “The CoRoT satellite in flight: description and
[46] C. Koen, “The analysis of indexed astronomical time series - xi. the
performance,” A&A, vol. 506, pp. 411–424, 2009.
statistics of oversampled white noise periodograms,” MNRAS, vol. 449,
[13] H. Rauer et al., “The PLATO 2.0 mission,” Experimental Astronomy,
no. 1, pp. 1098–1105, 2015.
vol. 38, pp. 249–330, 2014.
[47] C. Koen, “The analysis of indexed astronomical time series - xii. the
[14] D.A. Fischer et al., “Exoplanet Detection Techniques,” Protostars and
statistics of oversampled fourier spectra of noise plus a single sinusoid,”
Planets VI, pp. 715–737, 2014.
MNRAS, vol. 453, no. 2, pp. 1793–1798, 2015.
[15] M. Perryman, The exoplanet handbook, Cambridge Univ., 2011.
[48] M. Tuomi et al., “Signals embedded in the radial velocity noise. Periodic
[16] X. Dumusque et al., “An Earth-mass planet orbiting α Centauri B,”
variations in the τ Ceti velocities,” A&A, vol. 551, pp. A79, 2012.
[49] D. Donoho and J. Jin, “Higher criticism for detecting sparse heterogeNature, vol. 491, pp. 207–211, 2012.
neous mixtures,” Ann. Stat., 2004.
[17] A. Hatzes, “The Radial Velocity Detection of Earth-mass Planets in the
[50] Y.I. Ingster et al., “Detection boundary in sparse regression,” Electron.
Presence of Activity Noise: The Case of α Centauri Bb,” ApJ, vol. 770,
pp. 133, 2013.
J. Statist., vol. 4, pp. 1476–1526, 2010.
[18] V. Rajpaul et al., “Ghost in the time series: no planet for Alpha Cen
[51] G. Walther, “The Average Likelihood Ratio for Large-scale Multiple
B,” MNRAS, vol. 456, pp. L6–L10, 2016.
Testing and Detecting Sparse Mixtures,” IMS Collections : From
[19] L. Bigot et al., “The diameter of the CoRoT target HD 49933,” A&A,
Probability to Stat. and Back: High-Dimensional Models and Processes ,
vol. 534, no. 3, 2011.
vol. 9, pp. 317–326, 2012.
[20] R.A. Fisher, “Tests of Significance in Harmonic Analysis,” Proc. R. [52] A. Moscovich et al., “On the exact Berk-Jones statistics and their p-value
Soc. London, Ser. A, vol. 125, pp. 54–59, 1929.
calculation,” Electron. J. Stat., vol. 10, pp. 2329–2354, 2016.
[21] S.T. Chiu, “Detecting periodic components in a white gaussian time
[53] V. Gontscharuk et al., “Goodness of fit tests in terms of local levels
series,” J. R. Stat. Soc. Series B, vol. 51, no. 2, pp. 249–259, 1989.
with special emphasis on higher criticism tests,” Bernoulli, vol. 22, no.
[22] M. Shimshoni, “On fisher’s test of significance in harmonic analysis,”
3, pp. 1331–1363, 2016,.
Geophys. J. R. Astronom. Soc., pp. 373–377, 1971.
[54] P. Hall and J. Jin, “Innovated higher criticism for detecting sparse signals
[23] S. Kay, “Adaptive detection for unknown noise power spectral densities,”
in correlated noise,” Ann. Stat., vol. 38, pp. 1686–1732, 2010.
IEEE Trans. Signal Process, vol. 47, no. 1, pp. 10–21, 1999.
[55] S. Sulis, D. Mary, and L. Bigot, “Using hydrodynamical simulations
[24] T.H. Li, Time series with mixed spectra , CRC Press, 2014.
of stellar atmospheres for periodogram standardization: Application to
[25] P. Whittle, “The simultaneous estimation of a time series harmonic
exoplanet detection,” in IEEE ICASSP, March 2016, pp. 4428–4432.
[56] S.K. Gupta et al., “UPSO three channel fast photometer,” Bulletin
components and covariance structure,” Trabajos de Estadistica, vol. 3,
Astron. Soc. of India, vol. 29, pp. 479–486, 2001.
no. 1-2, pp. 43–57, 1952.
[26] M.S. Bartlett, “An introduction to stochastic processes,” Quart. J. R. [57] S. Sulis, D. Mary, and L. Bigot, “Overcoming the stellar noise
Meteorological Soc., vol. 81, no. 350, pp. 650, 1955.
barrier for the detection of telluric exoplanets: an approach based on
hydrodynamical simulations,” in EAS Publications Series, D. Mary
[27] A. Siegel, “Testing for periodicity in a time series,” J. Amer. Stat.
Assoc., vol. 75, no. 370, pp. 345–348, 1980.
et al., Eds., Sept. 2016, vol. 78-79, pp. 247–274.
N
2

−1−i


 [29]

14

[58] J.G. Proakis and D.G. Manolakis, Digital Signal Processing, PrenticeHall, 1996.
[59] M. Abramowitz et al., Spectral Analysis and Time Series, Dover
Publications, 1972.
[60] H.A. David and H.N. Nagaraja, Order Statistics, 3rd Ed., Wiley, 2003.
[61] S.M. Kay, Fundamentals of Statistical signal processing. Vol II :
Detection theory., Prentice-Hall, Inc, 1998.
[62] B. G. Quinn, “Testing for the presence of sinusoidal components,” J.
Appl. Probability, vol. 23, pp. 201–210, 1986.
[63] A. Schwarzenberg-Czerny, “The distribution of empirical periodograms:
Lomb–Scargle and PDM spectra,” MNRAS, vol. 301, pp. 831–840,
1998.
[64] T. Aittokallio et al., “Testing for Periodicity in Signals: An Application
to Detect Partial Upper Airway Obstruction during Sleep,,” J. Theoretical
Medicine, vol. 3, no. 4, pp. 231–245, 2001.
[65] J. Gutiérrez-Soto et al., “Low-amplitude variations detected by CoRoT
in the B8IIIe star HD 175869,” A&A, vol. 506, pp. 133–141, 2009.
[66] S. Aldor-Noiman et al., “The power to see: A new graphical test of
normality.,” Am. Stat., vol. 68, no. 4, pp. 318–318, 2013.
[67] D. Mary and A. Ferrari, “A non-asymptotic standardization of binomial
counts in higher criticism,” in Inform. Theory (ISIT), IEEE Int. Symp.,
June 2014, pp. 561–565.
[68] D.M. Kaplan and M. Goldman, “True equality (of pointwise sensitivity)
at last: a dirichlet alternative to Kolmogorov-Smirnov inference on
distributions.,” Tech. report, 2014.
[69] V. Gontscharuk et al., “The intermediates take it all: Asymptotics of
higher criticism statistics and a powerful alternative based on equal local
levels.,” Biom. J., vol. 57, no. 1, pp. 159–180, 2014.
[70] J. Li and D. Siegmund, “Higher criticism: p-values and criticism,” Ann.
Statist., vol. 43, no. 3, pp. 1323–1350, 2015.
[71] R.H. Berk and D.H. Jones, “Goodness-of-fit test statistics that dominate
the Kolmogorov statistics,” Z. Wahrscheinlichkeit. , vol. 47, pp. 47–59,
1979.
[72] A. Moscovich-Eiger and B. Nadler, “Fast calculation of boundary
crossing probabilities for Poisson processes,” ArXiv e-prints (V.3), 2015.
[73] J. Fan and Q. Yao, Nonlinear Time Series-Nonparametric and Parametric
Methods, Springer-Verlag New York, 2003.
[74] K. F. Turkman and A. M. Walker, “On the asymptotic distributions
of maxima of trigonometric polynomials with random coefficients,”
Advances in Applied Probability , vol. 16, no. 4, pp. 819–842, 1984.
[75] R.A Davis and T. Mikosch, “The maximum of the periodogram of a
non-gaussian sequence,” Ann. of Prob., vol. 27, pp. 522–536, 1999.
[76] D. Rife and R. Boorstyn, “Single tone parameter estimation from
discrete-time observations,” IEEE Trans. Inf. Theory, vol. 20, no. 5,
pp. 591–598, 1974.
[77] B.G. Quinn and P. J. Kootsookos, “Threshold behavior of the maximum
likelihood estimator of frequency,” IEEE Trans. Signal Process., vol. 42,
no. 11, pp. 3291–3294, 1994.
[78] H. Akaike, “Fitting autoregressive models for prediction,” Ann. Inst.
Stat. Math., vol. 21, no. 1, pp. 243–247, 1969.
[79] H. Akaike, “A new look at the statistical model identification,” IEEE
Trans. Automatic Control, vol. 19, no. 6, pp. 716–723, 1974.
[80] E. Parzen, “Multiple time series: determining the order of approximating
autoregressive schemes,” Tech. Report, , no. 23, pp. 716–723, 1975.
[81] E.J. Hannan and B.G. Quinn, “The determination of the order of an
autoregression,” J. R. Stat. Soc. Ser. B, vol. 41, pp. 190–195, 1979.
[82] J. Rissanen, “Universal coding, information prediction and estimation,”
IEEE Trans. Inf. Theory, , no. 30, pp. 629–636, 1984.
[83] A. Boardman et al., “A study on the optimum order of autoregressive
models for heart rate variability,” Physio. Meas., vol. 23, pp. 325, 2002.
[84] R.A. Garcia et al., “Global solar Doppler velocity determination with
the GOLF/SoHO instrument,” A&A, vol. 442, pp. 385–395, 2005.
[85] N. Meunier et al., “Variability of stellar granulation and convective
blueshift with spectral type and magnetic activity. I. K and G main
sequence stars,” ArXiv e-prints, 2016.
[86] R.V. Baluev, “Assessing the statistical significance of periodogram
peaks,” MNRAS, vol. 385, pp. 1279–1285, 2008.
[87] M. Süveges, “Extreme-value modelling for the significance assessment
of periodogram peaks,” MNRAS, vol. 440, no. 3, pp. 2099–2114, 2014.
[88] J.D. Scargle, “Studies in astronomical time series analysis. II - Statistical
aspects of spectral analysis of unevenly spaced data,” ApJ, vol. 263, pp.
835–853, 1982.
[89] G.L. Bretthorst, “Frequency Estimation And Generalized Lomb-Scargle
Periodograms.,” Stat. Challenges in Astronomy, pp. 309–329, 2003.
[90] T. Thong et al., “Lomb-wech periodogram for non-uniform sampling,”
Proc. 26th Annu. Int. Conf. IEEE EMBS, 2004.

[91] M. Zechmeister and M. Kürster, “The generalised Lomb-Scargle
periodogram. A new formalism for the floating-mean and Keplerian
periodograms,” A&A, vol. 496, pp. 577–584, 2009.
[92] R.V. Baluev, “Keplerian periodogram for Doppler exoplanet detection:
optimized computation and analytic significance thresholds,” MNRAS,
vol. 446, pp. 1478–1492, 2015.
[93] P.C. Gregory, “An apodized kepler periodogram for separating planetary
and stellar activity signals,” MNRAS, vol. 458, pp. 2604–2633, 2016.
[94] A.M. Zoubir, “Bootstrap: theory and applications,” in SPIE Conf., F.T.
Luk, Ed., 1993, vol. 2027, pp. 216–235.
[95] M. Süveges et al., “A comparative study of four significance measures
for periodicity detection in astronomical surveys,” MNRAS, vol. 450,
no. 2, pp. 2052–2066, 2015.
[96] H.S. Kasana, Complex variables : theory and applications. 2nd Edition ,
Prentice-Hall of India, 2005.

