Repository: StingraySoftware/stingray
Pull request #779 (2023-10-26 13:28:23)

Fixes to dynamical power spectrum + Dynamical Cross spectrum
+ Improves performance and consistency when using events (uses the internals already tested in `AveragePowerspectrum` to create light curves dynamically for each segment, rather than allocating the full light curve on the spot)
+ Resolves a bug when rebinning in time (it relied on dt, which was calculated in a weird manner)
+ Makes tests a little more understandable, I hope

BONUS: added Dynamical Cross spectrum, because why not


Patch:
diff --git a/docs/changes/779.bugfix.rst b/docs/changes/779.bugfix.rst
new file mode 100644
index 00000000..fd209edc
--- /dev/null
+++ b/docs/changes/779.bugfix.rst
@@ -0,0 +1,1 @@
+Various bug fixes in DynamicalPowerspectrum, on event loading and time rebinning
diff --git a/docs/notebooks b/docs/notebooks
index 2e24a896..60107493 160000
--- a/docs/notebooks
+++ b/docs/notebooks
@@ -1,1 +1,1 @@
-Subproject commit 2e24a8966070ffc4d1d9daf8b2afe28ae5cf0917
+Subproject commit 601074930f8f32899f0816f5bab4624c4d12592d
diff --git a/setup.cfg b/setup.cfg
index d9c274d6..342c3916 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:
 
 ;addopts = --disable-warnings
 
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 8a22b26b..a0cc3d77 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -15,6 +15,7 @@ from stingray.utils import rebin_data, rebin_data_log, simon
 
 from .base import StingrayObject
 from .events import EventList
+from .gti import cross_two_gtis, time_intervals_from_gtis
 from .lightcurve import Lightcurve
 from .utils import show_progress
 from .fourier import avg_cs_from_iterables, error_on_averaged_cross_spectrum
@@ -28,6 +29,7 @@ from scipy.special import factorial
 __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,253 @@ class AveragedCrossspectrum(Crossspectrum):
         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. 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. 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
+         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.
+
+    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.
+
+    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.
+
+    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 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, 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):
+            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 * sample_time:
+            raise ValueError("Length of the segment is too short to form a light curve!")
+
+        self.segment_size = segment_size
+        self.sample_time = sample_time
+        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.sample_time,
+            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 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
+        ----------
+        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.
+
+        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
+        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
+        ----------
+        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 Cross 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 da79909f..6fb037bf 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -7,17 +7,15 @@ import scipy
 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
-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
 
@@ -885,7 +883,7 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
         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*.
@@ -902,7 +900,7 @@ class DynamicalPowerspectrum(AveragedPowerspectrum):
     ----------
     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
@@ -921,6 +919,11 @@ class DynamicalPowerspectrum(AveragedPowerspectrum):
         use this and only give GTIs to the input object before making
         the power spectrum.
 
+    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.
+
     Attributes
     ----------
     segment_size: float
@@ -940,29 +943,37 @@ class DynamicalPowerspectrum(AveragedPowerspectrum):
     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. 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, 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):
+            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 * sample_time:
+            raise ValueError("Length of the segment is too short to form a light curve!")
 
-        if isinstance(lc, EventList):
-            lc = lc.to_lc(dt)
+        self.segment_size = segment_size
+        self.sample_time = sample_time
+        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,142 +989,22 @@ class DynamicalPowerspectrum(AveragedPowerspectrum):
             power spectrum.
         """
         avg = AveragedPowerspectrum(
-            lc, segment_size=self.segment_size, norm=self.norm, gti=self.gti, save_all=True
+            lc,
+            dt=self.sample_time,
+            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
-        )
-
-        tstart = lc.time[start_inds]
-        tend = lc.time[end_inds]
-
-        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
-
-    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
-            )
-            dynspec_new.append(bin_counts)
+        tstart, _ = time_intervals_from_gtis(current_gti, self.segment_size)
 
-        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
+        self.time = tstart + 0.5 * self.segment_size
+        self.df = avg.df
+        self.dt = self.segment_size
 
 
 def powerspectrum_from_time_array(


Test patch:
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index 32021b0d..d466b206 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -6,7 +6,7 @@ import matplotlib.pyplot as plt
 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,181 @@ class TestRoundTrip:
         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, 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.random.poisson(10, size=lc.counts.size)
+        ev1 = EventList()
+        ev1.simulate_times(lc)
+        ev2 = EventList()
+        ev2.simulate_times(lc)
+
+        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 > 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)
+
+    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)
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index 5b44516d..1a4cff65 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -151,7 +151,6 @@ class TestAveragedPowerspectrumEvents(object):
 
         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 @@ class TestAveragedPowerspectrum(object):
 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 @@ class TestDynamicalPowerspectrum(object):
         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
@@ -999,7 +997,7 @@ class TestDynamicalPowerspectrum(object):
             # 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):
@@ -1044,43 +1042,39 @@ class TestDynamicalPowerspectrum(object):
 
     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 @@ class TestDynamicalPowerspectrum(object):
         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)