Cross Spectra

This tutorial shows how to make and manipulate a cross spectrum of two light curves using Stingray.

import numpy as np
from stingray import Lightcurve, Crossspectrum, AveragedCrossspectrum

import matplotlib.pyplot as plt
import matplotlib.font_manager as font_manager
%matplotlib inline
font_prop = font_manager.FontProperties(size=16)

1. Create two light curves

There are two ways to make Lightcurve objects. We’ll show one way here. Check out “Lightcurve/Lightcurve tutorial.ipynb” for more examples.

Generate an array of relative timestamps that’s 8 seconds long, with dt = 0.03125 s, and make two signals in units of counts. The first is a sine wave with amplitude = 300 cts/s, frequency = 2 Hz, phase offset = 0 radians, and mean = 1000 cts/s. The second is a sine wave with amplitude = 200 cts/s, frequency = 2 Hz, phase offset = pi/4 radians, and mean = 900 cts/s. We then add Poisson noise to the light curves.

dt = 0.03125  # seconds
exposure = 8.  # seconds
times = np.arange(0, exposure, dt)  # seconds

signal_1 = 300 * np.sin(2.*np.pi*times/0.5) + 1000  # counts/s
signal_2 = 200 * np.sin(2.*np.pi*times/0.5 + np.pi/4) + 900  # counts/s
noisy_1 = np.random.poisson(signal_1*dt)  # counts
noisy_2 = np.random.poisson(signal_2*dt)  # counts

Now let’s turn noisy_1 and noisy_2 into Lightcurve objects.

lc1 = Lightcurve(times, noisy_1)
lc2 = Lightcurve(times, noisy_2)

Here we’re plotting them to see what they look like.

fig, ax = plt.subplots(1,1,figsize=(10,6))
ax.plot(lc1.time, lc1.counts, lw=2, color='blue')
ax.plot(lc1.time, lc2.counts, lw=2, color='red')
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.show()

2. Pass both of the light curves to the Crossspectrum class to create a Crossspectrum object.

The first Lightcurve passed is the channel of interest or interest band, and the second Lightcurve passed is the reference band. You can also specify the optional attribute norm if you wish to normalize the real part of the cross spectrum to squared fractional rms, Leahy, or squared absolute normalization. The default normalization is ‘frac’.

cs = Crossspectrum.from_lightcurve(lc1, lc2)
print(cs)

<stingray.crossspectrum.Crossspectrum object at 0x7f7aa3d518b0>

Note that, in principle, the Crossspectrum object could have been initialized directly as

ps = Crossspectrum(lc1, lc2, norm="leahy")

However, we recommend using the specific method for input light curve objects used above, for clarity. Equivalently, one can initialize a Crossspectrum object:

1. from EventList objects as

   bin_time = 0.1
   ps = Crossspectrum.from_events(events1, events2, dt=bin_time, norm="leahy")
   

   where the light curves, uniformly binned at 0.1 s, are created internally.

2. from numpy arrays of times, as

   bin_time = 0.1
   ps = Crossspectrum.from_events(times1, times2, dt=bin_time, gti=[[t0, t1], [t2, t3], ...], norm="leahy")
   

   where the light curves, uniformly binned at 0.1 s in this case, are created internally, and the good time intervals (time interval where the instrument was collecting data nominally) are passed by hand. Note that the frequencies of the cross spectrum will be expressed in inverse units as the input time arrays. If the times are expressed in seconds, frequencies will be in Hz; with times in days, frequencies will be in 1/d, and so on. We do not support units (e.g. astropy units) yet, so the user should pay attention to these details.

3. from an iterable of light curves

   ps = Crossspectrum.from_lc_iter(lc_iterable1, lc_iterable2, dt=bin_time, norm="leahy")
   

   where lc_iterableX is any iterable of Lightcurve objects (list, tuple, generator, etc.) and dt is the sampling time of the light curves. Note that this dt is needed because the iterables might be generators, in which case the light curves are lazy-loaded after a bunch of operations using dt have been done.

We can print the first five values in the arrays of the positive Fourier frequencies and the cross power. The cross power has a real and an imaginary component.

print(cs.freq[0:5])
print(cs.power[0:5])

[0.125 0.25  0.375 0.5   0.625]
[-3264.54599394-1077.46450232j  1066.6390401 -2783.16358879j
  3275.00416926 +196.64355198j -8345.12445869-6661.52326503j
  5916.3705245 +3602.05210672j]

Since the negative Fourier frequencies (and their associated cross powers) are discarded, the number of time bins per segment n is twice the length of freq and power.

print("Size of positive Fourier frequencies: %d" % len(cs.freq))
print("Number of data points per segment: %d" % cs.n)

Size of positive Fourier frequencies: 127
Number of data points per segment: 256

Properties

A Crossspectrum object has the following properties :

1. freq : Numpy array of mid-bin frequencies that the Fourier transform samples.

2. power : Numpy array of the cross spectrum (complex numbers).

3. df : The frequency resolution.

4. m : The number of cross spectra averaged together. For a Crossspectrum of a single segment, m=1.

5. n : The number of data points (time bins) in one segment of the light curves.

6. nphots1 : The total number of photons in the first (interest) light curve.

7. nphots2 : The total number of photons in the second (reference) light curve.

We can compute the amplitude of the cross spectrum, and plot it as a function of Fourier frequency. Notice how there’s a spike at our signal frequency of 2 Hz!

cs_amplitude = np.abs(cs.power)  # The mod square of the real and imaginary components

fig, ax1 = plt.subplots(1,1,figsize=(9,6), sharex=True)
ax1.plot(cs.freq, cs_amplitude, lw=2, color='blue')
ax1.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax1.set_ylabel("Cross spectral amplitude", fontproperties=font_prop)
ax1.set_yscale('log')
ax1.tick_params(axis='x', labelsize=16)
ax1.tick_params(axis='y', labelsize=16)
ax1.tick_params(which='major', width=1.5, length=7)
ax1.tick_params(which='minor', width=1.5, length=4)
for axis in ['top', 'bottom', 'left', 'right']:
    ax1.spines[axis].set_linewidth(1.5)
plt.show()

You’ll notice that the cross spectrum is a bit noisy. This is because we’re only using one segment of data. Let’s try averaging together multiple segments of data. # Averaged cross spectrum example You could use two long Lightcurves and have AveragedCrossspectrum chop them into specified segments, or give two lists of Lightcurves where each segment of Lightcurve is the same length. We’ll show the first way here. Remember to check the Lightcurve tutorial notebook for fancier ways of making light curves. ## 1. Create two long light curves. Generate an array of relative timestamps that’s 1600 seconds long, and two signals in count rate units, with the same properties as the previous example. We then add Poisson noise and turn them into Lightcurve objects.

long_dt = 0.03125  # seconds
long_exposure = 1600.  # seconds
long_times = np.arange(0, long_exposure, long_dt)  # seconds

# In count rate units here
long_signal_1 = 300 * np.sin(2.*np.pi*long_times/0.5) + 1000  # counts/s
long_signal_2 = 200 * np.sin(2.*np.pi*long_times/0.5 + np.pi/4) + 900  # counts/s

# Multiply by dt to get count units, then add Poisson noise
long_noisy_1 = np.random.poisson(long_signal_1*dt)  # counts
long_noisy_2 = np.random.poisson(long_signal_2*dt)  # counts

long_lc1 = Lightcurve(long_times, long_noisy_1)
long_lc2 = Lightcurve(long_times, long_noisy_2)

fig, ax = plt.subplots(1,1,figsize=(10,6))
ax.plot(long_lc1.time, long_lc1.counts, lw=2, color='blue')
ax.plot(long_lc1.time, long_lc2.counts, lw=2, color='red')
ax.set_xlim(0,20)
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.show()

2. Pass both light curves to the AveragedCrossspectrum class with a specified segment_size.

If the exposure (length) of the light curve cannot be divided by segment_size with a remainder of zero, the last incomplete segment is thrown out, to avoid signal artefacts. Here we’re using 8 second segments.

avg_cs = AveragedCrossspectrum.from_lightcurve(long_lc1, long_lc2, 8.)

200it [00:00, 12346.54it/s]

Note that also the AveragedCrossspectrum object could have been initialized using different input types:

1. from EventList objects as

   bin_time = 0.1
   ps = AveragedCrossspectrum.from_events(
       events1, events2, dt=bin_time, segment_size=segment_size,
       norm="leahy")
   

   (note, again, the necessity of the bin time)

2. from numpy arrays of times, as

   bin_time = 0.1
   ps = AveragedCrossspectrum.from_events(
       times1, times2, dt=bin_time, segment_size=segment_size,
       gti=[[t0, t1], [t2, t3], ...], norm="leahy")
   

   where the light curves, uniformly binned at 0.1 s in this case, are created internally, and the good time intervals (time interval where the instrument was collecting data nominally) are passed by hand. Note that the frequencies of the cross spectrum will be expressed in inverse units as the input time arrays. If the times are expressed in seconds, frequencies will be in Hz; with times in days, frequencies will be in 1/d, and so on. We do not support units (e.g. astropy units) yet, so the user should pay attention to these details.

3. from iterables of light curves

   ps = AveragedCrossspectrum.from_lc_iter(
       lc_iterable1, lc_iterable2, dt=bin_time, segment_size=segment_size,
       norm="leahy")
   

   where lc_iterableX is any iterable of Lightcurve objects (list, tuple, generator, etc.) and dt is the sampling time of the light curves. Note that this dt is needed because the iterables might be generators, in which case the light curves are lazy-loaded after a bunch of operations using dt have been done.

Again we can print the first five Fourier frequencies and first five cross spectral values, as well as the number of segments.

print(avg_cs.freq[0:5])
print(avg_cs.power[0:5])
print("\nNumber of segments: %d" % avg_cs.m)

[0.125 0.25  0.375 0.5   0.625]
[291.76338464-640.48290689j 182.72485752 -35.81942269j
 293.42490539+276.16187738j 771.98935476-595.89062793j
 361.32859119-101.50371039j]

Number of segments: 200

If m is less than 50 and you try to compute the coherence, a warning will pop up letting you know that your number of segments is significantly low, so the error on coherence might not follow the expected (Gaussian) statistical distributions.

test_cs = AveragedCrossspectrum.from_lightcurve(long_lc1, long_lc2, 40.)
print(test_cs.m)
coh, err = test_cs.coherence()

40it [00:00, 7645.47it/s]

40

Properties

An AveragedCrossspectrum object has the following properties, same as Crossspectrum :

1. freq : Numpy array of mid-bin frequencies that the Fourier transform samples.

2. power : Numpy array of the averaged cross spectrum (complex numbers).

3. df : The frequency resolution (in Hz).

4. m : The number of cross spectra averaged together, equal to the number of whole segments in a light curve.

5. n : The number of data points (time bins) in one segment of the light curves.

6. nphots1 : The total number of photons in the first (interest) light curve.

7. nphots2 : The total number of photons in the second (reference) light curve.

Let’s plot the amplitude of the averaged cross spectrum!

avg_cs_amplitude = np.abs(avg_cs.power)

fig, ax1 = plt.subplots(1,1,figsize=(9,6))
ax1.plot(avg_cs.freq, avg_cs_amplitude, lw=2, color='blue')
ax1.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax1.set_ylabel("Cross spectral amplitude", fontproperties=font_prop)
ax1.set_yscale('log')
ax1.tick_params(axis='x', labelsize=16)
ax1.tick_params(axis='y', labelsize=16)
ax1.tick_params(which='major', width=1.5, length=7)
ax1.tick_params(which='minor', width=1.5, length=4)
for axis in ['top', 'bottom', 'left', 'right']:
    ax1.spines[axis].set_linewidth(1.5)
plt.show()

Now we’ll show examples of all the things you can do with a Crossspectrum or AveragedCrossspectrum object using built-in stingray methods.

Normalizating the cross spectrum

The three kinds of normalization are: * leahy: Leahy normalization. Makes the Poisson noise level \(= 2\). See Leahy et al. 1983, ApJ, 266, 160L. * frac: Fractional rms-squared normalization, also known as rms normalization. Makes the Poisson noise level \(= 2 / \sqrt(meanrate_1\times meanrate_2)\). See Belloni & Hasinger 1990, A&A, 227, L33, and Miyamoto et al. 1992, ApJ, 391, L21.. This is the default. * abs: Absolute rms-squared normalization, also known as absolute normalization. Makes the Poisson noise level \(= 2 \times \sqrt(meanrate_1\times meanrate_2)\). See insert citation. * none: No normalization applied.

Note that these normalizations and the Poisson noise levels apply to the “cross power”, not the cross-spectral amplitude.

avg_cs_leahy = AveragedCrossspectrum.from_lightcurve(long_lc1, long_lc2, 8., norm='leahy')
avg_cs_frac = AveragedCrossspectrum.from_lightcurve(long_lc1, long_lc2, 8., norm='frac')
avg_cs_abs = AveragedCrossspectrum.from_lightcurve(long_lc1, long_lc2, 8., norm='abs')

200it [00:00, 15141.07it/s]
200it [00:00, 12807.43it/s]
200it [00:00, 13023.36it/s]

Here we plot the three normalized averaged cross spectra.

fig, [ax1, ax2, ax3] = plt.subplots(3,1,figsize=(6,12))
ax1.plot(avg_cs_leahy.freq, avg_cs_leahy.power, lw=2, color='black')
ax1.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax1.set_ylabel("Leahy cross-power", fontproperties=font_prop)
ax1.set_yscale('log')
ax1.tick_params(axis='x', labelsize=14)
ax1.tick_params(axis='y', labelsize=14)
ax1.tick_params(which='major', width=1.5, length=7)
ax1.tick_params(which='minor', width=1.5, length=4)
ax1.set_title("Leahy norm.", fontproperties=font_prop)

ax2.plot(avg_cs_frac.freq, avg_cs_frac.power, lw=2, color='black')
ax2.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax2.set_ylabel("rms cross-power", fontproperties=font_prop)
ax2.tick_params(axis='x', labelsize=14)
ax2.tick_params(axis='y', labelsize=14)
ax2.set_yscale('log')
ax2.tick_params(which='major', width=1.5, length=7)
ax2.tick_params(which='minor', width=1.5, length=4)
ax2.set_title("Fractional rms-squared norm.", fontproperties=font_prop)

ax3.plot(avg_cs_abs.freq, avg_cs_abs.power, lw=2, color='black')
ax3.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax3.set_ylabel("Absolute cross-power", fontproperties=font_prop)
ax3.tick_params(axis='x', labelsize=14)
ax3.tick_params(axis='y', labelsize=14)
ax3.set_yscale('log')
ax3.tick_params(which='major', width=1.5, length=7)
ax3.tick_params(which='minor', width=1.5, length=4)
ax3.set_title("Absolute rms-squared norm.", fontproperties=font_prop)

for axis in ['top', 'bottom', 'left', 'right']:
    ax1.spines[axis].set_linewidth(1.5)
    ax2.spines[axis].set_linewidth(1.5)
    ax3.spines[axis].set_linewidth(1.5)
plt.tight_layout()
plt.show()

Re-binning a cross spectrum in frequency

Typically, rebinning is done on an averaged, normalized cross spectrum. ## 1. We can linearly re-bin a cross spectrum (although this is not done much in practice)

print("DF before:", avg_cs.df)
# Both of the following ways are allowed syntax:
# lin_rb_cs = Crossspectrum.rebin(avg_cs, 0.25, method='mean')
lin_rb_cs = avg_cs.rebin(0.25, method='mean')
print("DF after:", lin_rb_cs.df)

DF before: 0.125
DF after: 0.25

2. And we can logarithmically/geometrically re-bin a cross spectrum

In this re-binning, each bin size is 1+f times larger than the previous bin size, where f is user-specified and normally in the range 0.01-0.1. The default value is f=0.01.

Logarithmic rebinning only keeps the real part of the cross spectum.

# Both of the following ways are allowed syntax:
# log_rb_cs, log_rb_freq, binning = Crossspectrum.rebin_log(avg_cs, f=0.02)
log_rb_cs = avg_cs.rebin_log(f=0.02)

Note that like rebin, rebin_log returns a Crossspectrum or AveragedCrossspectrum object (depending on the input object):

print(type(lin_rb_cs))

<class 'stingray.crossspectrum.AveragedCrossspectrum'>

Time lags / phase lags

1. Frequency-dependent lags

The lag-frequency spectrum shows the time lag between two light curves (usually non-overlapping broad energy bands) as a function of Fourier frequency. See Uttley et al. 2014, A&ARev, 22, 72 section 2.2.1.

In AveragedCrossspectrum, the second light curve is what is considered the reference in Uttley et al. and in most other spectral timing literature.

long_dt = 0.0015231682473469295763529  # seconds
long_exposure = 1600.  # seconds
long_times = np.arange(0, long_exposure, long_dt)  # seconds
frequency = 3.
phase_lag = np.pi / 3

# long_signal_1 = 300 * np.sin(2.*np.pi*long_times/0.5) + 100 * np.sin(2.*np.pi*long_times*5 + np.pi/6) + 1000
# long_signal_2 = 200 * np.sin(2.*np.pi*long_times/0.5 + np.pi/4) + 80 * np.sin(2.*np.pi*long_times*5) + 900

long_signal_1 = (300 * np.sin(2.*np.pi*long_times*frequency) + 1000) * dt
long_signal_2 = (200 * np.sin(2.*np.pi*long_times*frequency - phase_lag) + 900) * dt

long_lc1 = Lightcurve(long_times, np.random.normal(long_signal_1, 0.03))
long_lc2 = Lightcurve(long_times, np.random.normal(long_signal_2, 0.03))

# Note: the second light curve is what we use as a reference.
avg_cs = AveragedCrossspectrum.from_lightcurve(long_lc2, long_lc1, 53.)

fig, ax = plt.subplots(1,1,figsize=(10,6))
ax.plot(long_lc1.time, long_lc1.counts, lw=2, color='blue')
ax.plot(long_lc1.time, long_lc2.counts, lw=2, color='red')
ax.set_xlim(0,4)
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)
plt.show()

fig, ax = plt.subplots(1,1,figsize=(10,6))
ax.plot(avg_cs.freq, avg_cs.power, lw=2, color='blue')
plt.semilogy()
plt.show()

30it [00:00, 264.86it/s]

The time_lag method returns an np.ndarray with the time lag in seconds per positive Fourier frequency.

freq_lags, freq_lags_err = avg_cs.time_lag()
freq_plags, freq_plags_err = avg_cs.phase_lag()

# Expected time lag, given the input time lag
time_lag = phase_lag / (2. * np.pi * avg_cs.freq)

And this is a plot of the lag-frequency spectrum:

fig, ax = plt.subplots(1,1,figsize=(8,5))
ax.hlines(0, avg_cs.freq[0], avg_cs.freq[-1], color='black', linestyle='dashed', lw=2)
ax.errorbar(avg_cs.freq, freq_lags, yerr=freq_lags_err,fmt="o", lw=1, color='blue')
ax.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax.set_ylabel("Time lag (s)", fontproperties=font_prop)
ax.tick_params(axis='x', labelsize=14)
ax.tick_params(axis='y', labelsize=14)
ax.tick_params(which='major', width=1.5, length=7)
ax.tick_params(which='minor', width=1.5, length=4)
for axis in ['top', 'bottom', 'left', 'right']:
    ax.spines[axis].set_linewidth(1.5)
# plt.semilogx()
plt.axvline(frequency)
plt.xlim([2, 5])
plt.ylim([-0.05, 0.2])
plt.plot(avg_cs.freq, time_lag, label="Input time lag", lw=2, zorder=10)
plt.legend()
plt.show()

fig, ax = plt.subplots(1,1,figsize=(8,5))
ax.hlines(0, avg_cs.freq[0], avg_cs.freq[-1], color='black', linestyle='dashed', lw=2)
ax.errorbar(avg_cs.freq, freq_plags, yerr=freq_plags_err,fmt="o", lw=1, color='blue')
ax.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax.set_ylabel("Phase lag (rad)", fontproperties=font_prop)
ax.tick_params(axis='x', labelsize=14)
ax.tick_params(axis='y', labelsize=14)
ax.tick_params(which='major', width=1.5, length=7)
ax.tick_params(which='minor', width=1.5, length=4)
for axis in ['top', 'bottom', 'left', 'right']:
    ax.spines[axis].set_linewidth(1.5)
# plt.semilogx()
plt.axvline(frequency)
plt.xlim([2, 5])
plt.ylim([0, np.pi/ 2])
plt.axhline(phase_lag, label="Input phase lag", lw=2, zorder=10)
plt.legend()
plt.show()

2. Energy-dependent lags

The lag vs energy spectrum can be calculated using the LagEnergySpectrum from stingray.varenergy. Refer to the Spectral Timing documentation.

Coherence

Coherence is a Fourier-frequency-dependent measure of the linear correlation between time series measured simultaneously in two energy channels. See Vaughan and Nowak 1997, ApJ, 474, L43 and Uttley et al. 2014, A&ARev, 22, 72 section 2.1.3.

long_dt = 0.03125  # seconds
long_exposure = 1600.  # seconds
long_times = np.arange(0, long_exposure, long_dt)  # seconds

long_signal_1 = 300 * np.sin(2.*np.pi*long_times/0.5) + 1000
long_signal_2 = 200 * np.sin(2.*np.pi*long_times/0.5 + np.pi/4) + 900

long_noisy_1 = np.random.poisson(long_signal_1*dt)
long_noisy_2 = np.random.poisson(long_signal_2*dt)

long_lc1 = Lightcurve(long_times, long_noisy_1)
long_lc2 = Lightcurve(long_times, long_noisy_2)

avg_cs = AveragedCrossspectrum.from_lightcurve(long_lc1, long_lc2, 8.)

200it [00:00, 14681.05it/s]

The coherence method returns two np.ndarrays, of the coherence and uncertainty.

coh, err_coh = avg_cs.coherence()

The coherence and uncertainty have the same length as the positive Fourier frequencies.

print(len(coh) == len(avg_cs.freq))

True

And we can plot the coherence vs the frequency.

fig, ax = plt.subplots(1,1,figsize=(8,5))
# ax.hlines(0, avg_cs.freq[0], avg_cs.freq[-1], color='black', linestyle='dashed', lw=2)
ax.errorbar(avg_cs.freq, coh, yerr=err_coh, lw=2, color='blue')
ax.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax.set_ylabel("Coherence", fontproperties=font_prop)
ax.tick_params(axis='x', labelsize=14)
ax.tick_params(axis='y', labelsize=14)
ax.tick_params(which='major', width=1.5, length=7)
ax.tick_params(which='minor', width=1.5, length=4)
for axis in ['top', 'bottom', 'left', 'right']:
    ax.spines[axis].set_linewidth(1.5)
plt.show()