From 240c2411d6538e215c63062519bed9c46b9df487 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Thu, 26 Oct 2023 15:20:19 +0200
Subject: [PATCH 01/21] Fix bug with time resolution when rebinning

---
 stingray/powerspectrum.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index da79909ff..5b9969287 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -1082,7 +1082,7 @@ def rebin_time(self, dt_new, method="sum"):
         Parameters
         ----------
         dt_new: float
-            The new time resolution of  the dynamical power spectrum.
+            The new time resolution of the dynamical power spectrum.
             Must be larger than the time resolution of the old dynamical power
             spectrum!
 
@@ -1099,14 +1099,14 @@ def rebin_time(self, dt_new, method="sum"):
             New rebinned Dynamical Power Spectrum.
         """
         if dt_new < self.dt:
-            raise ValueError("New time resolution must be larger than " "old time resolution!")
+            raise ValueError("New time resolution must be larger than old time resolution!")
 
         new_dynspec_object = copy.deepcopy(self)
 
         dynspec_new = []
         for data in self.dyn_ps:
             time_new, bin_counts, bin_err, _ = utils.rebin_data(
-                self.time, data, dt_new, method=method
+                self.time, data, dt_new, method=method, dx=self.dt
             )
             dynspec_new.append(bin_counts)
 

From 48e0f6445e0c71d5c6cea24b6d4dce4399b64137 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Thu, 26 Oct 2023 15:25:30 +0200
Subject: [PATCH 02/21] Add changelog

---
 docs/changes/779.bugfix.rst | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 docs/changes/779.bugfix.rst

diff --git a/docs/changes/779.bugfix.rst b/docs/changes/779.bugfix.rst
new file mode 100644
index 000000000..fd209edcc
--- /dev/null
+++ b/docs/changes/779.bugfix.rst
@@ -0,0 +1 @@
+Various bug fixes in DynamicalPowerspectrum, on event loading and time rebinning

From f11273a2ba40424a9ea72d16b77f01dc1878fd83 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Thu, 26 Oct 2023 15:20:04 +0200
Subject: [PATCH 03/21] Fix initialization with event lists

---
 stingray/powerspectrum.py | 52 +++++++++++++++++----------------------
 1 file changed, 23 insertions(+), 29 deletions(-)

diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 5b9969287..0758988a9 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -950,19 +950,21 @@ class DynamicalPowerspectrum(AveragedPowerspectrum):
     def __init__(self, lc, segment_size, norm="frac", gti=None, dt=None):
         if isinstance(lc, EventList) and dt is None:
             raise ValueError("To pass an input event lists, please specify dt")
+        elif isinstance(lc, Lightcurve):
+            dt = lc.dt
+            if segment_size < 2 * lc.dt:
+                raise ValueError("Length of the segment is too short to form a " "light curve!")
+            elif segment_size > lc.tseg:
+                raise ValueError(
+                    "Length of the segment is too long to create "
+                    "any segments of the light curve!"
+                )
 
-        if isinstance(lc, EventList):
-            lc = lc.to_lc(dt)
+        self.segment_size = segment_size
+        self.input_dt = dt
+        self.gti = gti
+        self.norm = norm
 
-        if segment_size < 2 * lc.dt:
-            raise ValueError("Length of the segment is too short to form a " "light curve!")
-        elif segment_size > lc.tseg:
-            raise ValueError(
-                "Length of the segment is too long to create " "any segments of the light curve!"
-            )
-        AveragedPowerspectrum.__init__(
-            self, data=lc, segment_size=segment_size, norm=norm, gti=gti, dt=dt
-        )
         self._make_matrix(lc)
 
     def _make_matrix(self, lc):
@@ -978,33 +980,25 @@ def _make_matrix(self, lc):
             power spectrum.
         """
         avg = AveragedPowerspectrum(
-            lc, segment_size=self.segment_size, norm=self.norm, gti=self.gti, save_all=True
+            lc,
+            dt=self.input_dt,
+            segment_size=self.segment_size,
+            norm=self.norm,
+            gti=self.gti,
+            save_all=True,
         )
         self.dyn_ps = np.array(avg.cs_all).T
 
         self.freq = avg.freq
         current_gti = avg.gti
 
-        start_inds, end_inds = bin_intervals_from_gtis(
-            current_gti, self.segment_size, lc.time, dt=lc.dt
-        )
+        from .gti import time_intervals_from_gtis
 
-        tstart = lc.time[start_inds]
-        tend = lc.time[end_inds]
+        tstart, tend = time_intervals_from_gtis(current_gti, self.segment_size)
 
         self.time = tstart + 0.5 * (tend - tstart)
-
-        # Assign length of lightcurve as time resolution if only one value
-        if len(self.time) > 1:
-            self.dt = self.time[1] - self.time[0]
-        else:
-            self.dt = lc.n
-
-        # Assign biggest freq. resolution if only one value
-        if len(self.freq) > 1:
-            self.df = self.freq[1] - self.freq[0]
-        else:
-            self.df = 1 / lc.n
+        self.df = self.freq[1] - self.freq[0]
+        self.dt = self.segment_size
 
     def rebin_frequency(self, df_new, method="sum"):
         """

From 8cab3857769cf03eccb28bb5383bafb85d4766fc Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Thu, 26 Oct 2023 15:28:02 +0200
Subject: [PATCH 04/21] Fix check on input dt

---
 stingray/powerspectrum.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 0758988a9..84a3329f5 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -952,13 +952,13 @@ def __init__(self, lc, segment_size, norm="frac", gti=None, dt=None):
             raise ValueError("To pass an input event lists, please specify dt")
         elif isinstance(lc, Lightcurve):
             dt = lc.dt
-            if segment_size < 2 * lc.dt:
-                raise ValueError("Length of the segment is too short to form a " "light curve!")
-            elif segment_size > lc.tseg:
+            if segment_size > lc.tseg:
                 raise ValueError(
                     "Length of the segment is too long to create "
                     "any segments of the light curve!"
                 )
+        if segment_size < 2 * dt:
+            raise ValueError("Length of the segment is too short to form a light curve!")
 
         self.segment_size = segment_size
         self.input_dt = dt

From 545af310ec4fddc39be2c957f10c39e32585324a Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Thu, 26 Oct 2023 17:59:42 +0200
Subject: [PATCH 05/21] Fix bug when single frequency bin

---
 stingray/powerspectrum.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 84a3329f5..5f93a7c0b 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -988,7 +988,6 @@ def _make_matrix(self, lc):
             save_all=True,
         )
         self.dyn_ps = np.array(avg.cs_all).T
-
         self.freq = avg.freq
         current_gti = avg.gti
 
@@ -997,7 +996,10 @@ def _make_matrix(self, lc):
         tstart, tend = time_intervals_from_gtis(current_gti, self.segment_size)
 
         self.time = tstart + 0.5 * (tend - tstart)
-        self.df = self.freq[1] - self.freq[0]
+        if len(self.freq) < 2:
+            self.df = 1 / self.segment_size
+        else:
+            self.df = self.freq[1] - self.freq[0]
         self.dt = self.segment_size
 
     def rebin_frequency(self, df_new, method="sum"):

From 97480fa2aeefef697c5e8ca61251d18fceed220e Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Thu, 26 Oct 2023 17:59:52 +0200
Subject: [PATCH 06/21] Fix tests, making them more predictable

---
 stingray/tests/test_powerspectrum.py | 114 ++++++++++-----------------
 1 file changed, 40 insertions(+), 74 deletions(-)

diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index 5b44516d7..f1b5ace85 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -151,7 +151,6 @@ def test_method_norm(self, norm, err_dist):
 
         assert loc_pds.norm == renorm_pds.norm
         for attr in ["power", "unnorm_power", "power_err"]:
-            print(attr)
             loc = getattr(loc_pds, attr)
             renorm = getattr(renorm_pds, attr)
             assert np.allclose(loc, renorm, atol=0.5)
@@ -965,8 +964,8 @@ def test_leahy_correct_for_multiple(self, use_common_mean):
 class TestDynamicalPowerspectrum(object):
     def setup_class(cls):
         # generate timestamps
-        timestamps = np.linspace(1, 100, 10000)
-        dt = np.median(np.diff(timestamps))
+        timestamps = np.arange(0.005, 100.01, 0.01)
+        dt = 0.01
         freq = 25 + 1.2 * np.sin(2 * np.pi * timestamps / 130)
         # variability signal with drifiting frequency
         vari = 25 * np.sin(2 * np.pi * freq * timestamps)
@@ -975,9 +974,8 @@ def setup_class(cls):
         with warnings.catch_warnings():
             warnings.simplefilter("ignore", category=UserWarning)
 
-            lc = Lightcurve(
-                timestamps, signal, err_dist="poisson", dt=dt, gti=[[1 - dt / 2, 100 + dt / 2]]
-            )
+            lc = Lightcurve(timestamps, signal, err_dist="poisson", dt=dt, gti=[[0, 100]])
+
         cls.lc = lc
 
         # Simple lc to demonstrate rebinning of dyn ps
@@ -1044,43 +1042,39 @@ def test_size_of_trace_maximum(self):
 
     def test_rebin_small_dt(self):
         segment_size = 3
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore", category=UserWarning)
-            dps = DynamicalPowerspectrum(self.lc_test, segment_size=segment_size)
+        dps = DynamicalPowerspectrum(self.lc_test, segment_size=segment_size)
         with pytest.raises(ValueError):
             dps.rebin_time(dt_new=2.0)
 
     def test_rebin_small_df(self):
         segment_size = 3
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore", category=UserWarning)
-            dps = DynamicalPowerspectrum(self.lc, segment_size=segment_size)
+        dps = DynamicalPowerspectrum(self.lc, segment_size=segment_size)
         with pytest.raises(ValueError):
             dps.rebin_frequency(df_new=dps.df / 2.0)
 
     def test_rebin_time_default_method(self):
         segment_size = 3
-        dt_new = 4.0
-        rebin_time = np.array([2.0, 6.0, 10.0])
-        rebin_dps = np.array([[0.7962963, 1.16402116, 0.28571429]])
+        dt_new = 6.0
+        rebin_time = np.array([2.5, 8.5])
+        rebin_dps = np.array([[1.73611111, 0.81018519]])
         dps = DynamicalPowerspectrum(self.lc_test, segment_size=segment_size)
-        print(dps.dyn_ps)
         new_dps = dps.rebin_time(dt_new=dt_new)
         assert np.allclose(new_dps.time, rebin_time)
-        # assert np.allclose(new_dps.dyn_ps, rebin_dps)
+        assert np.allclose(new_dps.dyn_ps, rebin_dps)
         assert np.isclose(new_dps.dt, dt_new)
 
     def test_rebin_frequency_default_method(self):
         segment_size = 50
         df_new = 10.0
-        rebin_freq = np.array([5.01000198, 15.01000198, 25.01000198, 35.01000198, 45.01000198])
+        rebin_freq = np.array([5.01, 15.01, 25.01, 35.01])
         rebin_dps = np.array(
             [
-                [5.76369293e-06],
-                [7.07524761e-05],
-                [6.24846189e00],
-                [5.77470465e-05],
-                [1.76918128e-05],
+                [
+                    [1.71989342e-04, 6.42756881e-05],
+                    [7.54455204e-04, 2.14785049e-04],
+                    [6.24831554e00, 6.24984615e00],
+                    [6.71135792e-04, 7.42516599e-05],
+                ]
             ]
         )
         with warnings.catch_warnings():
@@ -1091,71 +1085,43 @@ def test_rebin_frequency_default_method(self):
         assert np.allclose(new_dps.dyn_ps, rebin_dps, atol=0.01)
         assert np.isclose(new_dps.df, df_new)
 
-    def test_rebin_time_mean_method(self):
+    @pytest.mark.parametrize("method", ["mean", "average"])
+    def test_rebin_time_mean_method(self, method):
         segment_size = 3
-        dt_new = 4.0
-        rebin_time = np.array([2.0, 6.0, 10.0])
-        rebin_dps = np.array([[0.59722222, 0.87301587, 0.21428571]])
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore", category=UserWarning)
-            dps = DynamicalPowerspectrum(self.lc_test, segment_size=segment_size)
-        new_dps = dps.rebin_time(dt_new=dt_new, method="mean")
-        assert np.allclose(new_dps.time, rebin_time)
-        # assert np.allclose(new_dps.dyn_ps, rebin_dps)
-        assert np.isclose(new_dps.dt, dt_new)
-
-    def test_rebin_frequency_mean_method(self):
-        segment_size = 50
-        df_new = 10.0
-        rebin_freq = np.array([5.01000198, 15.01000198, 25.01000198, 35.01000198, 45.01000198])
-        rebin_dps = np.array(
-            [
-                [1.15296690e-08],
-                [1.41532979e-07],
-                [1.24993989e-02],
-                [1.15516968e-07],
-                [3.53906336e-08],
-            ]
-        )
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore", category=UserWarning)
-            dps = DynamicalPowerspectrum(self.lc, segment_size=segment_size)
-        new_dps = dps.rebin_frequency(df_new=df_new, method="mean")
-        assert np.allclose(new_dps.freq, rebin_freq)
-        assert np.allclose(new_dps.dyn_ps, rebin_dps, atol=0.00001)
-        assert np.isclose(new_dps.df, df_new)
-
-    def test_rebin_time_average_method(self):
-        segment_size = 3
-        dt_new = 4.0
-        rebin_time = np.array([2.0, 6.0, 10.0])
-        rebin_dps = np.array([[0.59722222, 0.87301587, 0.21428571]])
+        dt_new = 6.0
+        rebin_time = np.array([2.5, 8.5])
+        rebin_dps = np.array([[1.73611111, 0.81018519]]) / 2
 
         with warnings.catch_warnings():
             warnings.simplefilter("ignore", category=UserWarning)
             dps = DynamicalPowerspectrum(self.lc_test, segment_size=segment_size)
-        new_dps = dps.rebin_time(dt_new=dt_new, method="average")
+        new_dps = dps.rebin_time(dt_new=dt_new, method=method)
         assert np.allclose(new_dps.time, rebin_time)
-        # assert np.allclose(new_dps.dyn_ps, rebin_dps)
+        assert np.allclose(new_dps.dyn_ps, rebin_dps)
         assert np.isclose(new_dps.dt, dt_new)
 
-    def test_rebin_frequency_average_method(self):
+    @pytest.mark.parametrize("method", ["mean", "average"])
+    def test_rebin_frequency_mean_method(self, method):
         segment_size = 50
         df_new = 10.0
-        rebin_freq = np.array([5.01000198, 15.01000198, 25.01000198, 35.01000198, 45.01000198])
-        rebin_dps = np.array(
-            [
-                [1.15296690e-08],
-                [1.41532979e-07],
-                [1.24993989e-02],
-                [1.15516968e-07],
-                [3.53906336e-08],
-            ]
+        rebin_freq = np.array([5.01, 15.01, 25.01, 35.01])
+        rebin_dps = (
+            np.array(
+                [
+                    [
+                        [1.71989342e-04, 6.42756881e-05],
+                        [7.54455204e-04, 2.14785049e-04],
+                        [6.24831554e00, 6.24984615e00],
+                        [6.71135792e-04, 7.42516599e-05],
+                    ]
+                ]
+            )
+            / 500
         )
         with warnings.catch_warnings():
             warnings.simplefilter("ignore", category=UserWarning)
             dps = DynamicalPowerspectrum(self.lc, segment_size=segment_size)
-        new_dps = dps.rebin_frequency(df_new=df_new, method="average")
+        new_dps = dps.rebin_frequency(df_new=df_new, method=method)
         assert np.allclose(new_dps.freq, rebin_freq)
         assert np.allclose(new_dps.dyn_ps, rebin_dps, atol=0.00001)
         assert np.isclose(new_dps.df, df_new)

From ba422d7920e2b52ea004d9bab930ae8efc54fa00 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Mon, 30 Oct 2023 10:19:57 +0100
Subject: [PATCH 07/21] Just use df and segment size, without reinventing the
 wheel

---
 stingray/powerspectrum.py | 15 +++++----------
 1 file changed, 5 insertions(+), 10 deletions(-)

diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 5f93a7c0b..dfac1571f 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -13,11 +13,11 @@
 from stingray.stats import pds_probability, amplitude_upper_limit
 
 from .events import EventList
-from .gti import cross_two_gtis
+from .gti import cross_two_gtis, time_intervals_from_gtis
+
 from .lightcurve import Lightcurve
 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
 
@@ -991,15 +991,10 @@ def _make_matrix(self, lc):
         self.freq = avg.freq
         current_gti = avg.gti
 
-        from .gti import time_intervals_from_gtis
-
-        tstart, tend = time_intervals_from_gtis(current_gti, self.segment_size)
+        tstart, _ = time_intervals_from_gtis(current_gti, self.segment_size)
 
-        self.time = tstart + 0.5 * (tend - tstart)
-        if len(self.freq) < 2:
-            self.df = 1 / self.segment_size
-        else:
-            self.df = self.freq[1] - self.freq[0]
+        self.time = tstart + 0.5 * self.segment_size
+        self.df = avg.df
         self.dt = self.segment_size
 
     def rebin_frequency(self, df_new, method="sum"):

From d6bd42de8086b4bbe598645f2cbca3bbd0b4cd80 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Mon, 30 Oct 2023 11:04:58 +0100
Subject: [PATCH 08/21] Stub dynamical cross spectrum because why not

---
 stingray/crossspectrum.py            | 236 ++++++++++++++++++++++++++-
 stingray/powerspectrum.py            | 115 +------------
 stingray/tests/test_crossspectrum.py | 166 ++++++++++++++++++-
 3 files changed, 401 insertions(+), 116 deletions(-)

diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 8a22b26b2..b1db2094b 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -15,6 +15,7 @@
 
 from .base import StingrayObject
 from .events import EventList
+from .gti import cross_two_gtis, time_intervals_from_gtis
 from .lightcurve import Lightcurve
 from .utils import show_progress
 from .fourier import avg_cs_from_iterables, error_on_averaged_cross_spectrum
@@ -28,6 +29,7 @@
 __all__ = [
     "Crossspectrum",
     "AveragedCrossspectrum",
+    "DynamicalCrossspectrum",
     "cospectra_pvalue",
     "normalize_crossspectrum",
     "time_lag",
@@ -377,7 +379,7 @@ def cospectra_pvalue(power, nspec):
     the data.
 
     Important: the underlying assumption that make this calculation valid
-    is that the powers in the power spectrum follow a Laplace distribution,
+    is that the powers in the cross spectrum follow a Laplace distribution,
     and this requires that:
 
     1. the co-spectrum is normalized according to [Leahy 1983]_
@@ -403,7 +405,7 @@ def cospectra_pvalue(power, nspec):
         the signal-to-noise ratio, i.e. makes the statistical distributions
         of the noise narrower, such that a smaller power might be very
         significant in averaged spectra even though it would not be in a single
-        power spectrum.
+        cross spectrum.
 
     Returns
     -------
@@ -1880,6 +1882,236 @@ def time_lag(self):
         return lag, lag_err
 
 
+class DynamicalCrossspectrum(AveragedCrossspectrum):
+    type = "crossspectrum"
+    """
+    Create a dynamical cross spectrum, also often called a *spectrogram*.
+
+    This class will divide two :class:`Lightcurve` objects into segments of
+    length ``segment_size``, create a cross spectrum for each segment and store
+    all powers in a matrix as a function of both time (using the mid-point of
+    each segment) and frequency.
+
+    This is often used to trace changes in period of a (quasi-)periodic signal
+    over time.
+
+    Parameters
+    ----------
+    data1 : :class:`stingray.Lightcurve` or :class:`stingray.EventList` object
+        The time series or event list from the subject band, or channel, for which
+        the dynamical cross spectrum is to be calculated.
+
+    data2 : :class:`stingray.Lightcurve` or :class:`stingray.EventList` object
+        The time series or event list from the reference band, or channel, of the dynamical
+        cross spectrum.
+
+    segment_size : float, default 1
+         Length of the segment of light curve, default value is 1 (in whatever
+         units the ``time`` array in the :class:`Lightcurve`` object uses).
+
+    norm: {"leahy" | "frac" | "abs" | "none" }, optional, default "frac"
+        The normaliation of the periodogram to be used.
+
+    Other Parameters
+    ----------------
+    gti: 2-d float array
+        ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]`` -- Good Time intervals.
+        This choice overrides the GTIs in the single light curves. Use with
+        care, especially 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 object before making
+        the cross spectrum.
+
+    Attributes
+    ----------
+    segment_size: float
+        The size of each segment to average. Note that if the total
+        duration of each input object in lc is not an integer multiple
+        of the ``segment_size``, then any fraction left-over at the end of the
+        time series will be lost.
+
+    dyn_ps : np.ndarray
+        The matrix of normalized squared absolute values of Fourier
+        amplitudes. The axis are given by the ``freq``
+        and ``time`` attributes.
+
+    norm: {``leahy`` | ``frac`` | ``abs`` | ``none``}
+        The normalization of the periodogram.
+
+    freq: numpy.ndarray
+        The array of mid-bin frequencies that the Fourier transform samples.
+
+    df: float
+        The frequency resolution.
+
+    dt: float
+        The time resolution.
+    """
+
+    def __init__(self, data1, data2, segment_size, norm="frac", gti=None, dt=None):
+        if isinstance(data1, EventList) and dt is None:
+            raise ValueError("To pass input event lists, please specify dt")
+        elif isinstance(data1, Lightcurve):
+            dt = data1.dt
+            if segment_size > data1.tseg:
+                raise ValueError(
+                    "Length of the segment is too long to create "
+                    "any segments of the light curve!"
+                )
+        if segment_size < 2 * dt:
+            raise ValueError("Length of the segment is too short to form a light curve!")
+
+        self.segment_size = segment_size
+        self.input_dt = dt
+        self.gti = gti
+        self.norm = norm
+
+        self._make_matrix(data1, data2)
+
+    def _make_matrix(self, data1, data2):
+        """
+        Create a matrix of powers for each time step and each frequency step.
+
+        Time increases with row index, frequency with column index.
+
+        Parameters
+        ----------
+        data1 : :class:`Lightcurve` or :class:`EventList`
+            The object providing the subject band/channel for the dynamical
+            cross spectrum.
+        data2 : :class:`Lightcurve` or :class:`EventList`
+            The object providing the reference band for the dynamical
+            cross spectrum.
+        """
+        avg = AveragedCrossspectrum(
+            data1,
+            data2,
+            dt=self.input_dt,
+            segment_size=self.segment_size,
+            norm=self.norm,
+            gti=self.gti,
+            save_all=True,
+        )
+        self.dyn_ps = np.array(avg.cs_all).T
+        self.freq = avg.freq
+        current_gti = avg.gti
+
+        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
+
+    def rebin_frequency(self, df_new, method="sum"):
+        """
+        Rebin the Dynamic Power Spectrum to a new frequency resolution.
+        Rebinning is an in-place operation, i.e. will replace the existing
+        ``dyn_ps`` attribute.
+
+        While the new resolution need not be an integer multiple of the
+        previous frequency resolution, be aware that if it is not, the last
+        bin will be cut off by the fraction left over by the integer division.
+
+        Parameters
+        ----------
+        df_new: float
+            The new frequency resolution of the dynamical power spectrum.
+            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
+            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(
+                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
+        return new_dynspec_object
+
+    def rebin_time(self, dt_new, method="sum"):
+        """
+        Rebin the Dynamic Power Spectrum to a new time resolution.
+        While the new resolution need not be an integer multiple of the
+        previous time resolution, be aware that if it is not, the last bin
+        will be cut off by the fraction left over by the integer division.
+
+        Parameters
+        ----------
+        dt_new: float
+            The new time resolution of the dynamical power spectrum.
+            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
+            should be summed or averaged.
+
+        Returns
+        -------
+        time_new: numpy.ndarray
+            Time axis with new rebinned time resolution.
+
+        dynspec_new: numpy.ndarray
+            New rebinned Dynamical Power Spectrum.
+        """
+        if dt_new < self.dt:
+            raise ValueError("New time resolution must be larger than old time resolution!")
+
+        new_dynspec_object = copy.deepcopy(self)
+
+        dynspec_new = []
+        for data in self.dyn_ps:
+            time_new, bin_counts, _, _ = 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
+        return new_dynspec_object
+
+    def trace_maximum(self, min_freq=None, max_freq=None):
+        """
+        Return the indices of the maximum powers in each segment
+        :class:`Powerspectrum` between specified frequencies.
+
+        Parameters
+        ----------
+        min_freq: float, default ``None``
+            The lower frequency bound.
+
+        max_freq: float, default ``None``
+            The upper frequency bound.
+
+        Returns
+        -------
+        max_positions : np.array
+            The array of indices of the maximum power in each segment having
+            frequency between ``min_freq`` and ``max_freq``.
+        """
+        if min_freq is None:
+            min_freq = np.min(self.freq)
+        if max_freq is None:
+            max_freq = np.max(self.freq)
+
+        max_positions = []
+        for ps in self.dyn_ps.T:
+            indices = np.logical_and(self.freq <= max_freq, min_freq <= self.freq)
+            max_power = np.max(ps[indices])
+            max_positions.append(np.where(ps == max_power)[0][0])
+
+        return np.array(max_positions)
+
+
 def crossspectrum_from_time_array(
     times1,
     times2,
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index dfac1571f..32ff8fa69 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -7,9 +7,7 @@
 import scipy.optimize
 import scipy.stats
 
-import stingray.utils as utils
-from stingray.crossspectrum import AveragedCrossspectrum, Crossspectrum
-from stingray.gti import bin_intervals_from_gtis, check_gtis
+from stingray.crossspectrum import AveragedCrossspectrum, Crossspectrum, DynamicalCrossspectrum
 from stingray.stats import pds_probability, amplitude_upper_limit
 
 from .events import EventList
@@ -885,7 +883,7 @@ def initial_checks(self, *args, **kwargs):
         return AveragedCrossspectrum.initial_checks(self, *args, **kwargs)
 
 
-class DynamicalPowerspectrum(AveragedPowerspectrum):
+class DynamicalPowerspectrum(DynamicalCrossspectrum):
     type = "powerspectrum"
     """
     Create a dynamical power spectrum, also often called a *spectrogram*.
@@ -997,115 +995,6 @@ def _make_matrix(self, lc):
         self.df = avg.df
         self.dt = self.segment_size
 
-    def rebin_frequency(self, df_new, method="sum"):
-        """
-        Rebin the Dynamic Power Spectrum to a new frequency resolution.
-        Rebinning is an in-place operation, i.e. will replace the existing
-        ``dyn_ps`` attribute.
-
-        While the new resolution need not be an integer multiple of the
-        previous frequency resolution, be aware that if it is not, the last
-        bin will be cut off by the fraction left over by the integer division.
-
-        Parameters
-        ----------
-        df_new: float
-            The new frequency resolution of the dynamical power spectrum.
-            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
-            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, _ = utils.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
-        return new_dynspec_object
-
-    def trace_maximum(self, min_freq=None, max_freq=None):
-        """
-        Return the indices of the maximum powers in each segment
-        :class:`Powerspectrum` between specified frequencies.
-
-        Parameters
-        ----------
-        min_freq: float, default ``None``
-            The lower frequency bound.
-
-        max_freq: float, default ``None``
-            The upper frequency bound.
-
-        Returns
-        -------
-        max_positions : np.array
-            The array of indices of the maximum power in each segment having
-            frequency between ``min_freq`` and ``max_freq``.
-        """
-        if min_freq is None:
-            min_freq = np.min(self.freq)
-        if max_freq is None:
-            max_freq = np.max(self.freq)
-
-        max_positions = []
-        for ps in self.dyn_ps.T:
-            indices = np.logical_and(self.freq <= max_freq, min_freq <= self.freq)
-            max_power = np.max(ps[indices])
-            max_positions.append(np.where(ps == max_power)[0][0])
-
-        return np.array(max_positions)
-
-    def rebin_time(self, dt_new, method="sum"):
-        """
-        Rebin the Dynamic Power Spectrum to a new time resolution.
-        While the new resolution need not be an integer multiple of the
-        previous time resolution, be aware that if it is not, the last bin
-        will be cut off by the fraction left over by the integer division.
-
-        Parameters
-        ----------
-        dt_new: float
-            The new time resolution of the dynamical power spectrum.
-            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
-            should be summed or averaged.
-
-        Returns
-        -------
-        time_new: numpy.ndarray
-            Time axis with new rebinned time resolution.
-
-        dynspec_new: numpy.ndarray
-            New rebinned Dynamical Power Spectrum.
-        """
-        if dt_new < self.dt:
-            raise ValueError("New time resolution must be larger than old time resolution!")
-
-        new_dynspec_object = copy.deepcopy(self)
-
-        dynspec_new = []
-        for data in self.dyn_ps:
-            time_new, bin_counts, bin_err, _ = utils.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
-        return new_dynspec_object
-
 
 def powerspectrum_from_time_array(
     times,
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index 32021b0da..0b95cc6b0 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -6,7 +6,7 @@
 import scipy.special
 from astropy.io import fits
 from stingray import Lightcurve
-from stingray import Crossspectrum, AveragedCrossspectrum
+from stingray import Crossspectrum, AveragedCrossspectrum, DynamicalCrossspectrum
 from stingray.crossspectrum import cospectra_pvalue
 from stingray.crossspectrum import normalize_crossspectrum, normalize_crossspectrum_gauss
 from stingray.crossspectrum import coherence, time_lag
@@ -1316,3 +1316,167 @@ def test_file_roundtrip(self, fmt):
         os.unlink(fname)
 
         self._check_equal(so, new_so)
+
+
+class TestDynamicalCrossspectrum(object):
+    def setup_class(cls):
+        # generate timestamps
+        timestamps = np.arange(0.005, 100.01, 0.01)
+        dt = 0.01
+        freq = 25 + 1.2 * np.sin(2 * np.pi * timestamps / 130)
+        # variability signal with drifiting frequency
+        vari = 25 * np.sin(2 * np.pi * freq * timestamps)
+        signal = vari + 50
+        # create a lightcurve
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=UserWarning)
+
+            lc = Lightcurve(timestamps, signal, err_dist="poisson", dt=dt, gti=[[0, 100]])
+
+        cls.lc = lc
+
+        # Simple lc to demonstrate rebinning of dyn ps
+        # Simple lc to demonstrate rebinning of dyn ps
+        test_times = np.arange(16)
+        test_counts = [2, 3, 1, 3, 1, 5, 2, 1, 4, 2, 2, 2, 3, 4, 1, 7]
+        cls.lc_test = Lightcurve(test_times, test_counts)
+
+    def test_with_short_seg_size(self):
+        with pytest.raises(ValueError):
+            dps = DynamicalCrossspectrum(self.lc, self.lc, segment_size=0)
+
+    def test_works_with_events(self):
+        lc = copy.deepcopy(self.lc)
+        lc.counts = np.floor(lc.counts)
+        ev = EventList.from_lc(lc)
+        dps = DynamicalCrossspectrum(lc, lc, segment_size=10)
+        with pytest.raises(ValueError):
+            # Without dt, it fails
+            _ = DynamicalCrossspectrum(ev, ev, segment_size=10)
+
+        dps_ev = DynamicalCrossspectrum(ev, ev, segment_size=10, dt=self.lc.dt)
+        assert np.allclose(dps.dyn_ps, dps_ev.dyn_ps)
+
+    def test_with_long_seg_size(self):
+        with pytest.raises(ValueError):
+            dps = DynamicalCrossspectrum(self.lc, self.lc, segment_size=1000)
+
+    def test_matrix(self):
+        dps = DynamicalCrossspectrum(self.lc, self.lc, segment_size=3)
+        nsegs = int(self.lc.tseg / dps.segment_size)
+        nfreq = int((1 / self.lc.dt) / (2 * (dps.freq[1] - dps.freq[0])) - (1 / self.lc.tseg))
+        assert dps.dyn_ps.shape == (nfreq, nsegs)
+
+    def test_trace_maximum_without_boundaries(self):
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=UserWarning)
+            dps = DynamicalCrossspectrum(self.lc, self.lc, segment_size=3)
+        max_pos = dps.trace_maximum()
+
+        assert np.max(dps.freq[max_pos]) <= 1 / self.lc.dt
+        assert np.min(dps.freq[max_pos]) >= 1 / dps.segment_size
+
+    def test_trace_maximum_with_boundaries(self):
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=UserWarning)
+            dps = DynamicalCrossspectrum(self.lc, self.lc, segment_size=3)
+        minfreq = 21
+        maxfreq = 24
+        max_pos = dps.trace_maximum(min_freq=minfreq, max_freq=maxfreq)
+
+        assert np.max(dps.freq[max_pos]) <= maxfreq
+        assert np.min(dps.freq[max_pos]) >= minfreq
+
+    def test_size_of_trace_maximum(self):
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=UserWarning)
+            dps = DynamicalCrossspectrum(self.lc, self.lc, segment_size=3)
+        max_pos = dps.trace_maximum()
+        nsegs = int(self.lc.tseg / dps.segment_size)
+        assert len(max_pos) == nsegs
+
+    def test_rebin_small_dt(self):
+        segment_size = 3
+        dps = DynamicalCrossspectrum(self.lc_test, self.lc_test, segment_size=segment_size)
+        with pytest.raises(ValueError):
+            dps.rebin_time(dt_new=2.0)
+
+    def test_rebin_small_df(self):
+        segment_size = 3
+        dps = DynamicalCrossspectrum(self.lc, self.lc, segment_size=segment_size)
+        with pytest.raises(ValueError):
+            dps.rebin_frequency(df_new=dps.df / 2.0)
+
+    def test_rebin_time_default_method(self):
+        segment_size = 3
+        dt_new = 6.0
+        rebin_time = np.array([2.5, 8.5])
+        rebin_dps = np.array([[1.73611111, 0.81018519]])
+        dps = DynamicalCrossspectrum(self.lc_test, self.lc_test, segment_size=segment_size)
+        new_dps = dps.rebin_time(dt_new=dt_new)
+        assert np.allclose(new_dps.time, rebin_time)
+        assert np.allclose(new_dps.dyn_ps, rebin_dps)
+        assert np.isclose(new_dps.dt, dt_new)
+
+    def test_rebin_frequency_default_method(self):
+        segment_size = 50
+        df_new = 10.0
+        rebin_freq = np.array([5.01, 15.01, 25.01, 35.01])
+        rebin_dps = np.array(
+            [
+                [
+                    [1.71989342e-04, 6.42756881e-05],
+                    [7.54455204e-04, 2.14785049e-04],
+                    [6.24831554e00, 6.24984615e00],
+                    [6.71135792e-04, 7.42516599e-05],
+                ]
+            ]
+        )
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=UserWarning)
+            dps = DynamicalCrossspectrum(self.lc, self.lc, segment_size=segment_size)
+        new_dps = dps.rebin_frequency(df_new=df_new)
+        assert np.allclose(new_dps.freq, rebin_freq)
+        assert np.allclose(new_dps.dyn_ps, rebin_dps, atol=0.01)
+        assert np.isclose(new_dps.df, df_new)
+
+    @pytest.mark.parametrize("method", ["mean", "average"])
+    def test_rebin_time_mean_method(self, method):
+        segment_size = 3
+        dt_new = 6.0
+        rebin_time = np.array([2.5, 8.5])
+        rebin_dps = np.array([[1.73611111, 0.81018519]]) / 2
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=UserWarning)
+            dps = DynamicalCrossspectrum(self.lc_test, self.lc_test, segment_size=segment_size)
+        new_dps = dps.rebin_time(dt_new=dt_new, method=method)
+        assert np.allclose(new_dps.time, rebin_time)
+        assert np.allclose(new_dps.dyn_ps, rebin_dps)
+        assert np.isclose(new_dps.dt, dt_new)
+
+    @pytest.mark.parametrize("method", ["mean", "average"])
+    def test_rebin_frequency_mean_method(self, method):
+        segment_size = 50
+        df_new = 10.0
+        rebin_freq = np.array([5.01, 15.01, 25.01, 35.01])
+        rebin_dps = (
+            np.array(
+                [
+                    [
+                        [1.71989342e-04, 6.42756881e-05],
+                        [7.54455204e-04, 2.14785049e-04],
+                        [6.24831554e00, 6.24984615e00],
+                        [6.71135792e-04, 7.42516599e-05],
+                    ]
+                ]
+            )
+            / 500
+        )
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=UserWarning)
+            dps = DynamicalCrossspectrum(self.lc, self.lc, segment_size=segment_size)
+        new_dps = dps.rebin_frequency(df_new=df_new, method=method)
+        assert np.allclose(new_dps.freq, rebin_freq)
+        assert np.allclose(new_dps.dyn_ps, rebin_dps, atol=0.00001)
+        assert np.isclose(new_dps.df, df_new)

From f2ad54c27073216af36d47aa250e6b280b1c901c Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Mon, 30 Oct 2023 12:19:09 +0100
Subject: [PATCH 09/21] Update notebooks

---
 docs/notebooks | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/notebooks b/docs/notebooks
index 2e24a8966..601074930 160000
--- a/docs/notebooks
+++ b/docs/notebooks
@@ -1 +1 @@
-Subproject commit 2e24a8966070ffc4d1d9daf8b2afe28ae5cf0917
+Subproject commit 601074930f8f32899f0816f5bab4624c4d12592d

From f28c8cedeb155b261d18f193263fa514f2a43124 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Tue, 14 Nov 2023 16:48:03 +0100
Subject: [PATCH 10/21] Add doc on dt

---
 stingray/crossspectrum.py | 10 ++++++++--
 stingray/powerspectrum.py |  7 ++++++-
 2 files changed, 14 insertions(+), 3 deletions(-)

diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index b1db2094b..d23db8443 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -1899,11 +1899,12 @@ class DynamicalCrossspectrum(AveragedCrossspectrum):
     ----------
     data1 : :class:`stingray.Lightcurve` or :class:`stingray.EventList` object
         The time series or event list from the subject band, or channel, for which
-        the dynamical cross spectrum is to be calculated.
+        the dynamical cross spectrum is to be calculated. If :class:`stingray.EventList`, ``dt``
+        must be specified as well.
 
     data2 : :class:`stingray.Lightcurve` or :class:`stingray.EventList` object
         The time series or event list from the reference band, or channel, of the dynamical
-        cross spectrum.
+        cross spectrum. If :class:`stingray.EventList`, ``dt`` must be specified as well.
 
     segment_size : float, default 1
          Length of the segment of light curve, default value is 1 (in whatever
@@ -1922,6 +1923,11 @@ class DynamicalCrossspectrum(AveragedCrossspectrum):
         use this and only give GTIs to the input object before making
         the cross spectrum.
 
+    dt: float
+        Compulsory for input :class:`stingray.EventList` data. The time resolution of the
+        lightcurve that is  created internally from the input event lists. Drives the
+        Nyquist frequency.
+
     Attributes
     ----------
     segment_size: float
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 32ff8fa69..91476a16c 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -900,7 +900,7 @@ class DynamicalPowerspectrum(DynamicalCrossspectrum):
     ----------
     lc : :class:`stingray.Lightcurve` or :class:`stingray.EventList` object
         The time series or event list of which the dynamical power spectrum is
-        to be calculated.
+        to be calculated. If :class:`stingray.EventList`, ``dt`` must be specified as well.
 
     segment_size : float, default 1
          Length of the segment of light curve, default value is 1 (in whatever
@@ -919,6 +919,11 @@ class DynamicalPowerspectrum(DynamicalCrossspectrum):
         use this and only give GTIs to the input object before making
         the power spectrum.
 
+    dt: float
+        Compulsory for input :class:`stingray.EventList` data. The time resolution of the
+        lightcurve that is created internally from the input event lists. Drives the
+        Nyquist frequency.
+
     Attributes
     ----------
     segment_size: float

From d32be51e9657b07a66b39b2e3e79429dbe783329 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Tue, 14 Nov 2023 16:56:44 +0100
Subject: [PATCH 11/21] Fix more docstrings [docs only]

---
 stingray/crossspectrum.py | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index d23db8443..7992c925f 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -1925,7 +1925,7 @@ class DynamicalCrossspectrum(AveragedCrossspectrum):
 
     dt: float
         Compulsory for input :class:`stingray.EventList` data. The time resolution of the
-        lightcurve that is  created internally from the input event lists. Drives the
+        lightcurve that is created internally from the input event lists. Drives the
         Nyquist frequency.
 
     Attributes
@@ -2014,9 +2014,9 @@ def rebin_frequency(self, df_new, method="sum"):
         Rebinning is an in-place operation, i.e. will replace the existing
         ``dyn_ps`` attribute.
 
-        While the new resolution need not be an integer multiple of the
-        previous frequency resolution, be aware that if it is not, the last
-        bin will be cut off by the fraction left over by the integer division.
+        While the new resolution does not need to be an integer of the previous frequency
+        resolution, be aware that if this is the case, the last frequency bin will be cut
+        off by the fraction left over by the integer division
 
         Parameters
         ----------
@@ -2045,9 +2045,10 @@ def rebin_frequency(self, df_new, method="sum"):
     def rebin_time(self, dt_new, method="sum"):
         """
         Rebin the Dynamic Power Spectrum to a new time resolution.
-        While the new resolution need not be an integer multiple of the
-        previous time resolution, be aware that if it is not, the last bin
-        will be cut off by the fraction left over by the integer division.
+
+        While the new resolution does not need to be an integer of the previous time
+        resolution, be aware that if this is the case, the last frequency bin will be cut
+        off by the fraction left over by the integer division
 
         Parameters
         ----------
@@ -2066,7 +2067,7 @@ def rebin_time(self, dt_new, method="sum"):
             Time axis with new rebinned time resolution.
 
         dynspec_new: numpy.ndarray
-            New rebinned Dynamical Power Spectrum.
+            New rebinned Dynamical Cross Spectrum.
         """
         if dt_new < self.dt:
             raise ValueError("New time resolution must be larger than old time resolution!")

From 6791e80732de0916be7ecd2d67c25a6a5a51238d Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Tue, 14 Nov 2023 17:11:23 +0100
Subject: [PATCH 12/21] let's call dt sample_time

---
 stingray/crossspectrum.py            | 23 ++++++++++++++---------
 stingray/powerspectrum.py            | 16 ++++++++--------
 stingray/tests/test_crossspectrum.py |  2 +-
 stingray/tests/test_powerspectrum.py |  2 +-
 4 files changed, 24 insertions(+), 19 deletions(-)

diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 7992c925f..13398efb0 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -1923,7 +1923,7 @@ class DynamicalCrossspectrum(AveragedCrossspectrum):
         use this and only give GTIs to the input object before making
         the cross spectrum.
 
-    dt: float
+    sample_time: float
         Compulsory for input :class:`stingray.EventList` data. The time resolution of the
         lightcurve that is created internally from the input event lists. Drives the
         Nyquist frequency.
@@ -1947,28 +1947,33 @@ class DynamicalCrossspectrum(AveragedCrossspectrum):
     freq: numpy.ndarray
         The array of mid-bin frequencies that the Fourier transform samples.
 
+    time: numpy.ndarray
+        The array of mid-point times of each interval used for the dynamical
+        power spectrum.
+
     df: float
         The frequency resolution.
 
     dt: float
-        The time resolution.
+        The time resolution of the dynamical spectrum. It is **not** the time resolution of the
+        input light curve.
     """
 
-    def __init__(self, data1, data2, segment_size, norm="frac", gti=None, dt=None):
-        if isinstance(data1, EventList) and dt is None:
-            raise ValueError("To pass input event lists, please specify dt")
+    def __init__(self, data1, data2, segment_size, norm="frac", gti=None, sample_time=None):
+        if isinstance(data1, EventList) and sample_time is None:
+            raise ValueError("To pass input event lists, please specify sample_time")
         elif isinstance(data1, Lightcurve):
-            dt = data1.dt
+            sample_time = data1.dt
             if segment_size > data1.tseg:
                 raise ValueError(
                     "Length of the segment is too long to create "
                     "any segments of the light curve!"
                 )
-        if segment_size < 2 * dt:
+        if segment_size < 2 * sample_time:
             raise ValueError("Length of the segment is too short to form a light curve!")
 
         self.segment_size = segment_size
-        self.input_dt = dt
+        self.sample_time = sample_time
         self.gti = gti
         self.norm = norm
 
@@ -1992,7 +1997,7 @@ def _make_matrix(self, data1, data2):
         avg = AveragedCrossspectrum(
             data1,
             data2,
-            dt=self.input_dt,
+            dt=self.sample_time,
             segment_size=self.segment_size,
             norm=self.norm,
             gti=self.gti,
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 91476a16c..71cc77d5a 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -919,7 +919,7 @@ class DynamicalPowerspectrum(DynamicalCrossspectrum):
         use this and only give GTIs to the input object before making
         the power spectrum.
 
-    dt: float
+    sample_time: float
         Compulsory for input :class:`stingray.EventList` data. The time resolution of the
         lightcurve that is created internally from the input event lists. Drives the
         Nyquist frequency.
@@ -950,21 +950,21 @@ class DynamicalPowerspectrum(DynamicalCrossspectrum):
         The time resolution.
     """
 
-    def __init__(self, lc, segment_size, norm="frac", gti=None, dt=None):
-        if isinstance(lc, EventList) and dt is None:
-            raise ValueError("To pass an input event lists, please specify dt")
+    def __init__(self, lc, segment_size, norm="frac", gti=None, sample_time=None):
+        if isinstance(lc, EventList) and sample_time is None:
+            raise ValueError("To pass an input event lists, please specify sample_time")
         elif isinstance(lc, Lightcurve):
-            dt = lc.dt
+            sample_time = lc.dt
             if segment_size > lc.tseg:
                 raise ValueError(
                     "Length of the segment is too long to create "
                     "any segments of the light curve!"
                 )
-        if segment_size < 2 * dt:
+        if segment_size < 2 * sample_time:
             raise ValueError("Length of the segment is too short to form a light curve!")
 
         self.segment_size = segment_size
-        self.input_dt = dt
+        self.sample_time = sample_time
         self.gti = gti
         self.norm = norm
 
@@ -984,7 +984,7 @@ def _make_matrix(self, lc):
         """
         avg = AveragedPowerspectrum(
             lc,
-            dt=self.input_dt,
+            dt=self.sample_time,
             segment_size=self.segment_size,
             norm=self.norm,
             gti=self.gti,
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index 0b95cc6b0..b45d31f3b 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -1354,7 +1354,7 @@ def test_works_with_events(self):
             # Without dt, it fails
             _ = DynamicalCrossspectrum(ev, ev, segment_size=10)
 
-        dps_ev = DynamicalCrossspectrum(ev, ev, segment_size=10, dt=self.lc.dt)
+        dps_ev = DynamicalCrossspectrum(ev, ev, segment_size=10, sample_time=self.lc.dt)
         assert np.allclose(dps.dyn_ps, dps_ev.dyn_ps)
 
     def test_with_long_seg_size(self):
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index f1b5ace85..1a4cff650 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -997,7 +997,7 @@ def test_works_with_events(self):
             # Without dt, it fails
             _ = DynamicalPowerspectrum(ev, segment_size=10)
 
-        dps_ev = DynamicalPowerspectrum(ev, segment_size=10, dt=self.lc.dt)
+        dps_ev = DynamicalPowerspectrum(ev, segment_size=10, sample_time=self.lc.dt)
         assert np.allclose(dps.dyn_ps, dps_ev.dyn_ps)
 
     def test_with_long_seg_size(self):

From eef15e68d6b17a24b838b1fdf67fae336e760e51 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Tue, 14 Nov 2023 17:15:27 +0100
Subject: [PATCH 13/21] Fix docstring [docs only]

---
 stingray/powerspectrum.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 71cc77d5a..ef82706f8 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -943,6 +943,10 @@ class DynamicalPowerspectrum(DynamicalCrossspectrum):
     freq: numpy.ndarray
         The array of mid-bin frequencies that the Fourier transform samples.
 
+    time: numpy.ndarray
+        The array of mid-point times of each interval used for the dynamical
+        power spectrum.
+
     df: float
         The frequency resolution.
 

From e6ef2efa42a4a614dc42c10e39a7efc9721be338 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Tue, 14 Nov 2023 17:40:43 +0100
Subject: [PATCH 14/21] Ignore more jax deprecation from dependencies

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index d9c274d6d..e49f43f66 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -107,7 +107,7 @@ filterwarnings =
     ignore:.*is a deprecated alias for:DeprecationWarning
     ignore:.*HIERARCH card will be created.*:
     ignore:.*FigureCanvasAgg is non-interactive.*:UserWarning
-    ignore:.*jax.* deprecated. Use jax.*instead:DeprecationWarning
+    ignore:.*jax.* deprecated:DeprecationWarning:haiku
 
 ;addopts = --disable-warnings
 

From 4004a4d327b82557c0aa22f8651dcc684e8bcd9c Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Tue, 14 Nov 2023 17:50:17 +0100
Subject: [PATCH 15/21] Ignore more jax deprecation from dependencies

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index e49f43f66..342c39168 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -107,7 +107,7 @@ filterwarnings =
     ignore:.*is a deprecated alias for:DeprecationWarning
     ignore:.*HIERARCH card will be created.*:
     ignore:.*FigureCanvasAgg is non-interactive.*:UserWarning
-    ignore:.*jax.* deprecated:DeprecationWarning:haiku
+    ignore:.*jax.* deprecated:DeprecationWarning:
 
 ;addopts = --disable-warnings
 

From a5cebc920d83ab2534f0740c0eb0c66e57844106 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Wed, 15 Nov 2023 16:04:02 +0100
Subject: [PATCH 16/21] Reiterate that dt is not lc.dt

---
 stingray/crossspectrum.py | 7 +++++--
 stingray/powerspectrum.py | 4 +++-
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 13398efb0..6cfbca49c 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -1956,7 +1956,8 @@ class DynamicalCrossspectrum(AveragedCrossspectrum):
 
     dt: float
         The time resolution of the dynamical spectrum. It is **not** the time resolution of the
-        input light curve.
+        input light curve. It is the integration time of each line of the dynamical power
+        spectrum (typically, an integer multiple of ``segment_size``).
     """
 
     def __init__(self, data1, data2, segment_size, norm="frac", gti=None, sample_time=None):
@@ -2060,7 +2061,9 @@ def rebin_time(self, dt_new, method="sum"):
         dt_new: float
             The new time resolution of the dynamical power spectrum.
             Must be larger than the time resolution of the old dynamical power
-            spectrum!
+            spectrum! Note: this 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``).
 
         method: {"sum" | "mean" | "average"}, optional, default "sum"
             This keyword argument sets whether the counts in the new bins
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index ef82706f8..6fb037bf1 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -951,7 +951,9 @@ class DynamicalPowerspectrum(DynamicalCrossspectrum):
         The frequency resolution.
 
     dt: float
-        The time resolution.
+        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``).
     """
 
     def __init__(self, lc, segment_size, norm="frac", gti=None, sample_time=None):

From e9b207f60665e35b75f80ea02d4d12fb7f7bea01 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Wed, 15 Nov 2023 16:08:43 +0100
Subject: [PATCH 17/21] Fix docstring

---
 stingray/crossspectrum.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 6cfbca49c..4a2f3021c 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -2053,7 +2053,7 @@ def rebin_time(self, dt_new, method="sum"):
         Rebin the Dynamic Power Spectrum to a new time resolution.
 
         While the new resolution does not need to be an integer of the previous time
-        resolution, be aware that if this is the case, the last frequency bin will be cut
+        resolution, be aware that if this is the case, the last time bin will be cut
         off by the fraction left over by the integer division
 
         Parameters

From 4a0d012046fdd447cebbd61d83864b92a475cb35 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Wed, 15 Nov 2023 18:02:52 +0100
Subject: [PATCH 18/21] Make message about time resolution more prominent

---
 stingray/crossspectrum.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 4a2f3021c..c87e83629 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -2052,6 +2052,10 @@ def rebin_time(self, dt_new, method="sum"):
         """
         Rebin the Dynamic Power Spectrum to a new time resolution.
 
+        Note: this 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``).
+
         While the new resolution does not need to be an integer of the previous time
         resolution, be aware that if this is the case, the last time bin will be cut
         off by the fraction left over by the integer division
@@ -2061,9 +2065,7 @@ def rebin_time(self, dt_new, method="sum"):
         dt_new: float
             The new time resolution of the dynamical power spectrum.
             Must be larger than the time resolution of the old dynamical power
-            spectrum! Note: this 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``).
+            spectrum!
 
         method: {"sum" | "mean" | "average"}, optional, default "sum"
             This keyword argument sets whether the counts in the new bins

From 78f59f36ef5da4070e8399d810c30fe6a71d6c45 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Wed, 15 Nov 2023 21:46:25 +0100
Subject: [PATCH 19/21] Make additional fix in docstring

---
 stingray/crossspectrum.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index c87e83629..a0cc3d77c 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -2052,8 +2052,8 @@ def rebin_time(self, dt_new, method="sum"):
         """
         Rebin the Dynamic Power Spectrum to a new time resolution.
 
-        Note: this is *not* the time resolution of the input light
-        curve! It is the integration time of each line of the dynamical power
+        Note: this is *not* changing the time resolution of the input light
+        curve! ``dt`` is the integration time of each line of the dynamical power
         spectrum (typically, an integer multiple of ``segment_size``).
 
         While the new resolution does not need to be an integer of the previous time

From e3d2bf30a672de530bef3e6f78bf7e0926537020 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Wed, 15 Nov 2023 21:51:02 +0100
Subject: [PATCH 20/21] Add additional test to show that dps is complex

---
 stingray/tests/test_crossspectrum.py | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index b45d31f3b..0a1aad166 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -1357,6 +1357,20 @@ def test_works_with_events(self):
         dps_ev = DynamicalCrossspectrum(ev, ev, segment_size=10, sample_time=self.lc.dt)
         assert np.allclose(dps.dyn_ps, dps_ev.dyn_ps)
 
+    def test_works_with_events_and_its_complex(self):
+        lc = copy.deepcopy(self.lc)
+        lc.counts = np.floor(lc.counts)
+        ev1 = EventList()
+        ev1.simulate_times(lc)
+        ev2 = EventList()
+        ev2.simulate_times(lc)
+
+        dps_ev = DynamicalCrossspectrum(ev1, ev1, segment_size=10, sample_time=self.lc.dt)
+        assert np.iscomplexobj(dps_ev.dyn_ps)
+        assert np.any(dps_ev.dyn_ps.imag > np.dps_ev.dyn_ps.real)
+        assert np.any(dps_ev.dyn_ps.imag < 0)
+        assert np.any(dps_ev.dyn_ps.imag > 0)
+
     def test_with_long_seg_size(self):
         with pytest.raises(ValueError):
             dps = DynamicalCrossspectrum(self.lc, self.lc, segment_size=1000)

From 36c5c0c84e7832a785ea250b354d3e72f6e36dba Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Wed, 15 Nov 2023 21:52:53 +0100
Subject: [PATCH 21/21] Fix test

---
 stingray/tests/test_crossspectrum.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index 0a1aad166..d466b2061 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -1359,15 +1359,15 @@ def test_works_with_events(self):
 
     def test_works_with_events_and_its_complex(self):
         lc = copy.deepcopy(self.lc)
-        lc.counts = np.floor(lc.counts)
+        lc.counts = np.random.poisson(10, size=lc.counts.size)
         ev1 = EventList()
         ev1.simulate_times(lc)
         ev2 = EventList()
         ev2.simulate_times(lc)
 
-        dps_ev = DynamicalCrossspectrum(ev1, ev1, segment_size=10, sample_time=self.lc.dt)
+        dps_ev = DynamicalCrossspectrum(ev1, ev2, segment_size=10, sample_time=self.lc.dt)
         assert np.iscomplexobj(dps_ev.dyn_ps)
-        assert np.any(dps_ev.dyn_ps.imag > np.dps_ev.dyn_ps.real)
+        assert np.any(dps_ev.dyn_ps.imag > dps_ev.dyn_ps.real)
         assert np.any(dps_ev.dyn_ps.imag < 0)
         assert np.any(dps_ev.dyn_ps.imag > 0)