Merge pull request #780 from StingraySoftware/power_colors

Introduce power colors

docs/changes/780.feature.rst

@@ -0,0 +1 @@
+Introduce power colors à la Heil et al. 2015

docs/index.rst

@@ -24,15 +24,16 @@ Current Capabilities
 1. Data handling and simulation
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-* loading event lists from fits files of a few missions (RXTE/PCA, NuSTAR/FPM, XMM-Newton/EPIC, NICER/XTI)
-* constructing light curves from event data, various operations on light curves (e.g. addition, subtraction, joining, and truncation)
+* loading event lists from fits files (and generally good handling of OGIP-compliant missions, like RXTE/PCA, NuSTAR/FPM, XMM-Newton/EPIC, NICER/XTI)
+* constructing light curves and time series from event data
+* various operations on time series (e.g. addition, subtraction, joining, and truncation)
 * simulating a light curve with a given power spectrum
 * simulating a light curve from another light curve and a 1-d (time) or 2-d (time-energy) impulse response
 * simulating an event list from a given light curve _and_ with a given energy spectrum
 * Good Time Interval operations
 * Filling gaps in light curves with statistically sound fake data
 
-2. Fourier methods
+1. Fourier methods
 ~~~~~~~~~~~~~~~~~~
 * power spectra and cross spectra in Leahy, rms normalization, absolute rms and no normalization
 * averaged power spectra and cross spectra
@@ -45,6 +46,7 @@ Current Capabilities
 * bispectra; *needs testing*
 * (Bayesian) quasi-periodic oscillation searches
 * Lomb-Scargle periodograms and cross spectra
+* Power Colors
 
 3. Other time series methods
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -62,7 +64,6 @@ Other future additions we are currently implementing are:
 * bicoherence
 * phase-resolved spectroscopy of quasi-periodic oscillations
 * Fourier-frequency-resolved spectroscopy
-* power colours
 * full HEASARC-compatible mission support
 * pulsar searches with :math:`H`-test
 * binary pulsar searches

stingray/crossspectrum.py

@@ -21,7 +21,7 @@
 from .fourier import avg_cs_from_iterables, error_on_averaged_cross_spectrum
 from .fourier import avg_cs_from_events, poisson_level
 from .fourier import fftfreq, fft, normalize_periodograms, raw_coherence
-from .fourier import get_flux_iterable_from_segments
+from .fourier import get_flux_iterable_from_segments, power_color
 
 from scipy.special import factorial
 
@@ -2028,6 +2028,9 @@ class DynamicalCrossspectrum(AveragedCrossspectrum):
         The time resolution of the dynamical spectrum. It is **not** the time resolution of the
         input light curve. It is the integration time of each line of the dynamical power
         spectrum (typically, an integer multiple of ``segment_size``).
+
+    m: int
+        The number of averaged powers in each spectral bin (initially 1, it changes after rebinning).
     """
 
     def __init__(self, data1, data2, segment_size, norm="frac", gti=None, sample_time=None):
@@ -2075,16 +2078,21 @@ def _make_matrix(self, data1, data2):
             save_all=True,
         )
         self.dyn_ps = np.array(avg.cs_all).T
+        conv = avg.cs_all / avg.unnorm_cs_all
+        self.unnorm_conversion = np.nanmean(conv)
         self.freq = avg.freq
         current_gti = avg.gti
+        self.nphots1 = avg.nphots1
+        self.nphots2 = avg.nphots2
 
         tstart, _ = time_intervals_from_gtis(current_gti, self.segment_size)
 
         self.time = tstart + 0.5 * self.segment_size
         self.df = avg.df
         self.dt = self.segment_size
+        self.m = 1
 
-    def rebin_frequency(self, df_new, method="sum"):
+    def rebin_frequency(self, df_new, method="average"):
         """
         Rebin the Dynamic Power Spectrum to a new frequency resolution.
         Rebinning is an in-place operation, i.e. will replace the existing
@@ -2101,24 +2109,26 @@ def rebin_frequency(self, df_new, method="sum"):
             Must be larger than the frequency resolution of the old dynamical
             power spectrum!
 
-        method: {"sum" | "mean" | "average"}, optional, default "sum"
-            This keyword argument sets whether the counts in the new bins
+        method: {"sum" | "mean" | "average"}, optional, default "average"
+            This keyword argument sets whether the powers in the new bins
             should be summed or averaged.
         """
         new_dynspec_object = copy.deepcopy(self)
         dynspec_new = []
         for data in self.dyn_ps.T:
-            freq_new, bin_counts, bin_err, _ = rebin_data(
+            freq_new, bin_counts, bin_err, step = rebin_data(
                 self.freq, data, dx_new=df_new, method=method
             )
             dynspec_new.append(bin_counts)
 
         new_dynspec_object.freq = freq_new
         new_dynspec_object.dyn_ps = np.array(dynspec_new).T
         new_dynspec_object.df = df_new
+        new_dynspec_object.m = int(step) * self.m
+
         return new_dynspec_object
 
-    def rebin_time(self, dt_new, method="sum"):
+    def rebin_time(self, dt_new, method="average"):
         """
         Rebin the Dynamic Power Spectrum to a new time resolution.
 
@@ -2137,8 +2147,8 @@ def rebin_time(self, dt_new, method="sum"):
             Must be larger than the time resolution of the old dynamical power
             spectrum!
 
-        method: {"sum" | "mean" | "average"}, optional, default "sum"
-            This keyword argument sets whether the counts in the new bins
+        method: {"sum" | "mean" | "average"}, optional, default "average"
+            This keyword argument sets whether the powers in the new bins
             should be summed or averaged.
 
         Returns
@@ -2156,14 +2166,77 @@ def rebin_time(self, dt_new, method="sum"):
 
         dynspec_new = []
         for data in self.dyn_ps:
-            time_new, bin_counts, _, _ = rebin_data(
+            time_new, bin_counts, _, step = rebin_data(
                 self.time, data, dt_new, method=method, dx=self.dt
             )
             dynspec_new.append(bin_counts)
 
         new_dynspec_object.time = time_new
         new_dynspec_object.dyn_ps = np.array(dynspec_new)
         new_dynspec_object.dt = dt_new
+        new_dynspec_object.m = int(step) * self.m
+        return new_dynspec_object
+
+    def rebin_by_n_intervals(self, n, method="average"):
+        """
+        Rebin the Dynamic Power Spectrum to a new time resolution, by summing contiguous intervals.
+
+        This is different from meth:`DynamicalPowerspectrum.rebin_time` in that it averages ``n``
+        consecutive intervals regardless of their distance in time. ``rebin_time`` will instead
+        average intervals that are separated at most by a time ``dt_new``.
+
+        Parameters
+        ----------
+        n: int
+            The number of intervals to be combined into one.
+
+        method: {"sum" | "mean" | "average"}, optional, default "average"
+            This keyword argument sets whether the powers in the new bins
+            should be summed or averaged.
+
+        Returns
+        -------
+        time_new: numpy.ndarray
+            Time axis with new rebinned time resolution.
+
+        dynspec_new: numpy.ndarray
+            New rebinned Dynamical Cross Spectrum.
+        """
+        if not np.issubdtype(type(n), np.integer):
+            warnings.warn("n must be an integer. Casting to int")
+
+            n = int(n)
+
+        if n == 1:
+            return copy.deepcopy(self)
+        if n < 1:
+            raise ValueError("n must be >= 1")
+
+        new_dynspec_object = copy.deepcopy(self)
+
+        dynspec_new = []
+        time_new = []
+        for i, data in enumerate(self.dyn_ps.T):
+            if i % n == 0:
+                count = 1
+                bin_counts = data
+                bin_times = self.time[i]
+            else:
+                count += 1
+                bin_counts += data
+                bin_times += self.time[i]
+            if count == n:
+                if method in ["mean", "average"]:
+                    bin_counts /= n
+                dynspec_new.append(bin_counts)
+                bin_times /= n
+                time_new.append(bin_times)
+
+        new_dynspec_object.time = time_new
+        new_dynspec_object.dyn_ps = np.array(dynspec_new).T
+        new_dynspec_object.dt *= n
+        new_dynspec_object.m = n * self.m
+
         return new_dynspec_object
 
     def trace_maximum(self, min_freq=None, max_freq=None):
@@ -2198,6 +2271,125 @@ def trace_maximum(self, min_freq=None, max_freq=None):
 
         return np.array(max_positions)
 
+    def power_colors(
+        self,
+        freq_edges=[1 / 256, 1 / 32, 0.25, 2.0, 16.0],
+        freqs_to_exclude=None,
+        poisson_power=0,
+    ):
+        """
+        Compute the power colors of the dynamical power spectrum.
+
+        See Heil et al. 2015, MNRAS, 448, 3348
+
+        Parameters
+        ----------
+        freq_edges: iterable
+            The edges of the frequency bins to be used for the power colors.
+
+        freqs_to_exclude : 1-d or 2-d iterable, optional, default None
+            The ranges of frequencies to exclude from the calculation of the power color.
+            For example, the frequencies containing strong QPOs.
+            A 1-d iterable should contain two values for the edges of a single range. (E.g.
+            ``[0.1, 0.2]``). ``[[0.1, 0.2], [3, 4]]`` will exclude the ranges 0.1-0.2 Hz and
+            3-4 Hz.
+
+        poisson_power : float or iterable, optional, default 0
+            The Poisson noise level of the power spectrum. If iterable, it should have the same
+            length as ``freq``. (This might apply to the case of a power spectrum with a
+            strong dead time distortion)
+
+        Returns
+        -------
+        pc0: np.ndarray
+        pc0_err: np.ndarray
+        pc1: np.ndarray
+        pc1_err: np.ndarray
+            The power colors for each spectrum and their respective errors
+        """
+        power_colors = []
+        if np.iscomplexobj(self.dyn_ps):
+            warnings.warn("When using power_colors, complex powers will be cast to real.")
+
+        for ps in self.dyn_ps.T:
+            power_colors.append(
+                power_color(
+                    self.freq,
+                    ps.real,
+                    freq_edges=freq_edges,
+                    freqs_to_exclude=freqs_to_exclude,
+                    df=self.df,
+                    poisson_power=poisson_power,
+                    m=self.m,
+                )
+            )
+
+        pc0, pc0_err, pc1, pc1_err = np.array(power_colors).T
+        return pc0, pc0_err, pc1, pc1_err
+
+    def compute_rms(self, min_freq, max_freq, poisson_noise_level=0):
+        """
+        Compute the fractional rms amplitude in the power spectrum
+        between two frequencies.
+
+        Parameters
+        ----------
+        min_freq: float
+            The lower frequency bound for the calculation.
+
+        max_freq: float
+            The upper frequency bound for the calculation.
+
+        Other parameters
+        ----------------
+        poisson_noise_level : float, default is None
+            This is the Poisson noise level of the PDS with same
+            normalization as the PDS.
+
+        Returns
+        -------
+        rms: float
+            The fractional rms amplitude contained between ``min_freq`` and
+            ``max_freq``.
+
+        rms_err: float
+            The error on the fractional rms amplitude.
+
+        """
+        from .fourier import rms_calculation
+
+        dyn_ps_unnorm = self.dyn_ps / self.unnorm_conversion
+        poisson_noise_unnrom = poisson_noise_level / self.unnorm_conversion
+        if not hasattr(self, "nphots"):
+            nphots = (self.nphots1 * self.nphots2) ** 0.5
+        else:
+            nphots = self.nphots
+        T = self.dt * self.m
+        M_freqs = self.m
+        K_freqs = 1
+        minind = self.freq.searchsorted(min_freq)
+        maxind = self.freq.searchsorted(max_freq)
+        min_freq = self.freq[minind]
+        freq_bins = maxind - minind
+        rmss = []
+        rms_errs = []
+        for unnorm_powers in dyn_ps_unnorm.T:
+            r, re = rms_calculation(
+                unnorm_powers[minind:maxind],
+                min_freq,
+                max_freq,
+                nphots,
+                T,
+                M_freqs,
+                K_freqs,
+                freq_bins,
+                poisson_noise_unnrom,
+                deadtime=0.0,
+            )
+            rmss.append(r)
+            rms_errs.append(re)
+        return rmss, rms_errs
+
 
 def crossspectrum_from_time_array(
     times1,