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()
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_[real_data]_12_0.png

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)
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_[real_data]_20_1.png
[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)
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_[real_data]_29_1.png
[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)
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_[real_data]_30_1.png

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)
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_[real_data]_37_1.png

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()
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_[real_data]_42_0.png

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()
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_[real_data]_44_0.png

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 ;)