Dynamical Power Spectra (on real data)¶
[1]:
%matplotlib inline
[2]:
# load auxiliary libraries
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
# import stingray
import stingray
plt.style.use('seaborn-talk')
All starts with a lightcurve..¶
Open the event file with astropy.io.fits
[3]:
f = fits.open('emr_cleaned.fits')
The time resolution is stored in the header of the first extension under the Keyword TIMEDEL
[4]:
dt = f[1].header['TIMEDEL']
The collumn TIME
of the first extension stores the time of each event
[5]:
toa = f[1].data['Time']
Let’s create a Lightcurve from the Events time of arrival witha a given time resolution
[6]:
lc = stingray.Lightcurve.make_lightcurve(toa=toa, dt=dt)
[7]:
lc.plot()

DynamicPowerspectrum¶
Let’s create a dynamic powerspectrum with the a segment size of 16s and the powers with a “leahy” normalization
[8]:
dynspec = stingray.DynamicalPowerspectrum(lc=lc, segment_size=16, norm='leahy')
The dyn_ps attribute stores the power matrix, each column corresponds to the powerspectrum of each segment of the light curve
[9]:
dynspec.dyn_ps
[9]:
array([[ 2.01901704e+00, 2.32485459e+00, 5.14704363e+00, ...,
9.76872866e-01, 9.49269045e-01, 4.60522187e+02],
[ 2.93960257e+00, 2.48892516e+00, 3.39280288e+00, ...,
6.23511732e+00, 4.27550837e+00, 1.06261843e+02],
[ 3.64619904e+00, 1.58266627e+00, 3.42614944e-01, ...,
1.16952148e+00, 3.54994270e+00, 4.56956463e+01],
...,
[ 1.69311108e+00, 5.18784072e-01, 1.57151667e+00, ...,
1.09923562e+00, 3.40274378e-01, 2.53108287e+00],
[ 2.95675687e-01, 2.47939959e+00, 2.84930818e+00, ...,
2.99674579e-01, 1.48585951e+00, 7.49068264e+00],
[ 8.84156884e-01, 1.65514790e+00, 4.17385519e-01, ...,
7.54942692e+00, 9.99801389e-01, 2.03835451e-01]])
To plot the DynamicalPowerspectrum matrix, we use the attributes time
and freq
to set the extend of the image axis. have a look at the documentation of matplotlib’s imshow()
.
[10]:
extent = min(dynspec.time), max(dynspec.time), max(dynspec.freq), min(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower left", aspect="auto", vmin=1.98, vmax=3.0,
interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)
[10]:
(700, 850)

[11]:
print("The dynamical powerspectrun has {} frequency bins and {} time bins".format(len(dynspec.freq), len(dynspec.time)))
The dynamical powerspectrun has 65535 frequency bins and 104 time bins
# Rebinning in Frequency
[12]:
print("The current frequency resolution is {}".format(dynspec.df))
The current frequency resolution is 0.0625
Let’s rebin to a frequency resolution of 2 Hz and using the average of the power
[13]:
dynspec.rebin_frequency(df_new=2.0, method="average")
[14]:
print("The new frequency resolution is {}".format(dynspec.df))
The new frequency resolution is 2.0
Let’s see how the Dynamical Powerspectrum looks now
[15]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=1.98, vmax=3.0,
interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(500, 1000)
[15]:
(500, 1000)

[16]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)
[16]:
(700, 850)

Rebin time¶
Let’s try to improve the visualization by rebinnin our matrix in the time axis
[17]:
print("The current time resolution is {}".format(dynspec.dt))
The current time resolution is 16.0
Let’s rebin to a time resolution of 64 s
[18]:
dynspec.rebin_time(dt_new=64.0, method="average")
[19]:
print("The new time resolution is {}".format(dynspec.dt))
The new time resolution is 64.0
[20]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)
[20]:
(700, 850)

Trace maximun¶
Let’s use the method trace_maximum()
to find the index of the maximum on each powerspectrum in a certain frequency range. For example, between 755 and 782Hz)
[21]:
tracing = dynspec.trace_maximum(min_freq=755, max_freq=782)
This is how the trace function looks like
[22]:
plt.plot(dynspec.time, dynspec.freq[tracing], color='red', alpha=1)
plt.show()

Let’s plot it on top of the dynamic spectrum
[23]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
interpolation="none", extent=extent, alpha=0.7)
plt.colorbar()
plt.ylim(740,800)
plt.plot(dynspec.time, dynspec.freq[tracing], color='red', lw=3, alpha=1)
plt.show()

The spike at 400 Hz is probably a statistical fluctutations, tracing by the maximum power can be dangerous!
We will implement better methods in the future, stay tunned ;)