LombScargle.jl

Documentation	Build Status	Code Coverage

Introduction

LombScargle.jl is a Julia package for a fast multi-threaded estimation of the frequency spectrum of a periodic signal with the Lomb–Scargle periodogram.

Another Julia package that provides tools to perform spectral analysis of signals is DSP.jl, but its methods require that the signal has been sampled at equally spaced times. Instead, the Lomb–Scargle periodogram enables you to analyze unevenly sampled data as well, which is a fairly common case in astronomy, a field where this periodogram is widely used.

The algorithms used in this package are reported in the following papers:

• Press, W. H., Rybicki, G. B. 1989, ApJ, 338, 277 (URL: http://dx.doi.org/10.1086/167197, Bibcode: http://adsabs.harvard.edu/abs/1989ApJ...338..277P)
• Townsend, R. H. D. 2010, ApJS, 191, 247 (URL: http://dx.doi.org/10.1088/0067-0049/191/2/247, Bibcode: http://adsabs.harvard.edu/abs/2010ApJS..191..247T)
• Zechmeister, M., Kürster, M. 2009, A&A, 496, 577 (URL: http://dx.doi.org/10.1051/0004-6361:200811296, Bibcode: http://adsabs.harvard.edu/abs/2009A%26A...496..577Z)

The package provides facilities to:

• compute the periodogram using different methods (with different speeds) and different normalizations. This is one of the fastest implementations of these methods available as free software. If Julia is run with more than one thread, computation is automatically multi-threaded, further speeding up calculations;
• access the frequency and period grid of the resulting periodogram, together with the power spectrum;
• find the maximum power in the periodogram and the frequency and period corresponding to the peak. All these queries can be restricted to a specified region, in order to search a local maximum, instead of the global one;
• calculate the probability that a peak arises from noise only (false-alarm probability) using analytic formulas, in order to assess the significance of the peak;
• perform bootstrap resamplings in order to compute the false-alarm probability with a statistical method;
• determine the best-fitting Lomb–Scargle model for the given data set at the given frequency.

All these features are thoroughly described in the full documentation, see below. Here we only give basic information.

Documentation

The complete manual of LombScargle.jl is available here. It has detailed explanation of all functions provided by the package and more examples than what you will find here, also with some plots.

Installation

The latest version of LombScargle.jl is available for Julia 1.0 and later versions, and can be installed with Julia built-in package manager. In a Julia session, after entering the package manager mode with ], run the command

pkg> add LombScargle

Older versions are also available for Julia 0.4-0.7.

Usage

After installing the package, you can start using it with

julia> 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, signal[, errors])

which returns a LombScargle.Periodogram. The only mandatory arguments are:

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

All these vectors must have the same length. The only optional argument is:

• errors: the uncertainties associated to each signal point. This vector must have the same length as times and signal.

Besides the two arguments introduced above, lombscargle has a number of other optional keywords in order to choose the right algorithm to use and tweak the periodogram. For the description of all these arguments see the complete manual.

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/latest/ for more details.

With the LombScargle.plan function you can pre-plan a periodogram and save time and memory for the actual computation of the periodogram. See the manual for details.

Examples

Here is an example of a noisy periodic signal (sin(π*t) + 1.5*cos(2π*t)) sampled at unevenly spaced times.

julia> using LombScargle

julia> ntimes = 1001
1001

# Observation times
julia> t = range(0.01, stop=10pi, length=ntimes)
0.01:0.03140592653589793:31.41592653589793

# Randomize times
julia> t += step(t)*rand(ntimes);

# The signal
julia> s = sinpi.(t) .+ 1.5 .* cospi.(2t) .+ rand(ntimes);

# Pre-plan the periodogram (see the documentation)
julia> plan = LombScargle.plan(t, s);

# Compute the periodogram
julia> pgram = lombscargle(plan)

You can plot the result, for example with Plots package. Use freqpower function to get the frequency grid and the power of the periodogram as a 2-tuple.

using Plots
plot(freqpower(pgram)...)

Signal with Uncertainties

The generalised Lomb–Scargle periodogram (used when the fit_mean optional keyword is true) is able to handle a signal with uncertainties, and they will be used as weights in the algorithm. The uncertainties can be passed either as the third optional argument errors to lombscargle or by providing this function with a signal vector of type Measurement (from Measurements.jl package).

using Measurements, Plots
ntimes = 1001
t = range(0.01, stop=10pi, length=ntimes)
s = sinpi.(2t)
errors = rand(0.1:1e-3:4.0, ntimes)
plot(freqpower(lombscargle(t, s, errors, maximum_frequency=1.5))...)
plot(freqpower(lombscargle(t, measurement(s, errors), maximum_frequency=1.5))...)

Performance

A pre-planned periodogram in LombScargle.jl computed in single thread mode with the fast method is more than 2 times faster than the implementation of the same algorithm provided by Astropy, and more than 4 times faster if 4 FFTW threads are used (on machines with at least 4 physical CPUs).

The following plot shows a comparison between the times needed to compute a periodogram for a signal with N datapoints using LombScargle.jl, with 1 or 4 FFTW threads (with flags = FFTW.MEASURE for better performance), and the single-threaded Astropy implementation. (Julia version: 1.6.0; LombScargle.jl version: 1.0.0; Python version: 3.8.6; Astropy version: 4.1. CPU: Intel(R) Core(TM) i7-4870HQ CPU @ 2.50GHz.)

Note that this comparison is unfair, as Astropy doesn’t support pre-planning a periodogram nor multi-threading, and it pads vectors for FFT to a length which is a power of 2, while by default LombScargle.jl uses length which are multiples of 2, 3, 5, 7. A non-planned periodogram in single thread mode in LombScargle.jl is still twice as fast as Astropy.

Development

The package is developed at https://github.com/JuliaAstro/LombScargle.jl. There you can submit bug reports, make suggestions, and propose pull requests.

History

The ChangeLog of the package is available in NEWS.md file in top directory.

License

The LombScargle.jl package is licensed under the BSD 3-clause "New" or "Revised" License. The original author is Mosè Giordano.

Acknowledgemets

This package adapts the implementation in Astropy of the the fast Lomb–Scargle method by Press & Rybicki (1989). We claim no endorsement nor promotion by the Astropy Team.