Lomb Scargle Power Spectra¶
This tutorial shows how to make and manipulate a Lomb Scargle power spectrum of two light curves using Stingray.
[1]:
from stingray.lightcurve import Lightcurve
from stingray.lombscargle import LombScarglePowerspectrum
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline
import matplotlib.font_manager as font_manager
%matplotlib inline
plt.style.use('seaborn-v0_8-talk')
font_prop = font_manager.FontProperties(size=16)
1. Create a light curve¶
There are two ways to make Lightcurve
objects. We’ll show one way here. Check out Lightcurve for more examples.
Make one with signals in units of counts. It is a sine wave with random normal noise, frequency of 3 and at random times and make its counts non-negative by subtracting its least value.
[2]:
rand = np.random.default_rng(42)
n = 100
t = np.sort(rand.random(n)) * 10
y = np.sin(2 * np.pi * 3.0 * t) + 0.1 * rand.standard_normal(n)
sub = np.min(y)
y -= sub
t0 = np.linspace(0, 10, 1000)
y0 = np.sin(2 * np.pi * 3.0 * t0) + 0.1 * rand.standard_normal(t0.size)
sub = np.min(y0)
y0 -= sub
spline = make_interp_spline(t, y)
Lets convert them into Lightcurve
objects
[3]:
lc = Lightcurve(t, y)
Let us plot them to see how they look
[4]:
fig, ax = plt.subplots(1,1,figsize=(10,6))
ax.scatter(lc.time, lc.counts, lw=2, color='blue',label='lc')
ax.plot(t0, y0, lw=2, color='red',label='source of lc')
ax.set_xlabel("Time (s)", fontproperties=font_prop)
ax.set_ylabel("Counts (cts)", fontproperties=font_prop)
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
ax.tick_params(which='major', width=1.5, length=7)
ax.tick_params(which='minor', width=1.5, length=4)
plt.legend()
plt.show()
2. Pass the light curve to the LombScarglePowerspectrum
class to create a LombScarglePowerspectrum
object.¶
You can also specify the optional attribute norm
if you wish to normalize the real part of the power spectrum to squared fractional rms, Leahy, or squared absolute normalization. The default normalization is ‘none’.
[5]:
lps = LombScarglePowerspectrum(
lc,
min_freq=0,
max_freq=None,
method="fast",
power_type="all",
norm="none",
)
We can print the first five values in the arrays of the positive Fourier frequencies and the power. The power has only real component, and imaginary component is zero.
[6]:
print(lps.freq[0:5])
print(lps.power[0:5])
[0.05163902 0.15491705 0.25819509 0.36147313 0.46475116]
[ 15.49526224+0.j 120.05686691+0.j 96.589673 +0.j 127.2231466 +0.j
30.42053746+0.j]
Parameters¶
data
: This parameter allows you to provide the light curve data to be Fourier-transformed. It can be either a`stingray.lightcurve.Lightcurve
<https://docs.stingray.science/en/stable/core.html#working-with-lightcurves>`__ or`stingray.events.EventList
<https://docs.stingray.science/en/stable/core.html#working-with-event-data>`__ object. It is optional, and the default value isNone
.norm
: Thenorm
parameter defines the normalization of the power spectrum. It accepts string values from the set {frac
,abs
,leahy
,none
}. The default normalization is set tonone
.power_type
: Thepower_type
parameter allows you to specify the type of power spectral power you want to compute. The options are:real
for the real part,absolute
for the magnitude, andall
to compute both real part and magnitude. The default isall
.fullspec
: This is a boolean parameter that determines whether to keep only the positive frequencies or include both positive and negative frequencies in the power spectrum. When set toFalse
(default), only positive frequencies are kept; when set toTrue
, both positive and negative frequencies are included.
Other Parameters¶
dt
: When constructing light curves using`stingray.events.EventList
<https://docs.stingray.science/en/stable/core.html#working-with-event-data>`__ objects, thedt
parameter represents the time resolution of the light curve. It is a float value that needs to be provided.skip_checks
: This is a boolean parameter that, when set toTrue
, skips initial checks for speed or other reasons. It’s useful when you have confidence in the inputs and want to improve processing speed.min_freq
: This parameter specifies the minimum frequency at which the Lomb-Scargle Fourier Transform should be computed.max_freq
: Similarly, themax_freq
parameter sets the maximum frequency for the Lomb-Scargle Fourier Transform.df
: Thedf
parameter, a float, represents the frequency resolution. It’s relevant when constructing light curves using`stingray.events.EventList
<https://docs.stingray.science/en/stable/core.html#working-with-event-data>`__ objects.method
: Themethod
parameter determines the method used by the Lomb-Scargle Fourier Transformation function. The allowed values arefast
andslow
, with the default beingfast
. Thefast
method uses the optimized Press and Rybicki O(n*log(n)) algorithm.oversampling
: This optional float parameter represents the interpolation oversampling factor. It is applicable when using the fast algorithm for the Lomb-Scargle Fourier Transform. The default value is 5.
Attributes¶
freq
: Thefreq
attribute is a numpy array that contains the mid-bin frequencies at which the Fourier transform samples the power spectrum.power
: Thepower
attribute is a numpy array that contains the normalized squared absolute values of Fourier amplitudes.power_err
: Thepower_err
attribute is a numpy array that provides the uncertainties associated with thepower
. The uncertainties are approximated using the formulapower_err = power / sqrt(m)
, wherem
is the number of power values averaged in each bin. For a single realization (m=1
), the error is equal to the power.df
: Thedf
attribute is a float that indicates the frequency resolution.m
: Them
attribute is an integer representing the number of averaged powers in each bin.n
: Then
attribute is an integer indicating the number of data points in the light curve.nphots
: Thenphots
attribute is a float representing the total number of photons in the light curve.
We can plot the power spectrum by using the plot function or manually taking the freq
and power
attributes
[7]:
fig, ax = plt.subplots(1,3,figsize=(15,6),sharey=True)
lps.plot(ax=ax[0])
ax[0].set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax[0].set_ylabel("Power", fontproperties=font_prop)
ax[1].plot(lps.freq, lps.power.real, lw=2, color='red')
ax[1].set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax[1].set_ylabel("Power(Real Component)", fontproperties=font_prop)
ax[2].plot(lps.freq, lps.power.imag, lw=2, color='blue')
ax[2].set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax[2].set_ylabel("Power(Imaginary Component)", fontproperties=font_prop)
[7]:
Text(0, 0.5, 'Power(Imaginary Component)')