Examples¶
In the following, the basic usage of the periodogram classes is demonstrated.
Fourier spectrum of evenly-sampled Poisson data¶
from __future__ import print_function, division
import numpy
import matplotlib.pylab as plt
# Import pyTiming
from PyAstronomy.pyTiming import pyPeriod
# Create some evenly sampled artificial data (Poisson noise)
time = numpy.arange(1000.)/10.
flux = numpy.random.poisson(1, len(time))
# Build the TimeSeries instance
lc = pyPeriod.TimeSeries(time, flux)
# Compute the Leahy-normalized Fourier transform,
# plot the time series, and check that the mean
# power level is 2 as expected.
fft = pyPeriod.Fourier(lc)
fig, ax = fft.plot()
print('Mean power level:', numpy.mean(fft.power))
plt.show()
Error-weighted (generalized) Lomb periodogram¶
from __future__ import print_function, division
import numpy
import matplotlib.pylab as plt
from PyAstronomy.pyTiming import pyPeriod
# Create unevenly saplted data with frequency=0.1,
# measurement error and Gaussian noise
time = numpy.arange(1000.) + numpy.random.normal(0., 0.1, 1000)
flux = 0.15 * numpy.sin(2. * numpy.pi * time / 10.)
# Add some noise
flux += numpy.random.normal(0, 1, time.size) * 0.5
error = numpy.ones(time.size) * 0.5
# Plot the light curve in top panel
plt.subplot(3, 1, 1)
plt.errorbar(time, flux, yerr=error)
# Build the TimeSeries instance
lc = pyPeriod.TimeSeries(time, flux, error)
# Compute and plot fast Lomb-Scargle periodogram,
# which does not take errors into account.
ls = pyPeriod.LombScargle(lc, ofac=1, hifac=1)
# Plot the Lomb-Scargle periodogram in middle panel
plt.subplot(3, 1, 2)
plt.plot(ls.freq, ls.power, 'r-')
# Compute the full error-weighted Lomb-Periodogram
# in 'Cumming' normalization and calculate the
# significance of the maximum peak.
clp = pyPeriod.Gls(lc, ofac=10, hifac=1, norm="Cumming")
maxPower = numpy.max(clp.power)
print("GLS maximum power: ", maxPower)
print("GLS statistics of maximum power peak: ", clp.stats(maxPower))
# Plot the generalized Lomb-Scargle periodogram in
# bottom panel.
plt.subplot(3, 1, 3)
plt.plot(clp.freq, clp.power)
# Show the results
plt.show()