Examples

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

julia> using LombScargle
julia> using Random
julia> Random.seed!(42); # make plots reproducible
julia> ntimes = 1001;
julia> t = range(0.01, 10pi, length = ntimes) # Observation times0.01:0.03140592653589793:31.41592653589793
julia> t += step(t)*rand(ntimes); # Randomize times
julia> s = sinpi.(t) .+ 1.5cospi.(2t) .+ rand(ntimes); # The signal
julia> plan = LombScargle.plan(t, s); # Pre-plan the periodogram
julia> pgram = lombscargle(plan) # Compute the periodogramLombScargle.Periodogram{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}}([0.0009002061856020596, 0.0008836113626863781, 0.0008495252318279351, 0.0007945083538758122, 0.0007052919678129175, 0.0005425637744200134, 0.0002821925162241195, 5.119517579364273e-5, 3.589360431210285e-5, 0.00019032755028405127  …  0.0011197360659888319, 0.0013600579636300514, 0.0009567336181250871, 0.00022198735604679484, 7.888602659701585e-5, 0.0010836848233666905, 0.0026394015757926293, 0.0034258931687385576, 0.0027428417988911587, 0.0011372880435844058], 0.0031823671601760094:0.006364734320352019:79.63237344908428, [0.015451270375901948, 0.051508007716585986, 0.08093296934439527, 0.10944493943658266, 0.15217510074016963, 0.18219939042833327, 0.21070470491330284, 0.25505298746255434, 0.2920490151715063, 0.2956190459595937  …  31.143038468604786, 31.164857936510924, 31.20507145945872, 31.233375392399598, 31.271122085204944, 31.314799568543375, 31.344062482974675, 31.380022287634237, 31.38983942755385, 31.438601072634807], :standard)

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

You can also plot the power vs the period, instead of the frequency, with periodpower:

using Plots
plot(periodpower(pgram)...)
xlabel!("Period")
ylabel!("Lomb–Scargle power")

Warning

If you do not fit for the mean of the signal (fit_mean=false keyword to lombscargle function) without centering the data (center_data=false) you can get inaccurate results. For example, spurious peaks at low frequencies can appear and the real peaks lose power:

plot(freqpower(lombscargle(t, s, fit_mean=false, center_data=false))...)

Tip

You can tune the frequency grid with appropriate keywords to lombscargle function. For example, in order to increase the sampling increase samples_per_peak, and set maximum_frequency to lower values in order to narrow the frequency range:

plot(freqpower(lombscargle(t, s, samples_per_peak=20, maximum_frequency=1.5))...)

If you simply want to use your own frequency grid, directly set the frequencies keyword:

plot(freqpower(lombscargle(t, s, frequencies=0.001:1e-3:1.5))...)

Signal with Uncertainties

The generalised Lomb–Scargle periodogram 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)
# Run one of the two following equivalent commands
plot(freqpower(lombscargle(t, s, errors, maximum_frequency=1.5))...)
plot(freqpower(lombscargle(t, measurement.(s, errors), maximum_frequency=1.5))...)

This is the plot of the power versus the period:

# Run one of the two following equivalent commands
plot(periodpower(lombscargle(t, s, errors, maximum_frequency=1.5))...)
plot(periodpower(lombscargle(t, measurement(s, errors), maximum_frequency=1.5))...)

We recall that the generalised Lomb–Scargle algorithm is used when the fit_mean optional keyword to lombscargle is true if no error is provided, instead it is always used if the signal has uncertainties.

Find Highest Power and Associated Frequencies and Periods

findmaxfreq function tells you the frequencies with the highest power in the periodogram (and you can get the period by taking its inverse):

julia> t = range(0, stop = 10, length = 1001);

julia> s = sinpi.(t);

julia> plan = LombScargle.plan(t, s); # Plan the periodogram

julia> p = lombscargle(plan);

julia> findmaxperiod(p) # Period with highest power
1-element Vector{Float64}:
 0.004987779939149084

julia> findmaxfreq(p) # Frequency with the highest power
1-element Vector{Float64}:
 200.49

This peak is at high frequencies, very far from the expected value of the period of 2. In order to find the real peak, you can either narrow the ranges in order to exclude higher armonics

julia> findmaxperiod(p, [1, 10]) # Limit the search to periods in [1, 10]
1-element Vector{Float64}:
 2.0408163265306123

julia> findmaxfreq(p, [0.1, 1]) # Limit the search to frequencies in [0.1, 1]
1-element Vector{Float64}:
 0.49

or pass the threshold argument to findmaxfreq or findmaxperiod. You can use findmaxpower to discover the highest power in the periodogram:

julia> findmaxpower(p)
0.9996235276303144

julia> findmaxperiod(p, 0.95)
10-element Vector{Float64}:
 2.0408163265306123
 1.9607843137254901
 0.010051261433309882
 0.010049241282283187
 0.009951238929246691
 0.009949258780220873
 0.005012782595618828
 0.005012280086211218
 0.004987779939149084
 0.004987282429804

julia> findmaxfreq(p, 0.95)
10-element Vector{Float64}:
   0.49
   0.51
  99.49
  99.51
 100.49
 100.51
 199.49
 199.51
 200.49
 200.51

The first peak is the real one, the other double peaks appear at higher armonics.

Significance of the Peaks

The significance of the peaks in the Lomb–Scargle periodogram can be assessed by measuring the False-Alarm Probability. Analytic expressions of this quantity and its inverse can be obtained with the fap and fapinv functions, respectively.

julia> t = linspace(0.01, 20, samples_per_peak = 10);

julia> s = sinpi.(e.*t).^2 .- cos.(5t).^4;

julia> plan = LombScargle.plan(t, s);

julia> p = lombscargle(plan);

# Find the false-alarm probability for the highest peak.
julia> fap(p, 0.3)
0.028198095962262748

Thus, a peak with power $0.3$ has a probability of $0.028$ that it is due to noise only. A quantity that is often used is the inverse of the false-alarm probability as well: what is the minimum power whose false-alarm probability is lower than the given probability? For example, if you want to know the minimum power for which the false-alarm probability is at most $0.01$ you can use:

julia> fapinv(p, 0.01)
0.3304696923786712

As we already noted, analytic expressions of the false-alarm probability and its inverse may not be reliable if your data does not satisfy specific assumptions. A better way to calculate this quantity is to use statistical methods. One of this is bootstrapping. In LombScargle.jl, you can use the function LombScargle.bootstrap to create a bootstrap sample and then you can calculate the false-alarm probability and its inverse using this sample.

Create a bootstrap sample with 10000 resamplings of the original data, re-using the same periodogram plan. The larger the better. This may take some minutes.

julia> b = LombScargle.bootstrap(10000, plan)

Calculate the false-alarm probability of a peak with power 0.3 using this bootstrap sample.

julia> fap(b, 0.3)
0.0209

Calculate the lowest power that has probability less than 0.01 in this bootstrap sample.

julia> fapinv(b, 0.01)
0.3268290388848437

If you query fapinv with a too low probability, the corresponding power cannot be determined and you will get NaN as result.

julia> fapinv(b, 1e-5)
NaN

If you want to find the power corresponding to a false-alarm probability of $\text{prob} = 10^{-5}$, you have to create a new bootstrap sample with $N$ resamplings so that $N\cdot\text{prob}$ can be rounded to an integer larger than or equal to one (for example $N = 10^{5}$).

Find the Best-Fitting Model

The LombScargle.model function can help you to test whether a certain frequency fits well your data.

using LombScargle
using Random
using Plots
Random.seed!(42)
t = range(0.01, stop = 10pi, length = 1000) # Observation times
s = sinpi.(t) .+ 1.2cospi.(t) .+ 0.3rand(length(t)) # The noisy signal
# Pick-up the best frequency
f = findmaxfreq(lombscargle(t, s, maximum_frequency=10, samples_per_peak=20))[1]
t_fit = range(0, stop = 1, length = 50)
s_fit = LombScargle.model(t, s, f, t_fit/f) # Determine the model
scatter(mod.(t.*f, 1), s, lab="Phased data", title="Best Lomb-Scargle frequency: $f")
plot!(t_fit, s_fit, lab="Best-fitting model", linewidth=4)

Tip

If there are more than one dominant frequency you may need to consider more models. This task may require some work and patience. Plot the periodogram in order to find the best frequencies.

using LombScargle
using Random
using Plots
Random.seed!(42)
t = range(0.01, stop = 5, length = 1000) # Observation times
s = sinpi.(2t) .+ 1.2cospi.(4t) .+ 0.3rand(length(t)) # Noisy signal
plan = LombScargle.plan(t, s, samples_per_peak=50)
p = lombscargle(plan)
# After plotting the periodogram, you discover
# that it has two prominent peaks around 1 and 2.
f1 = findmaxfreq(p, [0.8, 1.2])[1] # Get peak frequency around 1
f2 = findmaxfreq(p, [1.8, 2.2])[1] # Get peak frequency around 2
fit1 = LombScargle.model(t, s, f1) # Determine the first model
fit2 = LombScargle.model(t, s, f2) # Determine the second model
scatter(t, s, lab="Data", title="Best-fitting Lomb-Scargle model")
plot!(t, fit1 + fit2, lab="Best-fitting model", linewidth=4)