Merge pull request #693 from mgullik/rebinning

Correct rms calculation of rebinned PSD

stingray/crossspectrum.py

@@ -551,6 +551,9 @@ class Crossspectrum(StingrayObject):
         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
 
@@ -601,6 +604,7 @@ def __init__(
         self.dt = dt
         norm = norm.lower()
         self.norm = norm
+        self.k = 1
 
         if not good_input:
             return self._initialize_empty()
@@ -1164,6 +1168,7 @@ def rebin_log(self, f=0.01):
         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
@@ -1704,6 +1709,7 @@ def _initialize_empty(self):
         self.m = 1
         self.n = None
         self.fullspec = None
+        self.k = 1
         return
 
 class AveragedCrossspectrum(Crossspectrum):

stingray/fourier.py

@@ -416,7 +416,7 @@ def normalize_periodograms(unnorm_power, dt, n_bin, mean_flux=None, n_ph=None,
     """
     Wrapper around all the normalize_NORM methods.
 
-    Normalize the real part of the cross spectrum to Leahy, absolute rms^2,
+    Normalize the cross-spectrum or the power-spectrum to Leahy, absolute rms^2,
     fractional rms^2 normalization, or not at all.
 
     Parameters
@@ -490,6 +490,88 @@ def normalize_periodograms(unnorm_power, dt, n_bin, mean_flux=None, n_ph=None,
     raise ValueError("Unrecognized power type")
 
 
+def unnormalize_periodograms(norm_power, dt, n_bin, n_ph,
+                           variance=None, background_flux=0., norm=None,
+                           power_type="all"):
+    """
+    Wrapper around all the normalize_NORM methods.
+
+    Unnormalize the power of the cross-spectrum to Leahy, absolute rms^2,
+    fractional rms^2 normalization, or not at all.
+
+    Parameters
+    ----------
+    norm_power: numpy.ndarray
+        The normalized cross-spectrum or poisson noise
+
+    dt: float
+        The sampling time of the light curve
+
+    n_bin: int
+        The number of bins in the light curve
+
+    Other parameters
+    ----------------
+    mean_flux: float
+        The mean of the light curve used to calculate the powers
+        (If a cross spectrum, the geometrical mean of the light
+        curves in the two channels). Only relevant for "frac" normalization
+
+    n_ph: int or float
+        The number of counts in the light curve used to calculate
+        the unnormalized periodogram. Only relevant for Leahy normalization.
+
+    variance: float
+        The average variance of the measurements in light curve (if a cross
+        spectrum,  the geometrical mean of the variances in the two channels).
+        **NOT** the variance of the light curve, but of each flux measurement
+        (square of light curve error bar)! Only relevant for the Leahy
+        normalization of non-Poissonian data.
+
+    norm : str
+        One of ``leahy`` (Leahy+83), ``frac`` (fractional rms), ``abs``
+        (absolute rms),
+
+    power_type : str
+        One of ``real`` (real part), ``all`` (all complex powers), ``abs``
+        (absolute value)
+
+    background_flux : float, default 0
+        The background flux, in the same units as `mean_flux`.
+
+    Returns
+    -------
+    power: numpy.nd.array
+        The normalized co-spectrum (real part of the cross spectrum). For
+        'none' normalization, imaginary part is returned as well.
+    """
+
+    if norm == "leahy" and variance is not None:
+        unnorm_power = norm_power * (variance * n_ph) / 2.
+    elif norm == "leahy":
+        unnorm_power = norm_power * n_ph / 2.
+    elif norm == "frac":
+        if background_flux > 0:
+            unnorm_power = norm_power * ((n_ph/n_bin - background_flux) ** 2 *
+                                          n_bin) / (2. * dt)
+        else:
+            unnorm_power = norm_power * (n_ph**2/n_bin) / (2. * dt)  
+    elif norm == "abs":
+        unnorm_power = norm_power * dt * n_bin / 2.
+    elif norm == "none":
+        unnorm_power = norm_power
+    else:
+        raise ValueError("Unknown value for the norm")
+
+    if power_type == "all":
+        return unnorm_power
+    if power_type == "real":
+        return unnorm_power.real
+    if power_type in ["abs", "absolute"]:
+        return np.abs(unnorm_power)
+    raise ValueError("Unrecognized power type")
+
+
 def bias_term(power1, power2, power1_noise, power2_noise, n_ave,
               intrinsic_coherence=1.0):
     """
@@ -662,6 +744,77 @@ def estimate_intrinsic_coherence(cross_power, power1, power2, power1_noise,
     return new_coherence
 
 
+def rms_calculation(unnorm_powers, min_freq, max_freq, nphots, 
+    T, M_freqs, K_freqs, freq_bins, poisson_noise_unnrom, deadtime=0.):
+    """
+    Compute the fractional rms amplitude in the given power or cross spectrum
+    
+    NOTE: all array quantities are already in the correct energy range
+
+    Parameters
+    ----------
+    unnrom_powers: array of float 
+        unnormalised power or cross spectrum, the array has already been 
+        filtered for the given frequency range
+
+    min_freq: float
+        The lower frequency bound for the calculation (from the freq grid).
+
+    max_freq: float
+        The upper frequency bound for the calculation (from the freq grid).
+    
+    nphots: float 
+        Number of photons for the full power or cross spectrum
+    
+    T: float 
+        Time length of the light curve 
+    
+    M_freq: scalar or array of float 
+        If scalar, it is the number of segments in the AveragedCrossspectrum
+        If array, it is the number of segments times the rebinning sample 
+        in the given frequency range.
+
+    K_freq: scalar or array of float 
+        If scalar, the power or cross spectrum is not rebinned (K_freq = 1) 
+        If array,  the power or cross spectrum is rebinned and it is the 
+        rebinned sample in the given frequency range.
+
+    freq_bins: integer 
+        if the cross or power spectrum is rebinned freq_bins = 1, 
+        if it NOT rebinned freq_bins is the number of frequency bins 
+        in the given frequency range.
+
+    poisson_noise_unnrom : float
+        This is the Poisson noise level unnormalised.
+
+    Other parameters
+    ----------------
+    deadtime: float 
+        Deadtime of the instrument 
+
+    Returns
+    -------
+    rms: float
+        The fractional rms amplitude contained between ``min_freq`` and
+        ``max_freq``.
+
+    rms_err: float
+        The error on the fractional rms amplitude.
+
+    """
+    rms_squared = np.sum((unnorm_powers - poisson_noise_unnrom) * 1/T * K_freqs) \
+        * 2 * T / nphots**2
+    rms = np.sqrt(rms_squared)
+
+    rms_noise_squared = poisson_noise_unnrom * (max_freq - min_freq) \
+        * 2 * T / nphots**2 #rms of the noise    
+    rms_err_squared = (2 * rms_squared * rms_noise_squared + rms_noise_squared**2) / \
+        (2 * np.sum(M_freqs) * freq_bins * rms_squared)
+    rms_err = np.sqrt(rms_err_squared)
+
+    return rms, rms_err
+
+
 def error_on_averaged_cross_spectrum(cross_power, seg_power, ref_power, n_ave,
                                      seg_power_noise, ref_power_noise,
                                      common_ref=False):

stingray/multitaper.py

@@ -102,6 +102,9 @@ class Multitaper(Powerspectrum):
 
     n: int
         The number of data points in the light curve
+    
+    k: array of int
+        The rebinning scheme if the object has been rebinned otherwise is set to 1. 
 
     nphots: float
         The total number of photons in the light curve
@@ -158,6 +161,7 @@ def __init__(self, data=None, norm="frac", gti=None, dt=None, lc=None,
             self.m = 1
             self.n = None
             self.nphots = None
+            self.k = 1
             self.jk_var_deg_freedom = None
             return
         elif not isinstance(data, EventList):
@@ -169,6 +173,7 @@ def __init__(self, data=None, norm="frac", gti=None, dt=None, lc=None,
         self.lc = lc
         self.power_type = 'real'
         self.fullspec = False
+        self.k = 1
 
         self._make_multitaper_periodogram(lc, NW=NW, adaptive=adaptive,
                                           jackknife=jackknife, low_bias=low_bias,

stingray/powerspectrum.py

@@ -17,10 +17,11 @@
 from .events import EventList
 from .gti import cross_two_gtis
 from .lightcurve import Lightcurve
-from .fourier import avg_pds_from_iterable
+from .fourier import avg_pds_from_iterable, unnormalize_periodograms
 from .fourier import avg_pds_from_events
 from .fourier import fftfreq, fft
 from .fourier import get_flux_iterable_from_segments
+from .fourier import rms_calculation, poisson_level
 
 try:
     from tqdm import tqdm as show_progress
@@ -166,7 +167,8 @@ def rebin(self, df=None, f=None, method="mean"):
 
         return bin_ps
 
-    def compute_rms(self, min_freq, max_freq, white_noise_offset=2.):
+    def compute_rms(self, min_freq, max_freq, poisson_noise_level=None, 
+        white_noise_offset=None, deadtime=0.):
         """
         Compute the fractional rms amplitude in the power spectrum
         between two frequencies.
@@ -181,7 +183,17 @@ def compute_rms(self, min_freq, max_freq, white_noise_offset=2.):
 
         Other parameters
         ----------------
-        white_noise_offset : float, default 0
+        poisson_noise_level : float, default is None
+            This is the Poisson noise level of the PDS with same 
+            normalization as the PDS. If poissoin_noise_level is None, 
+            the Poisson noise is calculated in the idealcase 
+            e.g. 2./<countrate> for fractional rms normalisation 
+            Dead time and other instrumental effects can alter it. 
+            The user can fit the Poisson noise level outside 
+            this function using the same normalisation of the PDS 
+            and it will get subtracted from powers here.
+
+        white_noise_offset : float, default None
             This is the white noise level, in Leahy normalization. In the ideal
             case, this is 2. Dead time and other instrumental effects can alter
             it. The user can fit the white noise level outside this function
@@ -199,19 +211,46 @@ def compute_rms(self, min_freq, max_freq, white_noise_offset=2.):
         """
         minind = self.freq.searchsorted(min_freq)
         maxind = self.freq.searchsorted(max_freq)
-        powers = self.power[minind:maxind]
         nphots = self.nphots
-
-        if self.norm.lower() == 'leahy':
-            powers_leahy = powers.copy()
+        # distinguish the rebinned and non-rebinned case
+        if isinstance(self.m, Iterable):
+            M_freq = self.m[minind:maxind]
+            K_freq = self.k[minind:maxind]
+            freq_bins = 1
         else:
-            powers_leahy = \
-                self.unnorm_power[minind:maxind].real * 2 / nphots
+            M_freq = self.m
+            K_freq = self.k
+            freq_bins = maxind - minind
+        T = self.dt * self.n 
+
+        if  white_noise_offset is not None:
+            powers = self.power[minind:maxind]
+            warnings.warn(
+            "the option white_noise_offset now deprecated and will be "
+            "removed in the next major release. The routine"
+            "is correct only with non-rebinned power-spectra.",
+            DeprecationWarning)
+
+            if self.norm.lower() == 'leahy':
+                powers_leahy = powers.copy()
+            else:
+                powers_leahy = \
+                    self.unnorm_power[minind:maxind].real * 2 / nphots
 
-        rms = np.sqrt(np.sum(powers_leahy - white_noise_offset) / nphots)
-        rms_err = self._rms_error(powers_leahy)
+            rms = np.sqrt(np.sum(powers_leahy - white_noise_offset) / nphots)
+            rms_err = self._rms_error(powers_leahy)
+            return rms, rms_err
+
+        else:
+            if poisson_noise_level is None: 
+                poisson_noise_unnorm = poisson_level('none', n_ph=self.nphots) 
+            else:   
+                poisson_noise_unnorm = unnormalize_periodograms(poisson_noise_level, 
+                self.dt, self.n, self.nphots, norm=self.norm)
+            return rms_calculation(self.unnorm_power[minind:maxind], self.freq[minind], 
+                self.freq[maxind], self.nphots, T, M_freq, K_freq, freq_bins, 
+                poisson_noise_unnorm, deadtime)
 
-        return rms, rms_err
 
     def _rms_error(self, powers):
         """
@@ -650,6 +689,7 @@ def _initialize_empty(self):
         self.nphots1 = None
         self.m = 1
         self.n = None
+        self.k = 1
         return
 
 
@@ -776,6 +816,7 @@ def __init__(self, data=None, segment_size=None, norm="frac", gti=None,
         self.save_all = save_all
         self.segment_size = segment_size
         self.show_progress = not silent
+        self.k = 1
 
         if not good_input:
             return self._initialize_empty()

stingray/tests/test_crossspectrum.py

@@ -395,9 +395,9 @@ def test_coherence_is_one_on_single_interval(self):
             cs = Crossspectrum(lc1, lc2)
             coh = cs.coherence()
 
-        assert len(coh) == 2
+        assert np.isclose(len(coh), 2, rtol = 0.001)
         # The raw coherence of a single interval is 1 by definition
-        assert np.abs(np.mean(coh)) == 1
+        assert np.isclose(np.abs(np.mean(coh)), 1, rtol = 0.001)
 
     def test_high_coherence(self):
         t = np.arange(1280)
@@ -1176,8 +1176,8 @@ def test_rebin_factor_log_rebins_all_attrs(self, norm):
             assert hasattr(new_cs.pds2, attr) and getattr(new_cs.pds2, attr).size == N
 
         for attr in cs1.meta_attrs():
-            if attr not in ["df", "gti", "m"]:
-                assert getattr(cs1, attr) == getattr(new_cs, attr)
+            if attr not in ["df", "gti", "m", "k"]:
+                assert np.all(getattr(cs1, attr) == getattr(new_cs, attr))
 
     def test_rebin(self):
         with warnings.catch_warnings(record=True) as w:
@@ -1289,9 +1289,8 @@ def test_coherence_computes_correctly(self):
         with pytest.warns(DeprecationWarning):
             coh = coherence(self.lc1, self.lc2)
 
-        assert len(coh) == 2
-        assert np.abs(np.mean(coh)) == 1
-
+        assert np.isclose(len(coh), 2, rtol = 0.001)
+        assert np.isclose(np.abs(np.mean(coh)), 1, rtol = 0.001)
 
 class TestTimelagFunction(object):