Usage

After installing the package, you can start using it with

using LombScargle

The module defines a new LombScargle.Periodogram data type, which, however, is not exported because you will most probably not need to directly manipulate such objects. This data type holds both the frequency and the power vectors of the periodogram.

The main function provided by the package is lombscargle:

lombscargle(times::AbstractVector{<:Real}, signal::AbstractVector{<:Real},
            [errors::AbstractVector{<:Real}]; keywords...)

lombscargle returns a LombScargle.Periodogram. The only two mandatory arguments are:

• times: the vector of observation times
• signal: the vector of observations associated with times

The optional argument is:

• errors: the uncertainties associated to each signal point.

All these vectors must have the same length.

If the signal has uncertainties, the signal vector can also be a vector of Measurement objects (from Measurements.jl package), in which case you need not to pass a separate errors vector for the uncertainties of the signal. You can create arrays of Measurement objects with the measurement function, see Measurements.jl manual at https://juliaphysics.github.io/Measurements.jl/stable for more details. The generalised Lomb–Scargle periodogram by Zechmeister and Kürster (2009) is always used when the signal has uncertainties, because the original Lomb–Scargle algorithm cannot handle them.

Planning the Periodogram

In a manner similar to planning Fourier transforms with FFTW, it is possible to speed-up computation of the Lomb–Scargle periodogram by pre-planning it with LombScargle.plan function. It has the same syntax as lombscargle, which in the base case is:

LombScargle.plan(times::AbstractVector{<:Real}, signal::AbstractVector{<:Real},
                 [errors::AbstractVector{<:Real}];
                 normalization::Symbol=:standard,
                 noise_level::Real=1,
                 center_data::Bool=true,
                 fit_mean::Bool=true,
                 fast::Bool=true,
                 flags::Integer=FFTW.ESTIMATE,
                 timelimit::Real=Inf,
                 oversampling::Integer=5,
                 padding_factors::Vector{Int}=[2],
                 Mfft::Integer=4,
                 samples_per_peak::Integer=5,
                 nyquist_factor::Integer=5,
                 minimum_frequency::Real=NaN,
                 maximum_frequency::Real=NaN,
                 frequencies::AbstractVector{Real}=
                 autofrequency(times,
                               samples_per_peak=samples_per_peak,
                               nyquist_factor=nyquist_factor,
                               minimum_frequency=minimum_frequency,
                               maximum_frequency=maximum_frequency))

autofrequency(times::AbstractVector{Real};
              samples_per_peak::Integer=5,
              nyquist_factor::Integer=5,
              minimum_frequency::Real=NaN,
              maximum_frequency::Real=NaN)

LombScargle.plan takes all the same argument as lombscargle shown above and returns a LombScargle.PeriodogramPlan object after having pre-computed certain quantities needed afterwards, and pre-allocated the memory for the periodogram. It is highly suggested to plan a periodogram before actually computing it, especially for the fast method. Once you plan a periodogram, you can pass the LombScargle.PeriodogramPlan to lombscargle as the only argument.

lombscargle(plan::PeriodogramPlan)

Planning the periodogram has a twofold advantage. First of all, the planning stage is type-unstable, because the type of the plan depends on the value of input parameters, and not on their types. Thus, separating the planning (inherently inefficient) from the actual computation of the periodogram (completely type-stable) makes overall computation faster than directly calling lombscargle. Secondly, the LombScargle.PeriodogramPlan bears the time vector, but the quantities that are pre-computed in planning stage do not actually depend on it. This is particularly useful if you want to calculate the False-Alarm Probability via bootstrapping with LombScargle.bootstrap: the vector time is randomly shuffled, but pre-computed quantities will remain the same, saving both time and memory in each iteration. In addition, you ensure that you will use the same options you used to compute the periodogram.

Fast Algorithm

When the frequency grid is evenly spaced, you can compute an approximate generalised Lomb–Scargle periodogram using a fast algorithm proposed by Press and Rybicki (1989) that greatly speeds up calculations, as it scales as $O[N \log(M)]$ for $N$ data points and $M$ frequencies. For comparison, the true Lomb–Scargle periodogram has complexity $O[NM]$. The larger the number of datapoints, the more accurate the approximation.

The only prerequisite in order to be able to employ this fast method is to provide a frequencies vector as an AbstractRange object, which ensures that the frequency grid is perfectly evenly spaced. This is the default, since LombScargle.autofrequency returns an AbstractRange object.

Since this fast method is accurate only for large datasets, it is enabled by default only if the number of output frequencies is larger than 200. You can override the default choice of using this method by setting the fast keyword to true or false. We recall that in any case, the frequencies vector must be a Range in order to use this method.

To summarize, provided that frequencies vector is an AbstractRange object, you can use the fast method:

• by default if the length of the output frequency grid is larger than 200 points
• in any case with the fast=true keyword

Setting fast=false always ensures you that this method will not be used, instead fast=true actually enables it only if frequencies is an AbstractRange.

Normalization

By default, the periodogram $p(f)$ is normalized so that it has values in the range $0 \leq p(f) \leq 1$, with $p = 0$ indicating no improvement of the fit and $p = 1$ a "perfect" fit (100% reduction of $\chi^2$ or $\chi^2 = 0$). This is the normalization suggested by Lomb (1976) and Zechmeister and Kürster (2009), and corresponds to the :standard normalization in lombscargle function. Zechmeister and Kürster (2009) wrote the formula for the power of the periodogram at frequency $f$ as

\[p(f) = \frac{1}{YY} \left[ \frac{YC^2_{\tau}}{CC_{\tau}} + \frac{YS^2_{\tau}}{SS_{\tau}} \right]\]

See the paper for details. The other normalizations for periodograms $P(f)$ are calculated from this one. In what follows, $N$ is the number of observations.

• :model:

  \[P(f) = \frac{p(f)}{1 - p(f)}\]

• :log:

  \[P(f) = -\log(1 - p(f))\]

• :psd:

  \[P(f) = \frac{W}{2}\left[\frac{YC^2_{\tau}}{CC_{\tau}} + \frac{YS^2_{\tau}}{SS_{\tau}}\right] = p(f) \frac{W*YY}{2}\]

  where W is the sum of the inverse of the individual errors, $W = \sum \frac{1}{\sigma_{i}}$, as given in Zechmeister and Kürster (2009).

• :Scargle:

  \[P(f) = \frac{p(f)}{\text{noise level}}\]

  This normalization can be used when you know the noise level (expected from the a priori known noise variance or population variance), but this isn't usually the case. See Scargle (1982)

• :HorneBaliunas:

  \[P(f) = \frac{N - 1}{2} p(f)\]

  This is like the :Scargle normalization, where the noise has been estimated for Gaussian noise to be $(N - 1)/2$. See Horne and Baliunas (1986)

• If the data contains a signal or if errors are under- or overestimated or if intrinsic variability is present, then $(N-1)/2$ may not be a good uncorrelated estimator for the noise level. Cumming et al. (1999) suggested to estimate the noise level a posteriori with the residuals of the best fit and normalised the periodogram as:

  \[P(f) = \frac{N - 3}{2} \frac{p(f)}{1 - p(f_{\text{best}})}\]

  This is the :Cumming normalization option

Access Frequency Grid and Power Spectrum of the Periodogram

lombscargle returns a LombScargle.Periodogram object, but you most probably want to use the frequency grid and the power spectrum. You can access these vectors with freq and power functions, just like in DSP.jl package. If you want to get the 2-tuple (freq(p), power(p)) use the freqpower function.

power(p::Periodogram)

freq(p::Periodogram)

freqpower(p::Periodogram)

Access Period Grid

The following utilities are the analogs of freq and freqpower, but relative to the periods instead of the frequencies. Thus period(p) returns the vector of periods in the periodogram, that is 1./freq(p), and periodpower(p) gives you the 2-tuple (period(p), power(p)).

period(p::Periodogram)

periodpower(p::Periodogram)

findmaxpower, findmaxfreq, and findmaxperiod Functions

Once you compute the periodogram, you usually want to know which are the frequencies or periods with highest power. To do this, you can use the findmaxfreq and findmaxperiod functions.

findmaxpower(p::Periodogram)

findmaxfreq(p::Periodogram, [interval::AbstractVector{Real}],
            threshold::Real=findmaxpower(p))

findmaxperiod(p::Periodogram, [interval::AbstractVector{Real}],
              threshold::Real=findmaxpower(p))

False-Alarm Probability

Noise in the data produce fluctuations in the periodogram that will present several local peaks, but not all of them related to real periodicities. The significance of the peaks can be tested by calculating the probability that its power can arise purely from noise. The higher the value of the power, the lower will be this probability.

The probability $\Pr(p > p_0)$ that the periodogram power $p$ can exceed the value $p_0$ can be calculated with the prob function, whose first argument is the periodogram and the second one is the $p_0$ value. The function probinv is its inverse: it takes the probability as second argument and returns the corresponding $p_0$ value.

prob(P::Periodogram, pow::Real)

probinv(P::Periodogram, prob::Real)

LombScargle.M(P::Periodogram)

fap(P::Periodogram, pow::Real)

fapinv(P::Periodogram, prob::Real)

Here are the probability functions for each normalization supported by LombScargle.jl:

• :standard ($p \in [0, 1]$):

  \[\Pr(p > p_0) = (1 - p_0)^{(N - 3)/2}\]

• :Scargle ($p \in [0, \infty)$):

  \[\Pr(p > p_0) = \exp(-p_0)\]

• :HorneBaliunas ($p \in [0, (N - 1)/2]$):

  \[\Pr(p > p_0) = \left(1 - \frac{2p_0}{N - 1}\right)^{(N - 3)/2}\]

• :Cumming ($p \in [0, \infty)$):

  \[\Pr(p > p_0) = \left(1 + \frac{2p_0}{N - 3}\right)^{-(N - 3)/2}\]

As explained by Sturrock and Scargle (2010), «the term "false-alarm probability denotes the probability that at least one out of $M$ independent power values in a prescribed search band of a power spectrum computed from a white-noise time series is expected to be as large as or larger than a given value». LombScargle.jl provides the fap function to calculate the false-alarm probability (FAP) of a given power in a periodogram. Its first argument is the periodogram, the second one is the value $p_0$ of the power of which you want to calculate the FAP. The function fap uses the formula

\[\mathrm{FAP} = 1 - (1 - \Pr(p > p_0))^M\]

where $M$ is the number of independent frequencies estimated with $M = T \cdot \Delta f$, being $T$ the duration of the observations and $\Delta f$ the width of the frequency range in which the periodogram has been calculated (see Cumming (2004)). The function fapinv is the inverse of fap: it takes as second argument the value of the FAP and returns the corresponding value $p_0$ of the power.

The detection threshold $p_0$ is the periodogram power corresponding to some (small) value of $\mathrm{FAP}$, i.e. the value of $p$ exceeded due to noise alone in only a small fraction $\mathrm{FAP}$ of trials. An observed power larger than $p_0$ indicates that a signal is likely present (see Cumming (2004)).

Bootstrapping

One of the possible and simplest statistical methods that you can use to measure the false-alarm probability and its inverse is bootstrapping (see section 4.2.2 of Murdoch et al. (1993)).

The recipe of the bootstrap method is very simple to implement:

• repeat the Lomb–Scargle analysis a large number $N$ of times on the original data, but with the signal (and errors, if present) vector randomly shuffled. As an alternative, shuffle only the time vector;
• out of all these simulations, store the powers of the highest peaks;
• in order to estimate the false-alarm probability of a given power, count how many times the highest peak of the simulations exceeds that power, as a fraction of $N$. If you instead want to find the inverse of the false-alarm probability $\text{prob}$, looks for the $N\cdot\text{prob}$-th element of the highest peaks vector sorted in descending order.

Remember to pass to lombscargle function the same options, if any, you used to compute the Lomb–Scargle periodogram before.

LombScargle.jl provides simple methods to perform such analysis. The LombScargle.bootstrap function allows you to create a bootstrap sample with N permutations of the original data.

LombScargle.bootstrap(N::Integer,
                      times::AbstractVector{Real},
                      signal::AbstractVector{Real},
                      errors::AbstractVector{Real}=ones(signal); ...)

LombScargle.bootstrap(N::Integer, plan::PeriodogramPlan)

The false-alarm probability and its inverse can be calculated with fap and fapinv functions respectively. Their syntax is the same as the methods introduced above, but with a LombScargle.Bootstrap object as first argument, instead of the LombScargle.Periodogram one.

fap(b::Bootstrap, power::Real)

fapinv(b::Bootstrap, prob::Real)

LombScargle.model Function

For each frequency $f$ (and hence for the corresponding angular frequency $\omega = 2\pi f$) the Lomb–Scargle algorithm looks for the sinusoidal function of the type

\[a_f\cos(\omega t) + b_f\sin(\omega t) + c_f\]

that best fits the data. In the original Lomb–Scargle algorithm the offset $c$ is null (see Lomb (1976)). In order to find the best-fitting coefficients $a_f$, $b_f$, and $c_f$ for the given frequency $f$, without actually performing the periodogram, you can solve the linear system $\mathbf{A} \mathbf{x} = \mathbf{y}$, where $\mathbf{A}$ is the matrix

\[\begin{bmatrix} \cos(\omega t) & \sin(\omega t) & 1 \end{bmatrix} = \begin{bmatrix} \cos(\omega t_1) & \sin(\omega t_1) & 1 \\ \vdots & \vdots & \vdots \\ \cos(\omega t_n) & \sin(\omega t_n) & 1 \end{bmatrix}\]

$t = [t_1, \dots, t_n]^\mathsf{T}$ is the column vector of observation times, $\mathbf{x}$ is the column vector with the unknown coefficients

\[\begin{bmatrix} a_f \\ b_f \\ c_f \end{bmatrix}\]

and $\mathbf{y}$ is the column vector of the signal. The solution of the matrix gives the wanted coefficients.

This is what the LombScargle.model function does in order to return the best fitting Lomb–Scargle model for the given signal at the given frequency.

LombScargle.model(times::AbstractVector{Real},
                  signal::AbstractVector{R2},
                  [errors::AbstractVector{R3},]
                  frequency::Real,
                  [times_fit::AbstractVector{R4}];
                  center_data::Bool=true,
                  fit_mean::Bool=true)