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

/Users/danielahuppenkothen/work/repositories/stingraysoftware/stingray/stingray/utils.py:21: UserWarning: Numba not installed. Faking it
warnings.warn("Numba not installed. Faking it")


## 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 ‘none’.

:

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

<stingray.crossspectrum.Crossspectrum object at 0x109779be0>

/Users/danielahuppenkothen/work/repositories/stingraysoftware/stingray/stingray/utils.py:74: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
warnings.warn("SIMON says: {0}".format(message), **kwargs)


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]
[ 1824.68118955+6010.74895924j -6353.40667544-2232.44448393j
-2341.92791487-5291.46491608j  5302.35901464-2806.81031101j
5172.60185652-3493.8451835j ]


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(long_lc1, long_lc2, 8.)

/Users/danielahuppenkothen/work/repositories/stingraysoftware/stingray/stingray/utils.py:74: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
warnings.warn("SIMON says: {0}".format(message), **kwargs)


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]
[ 105.15447941-102.47106937j -179.55852266-642.02444467j
-92.33657871 +60.11303608j -233.43523331-995.41309285j
78.42786218 -35.42395344j]

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(long_lc1, long_lc2, 40.)
print(test_cs.m)
coh, err = test_cs.coherence()

/Users/danielahuppenkothen/work/repositories/stingraysoftware/stingray/stingray/utils.py:74: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
warnings.warn("SIMON says: {0}".format(message), **kwargs)

40

/Users/danielahuppenkothen/work/repositories/stingraysoftware/stingray/stingray/utils.py:74: UserWarning: SIMON says: Number of segments used in averaging is significantly low. The result might not follow the expected statistical distributions.
warnings.warn("SIMON says: {0}".format(message), **kwargs)


# 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. * 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. This is the default.

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

:

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

/Users/danielahuppenkothen/work/repositories/stingraysoftware/stingray/stingray/utils.py:74: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
warnings.warn("SIMON says: {0}".format(message), **kwargs)


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.

:

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) + 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*3) + 1000) * dt
long_signal_2 = (200 * np.sin(2.*np.pi*long_times*3 + np.pi/3) + 900) * dt

long_lc1 = Lightcurve(long_times, long_signal_1)
long_lc2 = Lightcurve(long_times, long_signal_2)

avg_cs = AveragedCrossspectrum(long_lc1, long_lc2, 8.)

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.show()

/Users/danielahuppenkothen/work/repositories/stingraysoftware/stingray/stingray/utils.py:74: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
warnings.warn("SIMON says: {0}".format(message), **kwargs) /opt/local/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/numpy/core/numeric.py:531: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order) 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()

unnorm_power: [  3.51156230e-18 -5.84256529e-18j   2.60183217e-18 -3.93443595e-18j
4.08649368e-19 -4.19889793e-19j   2.36261362e-19 -1.57353749e-19j
6.02258895e-19 +5.16342834e-19j   3.46071523e-19 +2.39133138e-19j
3.88905078e-19 +3.49128137e-19j   3.76555157e-19 -1.04613037e-19j
8.57453958e-19 -7.16777928e-19j   4.74163606e-18 +6.22998684e-18j
5.55686678e-19 -4.93977389e-19j   1.96285729e-18 -3.29062895e-18j
1.35465281e-17 -2.28012208e-17j   7.34083227e-19 -1.04162844e-18j
2.70601379e-19 -3.31166964e-19j   2.36659779e-19 -2.79068969e-19j
1.38438881e-19 -3.88588884e-20j   2.09097473e-19 +1.22522937e-19j
1.40497006e-19 +3.01781468e-20j   2.53759710e-19 +2.83428365e-19j
1.42151435e-17 +2.39376102e-17j   2.89479562e-19 +3.55108965e-19j
1.99837536e-19 +2.12747336e-19j   9.60000000e+07 -1.66276878e+08j
8.71832062e-20 -2.52008051e-20j   6.84840390e-20 -2.14862443e-20j
1.07345188e-19 -1.37523646e-19j   6.88536095e-20 -1.14344617e-20j
1.57832102e-19 +1.05298764e-19j   9.17950252e-20 -4.24881813e-20j
1.66012115e-19 +4.27258685e-20j   1.30027117e-18 +1.51227952e-18j
5.00279697e-19 +2.57625930e-19j   6.72832457e-19 -8.03750359e-19j
1.29265914e-17 -2.27832436e-17j   2.31334358e-18 -3.00431708e-18j
9.76635813e-19 +3.30669917e-20j   5.81726520e-19 -9.36203498e-19j
7.88415887e-19 -1.06221780e-18j   3.70700018e-19 -2.72935648e-19j
1.90091781e-19 -1.34239411e-19j   1.57104830e-19 -8.10779819e-20j
4.15999461e-19 +1.68123908e-19j   9.36589175e-19 +1.07747658e-18j
4.03826883e-19 -4.13456354e-19j   2.40768556e-18 -4.05374087e-18j
3.71937530e-18 -5.94013540e-18j   4.87589707e-19 -5.57306638e-19j
2.83806434e-19 -5.82049357e-20j   3.89471163e-19 +1.43466979e-19j
9.38678814e-19 +1.24743600e-18j   2.23041969e-18 +3.77656911e-18j
6.92637827e-19 -4.28709692e-19j   2.93570859e-19 -2.88708689e-19j
7.29208394e-19 +1.24230061e-19j   2.29538858e-19 -2.52901513e-19j
2.80334115e-19 -3.58322764e-19j   4.27406930e-18 -7.13382265e-18j
1.97070335e-19 -2.22893692e-19j   1.10498272e-19 -1.03413708e-20j
1.39194731e-19 -8.73750345e-20j   9.21053463e-20 -5.07763852e-21j
1.99515929e-19 +1.99494842e-19j   1.32718699e-19 +1.01582536e-19j
3.50077510e-19 +3.86965968e-19j   9.95095373e-18 +1.65238467e-17j
5.72489059e-19 +7.80740130e-19j   2.96319586e-19 -9.89700691e-20j
1.37786640e-17 -2.38878851e-17j   3.05178081e-19 -3.68359591e-19j
2.42261593e-19 -2.05908439e-19j   1.17616106e-19 -9.65815734e-20j
8.01555479e-20 -1.33438480e-20j   1.94292533e-19 +1.97990901e-19j
6.89528984e-20 +3.01360432e-20j   7.33355585e-20 -3.86762730e-21j
2.25276850e-19 +3.50311664e-21j   1.51700398e-19 +1.26606190e-19j
1.45809708e-19 -5.47471092e-20j   1.16906974e-18 -1.61212904e-18j
1.20606568e-18 +7.54324248e-19j   3.88117315e-19 +5.50058656e-20j
1.56063042e-19 -8.36770346e-20j   3.06934112e-19 -2.49069019e-19j
6.15941081e-19 -5.51899503e-19j   1.59087246e-19 +1.25386968e-19j
8.38132281e-20 -8.53141653e-21j   1.61531773e-19 -4.02963098e-20j
3.80182166e-19 +5.30779147e-19j   9.49269102e-20 -2.09335914e-20j
2.65989813e-19 -3.92960759e-19j   8.59370550e-19 -1.25143301e-18j
1.35628901e-19 -7.76498659e-20j   1.26499316e-19 -2.31562116e-20j
3.15697566e-19 +2.07501418e-19j   2.93905657e-19 +2.78045784e-19j
1.30087468e-18 +2.03166513e-18j   3.56389874e-19 -2.35271359e-20j
8.17631916e-19 -1.32498365e-18j   2.66631256e-18 -3.23949452e-18j
2.56309783e-19 -2.89847628e-19j   1.02399399e-19 -9.41152715e-20j
4.32334848e-19 -6.23760411e-19j   7.67444673e-20 -1.33350113e-20j
7.54584646e-20 -1.00486173e-20j   8.40135882e-20 -5.64591672e-20j
7.37394299e-20 +1.99426657e-20j   4.22454107e-19 +5.70392427e-19j
8.81135898e-20 +2.78749443e-21j   2.59074497e-19 +2.87633700e-19j
1.52056301e-18 +1.96897478e-18j   2.56570835e-19 +2.21945494e-19j
4.78498190e-19 -1.10506973e-19j   9.55112397e-18 -1.66594160e-17j
6.06079437e-19 -8.33038601e-19j   2.16824128e-19 -1.64072102e-19j
1.64419616e-19 -1.92180213e-19j   1.00688482e-19 +5.88060237e-21j
4.27836181e-19 +5.02867773e-19j   9.83835409e-20 +5.14334458e-20j
8.08348323e-20 +1.34161082e-20j   2.68025800e-19 -7.74202404e-20j
1.86717882e-19 +1.67990114e-19j   9.49143760e-20 +6.78328565e-20j
1.71274886e-19 -1.18928604e-19j   3.51008234e-19 +1.63633300e-19j
9.38196620e-19 +1.28452493e-18j]
ph_lag: [-1.02961915 -0.98652014 -0.79896389 -0.58755172  0.70874201  0.60465566
0.73155439 -0.27098231 -0.6962732   0.92022969 -0.72667619 -1.03295497
-1.03471589 -0.95689239 -0.88570639 -0.86744514 -0.27365163  0.53003265
0.21158084  0.84057166  1.03491753  0.88686448  0.816678   -1.04719755
-0.28138618 -0.30401497 -0.90802317 -0.16456727  0.58834187 -0.43349612
0.25189923  0.86063527  0.47554685 -0.87382968 -1.05471302 -0.91461355
0.03384513 -1.01481334 -0.93228297 -0.63465645 -0.61486298 -0.47642529
0.38407414  0.85523614 -0.79717991 -1.0348434  -1.01137537 -0.8520208
-0.20228183  0.35293975  0.92570572  1.03732159 -0.55423849 -0.77704813
0.16874285 -0.83378636 -0.9069101  -1.03101862 -0.84681029 -0.09331672
-0.56055143 -0.05507285  0.78534532  0.65328278  0.83540566  1.02874773
0.93809222 -0.32234839 -1.04760621 -0.87892992 -0.70446072 -0.68751052
-0.16496162  0.79482568  0.41203485 -0.05268996  0.01554902  0.69547384
-0.35918229 -0.94337585  0.55891716  0.14078723 -0.49216675 -0.68169817
-0.73061543  0.66748258 -0.10144142 -0.24447382  0.94923293 -0.21704935
-0.97575182 -0.96904617 -0.51996632 -0.18104956  0.58147534  0.75767586
1.00127051 -0.06591952 -1.01790656 -0.8821491  -0.84672808 -0.74326784
-0.96470922 -0.1720409  -0.13238862 -0.59170252  0.26412915  0.9333141
0.03162469  0.8375892   0.91319955  0.71316427 -0.22696614 -1.0502279
-0.94181698 -0.6477791  -0.86308951  0.05833765  0.8658423   0.48170907
0.16447017 -0.28119971  0.73264943  0.62050668 -0.60693925  0.43622811
0.93996803]


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, 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.show() ## 2. Energy-dependent lags¶

:

## In development


# 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(long_lc1, long_lc2, 8.)

/Users/danielahuppenkothen/work/repositories/stingraysoftware/stingray/stingray/utils.py:74: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
warnings.warn("SIMON says: {0}".format(message), **kwargs)


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, 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() 