The Generalized Lomb-Scargle Periodogram (GLS)¶

The GLS class provides an implementation of the Generalized Lomb-Scargle periodogram as described by Zechmeister & Kuerster 2009 (A&A 496, 577) and earlier by Ferraz-Mello 1981 (AJ 86, 691). Compared to the periodogram presented, e.g., by Lomb 1976, the GLS takes into account measurement errors and a constant term in the fit of the trigonometric function.

The periodogram calculation can be controlled by a number of keyword argument. Amongst others, the boundaries of the frequency range and its sampling can be adjusted or specified explicitly, various conventions for the periodogram normalization can be adopted, and the classical Lomb-Scargle periodogram can be obtained.

In the following, we provide a number of examples, presenting some of the functionality.

Example: Simple periodogram with default options¶

This example shows a basic application of the GLS class, using only default options and an artificial input data set. No errors are specified.

from __future__ import print_function, division
import numpy as np
import matplotlib.pylab as plt
from PyAstronomy.pyTiming import pyPeriod

# Create some evenly sampled data including a periodic term.
# Number of data points, input frequency, amplitude of input signal,
# and standard deviation of noise.
N = 1000
f = 0.1
A = 0.15
sig = 0.2

time = np.arange(float(N))
flux = A * np.sin(2. * np.pi*time*f)
# Adding the noise
flux += np.random.normal(0, sig, time.size)

# Compute the GLS periodogram with default options.
clp = pyPeriod.Gls((time, flux))

# Print helpful information to screen
clp.info()

# and plot power vs. frequency.
plt.xlabel("Frequency")
plt.ylabel("Power")
plt.plot(clp.freq, clp.power, 'b.-')
plt.show()

Example: Adding an error column and FAP levels¶

In this example, an error column will be added. Additionally, we will calculate and show a number of False-Alarm-Probability levels to help distinguish between significant and spurious peaks.

from __future__ import print_function, division
import numpy as np
import matplotlib.pylab as plt
from PyAstronomy.pyTiming import pyPeriod

# Create some evenly sampled data including a periodic term.
# Number of data points, input frequency, amplitude of input signal,
# and standard deviation of noise.
N = 1000
f = 0.1
A = 0.05
sig = 0.2

time = np.arange(float(N))
flux = A * np.sin(2. * np.pi*time*f)
# Adding the noise
flux += np.random.normal(0, sig, time.size)
# Adding an error column
err = np.ones(N)*sig

# Compute the GLS periodogram with default options.
# Choose Zechmeister-Kuerster normalization explicitly
clp = pyPeriod.Gls((time, flux, err), norm="ZK")

# Print helpful information to screen
clp.info()

# Define FAP levels of 10%, 5%, and 1%
fapLevels = np.array([0.1, 0.05, 0.01])
# Obtain the associated power thresholds
plevels = clp.powerLevel(fapLevels)

# and plot power vs. frequency.
plt.xlabel("Frequency")
plt.ylabel("Power")
plt.plot(clp.freq, clp.power, 'b.-')
# Add the FAP levels to the plot
for i in range(len(fapLevels)):
    plt.plot([min(clp.freq), max(clp.freq)], [plevels[i]]*2, '--',
             label="FAP = %4.1f%%" % (fapLevels[i]*100))
plt.legend()

plt.show()

Example: Getting highest-power (best-fit) sine curve¶

The best-fit sine curve is extracted from the periodogram class.

from __future__ import print_function, division
import numpy as np
import matplotlib.pylab as plt
from PyAstronomy.pyTiming import pyPeriod

# Create some evenly sampled data including a periodic term.
# Number of data points, input frequency, amplitude of input signal,
# and standard deviation of noise.
N = 1000
f = 0.171
A = 0.5
sig = 0.2

time = np.arange(float(N))
flux = A * np.sin(2. * np.pi*time*f)
# Adding the noise
flux += np.random.normal(0, sig, time.size)
# Adding an error column
err = np.ones(N)*sig

# Compute the GLS periodogram with default options.
# Choose Zechmeister-Kuerster normalization explicitly
clp = pyPeriod.Gls((time, flux, err), norm="ZK")

# Get index associated with highest power
ifmax = np.argmax(clp.power)
# and highest power and associated frequency
pmax = clp.power[ifmax]
fmax = clp.freq[ifmax]
# Convert frequency into period
hpp = 1./fmax
print("Highest-power period: ", hpp)

# Calculate sine wave associated with 'best-fit' frequency
bestSine = clp.sinmod(time)

plt.subplot(2, 1, 1)
plt.title("Data and sine asscoiated with highest-power frequency")
plt.plot(time, flux, 'b.')
plt.plot(time, bestSine, 'r--')
plt.subplot(2, 1, 2)
plt.title("Folded data")
plt.plot(time/hpp-np.floor(time/hpp), flux, 'b.')
plt.show()

API documentation¶

class PyAstronomy.pyTiming.pyPeriod.Gls(lc, fbeg=None, fend=None, Pbeg=None, Pend=None, ofac=10, hifac=1, freq=None, norm='ZK', ls=False, fast=False, verbose=False, **kwargs)¶

Compute the Generalized Lomb-Scargle (GLS) periodogram.

The Gls class computes the error-weighted Lomb-Scargle periodogram as developed by [ZK09] using various possible normalizations.

The constructor of Gls takes a TimeSeries instance (i.e., a light curve) as first argument. The constructor allows to pass keywords to adjust the freq array, which will be used to calculate the periodogram.

The main result of the calculation, i.e., the power, are stored in the class property power.

Parameters

lctuple or list or TimeSeries object: The light curve data either in the form of a TimeSeries object (or any object providing the attributes time, flux, and error) or a tuple or list providing time as first element, flux as second element, and optionally, the error as third element.
fbeg, fendfloat, optional: The beginning and end frequencies for the periodogram (inverse units of time axis).
Pbeg, Pendfloat, optional: The beginning and end periods for the periodogram (same units as for time axis).
ofacint: Oversampling factor of frequency grid (default=10).
hifacfloat: Maximum frequency freq = hifac * (average Nyquist frequency) (default=1).
freqarray, optional: Contains the frequencies at which to calculate the periodogram. If given, fast and verbose option are not available. If not given, a frequency array will be automatically generated.
normstring, optional: The normalization; either of “ZK”, “Scargle”, “HorneBaliunas”, “Cumming”, “wrms”, “chisq”. The default is unity (“ZK”).
lsboolean, optional: If True, the conventional Lomb-Scargle periodogram will be computed (default is False).
fastboolean, optional: If True, recursive relations for trigonometric functions will be used leading to faster evaluation (default is False).
verboseboolean, optional: Set True to obtain some statistical output (default is False).

Examples

Create 1000 unevenly sampled data points with frequency=0.1, measurement error and Gaussian noise

>>> time = np.random.uniform(54000., 56000., 1000)
>>> flux = 0.15 * np.sin(2. * np.pi * time / 10.)

Add some noise

>>> error = 0.5 * np.ones(time.size)
>>> flux += np.random.normal(0, error)

Compute the full error-weighted Lomb-Periodogram in ‘ZK’ normalization and calculate the significance of the maximum peak.

>>> gls = Gls((time, flux, error), verbose=True)

>>> maxPower = gls.pmax
>>> print("GLS maximum power: ", maxPower)
>>> print("GLS statistics of maximum power peak: ", gls.stats(maxPower))
>>> gls.plot(block=True)

Attributes

powerarray: The normalized power of the GLS.
freqarray: The frequency array.
ofacint: The oversampling factor of frequency grid.
hifacfloat: The maximum frequency.
tarray: The abscissa data values.
yarray: The ordinate data values.
e_yarray: The errors of the data values.
normstring, {‘ZK’, ‘Scargle’, ‘HorneBaliunas’, ‘Cumming’, ‘wrms’, ‘chisq’}: The used normalization.

Methods

`FAP`(Pn)	Obtain the false-alarm probability (FAP).
`info`([noprint])	Prints some basic statistical output screen.
`plot`([block, period])	Create a plot.
`pnorm`([norm])	Assign or modify normalization (can be done afterwards).
`powerLevel`(FAPlevel)	Power threshold for FAP level.
`prob`(Pn)	Probability of obtaining the given power.
`probInv`(Prob)	Calculate minimum power for given probability.
`sinmod`(t)	Calculate best-fit sine curve.
`stats`(Pn)	Obtain basic statistics for power threshold.
`toFile`(ofile[, header])	Write periodogram to file.

FAP(Pn)¶

Obtain the false-alarm probability (FAP).

The FAP 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 as large as or larger than the threshold, Pn. It is assessed through

\[FAP(Pn) = 1 - (1-Prob(P>Pn))^M \; ,\]

where “Prob(P>Pn)” depends on the type of periodogram and normalization and is calculated by using the prob method; M is the number of independent power values and is computed internally.

Parameters

Pnfloat: Power threshold.

Returns

FAPfloat: False alarm probability.

info(noprint=False)¶

Prints some basic statistical output screen.

Parameters

noprintboolean: If True, printing is suppressed. Default is False.

Returns

Informationdictionary: A dictionary with the printed information.

plot(block=False, period=False)¶: Create a plot.

pnorm(norm='ZK')¶

Assign or modify normalization (can be done afterwards).

Parameters

normstring, optional: The normalization to be used (default is ‘ZK’).

Examples

>>> gls.pnorm('wrms')

powerLevel(FAPlevel)¶

Power threshold for FAP level.

Parameters

FAPlevelfloat or array_like: “False Alarm Probability” threshold

Returns

Thresholdfloat or array: The power threshold pertaining to a specified false-alarm probability (FAP). Powers exceeding this threshold have FAPs smaller than FAPlevel.

prob(Pn)¶

Probability of obtaining the given power.

Calculate the probability to obtain a power higher than Pn from the noise, which is assumed to be Gaussian.

Note

Normalization (see [ZK09] for further details).

Scargle:

\[exp(-Pn)\]

HorneBaliunas:

\[\left(1 - 2 \times \frac{Pn}{N-1} \right)^{(N-3)/2}\]

Cumming:

\[\left(1+2\times \frac{Pn}{N-3}\right)^{-(N-3)/2}\]

Parameters

Pnfloat: Power threshold.

Returns

Probabilityfloat: The probability to obtain a power equal or higher than the threshold from the noise.

probInv(Prob)¶

Calculate minimum power for given probability.

This function is the inverse of Prob(Pn). Returns the minimum power for a given probability threshold Prob.

Parameters

Probfloat: Probability threshold.

Returns

Power thresholdfloat: The minimum power for the given false-alarm probability threshold.

sinmod(t)¶

Calculate best-fit sine curve.

The parameters of the best-fit sine curve can be accessed via the dictionary attribute hpstat. Specifically, “amp” holds the amplitude, “fbest” the best-fit frequency, “T0” the reference time (i.e., time of zero phase), and “offset” holds the additive offset of the sine wave.

Parameters

tarray: Time array at which to calculate the sine.

Returns

Sine curvearray: The best-fit sine curve (i.e., that for which the power is maximal).

stats(Pn)¶

Obtain basic statistics for power threshold.

Parameters

Pnfloat: Power threshold.

Returns

Statisticsdictionary: A dictionary containing {‘Pn’: Pn, ‘Prob’: Prob(Pn) , ‘FAP’: FAP(Pn)} for the specified power threshold, Pn.

toFile(ofile, header=True)¶

Write periodogram to file.

The output file is a standard text file with two columns, viz., frequency and power.

Parameters

ofilestring: Name of the output file.

The Generalized Lomb-Scargle Periodogram (GLS)¶

Example: Simple periodogram with default options¶

Example: Adding an error column and FAP levels¶

Example: Getting highest-power (best-fit) sine curve¶

API documentation¶

Table of Contents

Previous topic

Next topic

This Page