This tutorial is intended to demonstrate bispectrum Analysis on Lightcurve data.
Bispectrum is an example of a Higher Order Spectra(HOS) and contains more information that simple Powerspectrum or non-ploy spectra.
For detailed information on Bispectra visit : https://arxiv.org/pdf/1308.3150.pdf
In Stingray, Bispectrum can be created from a Lightcurve(For more information on Lightcurve, visit Lightcurve Notebook).
First we import relevant classes.
In [2]:
from stingray import lightcurve
import numpy as np
from stingray.bispectrum import Bispectrum
import matplotlib.pyplot as plt
%matplotlib inline
Lightcurve Object can be created from an array of time stamps and an array of counts. Creating a simple lightcurve to demonstrate Bispectrum.
In [3]:
times = np.arange(1,11)
counts = np.array([2, 1, 3, 4, 2, 5, 1, 0, 2, 3])
lc = lightcurve.Lightcurve(times,counts)
lc.counts
Out[3]:
In [4]:
lc.plot(labels=['times','counts'])
A Bispectrum Object takes 4 parameter.
lc : The light curve (lc).maxlag : Maximum lag on both positive and negative sides of 3rd order cumulant (Similar to lags in correlation).window : Specifies the type of window to apply as as stringscale : 'biased' or 'unbiased' for normalizationArguments 2 and 3 are optional. If maxlag is not specified, it is set to no. of observations in lightcurve divided by 2. i.e lc.n/2 .
In [5]:
bs = Bispectrum(lc)
Different attribute values can be observed by calling relevant properties. Most common are:
In [6]:
bs.freq
Out[6]:
In [7]:
bs.lags
Out[7]:
In [8]:
bs.cum3
Out[8]:
In [9]:
bs.bispec_mag
Out[9]:
In [10]:
bs.bispec_phase
Out[10]:
Bispectrum in stingray also provides functionality for contour plots of:
In [11]:
p = bs.plot_cum3()
p.show()
In [12]:
p = bs.plot_mag()
p.show()
In [13]:
p = bs.plot_phase()
p.show()
In [14]:
dt = 0.0001 # seconds
freq = 1 #Hz
exposure = 50. # seconds
times = np.arange(0, exposure, dt) # seconds
signal = 300 * np.sin(2.*np.pi*freq*times/0.5) + 1000 # counts/s
noisy = np.random.poisson(signal*dt) # counts
lc = lightcurve.Lightcurve(times,noisy)
In [15]:
lc.n
Out[15]:
In [16]:
lc.plot()
In this example, 'unbiased' scaled Bispectrum is calculated.
In [17]:
bs = Bispectrum(lc, maxlag=25, scale='unbiased')
In [18]:
bs.freq[:5]
Out[18]:
In [19]:
bs.lags[-5:]
Out[19]:
In [20]:
bs.n
Out[20]:
In [21]:
bs.cum3[0]
Out[21]:
In [22]:
bs.bispec_mag[1]
Out[22]:
In [23]:
bs.bispec_phase[1]
Out[23]:
In [24]:
p = bs.plot_cum3()
p.show()
In [25]:
p = bs.plot_mag()
p.show()
In [26]:
p = bs.plot_phase()
p.show()
Bispectrum in Stingray now supports 2D windows to apply before calculating Bispectrum.
Windows currently available in Stingray include:
Windows are available in stingray.utils package and can be used by calling create_window function.
Now, we demonstrate Bispectrum with windows applied. By default, now window is applied.
In [29]:
window = 'uniform'
bs = Bispectrum(lc,maxlag=25,window = window, scale ='unbiased')
In [30]:
bs.window_name
Out[30]:
In [32]:
cont = plt.contourf(bs.lags, bs.lags, bs.window, 100, cmap=plt.cm.Spectral_r)
plt.colorbar(cont)
plt.title('2D Uniform window')
Out[32]:
In [34]:
mag_plot = bs.plot_mag()
mag_plot.show()
In [35]:
phase_plot = bs.plot_phase()
phase_plot.show()
Now, let us try some more window functions.
In [36]:
bs = Bispectrum(lc, maxlag=25,window = 'hamming',scale='biased')
In [37]:
bs.window_name
Out[37]:
In [38]:
cont = plt.contourf(bs.lags, bs.lags, bs.window, 100, cmap=plt.cm.Spectral_r)
plt.colorbar(cont)
plt.title('2D Hamming window')
Out[38]:
In [39]:
mag_plot = bs.plot_mag()
mag_plot.show()
In [40]:
phase_plot = bs.plot_phase()
phase_plot.show()
In [45]:
bs = Bispectrum(lc, maxlag = 25, window='triangular',scale='unbiased')
In [46]:
bs.window_name
Out[46]:
In [47]:
cont = plt.contourf(bs.lags, bs.lags, bs.window, 100, cmap=plt.cm.Spectral_r)
plt.colorbar(cont)
plt.title('2D Flat Top window')
Out[47]:
In [48]:
bs.plot_mag().show()
In [52]:
bs.plot_phase().show()