Merge pull request #672 from StingraySoftware/spectral_timing_fixes

Spectral timing fixes
StingraySoftware · Oct 1, 2022 · 91d7af6 · 91d7af6
2 parents 56e8ab1 + 93d785a
commit 91d7af6
Show file tree

Hide file tree

Showing 7 changed files with 66 additions and 40 deletions.
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -1,6 +1,21 @@
 Changelog
 =========
 
+v1.1 (unreleased)
+-----------------
+Bug fixes
+^^^^^^^^^
+- IMPORTANT: Fixed sign of time lags, which were calculated using the interest band as the reference.
+- Fixed an issue when the fractional exposure in FITS light curves is slightly >1 (as sometimes happens in NICER data)
+
+New
+^^^
+- Implemented the ``bexvar`` variability estimation method for light curves.
+
+Improvements
+^^^^^^^^^^^^
+- A less confusing default value of segment_size in Z searches
+
 v1.0 (2022-03-29)
 ---------------------
 TL,DR: these things will break your code with v1.0:

diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
@@ -503,9 +503,9 @@ class Crossspectrum(StingrayObject):
     ----------------
     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 
+        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
@@ -928,7 +928,7 @@ def _fourier_cross(self, lc1, lc2, fullspec=False):
         fourier_2 = fft(lc2.counts)  # do Fourier transform 2
 
         freqs = fftfreq(lc1.n, lc1.dt)
-        cross = np.multiply(fourier_1, np.conj(fourier_2))
+        cross = np.multiply(fourier_2, np.conj(fourier_1))
 
         if fullspec is True:
             return freqs, cross
@@ -1722,13 +1722,13 @@ class AveragedCrossspectrum(Crossspectrum):
     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 
+        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 
+        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
@@ -1745,9 +1745,9 @@ class AveragedCrossspectrum(Crossspectrum):
     ----------------
     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 
+        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
@@ -1794,16 +1794,16 @@ class AveragedCrossspectrum(Crossspectrum):
         after the average.
 
     legacy: bool
-        Use the legacy machinery of `AveragedCrossspectrum`. This might be 
-        useful to compare with old results, and is also needed to use light 
-        curve lists as an input, to conserve the spectra of each segment, or 
+        Use the legacy machinery of `AveragedCrossspectrum`. This might be
+        useful to compare with old results, and is also needed to use light
+        curve lists as an input, to conserve the spectra of each segment, or
         to use the large_data option.
 
     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 
+        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
@@ -1818,8 +1818,8 @@ class AveragedCrossspectrum(Crossspectrum):
         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 
+        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

diff --git a/stingray/fourier.py b/stingray/fourier.py
@@ -526,7 +526,7 @@ def bias_term(power1, power2, power1_noise, power2_noise, n_ave,
     bias : float `np.array`, same shape as ``power1`` and ``power2``
         The bias term
     """
-    if n_ave > 500:
+    if (isinstance(n_ave, Iterable) and np.all(n_ave > 500)) or (not isinstance(n_ave, Iterable) and n_ave > 500):
         return 0. * power1
     bsq = power1 * power2 - intrinsic_coherence * (power1 - power1_noise) * \
           (power2 - power2_noise)
@@ -571,6 +571,7 @@ def raw_coherence(cross_power, power1, power2, power1_noise, power2_noise,
     if isinstance(num, Iterable):
         num[num < 0] = (cross_power * np.conj(cross_power)).real[num < 0]
     elif num < 0:
+        warnings.warn("Negative numerator in raw_coherence calculation. Setting bias term to 0")
         num = (cross_power * np.conj(cross_power)).real
     den = power1 * power2
     return num / den
@@ -663,7 +664,7 @@ def estimate_intrinsic_coherence(cross_power, power1, power2, power1_noise,
 
 def error_on_averaged_cross_spectrum(cross_power, seg_power, ref_power, n_ave,
                                      seg_power_noise, ref_power_noise,
-                                     common_ref="False"):
+                                     common_ref=False):
     """
     Error on cross spectral quantities, From Ingram 2019.
 
@@ -702,7 +703,7 @@ def error_on_averaged_cross_spectrum(cross_power, seg_power, ref_power, n_ave,
         Error on the modulus of the cross spectrum
 
     """
-    if n_ave < 30:
+    if (not isinstance(n_ave, Iterable) and n_ave < 30) or (isinstance(n_ave, Iterable) and np.any(n_ave) < 30):
         warnings.warn("n_ave is below 30. Please note that the error bars "
                       "on the quantities derived from the cross spectrum "
                       "are only reliable for a large number of averaged "
@@ -1160,7 +1161,7 @@ def avg_cs_from_iterables_quick(
         ft2 = fft(flux2)
 
         # Calculate the unnormalized cross spectrum
-        unnorm_power = ft1 * ft2.conj()
+        unnorm_power = ft1.conj() * ft2
 
         # Accumulate the sum to calculate the total mean of the lc
         sum_of_photons1 += n_ph1
@@ -1358,7 +1359,7 @@ def local_show_progress(a):
         n_ph = np.sqrt(n_ph1 * n_ph2)
 
         # Calculate the unnormalized cross spectrum
-        unnorm_power = ft1 * ft2.conj()
+        unnorm_power = ft1.conj() * ft2
         unnorm_pd1 = unnorm_pd2 = 0
 
         # If requested, calculate the auxiliary PDSs

diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
@@ -1219,6 +1219,9 @@ def test_timelag(self):
                                   test_lc1.counts,
                                   err_dist=test_lc1.err_dist,
                                   dt=dt)
+            # The second light curve is delayed by two bins.
+            # The time lag should be -2 * dt, because this will
+            # become the reference band in AveragedCrossspectrum
             test_lc2 = Lightcurve(test_lc1.time,
                                   np.array(np.roll(test_lc1.counts, 2)),
                                   err_dist=test_lc1.err_dist,
@@ -1230,7 +1233,9 @@ def test_timelag(self):
 
             time_lag, time_lag_err = cs.time_lag()
 
-        assert np.all(np.abs(time_lag[:6] - 0.1) < 3 * time_lag_err[:6])
+        # The actual measured time lag will be half that for AveragedCrosspectrum
+        measured_lag = -dt
+        assert np.all(np.abs(time_lag[:6] - measured_lag) < 3 * time_lag_err[:6])
 
     def test_errorbars_legacy(self):
         time = np.arange(10000) * 0.1

diff --git a/stingray/tests/test_fourier.py b/stingray/tests/test_fourier.py
@@ -64,7 +64,7 @@ def setup_class(cls):
         good = (freq > 0) & (freq < 0.1)
         ft1, ft2 = ft1[good], ft2[good]
         cls.cross = normalize_periodograms(
-            ft1 * ft2.conj(), dt, cls.N, mean, norm="abs", power_type="all")
+            ft1.conj() * ft2, dt, cls.N, mean, norm="abs", power_type="all")
         cls.pds1 = normalize_periodograms(
             ft1 * ft1.conj(), dt, cls.N, mean, norm="abs", power_type="real")
         cls.pds2 = normalize_periodograms(

diff --git a/stingray/tests/test_varenergyspectrum.py b/stingray/tests/test_varenergyspectrum.py
@@ -358,35 +358,36 @@ def setup_class(cls):
         data = np.load(os.path.join(datadir, "sample_variable_lc.npy"))
         flux = data
         times = np.arange(data.size) * dt
-        maxfreq = 0.25 / cls.time_lag
+        maxfreq = 0.15 / cls.time_lag
         roll_amount = int(cls.time_lag // dt)
         good = slice(roll_amount, roll_amount + int(200 // dt))
 
+        # When rolling, a positive delay is introduced
         rolled_flux = np.array(np.roll(flux, roll_amount))
         times, flux, rolled_flux = times[good], flux[good], rolled_flux[good]
 
         length = times[-1] - times[0]
-        test_lc1 = Lightcurve(times, flux, err_dist="gauss", dt=dt, skip_checks=True)
-        test_lc2 = Lightcurve(test_lc1.time, rolled_flux, err_dist=test_lc1.err_dist, dt=dt)
-        test_ev1, test_ev2 = EventList(), EventList()
-        test_ev1.simulate_times(test_lc1)
-        test_ev2.simulate_times(test_lc2)
+        test_ref = Lightcurve(times, flux, err_dist="gauss", dt=dt, skip_checks=True)
+        test_sub = Lightcurve(test_ref.time, rolled_flux, err_dist=test_ref.err_dist, dt=dt)
+        test_ref_ev, test_sub_ev = EventList(), EventList()
+        test_ref_ev.simulate_times(test_ref)
+        test_sub_ev.simulate_times(test_sub)
 
-        test_ev1.energy = np.random.uniform(0.3, 9, len(test_ev1.time))
-        test_ev2.energy = np.random.uniform(9, 12, len(test_ev2.time))
+        test_sub_ev.energy = np.random.uniform(0.3, 9, len(test_sub_ev.time))
+        test_ref_ev.energy = np.random.uniform(9, 12, len(test_ref_ev.time))
 
         cls.lag = LagEnergySpectrum(
-            test_ev1,
-            freq_interval=[0, maxfreq],
+            test_sub_ev,
+            freq_interval=[maxfreq / 2, maxfreq],
             energy_spec=(0.3, 9, 1, "lin"),
             ref_band=[9, 12],
             bin_time=dt / 2,
             segment_size=length,
-            events2=test_ev2,
+            events2=test_ref_ev,
         )
 
         # Make single event list
-        test_ev = test_ev1.join(test_ev2)
+        test_ev = test_sub_ev.join(test_ref_ev)
 
         cls.lag_same = LagEnergySpectrum(
             test_ev,

diff --git a/stingray/varenergyspectrum.py b/stingray/varenergyspectrum.py
@@ -875,9 +875,13 @@ def _spectrum_function(self):
                 Cmean, mean_sub_power, mean_ref_power, m_tot, sub_power_noise, ref_power_noise, common_ref=common_ref
             )
 
-            lag = np.mean((np.angle(cross[good]) / (2 * np.pi * freq[good])))
+            # The frequency of these lags is measured from the *weighted* mean of the frequencies
+            # in the cross spectrum. The weight is just the absolute value of the CS
+            csabs = np.abs(cross[good])
+            fmean = np.sum(freq[good] * csabs) / np.sum(csabs)
+            lag = np.angle(Cmean) / (2 * np.pi * fmean)
 
-            lag_e = phi_e / (2 * np.pi * f)
+            lag_e = phi_e / (2 * np.pi * fmean)
             self.spectrum[i] = lag
             self.spectrum_error[i] = lag_e