Documentation |
Build Status |
Code Coverage |
---|---|---|

`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.

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.

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.

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.

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)...)
```

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))...)
```

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.

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

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

The `LombScargle.jl`

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

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.

07/11/2016

3 months ago

258 commits