+import copy
+import warnings
+from collections.abc import Iterable, Iterator, Generator
+
+import numpy as np
+import scipy
+import scipy.optimize
+import scipy.stats
+import matplotlib.pyplot as plt
+
+from stingray.exceptions import StingrayError
+from stingray.utils import rebin_data, rebin_data_log, simon
+
+from .base import StingrayObject
+from .events import EventList
+from .gti import cross_two_gtis, time_intervals_from_gtis
+from .lightcurve import Lightcurve
+from .fourier import avg_cs_from_iterables, error_on_averaged_cross_spectrum
+from .fourier import avg_cs_from_timeseries, poisson_level
+from .fourier import normalize_periodograms, raw_coherence
+from .fourier import get_flux_iterable_from_segments
+from .fourier import get_rms_from_unnorm_periodogram
+from .power_colors import power_color
+
+from scipy.special import factorial
+
+
+__all__ = [
+ "Crossspectrum",
+ "AveragedCrossspectrum",
+ "DynamicalCrossspectrum",
+ "cospectra_pvalue",
+ "normalize_crossspectrum",
+ "time_lag",
+ "coherence",
+ "get_flux_generator",
+]
+
+
+def get_flux_generator(data, segment_size, dt=None):
+ """Get a flux generator from different segments of a data object
+
+ It is just a wrapper around
+ ``stingray.fourier.get_flux_iterable_from_segments``, providing
+ this method with the information it needs to create the iterables,
+ starting from an event list or a light curve.
+
+ Only accepts `Lightcurve`s and `EventList`s.
+
+ Parameters
+ ----------
+ data : `Lightcurve` or `EventList`
+ Input data
+ segment_size : float
+ Segment size in seconds
+
+ Other parameters
+ ----------------
+ dt : float, default None
+ Sampling time of the output flux iterables. Required if input data
+ is an event list, otherwise the light curve sampling time is selected.
+
+ Returns
+ -------
+ flux_iterable : ``generator``
+ Generator of flux arrays.
+
+ Examples
+ --------
+ >>> mean = 256
+ >>> length = 128
+ >>> times = np.sort(np.random.uniform(0, length, int(mean * length)))
+ >>> events = EventList(time=times, gti=[[0, length]])
+ >>> dt = 0.125
+ >>> segment_size = 4
+
+ Create a light curve
+ >>> lc = events.to_lc(dt=dt)
+
+ Create a light curve with a different error distribution
+ >>> lc_renorm = copy.deepcopy(lc)
+ >>> lc_renorm.counts = lc.counts / mean
+ >>> lc_renorm.counts_err = lc.counts_err / mean
+ >>> lc_renorm.err_dist = "gauss"
+
+ Create an iterable from events, forgetting ``dt``. Should fail
+ >>> get_flux_generator(events, segment_size, dt=None)
+ Traceback (most recent call last):
+ ...
+ ValueError: If data is an EventList, you need to specify...
+
+ Create an iterable from events
+ >>> iter_ev = get_flux_generator(events, segment_size, dt=dt)
+
+ Create an iterable from the light curve
+ >>> iter_lc = get_flux_generator(lc, segment_size, dt=dt)
+
+ Create an iterable from the non-poisson light curve
+ >>> iter_lc_nonpois = get_flux_generator(lc_renorm, segment_size, dt=dt)
+
+ Verify that they are equivalent
+ >>> for l1, l2 in zip(iter_ev, iter_lc): assert np.allclose(l1, l2)
+
+ Note that the iterable for non-Poissonian light curves also returns the uncertainty
+ >>> for l1, (l2, l2e) in zip(iter_lc, iter_lc_nonpois): assert np.allclose(l1, l2 * mean)
+
+ """
+ times = data.time
+ gti = data.gti
+
+ counts = err = None
+ if isinstance(data, Lightcurve):
+ counts = data.counts
+ N = counts.size
+ if data.err_dist.lower() != "poisson":
+ err = data.counts_err
+ elif isinstance(data, EventList):
+ if dt is None:
+ raise ValueError("If data is an EventList, you need to specify the bin time dt")
+ N = int(np.rint(segment_size / dt))
+
+ flux_iterable = get_flux_iterable_from_segments(
+ times, gti, segment_size, N, fluxes=counts, errors=err
+ )
+ return flux_iterable
+
+
+
+
[docs]
+
def coherence(lc1, lc2):
+
"""
+
Estimate coherence function of two light curves.
+
For details on the definition of the coherence, see Vaughan and Nowak,
+
1996 [#]_.
+
+
Parameters
+
----------
+
lc1: :class:`stingray.Lightcurve` object
+
The first light curve data for the channel of interest.
+
lc2: :class:`stingray.Lightcurve` object
+
The light curve data for reference band
+
+
Returns
+
-------
+
coh : ``np.ndarray``
+
The array of coherence versus frequency
+
+
References
+
----------
+
.. [#] https://iopscience.iop.org/article/10.1086/310430
+
"""
+
+
warnings.warn(
+
"The coherence function, as implemented, does not work as expected. "
+
"Please use the coherence function of AveragedCrossspectrum, with the "
+
"correct parameters.",
+
DeprecationWarning,
+
)
+
if not isinstance(lc1, Lightcurve):
+
raise TypeError("lc1 must be a lightcurve.Lightcurve object")
+
+
if not isinstance(lc2, Lightcurve):
+
raise TypeError("lc2 must be a lightcurve.Lightcurve object")
+
+
cs = Crossspectrum(lc1, lc2, norm="none")
+
+
return cs.coherence()
+
[docs]
+
class Crossspectrum(StingrayObject):
+
main_array_attr = "freq"
+
type = "crossspectrum"
+
+
"""
+
Make a cross spectrum from a (binned) light curve.
+
You can also make an empty :class:`Crossspectrum` object to populate with your
+
own Fourier-transformed data (this can sometimes be useful when making
+
binned power spectra). Stingray uses the scipy.fft standards for the sign
+
of the Nyquist frequency.
+
+
Parameters
+
----------
+
data1: :class:`stingray.Lightcurve` or :class:`stingray.events.EventList`, optional, default ``None``
+
The dataset for the first channel/band of interest.
+
+
data2: :class:`stingray.Lightcurve` or :class:`stingray.events.EventList`, optional, default ``None``
+
The dataset for the second, or "reference", band.
+
+
norm: {``frac``, ``abs``, ``leahy``, ``none``}, default ``none``
+
The normalization of the (real part of the) cross spectrum.
+
+
power_type: string, optional, default ``real``
+
Parameter to choose among complete, real part and magnitude of the cross spectrum.
+
+
fullspec: boolean, optional, default ``False``
+
If False, keep only the positive frequencies, or if True, keep all of them .
+
+
Other Parameters
+
----------------
+
gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+
Good Time intervals. Defaults to the common GTIs from the two input
+
objects. Could throw errors if these GTIs have overlaps with the input
+
`Lightcurve` GTIs! If you're getting errors regarding your GTIs, don't
+
use this and only give GTIs to the `Lightcurve` objects before making
+
the cross spectrum.
+
+
lc1: :class:`stingray.Lightcurve`object OR iterable of :class:`stingray.Lightcurve` objects
+
For backwards compatibility only. Like ``data1``, but no
+
:class:`stingray.events.EventList` objects allowed
+
+
lc2: :class:`stingray.Lightcurve`object OR iterable of :class:`stingray.Lightcurve` objects
+
For backwards compatibility only. Like ``data2``, but no
+
:class:`stingray.events.EventList` objects allowed
+
+
dt: float
+
The time resolution of the light curve. Only needed when constructing
+
light curves in the case where ``data1``, ``data2`` are
+
:class:`EventList` objects
+
+
skip_checks: bool
+
Skip initial checks, for speed or other reasons (you need to trust your
+
inputs!)
+
+
+
Attributes
+
----------
+
freq: numpy.ndarray
+
The array of mid-bin frequencies that the Fourier transform samples
+
+
power: numpy.ndarray
+
The array of cross spectra (complex numbers)
+
+
power_err: numpy.ndarray
+
The uncertainties of ``power``.
+
An approximation for each bin given by ``power_err= power/sqrt(m)``.
+
Where ``m`` is the number of power averaged in each bin (by frequency
+
binning, or averaging more than one spectra). Note that for a single
+
realization (``m=1``) the error is equal to the power.
+
+
df: float
+
The frequency resolution
+
+
m: int
+
The number of averaged cross-spectra amplitudes in each bin.
+
+
n: int
+
The number of data points/time bins in one segment of the light
+
curves.
+
+
k: array of int
+
The rebinning scheme if the object has been rebinned otherwise is set to 1.
+
+
nphots1: float
+
The total number of photons in light curve 1
+
+
nphots2: float
+
The total number of photons in light curve 2
+
+
"""
+
+
def __init__(
+
self,
+
data1=None,
+
data2=None,
+
norm="frac",
+
gti=None,
+
lc1=None,
+
lc2=None,
+
power_type="all",
+
dt=None,
+
fullspec=False,
+
skip_checks=False,
+
save_all=False,
+
):
+
self._type = None
+
# for backwards compatibility
+
if data1 is None:
+
data1 = lc1
+
if data2 is None:
+
data2 = lc2
+
+
empty = data1 is None and data2 is None
+
good_input = not empty
+
+
if not skip_checks:
+
good_input = self.initial_checks(
+
data1=data1,
+
data2=data2,
+
norm=norm,
+
gti=gti,
+
lc1=lc1,
+
lc2=lc2,
+
power_type=power_type,
+
dt=dt,
+
fullspec=fullspec,
+
)
+
+
self.dt = dt
+
norm = norm.lower()
+
self.norm = norm
+
self.k = 1
+
+
if empty or not good_input:
+
return self._initialize_empty()
+
+
return self._initialize_from_any_input(
+
data1,
+
data2,
+
dt=dt,
+
norm=norm,
+
power_type=power_type,
+
fullspec=fullspec,
+
gti=gti,
+
save_all=save_all,
+
)
+
+
+
[docs]
+
def initial_checks(
+
self,
+
data1=None,
+
data2=None,
+
norm="frac",
+
gti=None,
+
lc1=None,
+
lc2=None,
+
segment_size=None,
+
power_type="real",
+
dt=None,
+
fullspec=False,
+
):
+
"""Run initial checks on the input.
+
+
Returns True if checks are passed, False if they are not.
+
+
Raises various errors for different bad inputs
+
+
Examples
+
--------
+
>>> times = np.arange(0, 10)
+
>>> counts = np.random.poisson(100, 10)
+
>>> lc1 = Lightcurve(times, counts, skip_checks=True)
+
>>> lc2 = Lightcurve(times, counts, skip_checks=True)
+
>>> ev1 = EventList(times)
+
>>> ev2 = EventList(times)
+
>>> c = Crossspectrum()
+
>>> ac = AveragedCrossspectrum()
+
+
If norm is not a string, raise a TypeError
+
>>> Crossspectrum.initial_checks(c, norm=1)
+
Traceback (most recent call last):
+
...
+
TypeError: norm must be a string...
+
+
If ``norm`` is not one of the valid norms, raise a ValueError
+
>>> Crossspectrum.initial_checks(c, norm="blabla")
+
Traceback (most recent call last):
+
...
+
ValueError: norm must be 'frac'...
+
+
If ``power_type`` is not one of the valid norms, raise a ValueError
+
>>> Crossspectrum.initial_checks(c, power_type="blabla")
+
Traceback (most recent call last):
+
...
+
ValueError: `power_type` not recognized!
+
+
If the user passes only one light curve, raise a ValueError
+
+
>>> Crossspectrum.initial_checks(c, data1=lc1, data2=None)
+
Traceback (most recent call last):
+
...
+
ValueError: You can't do a cross spectrum...
+
+
If the user passes an event list without dt, raise a ValueError
+
+
>>> Crossspectrum.initial_checks(c, data1=ev1, data2=ev2, dt=None)
+
Traceback (most recent call last):
+
...
+
ValueError: If using event lists, please specify...
+
"""
+
if isinstance(norm, str) is False:
+
raise TypeError("norm must be a string")
+
+
if norm.lower() not in ["frac", "abs", "leahy", "none"]:
+
raise ValueError("norm must be 'frac', 'abs', 'leahy', or 'none'!")
+
+
if power_type not in ["all", "absolute", "real"]:
+
raise ValueError("`power_type` not recognized!")
+
+
# check if input data is a Lightcurve object, if not make one or
+
# make an empty Crossspectrum object if lc1 == ``None`` or lc2 == ``None``
+
+
if lc1 is not None or lc2 is not None:
+
warnings.warn(
+
"The lcN keywords are now deprecated. Use dataN instead", DeprecationWarning
+
)
+
+
if data1 is None or data2 is None:
+
if data1 is not None or data2 is not None:
+
raise ValueError("You can't do a cross spectrum with just one light curve!")
+
else:
+
return False
+
+
dt_is_invalid = (dt is None) or (dt <= np.finfo(float).resolution)
+
+
if segment_size is None:
+
# checks to be run for non-averaged spectra
+
if gti is not None and len(gti) > 1:
+
raise TypeError("Non-averaged cross spectra need a single GTI")
+
+
if type(data1) != type(data2):
+
raise TypeError("Input data have to be of the same kind")
+
+
if isinstance(data1, EventList):
+
if dt_is_invalid:
+
raise ValueError(
+
"If using event lists, please specify the bin time to generate lightcurves."
+
)
+
elif isinstance(data1, Lightcurve):
+
if data1.err_dist.lower() != data2.err_dist.lower():
+
simon(
+
"Your lightcurves have different statistics."
+
"The errors in the Crossspectrum will be incorrect."
+
)
+
+
# If dt differs slightly, its propagated error must not be more than
+
# 1/100th of the bin
+
if not np.isclose(data1.dt, data2.dt, rtol=0.01 * data1.dt / data1.tseg):
+
raise StingrayError("Light curves do not have same time binning dt.")
+
+
if data1.tseg != data2.tseg:
+
simon(
+
"Lightcurves do not have same tseg. This means that the data"
+
"from the two channels are not completely in sync. This "
+
"might or might not be an issue. Keep an eye on it."
+
)
+
elif isinstance(data1, (list, tuple)):
+
if not isinstance(data1[0], Lightcurve) or not isinstance(data2[0], Lightcurve):
+
raise TypeError("Inputs lists have to contain light curve objects")
+
+
if data1[0].err_dist.lower() != data2[0].err_dist.lower():
+
simon(
+
"Your lightcurves have different statistics."
+
"The errors in the Crossspectrum will be incorrect."
+
)
+
elif isinstance(data1, (Generator, Iterator)):
+
pass
+
else:
+
raise TypeError("Input data are invalid")
+
+
return True
+
+
+
+
[docs]
+
def rebin(self, df=None, f=None, method="mean"):
+
"""
+
Rebin the cross spectrum to a new frequency resolution ``df``.
+
+
Parameters
+
----------
+
df: float
+
The new frequency resolution
+
+
Other Parameters
+
----------------
+
f: float
+
the rebin factor. If specified, it substitutes df with ``f*self.df``
+
+
Returns
+
-------
+
bin_cs = :class:`Crossspectrum` (or one of its subclasses) object
+
The newly binned cross spectrum or power spectrum.
+
Note: this object will be of the same type as the object
+
that called this method. For example, if this method is called
+
from :class:`AveragedPowerspectrum`, it will return an object of class
+
:class:`AveragedPowerspectrum`, too.
+
"""
+
+
if f is None and df is None:
+
raise ValueError("You need to specify at least one between f and df")
+
elif f is not None:
+
df = f * self.df
+
+
# rebin cross spectrum to new resolution
+
binfreq, bincs, binerr, step_size = rebin_data(
+
self.freq, self.power, df, self.power_err, method=method, dx=self.df
+
)
+
# make an empty cross spectrum object
+
# note: syntax deliberate to work with subclass Powerspectrum
+
bin_cs = copy.copy(self)
+
+
# store the binned periodogram in the new object
+
bin_cs.freq = binfreq
+
bin_cs.power = bincs
+
bin_cs.df = df
+
bin_cs.power_err = binerr
+
+
if hasattr(self, "unnorm_power") and self.unnorm_power is not None:
+
unnorm_power_err = None
+
if hasattr(self, "unnorm_power_err") and self.unnorm_power_err is not None:
+
unnorm_power_err = self.unnorm_power_err
+
+
_, binpower_unnorm, binpower_err_unnorm, _ = rebin_data(
+
self.freq, self.unnorm_power, df, dx=self.df, yerr=unnorm_power_err, method=method
+
)
+
+
if hasattr(self, "unnorm_power_err") and self.unnorm_power_err is not None:
+
bin_cs.unnorm_power_err = binpower_err_unnorm
+
+
bin_cs.unnorm_power = binpower_unnorm
+
+
if hasattr(self, "cs_all"):
+
cs_all = []
+
for c in self.cs_all:
+
cs_all.append(
+
rebin_data(self.freq, c, dx_new=df, yerr=None, method=method, dx=self.df)[1]
+
)
+
bin_cs.cs_all = cs_all
+
if hasattr(self, "pds1"):
+
bin_cs.pds1 = self.pds1.rebin(df=df, f=f, method=method)
+
if hasattr(self, "pds2"):
+
bin_cs.pds2 = self.pds2.rebin(df=df, f=f, method=method)
+
+
bin_cs.m = np.rint(step_size * self.m)
+
+
return bin_cs
+
+
+
+
[docs]
+
def to_norm(self, norm, inplace=False):
+
"""Convert Cross spectrum to new normalization.
+
+
Parameters
+
----------
+
norm : str
+
The new normalization of the spectrum
+
+
Other parameters
+
----------------
+
inplace: bool, default False
+
If True, change the current instance. Otherwise, return a new one
+
+
Returns
+
-------
+
new_spec : object, same class as input
+
The new, normalized, spectrum.
+
"""
+
if norm == self.norm:
+
return copy.deepcopy(self)
+
+
variance1 = variance2 = variance = None
+
if self.type == "powerspectrum":
+
# This is the case for Powerspectrum
+
mean = mean1 = mean2 = self.nphots / self.n
+
if hasattr(self, "err_dist") and self.err_dist != "poisson":
+
variance = self.variance
+
nph = self.nphots
+
else:
+
nph = np.sqrt(self.nphots1 * self.nphots2)
+
mean1 = self.nphots1 / self.n
+
mean2 = self.nphots2 / self.n
+
mean = nph / self.n
+
if hasattr(self, "err_dist") and self.err_dist != "poisson":
+
variance1 = self.variance1
+
variance2 = self.variance2
+
variance = np.sqrt(self.variance1 * self.variance2)
+
+
if inplace:
+
new_spec = self
+
else:
+
new_spec = copy.deepcopy(self)
+
+
power_type = "all"
+
if hasattr(self, "power_type"):
+
power_type = self.power_type
+
+
for attr in ["power", "power_err"]:
+
unnorm_attr = "unnorm_" + attr
+
if not hasattr(self, unnorm_attr) or getattr(self, unnorm_attr) is None:
+
continue
+
power = normalize_periodograms(
+
getattr(self, unnorm_attr),
+
self.dt,
+
self.n,
+
mean,
+
n_ph=nph,
+
variance=variance,
+
norm=norm,
+
power_type=power_type,
+
)
+
setattr(new_spec, attr, power)
+
new_spec.norm = norm
+
if hasattr(self, "pds1"):
+
p1 = normalize_periodograms(
+
getattr(self.pds1, unnorm_attr),
+
self.dt,
+
self.n,
+
mean1,
+
n_ph=self.nphots1,
+
variance=variance1,
+
norm=norm,
+
power_type=power_type,
+
)
+
setattr(new_spec.pds1, attr, p1)
+
p2 = normalize_periodograms(
+
getattr(self.pds2, unnorm_attr),
+
self.dt,
+
self.n,
+
mean2,
+
n_ph=self.nphots2,
+
variance=variance2,
+
norm=norm,
+
power_type=power_type,
+
)
+
setattr(new_spec.pds2, attr, p2)
+
new_spec.pds1.norm = new_spec.pds2.norm = norm
+
+
return new_spec
+
+
+
def _normalize_crossspectrum(self, unnorm_power):
+
"""
+
Normalize the real part of the cross spectrum to Leahy, absolute rms^2,
+
fractional rms^2 normalization, or not at all.
+
+
Parameters
+
----------
+
unnorm_power: numpy.ndarray
+
The unnormalized cross spectrum.
+
+
Returns
+
-------
+
power: numpy.nd.array
+
The normalized co-spectrum (real part of the cross spectrum). For
+
'none' normalization, imaginary part is returned as well.
+
"""
+
+
nph = np.sqrt(self.nphots1 * self.nphots2)
+
mean = nph / self.n
+
variance = None
+
if self.err_dist != "poisson":
+
variance = np.sqrt(self.variance1 * self.variance2)
+
return normalize_periodograms(
+
unnorm_power,
+
self.dt,
+
self.n,
+
mean,
+
n_ph=nph,
+
variance=variance,
+
norm=self.norm,
+
power_type=self.power_type,
+
)
+
+
+
[docs]
+
def rebin_log(self, f=0.01):
+
"""
+
Logarithmic rebin of the periodogram.
+
The new frequency depends on the previous frequency
+
modified by a factor f:
+
+
.. math::
+
+
d\\nu_j = d\\nu_{j-1} (1+f)
+
+
Parameters
+
----------
+
f: float, optional, default ``0.01``
+
parameter that steers the frequency resolution
+
+
+
Returns
+
-------
+
new_spec : :class:`Crossspectrum` (or one of its subclasses) object
+
The newly binned cross spectrum or power spectrum.
+
Note: this object will be of the same type as the object
+
that called this method. For example, if this method is called
+
from :class:`AveragedPowerspectrum`, it will return an object of class
+
"""
+
+
binfreq, binpower, binpower_err, nsamples = rebin_data_log(
+
self.freq, self.power, f, y_err=self.power_err, dx=self.df
+
)
+
+
new_spec = copy.copy(self)
+
new_spec.freq = binfreq
+
new_spec.power = binpower
+
new_spec.power_err = binpower_err
+
new_spec.m = nsamples * self.m
+
new_spec.dt = self.dt
+
new_spec.k = nsamples
+
+
if hasattr(self, "unnorm_power") and self.unnorm_power is not None:
+
unnorm_power_err = None
+
if hasattr(self, "unnorm_power_err") and self.unnorm_power_err is not None:
+
unnorm_power_err = self.unnorm_power_err
+
_, binpower_unnorm, binpower_err_unnorm, _ = rebin_data_log(
+
self.freq, self.unnorm_power, f, dx=self.df, y_err=unnorm_power_err
+
)
+
+
new_spec.unnorm_power = binpower_unnorm
+
if hasattr(self, "unnorm_power_err") and self.unnorm_power_err is not None:
+
new_spec.unnorm_power_err = binpower_err_unnorm
+
+
if hasattr(self, "pds1"):
+
new_spec.pds1 = self.pds1.rebin_log(f)
+
if hasattr(self, "pds2"):
+
new_spec.pds2 = self.pds2.rebin_log(f)
+
+
if hasattr(self, "cs_all"):
+
cs_all = []
+
+
for c in self.cs_all:
+
cs_all.append(rebin_data_log(self.freq, c, f, dx=self.df)[1])
+
new_spec.cs_all = cs_all
+
+
return new_spec
+
+
+
+
[docs]
+
def coherence(self):
+
"""Compute Coherence function of the cross spectrum.
+
+
Coherence is defined in Vaughan and Nowak, 1996 [#]_.
+
It is a Fourier frequency dependent measure of the linear correlation
+
between time series measured simultaneously in two energy channels.
+
+
Returns
+
-------
+
coh : numpy.ndarray
+
Coherence function
+
+
References
+
----------
+
.. [#] https://iopscience.iop.org/article/10.1086/310430
+
"""
+
# this computes the averaged power spectrum, but using the
+
# cross spectrum code to avoid circular imports
+
+
return raw_coherence(
+
self.unnorm_power, self.pds1.unnorm_power, self.pds2.unnorm_power, 0, 0, self.n
+
)
+
+
+
+
[docs]
+
def phase_lag(self):
+
"""Calculate the fourier phase lag of the cross spectrum.
+
+
This is defined as the argument of the complex cross spectrum, and gives
+
the delay at all frequencies, in cycles, of one input light curve with respect
+
to the other.
+
"""
+
return np.angle(self.unnorm_power)
+
+
+
+
[docs]
+
def time_lag(self):
+
r"""Calculate the fourier time lag of the cross spectrum.
+
The time lag is calculated by taking the phase lag :math:`\phi` and
+
+
..math::
+
+
\tau = \frac{\phi}{\two pi \nu}
+
+
where :math:`\nu` is the center of the frequency bins.
+
"""
+
if self.__class__ in [Crossspectrum, AveragedCrossspectrum]:
+
ph_lag = self.phase_lag()
+
+
return ph_lag / (2 * np.pi * self.freq)
+
else:
+
raise AttributeError("Object has no attribute named 'time_lag' !")
+
+
+
+
[docs]
+
def plot(
+
self, labels=None, axis=None, title=None, marker="-", save=False, filename=None, ax=None
+
):
+
"""
+
Plot the amplitude of the cross spectrum vs. the frequency using ``matplotlib``.
+
+
Parameters
+
----------
+
labels : iterable, default ``None``
+
A list of tuple with ``xlabel`` and ``ylabel`` as strings.
+
+
axis : list, tuple, string, default ``None``
+
Parameter to set axis properties of the ``matplotlib`` figure. For example
+
it can be a list like ``[xmin, xmax, ymin, ymax]`` or any other
+
acceptable argument for the``matplotlib.pyplot.axis()`` method.
+
+
title : str, default ``None``
+
The title of the plot.
+
+
marker : str, default '-'
+
Line style and color of the plot. Line styles and colors are
+
combined in a single format string, as in ``'bo'`` for blue
+
circles. See ``matplotlib.pyplot.plot`` for more options.
+
+
save : boolean, optional, default ``False``
+
If ``True``, save the figure with specified filename.
+
+
filename : str
+
File name of the image to save. Depends on the boolean ``save``.
+
+
ax : ``matplotlib.Axes`` object
+
An axes object to fill with the cross correlation plot.
+
"""
+
+
if ax is None:
+
fig = plt.figure("crossspectrum")
+
ax = fig.add_subplot(1, 1, 1)
+
+
ax2 = None
+
if np.any(np.iscomplex(self.power)):
+
ax.plot(self.freq, np.abs(self.power), marker, color="k", label="Amplitude")
+
+
ax2 = ax.twinx()
+
ax2.tick_params("y", colors="b")
+
ax2.plot(
+
self.freq, self.power.imag, marker, color="b", alpha=0.5, label="Imaginary Part"
+
)
+
+
ax.plot(self.freq, self.power.real, marker, color="r", alpha=0.5, label="Real Part")
+
+
lines, line_labels = ax.get_legend_handles_labels()
+
lines2, line_labels2 = ax2.get_legend_handles_labels()
+
lines = lines + lines2
+
line_labels = line_labels + line_labels2
+
+
else:
+
ax.plot(self.freq, np.abs(self.power), marker, color="b")
+
lines, line_labels = ax.get_legend_handles_labels()
+
+
xlabel = "Frequency (Hz)"
+
ylabel = f"Power ({self.norm})"
+
+
if labels is not None:
+
try:
+
xlabel = labels[0]
+
ylabel = labels[1]
+
+
except IndexError:
+
simon("``labels`` must have two labels for x and y axes.")
+
# Not raising here because in case of len(labels)==1, only
+
# x-axis will be labelled.
+
+
ax.set_xlabel(xlabel)
+
if ax2 is not None:
+
ax.set_ylabel(ylabel + "-Real")
+
ax2.set_ylabel(ylabel + "-Imaginary")
+
else:
+
ax.set_ylabel(ylabel)
+
+
ax.legend(lines, line_labels, loc="best")
+
+
if axis is not None:
+
ax.set_xlim(axis[0:2])
+
ax.set_ylim(axis[2:4])
+
if ax2 is not None:
+
ax2.set_ylim(axis[2:4])
+
if title is not None:
+
ax.set_title(title)
+
+
if save:
+
if filename is None:
+
plt.gcf().savefig("spec.png")
+
else:
+
plt.gcf().savefig(filename)
+
+
return ax
+
+
+
+
[docs]
+
def classical_significances(self, threshold=1, trial_correction=False):
+
"""
+
Compute the classical significances for the powers in the power
+
spectrum, assuming an underlying noise distribution that follows a
+
chi-square distributions with 2M degrees of freedom, where M is the
+
number of powers averaged in each bin.
+
+
Note that this function will *only* produce correct results when the
+
following underlying assumptions are fulfilled:
+
+
1. The power spectrum is Leahy-normalized
+
2. There is no source of variability in the data other than the
+
periodic signal to be determined with this method. This is important!
+
If there are other sources of (aperiodic) variability in the data, this
+
method will *not* produce correct results, but instead produce a large
+
number of spurious false positive detections!
+
3. There are no significant instrumental effects changing the
+
statistical distribution of the powers (e.g. pile-up or dead time)
+
+
By default, the method produces ``(index,p-values)`` for all powers in
+
the power spectrum, where index is the numerical index of the power in
+
question. If a ``threshold`` is set, then only powers with p-values
+
*below* that threshold with their respective indices. If
+
``trial_correction`` is set to ``True``, then the threshold will be corrected
+
for the number of trials (frequencies) in the power spectrum before
+
being used.
+
+
Parameters
+
----------
+
threshold : float, optional, default ``1``
+
The threshold to be used when reporting p-values of potentially
+
significant powers. Must be between 0 and 1.
+
Default is ``1`` (all p-values will be reported).
+
+
trial_correction : bool, optional, default ``False``
+
A Boolean flag that sets whether the ``threshold`` will be corrected
+
by the number of frequencies before being applied. This decreases
+
the ``threshold`` (p-values need to be lower to count as significant).
+
Default is ``False`` (report all powers) though for any application
+
where `threshold`` is set to something meaningful, this should also
+
be applied!
+
+
Returns
+
-------
+
pvals : iterable
+
A list of ``(index, p-value)`` tuples for all powers that have p-values
+
lower than the threshold specified in ``threshold``.
+
+
"""
+
if not self.norm == "leahy":
+
raise ValueError("This method only works on Leahy-normalized power spectra!")
+
+
if np.size(self.m) == 1:
+
# calculate p-values for all powers
+
# leave out zeroth power since it just encodes the number of photons!
+
pv = np.array([cospectra_pvalue(power, self.m) for power in self.power])
+
else:
+
pv = np.array([cospectra_pvalue(power, m) for power, m in zip(self.power, self.m)])
+
+
# if trial correction is used, then correct the threshold for
+
# the number of powers in the power spectrum
+
if trial_correction:
+
threshold /= self.power.shape[0]
+
+
# need to add 1 to the indices to make up for the fact that
+
# we left out the first power above!
+
indices = np.where(pv < threshold)[0]
+
+
pvals = np.vstack([pv[indices], indices])
+
+
return pvals
+
+
+
+
[docs]
+
@staticmethod
+
def from_time_array(
+
times1,
+
times2,
+
dt,
+
segment_size=None,
+
gti=None,
+
norm="none",
+
power_type="all",
+
silent=False,
+
fullspec=False,
+
use_common_mean=True,
+
):
+
"""Calculate AveragedCrossspectrum from two arrays of event times.
+
+
Parameters
+
----------
+
times1 : `np.array`
+
Event arrival times of channel 1
+
times2 : `np.array`
+
Event arrival times of channel 2
+
dt : float
+
The time resolution of the intermediate light curves
+
(sets the Nyquist frequency)
+
+
Other parameters
+
----------------
+
segment_size : float
+
The length, in seconds, of the light curve segments that will be
+
averaged. Only relevant (and required) for `AveragedCrossspectrum`.
+
gti : [[gti0, gti1], ...]
+
Good Time intervals. Defaults to the common GTIs from the two input
+
objects. Could throw errors if these GTIs have overlaps with the
+
input object GTIs! If you're getting errors regarding your GTIs,
+
don't use this and only give GTIs to the input objects before
+
making the cross spectrum.
+
norm : str, default "frac"
+
The normalization of the periodogram. "abs" is absolute rms, "frac" is
+
fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
+
unnormalized periodogram
+
use_common_mean : bool, default True
+
The mean of the light curve can be estimated in each interval, or on
+
the full light curve. This gives different results (Alston+2013).
+
Here we assume the mean is calculated on the full light curve, but
+
the user can set ``use_common_mean`` to False to calculate it on a
+
per-segment basis.
+
fullspec : bool, default False
+
Return the full periodogram, including negative frequencies
+
silent : bool, default False
+
Silence the progress bars
+
power_type : str, default 'all'
+
If 'all', give complex powers. If 'abs', the absolute value; if 'real',
+
the real part
+
"""
+
+
return crossspectrum_from_time_array(
+
times1,
+
times2,
+
dt,
+
segment_size=segment_size,
+
gti=gti,
+
norm=norm,
+
power_type=power_type,
+
silent=silent,
+
fullspec=fullspec,
+
use_common_mean=use_common_mean,
+
)
+
+
+
+
[docs]
+
@staticmethod
+
def from_events(
+
events1,
+
events2,
+
dt,
+
segment_size=None,
+
norm="none",
+
power_type="all",
+
silent=False,
+
fullspec=False,
+
use_common_mean=True,
+
gti=None,
+
):
+
"""Calculate AveragedCrossspectrum from two event lists
+
+
Parameters
+
----------
+
events1 : `stingray.EventList`
+
Events from channel 1
+
events2 : `stingray.EventList`
+
Events from channel 2
+
dt : float
+
The time resolution of the intermediate light curves
+
(sets the Nyquist frequency)
+
+
Other parameters
+
----------------
+
segment_size : float
+
The length, in seconds, of the light curve segments that will be averaged.
+
Only relevant (and required) for AveragedCrossspectrum
+
norm : str, default "frac"
+
The normalization of the periodogram. "abs" is absolute rms, "frac" is
+
fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
+
unnormalized periodogram
+
use_common_mean : bool, default True
+
The mean of the light curve can be estimated in each interval, or on
+
the full light curve. This gives different results (Alston+2013).
+
Here we assume the mean is calculated on the full light curve, but
+
the user can set ``use_common_mean`` to False to calculate it on a
+
per-segment basis.
+
fullspec : bool, default False
+
Return the full periodogram, including negative frequencies
+
silent : bool, default False
+
Silence the progress bars
+
power_type : str, default 'all'
+
If 'all', give complex powers. If 'abs', the absolute value; if 'real',
+
the real part
+
gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+
Good Time intervals. Defaults to the common GTIs from the two input
+
objects. Could throw errors if these GTIs have overlaps with the
+
input object GTIs! If you're getting errors regarding your GTIs,
+
don't use this and only give GTIs to the input objects before
+
making the cross spectrum.
+
"""
+
+
return crossspectrum_from_events(
+
events1,
+
events2,
+
dt,
+
segment_size=segment_size,
+
norm=norm,
+
power_type=power_type,
+
silent=silent,
+
fullspec=fullspec,
+
use_common_mean=use_common_mean,
+
gti=gti,
+
)
+
+
+
+
[docs]
+
@staticmethod
+
def from_lightcurve(
+
lc1,
+
lc2,
+
segment_size=None,
+
norm="none",
+
power_type="all",
+
silent=False,
+
fullspec=False,
+
use_common_mean=True,
+
gti=None,
+
):
+
"""Calculate AveragedCrossspectrum from two light curves
+
+
Parameters
+
----------
+
lc1 : `stingray.Lightcurve`
+
Light curve from channel 1
+
lc2 : `stingray.Lightcurve`
+
Light curve from channel 2
+
+
Other parameters
+
----------------
+
segment_size : float
+
The length, in seconds, of the light curve segments that will be averaged.
+
Only relevant (and required) for AveragedCrossspectrum
+
norm : str, default "frac"
+
The normalization of the periodogram. "abs" is absolute rms, "frac" is
+
fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
+
unnormalized periodogram
+
use_common_mean : bool, default True
+
The mean of the light curve can be estimated in each interval, or on
+
the full light curve. This gives different results (Alston+2013).
+
Here we assume the mean is calculated on the full light curve, but
+
the user can set ``use_common_mean`` to False to calculate it on a
+
per-segment basis.
+
fullspec : bool, default False
+
Return the full periodogram, including negative frequencies
+
silent : bool, default False
+
Silence the progress bars
+
power_type : str, default 'all'
+
If 'all', give complex powers. If 'abs', the absolute value; if 'real',
+
the real part
+
gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+
Good Time intervals. Defaults to the common GTIs from the two input
+
objects. Could throw errors if these GTIs have overlaps with the
+
input object GTIs! If you're getting errors regarding your GTIs,
+
don't use this and only give GTIs to the input objects before
+
making the cross spectrum.
+
"""
+
return crossspectrum_from_lightcurve(
+
lc1,
+
lc2,
+
segment_size=segment_size,
+
norm=norm,
+
power_type=power_type,
+
silent=silent,
+
fullspec=fullspec,
+
use_common_mean=use_common_mean,
+
gti=gti,
+
)
+
+
+
+
[docs]
+
@staticmethod
+
def from_stingray_timeseries(
+
ts1,
+
ts2,
+
flux_attr,
+
error_flux_attr=None,
+
segment_size=None,
+
norm="none",
+
power_type="all",
+
silent=False,
+
fullspec=False,
+
use_common_mean=True,
+
gti=None,
+
):
+
"""Calculate AveragedCrossspectrum from two light curves
+
+
Parameters
+
----------
+
ts1 : `stingray.Timeseries`
+
Time series from channel 1
+
ts2 : `stingray.Timeseries`
+
Time series from channel 2
+
flux_attr : `str`
+
What attribute of the time series will be used.
+
+
Other parameters
+
----------------
+
error_flux_attr : `str`
+
What attribute of the time series will be used as error bar.
+
segment_size : float
+
The length, in seconds, of the light curve segments that will be averaged.
+
Only relevant (and required) for AveragedCrossspectrum
+
norm : str, default "frac"
+
The normalization of the periodogram. "abs" is absolute rms, "frac" is
+
fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
+
unnormalized periodogram
+
use_common_mean : bool, default True
+
The mean of the light curve can be estimated in each interval, or on
+
the full light curve. This gives different results (Alston+2013).
+
Here we assume the mean is calculated on the full light curve, but
+
the user can set ``use_common_mean`` to False to calculate it on a
+
per-segment basis.
+
fullspec : bool, default False
+
Return the full periodogram, including negative frequencies
+
silent : bool, default False
+
Silence the progress bars
+
power_type : str, default 'all'
+
If 'all', give complex powers. If 'abs', the absolute value; if 'real',
+
the real part
+
gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+
Good Time intervals. Defaults to the common GTIs from the two input
+
objects. Could throw errors if these GTIs have overlaps with the
+
input object GTIs! If you're getting errors regarding your GTIs,
+
don't use this and only give GTIs to the input objects before
+
making the cross spectrum.
+
"""
+
return crossspectrum_from_timeseries(
+
ts1,
+
ts2,
+
flux_attr=flux_attr,
+
error_flux_attr=error_flux_attr,
+
segment_size=segment_size,
+
norm=norm,
+
power_type=power_type,
+
silent=silent,
+
fullspec=fullspec,
+
use_common_mean=use_common_mean,
+
gti=gti,
+
)
+
+
+
+
[docs]
+
@staticmethod
+
def from_lc_iterable(
+
iter_lc1,
+
iter_lc2,
+
dt,
+
segment_size,
+
norm="none",
+
power_type="all",
+
silent=False,
+
fullspec=False,
+
use_common_mean=True,
+
gti=None,
+
):
+
"""Calculate AveragedCrossspectrum from two light curves
+
+
Parameters
+
----------
+
iter_lc1 : iterable of `stingray.Lightcurve` objects or `np.array`
+
Light curves from channel 1. If arrays, use them as counts
+
iter_lc1 : iterable of `stingray.Lightcurve` objects or `np.array`
+
Light curves from channel 2. If arrays, use them as counts
+
dt : float
+
The time resolution of the light curves
+
(sets the Nyquist frequency)
+
+
Other parameters
+
----------------
+
segment_size : float
+
The length, in seconds, of the light curve segments that will be averaged.
+
Only relevant (and required) for AveragedCrossspectrum
+
norm : str, default "frac"
+
The normalization of the periodogram. "abs" is absolute rms, "frac" is
+
fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
+
unnormalized periodogram
+
use_common_mean : bool, default True
+
The mean of the light curve can be estimated in each interval, or on
+
the full light curve. This gives different results (Alston+2013).
+
Here we assume the mean is calculated on the full light curve, but
+
the user can set ``use_common_mean`` to False to calculate it on a
+
per-segment basis.
+
fullspec : bool, default False
+
Return the full periodogram, including negative frequencies
+
silent : bool, default False
+
Silence the progress bars
+
power_type : str, default 'all'
+
If 'all', give complex powers. If 'abs', the absolute value; if 'real',
+
the real part
+
gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+
Good Time intervals. Defaults to the common GTIs from the two input
+
objects. Could throw errors if these GTIs have overlaps with the
+
input object GTIs! If you're getting errors regarding your GTIs,
+
don't use this and only give GTIs to the input objects before
+
making the cross spectrum.
+
save_all : bool, default False
+
If True, save the cross spectrum of each segment in the ``cs_all``
+
attribute of the output :class:`Crossspectrum` object.
+
"""
+
+
return crossspectrum_from_lc_iterable(
+
iter_lc1,
+
iter_lc2,
+
dt,
+
segment_size,
+
norm=norm,
+
power_type=power_type,
+
silent=silent,
+
fullspec=fullspec,
+
use_common_mean=use_common_mean,
+
gti=gti,
+
)
+
+
+
def _initialize_from_any_input(
+
self,
+
data1,
+
data2,
+
dt=None,
+
segment_size=None,
+
norm="frac",
+
power_type="all",
+
silent=False,
+
fullspec=False,
+
gti=None,
+
use_common_mean=True,
+
save_all=False,
+
):
+
"""Initialize the class, trying to understand the input types.
+
+
The input arguments are the same as ``__init__()``. Based on the type
+
of ``data1``, this method will call the appropriate
+
``crossspectrum_from_XXXX`` function, and initialize ``self`` with
+
the correct attributes.
+
"""
+
if segment_size is None and isinstance(data1, StingrayObject):
+
common_gti = cross_two_gtis(data1.gti, data2.gti)
+
data1.gti = common_gti
+
data2.gti = common_gti
+
data1 = data1.apply_gtis(inplace=False)
+
data2 = data2.apply_gtis(inplace=False)
+
+
data1_is_binned = (
+
"counts" in data1.array_attrs() or "_counts" in data1.internal_array_attrs()
+
)
+
data2_is_binned = (
+
"counts" in data2.array_attrs() or "_counts" in data2.internal_array_attrs()
+
)
+
+
if data1_is_binned and data2_is_binned:
+
assert data2.time.size == data1.time.size and np.allclose(
+
data2.time - data1.time, 0
+
), "Time arrays are not the same"
+
elif data1_is_binned or data2_is_binned:
+
raise ValueError("Please use input data of the same kind")
+
+
if isinstance(data1, EventList):
+
spec = crossspectrum_from_events(
+
data1,
+
data2,
+
dt,
+
segment_size,
+
norm=norm,
+
power_type=power_type,
+
silent=silent,
+
fullspec=fullspec,
+
use_common_mean=use_common_mean,
+
gti=gti,
+
save_all=save_all,
+
)
+
elif isinstance(data1, Lightcurve):
+
spec = crossspectrum_from_lightcurve(
+
data1,
+
data2,
+
segment_size,
+
norm=norm,
+
power_type=power_type,
+
silent=silent,
+
fullspec=fullspec,
+
use_common_mean=use_common_mean,
+
gti=gti,
+
save_all=save_all,
+
)
+
spec.lc1 = data1
+
spec.lc2 = data2
+
elif isinstance(data1, (tuple, list)):
+
dt = data1[0].dt
+
# This is a list of light curves.
+
spec = crossspectrum_from_lc_iterable(
+
data1,
+
data2,
+
dt,
+
segment_size,
+
norm=norm,
+
power_type=power_type,
+
silent=silent,
+
fullspec=fullspec,
+
gti=gti,
+
use_common_mean=use_common_mean,
+
save_all=save_all,
+
)
+
else: # pragma: no cover
+
raise TypeError(f"Bad inputs to Crosssspectrum: {type(data1)}")
+
+
for key, val in spec.__dict__.items():
+
setattr(self, key, val)
+
return
+
+
def _initialize_empty(self):
+
"""Set all attributes to None."""
+
self.freq = None
+
self.power = None
+
self.power_err = None
+
self.unnorm_power = None
+
self.unnorm_power_err = None
+
self.df = None
+
self.dt = None
+
self.nphots1 = None
+
self.nphots2 = None
+
self.m = 1
+
self.n = None
+
self.fullspec = None
+
self.k = 1
+
return
+
+
+
[docs]
+
def deadtime_correct(
+
self, dead_time, rate, background_rate=0, limit_k=200, n_approx=None, paralyzable=False
+
):
+
"""
+
Correct the power spectrum for dead time effects.
+
+
This correction is based on the formula given in Zhang et al. 2015, assuming
+
a constant dead time for all events.
+
For more advanced dead time corrections, see the FAD method from `stingray.deadtime.fad`
+
+
Parameters
+
----------
+
dead_time: float
+
The dead time of the detector.
+
rate : float
+
Detected source count rate
+
+
Other Parameters
+
----------------
+
background_rate : float, default 0
+
Detected background count rate. This is important to estimate when deadtime is given by the
+
combination of the source counts and background counts (e.g. in an imaging X-ray detector).
+
paralyzable: bool, default False
+
If True, the dead time correction is done assuming a paralyzable
+
dead time. If False, the correction is done assuming a non-paralyzable
+
(more common) dead time.
+
limit_k : int, default 200
+
Limit to this value the number of terms in the inner loops of
+
calculations. Check the plots returned by the `check_B` and
+
`check_A` functions to test that this number is adequate.
+
n_approx : int, default None
+
Number of bins to calculate the model power spectrum. If None, it will use the size of
+
the input frequency array. Relatively simple models (e.g., low count rates compared to
+
dead time) can use a smaller number of bins to speed up the calculation, and the final
+
power values will be interpolated.
+
+
Returns
+
-------
+
spectrum: :class:`Crossspectrum` or derivative.
+
The dead-time corrected spectrum.
+
"""
+
# I put it here to avoid circular imports
+
from stingray.deadtime import non_paralyzable_dead_time_model
+
+
if paralyzable:
+
raise NotImplementedError("Paralyzable dead time correction is not implemented yet.")
+
+
model = non_paralyzable_dead_time_model(
+
self.freq,
+
dead_time,
+
rate,
+
bin_time=self.dt,
+
limit_k=limit_k,
+
background_rate=background_rate,
+
n_approx=n_approx,
+
)
+
correction = 2 / model
+
new_spec = copy.deepcopy(self)
+
new_spec.power *= correction
+
+
# Now correct internal attributes for both the spectrum and its possible
+
# sub-spectra (PDSs from single channels, dynamical spectra, etc.)
+
objects = [new_spec]
+
if hasattr(new_spec, "pds1"):
+
objects += [new_spec.pds1, new_spec.pds2]
+
+
for obj in objects:
+
for attr in ["unnorm_power", "power_err", "unnorm_power_err"]:
+
if hasattr(obj, attr):
+
setattr(obj, attr, getattr(obj, attr) * correction)
+
+
for attr in ["cs_all", "unnorm_cs_all"]:
+
if hasattr(new_spec, attr):
+
dynsp = getattr(new_spec, attr)
+
for i, power in enumerate(dynsp):
+
dynsp[i] = power * correction
+
+
return new_spec
+
+
[docs]
+
class AveragedCrossspectrum(Crossspectrum):
+
type = "crossspectrum"
+
"""
+
Make an averaged cross spectrum from a light curve by segmenting two
+
light curves, Fourier-transforming each segment and then averaging the
+
resulting cross spectra.
+
+
Parameters
+
----------
+
data1: :class:`stingray.Lightcurve`object OR iterable of :class:`stingray.Lightcurve` objects OR :class:`stingray.EventList` object
+
A light curve from which to compute the cross spectrum. In some cases,
+
this would be the light curve of the wavelength/energy/frequency band
+
of interest.
+
+
data2: :class:`stingray.Lightcurve`object OR iterable of :class:`stingray.Lightcurve` objects OR :class:`stingray.EventList` object
+
A second light curve to use in the cross spectrum. In some cases, this
+
would be the wavelength/energy/frequency reference band to compare the
+
band of interest with.
+
+
segment_size: float
+
The size of each segment to average. Note that if the total duration of
+
each :class:`Lightcurve` object in ``lc1`` or ``lc2`` is not an
+
integer multiple of the ``segment_size``, then any fraction left-over
+
at the end of the time series will be lost. Otherwise you introduce
+
artifacts.
+
+
norm: {``frac``, ``abs``, ``leahy``, ``none``}, default ``none``
+
The normalization of the (real part of the) cross spectrum.
+
+
Other Parameters
+
----------------
+
gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+
Good Time intervals. Defaults to the common GTIs from the two input
+
objects. Could throw errors if these GTIs have overlaps with the
+
input object GTIs! If you're getting errors regarding your GTIs,
+
don't use this and only give GTIs to the input objects before
+
making the cross spectrum.
+
+
dt : float
+
The time resolution of the light curve. Only needed when constructing
+
light curves in the case where data1 or data2 are of :class:EventList
+
+
power_type: string, optional, default ``all``
+
Parameter to choose among complete, real part and magnitude of
+
the cross spectrum.
+
+
silent : bool, default False
+
Do not show a progress bar when generating an averaged cross spectrum.
+
Useful for the batch execution of many spectra
+
+
lc1: :class:`stingray.Lightcurve`object OR iterable of :class:`stingray.Lightcurve` objects
+
For backwards compatibility only. Like ``data1``, but no
+
:class:`stingray.events.EventList` objects allowed
+
+
lc2: :class:`stingray.Lightcurve`object OR iterable of :class:`stingray.Lightcurve` objects
+
For backwards compatibility only. Like ``data2``, but no
+
:class:`stingray.events.EventList` objects allowed
+
+
fullspec: boolean, optional, default ``False``
+
If True, return the full array of frequencies, otherwise return just the
+
positive frequencies.
+
+
save_all : bool, default False
+
Save all intermediate PDSs used for the final average. Use with care.
+
This is likely to fill up your RAM on medium-sized datasets, and to
+
slow down the computation when rebinning.
+
+
skip_checks: bool
+
Skip initial checks, for speed or other reasons (you need to trust your
+
inputs!)
+
+
use_common_mean: bool
+
Averaged cross spectra are normalized in two possible ways: one is by normalizing
+
each of the single spectra that get averaged, the other is by normalizing after the
+
averaging. If `use_common_mean` is selected, the spectrum will be normalized
+
after the average.
+
+
gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+
Good Time intervals. Defaults to the common GTIs from the two input
+
objects. Could throw errors if these GTIs have overlaps with the
+
input object GTIs! If you're getting errors regarding your GTIs,
+
don't use this and only give GTIs to the input objects before
+
making the cross spectrum.
+
+
Attributes
+
----------
+
freq: numpy.ndarray
+
The array of mid-bin frequencies that the Fourier transform samples.
+
+
power: numpy.ndarray
+
The array of cross spectra.
+
+
power_err: numpy.ndarray
+
The uncertainties of ``power``.
+
An approximation for each bin given by ``power_err= power/sqrt(m)``.
+
Where ``m`` is the number of power averaged in each bin (by frequency
+
binning, or averaging power spectra of segments of a light curve).
+
Note that for a single realization (``m=1``) the error is equal to the
+
power.
+
+
df: float
+
The frequency resolution.
+
+
m: int
+
The number of averaged cross spectra.
+
+
n: int
+
The number of time bins per segment of light curve.
+
+
nphots1: float
+
The total number of photons in the first (interest) light curve.
+
+
nphots2: float
+
The total number of photons in the second (reference) light curve.
+
+
gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+
Good Time intervals.
+
"""
+
+
def __init__(
+
self,
+
data1=None,
+
data2=None,
+
segment_size=None,
+
norm="frac",
+
gti=None,
+
power_type="all",
+
silent=False,
+
lc1=None,
+
lc2=None,
+
dt=None,
+
fullspec=False,
+
save_all=False,
+
use_common_mean=True,
+
skip_checks=False,
+
):
+
self._type = None
+
# for backwards compatibility
+
if data1 is None:
+
data1 = lc1
+
if data2 is None:
+
data2 = lc2
+
+
empty = data1 is None and data2 is None
+
good_input = not empty
+
+
if not skip_checks:
+
good_input = self.initial_checks(
+
data1=data1,
+
data2=data2,
+
norm=norm,
+
gti=gti,
+
lc1=lc1,
+
lc2=lc2,
+
power_type=power_type,
+
dt=dt,
+
fullspec=fullspec,
+
segment_size=segment_size,
+
)
+
norm = norm.lower()
+
self.norm = norm
+
self.dt = dt
+
self.save_all = save_all
+
self.segment_size = segment_size
+
self.show_progress = not silent
+
+
if empty or not good_input:
+
return self._initialize_empty()
+
+
if isinstance(data1, Generator):
+
warnings.warn(
+
"The averaged Cross spectrum from a generator of "
+
"light curves pre-allocates the full list of light "
+
"curves, losing all advantage of lazy loading. If it "
+
"is important for you, use the "
+
"AveragedCrossspectrum.from_lc_iterable static "
+
"method, specifying the sampling time `dt`."
+
)
+
data1 = list(data1)
+
data2 = list(data2)
+
+
return self._initialize_from_any_input(
+
data1,
+
data2,
+
dt=dt,
+
segment_size=segment_size,
+
gti=gti,
+
norm=norm,
+
power_type=power_type,
+
silent=silent,
+
fullspec=fullspec,
+
use_common_mean=use_common_mean,
+
save_all=save_all,
+
)
+
+
+
[docs]
+
def initial_checks(self, data1, segment_size=None, **kwargs):
+
"""
+
+
Examples
+
--------
+
>>> times = np.arange(0, 10)
+
>>> ev1 = EventList(times)
+
>>> ev2 = EventList(times)
+
>>> ac = AveragedCrossspectrum()
+
+
If AveragedCrossspectrum, you need ``segment_size``
+
>>> AveragedCrossspectrum.initial_checks(ac, data1=ev1, data2=ev2, dt=1)
+
Traceback (most recent call last):
+
...
+
ValueError: segment_size must be specified...
+
+
And it needs to be finite!
+
>>> AveragedCrossspectrum.initial_checks(ac, data1=ev1, data2=ev2, dt=1., segment_size=np.nan)
+
Traceback (most recent call last):
+
...
+
ValueError: segment_size must be finite!
+
"""
+
good = Crossspectrum.initial_checks(self, data1, segment_size=segment_size, **kwargs)
+
if not good:
+
return False
+
if isinstance(self, AveragedCrossspectrum) and segment_size is None and data1 is not None:
+
raise ValueError("segment_size must be specified")
+
+
if (
+
isinstance(self, AveragedCrossspectrum)
+
and segment_size is not None
+
and not np.isfinite(segment_size)
+
):
+
raise ValueError("segment_size must be finite!")
+
return True
+
+
+
+
[docs]
+
def coherence(self):
+
"""Averaged Coherence function.
+
+
+
Coherence is defined in Vaughan and Nowak, 1996 [#]_.
+
It is a Fourier frequency dependent measure of the linear correlation
+
between time series measured simultaneously in two energy channels.
+
+
Compute an averaged Coherence function of cross spectrum by computing
+
coherence function of each segment and averaging them. The return type
+
is a tuple with first element as the coherence function and the second
+
element as the corresponding uncertainty associated with it.
+
+
Note : The uncertainty in coherence function is strictly valid for Gaussian \
+
statistics only.
+
+
Returns
+
-------
+
(coh, uncertainty) : tuple of np.ndarray
+
Tuple comprising the coherence function and uncertainty.
+
+
References
+
----------
+
.. [#] https://iopscience.iop.org/article/10.1086/310430
+
"""
+
if np.any(self.m < 50):
+
simon(
+
"Number of segments used in averaging is "
+
"significantly low. The result might not follow the "
+
"expected statistical distributions."
+
)
+
c = self.unnorm_power
+
p1 = self.pds1.unnorm_power
+
p2 = self.pds2.unnorm_power
+
+
meanrate1 = self.nphots1 / self.n / self.dt
+
meanrate2 = self.nphots2 / self.n / self.dt
+
+
P1noise = poisson_level(norm="none", meanrate=meanrate1, n_ph=self.nphots1)
+
P2noise = poisson_level(norm="none", meanrate=meanrate2, n_ph=self.nphots2)
+
+
coh = raw_coherence(c, p1, p2, P1noise, P2noise, self.n)
+
+
# Calculate uncertainty
+
uncertainty = (2**0.5 * coh * (1 - coh)) / (np.sqrt(coh) * self.m**0.5)
+
+
uncertainty[coh == 0] = 0.0
+
+
return (coh, uncertainty)
+
+
+
+
[docs]
+
def phase_lag(self):
+
"""Return the fourier phase lag of the cross spectrum."""
+
lag = np.angle(self.unnorm_power)
+
coh, uncert = self.coherence()
+
+
dum = (1.0 - coh) / (2.0 * coh)
+
+
dum[coh == 0] = 0.0
+
+
lag_err = np.sqrt(dum / self.m)
+
return lag, lag_err
+
+
+
+
[docs]
+
def time_lag(self):
+
"""Calculate time lag and uncertainty.
+
+
Equation from Bendat & Piersol, 2011 [bendat-2011]__.
+
+
Returns
+
-------
+
lag : np.ndarray
+
The time lag
+
+
lag_err : np.ndarray
+
The uncertainty in the time lag
+
"""
+
ph_lag, ph_lag_err = self.phase_lag()
+
+
lag = ph_lag / (2 * np.pi * self.freq)
+
lag_err = ph_lag_err / (2 * np.pi * self.freq)
+
+
return lag, lag_err
+