Repository: StingraySoftware/stingray
Pull request #193 (2017-05-23 08:14:02)

[WIP] Dynamical Powerspectrum
This commit picks up the changes of #132 and try to implement the suggestions of @dhuppenkothen 

Please, review away the current implementation.

Showcase notebook: https://gist.github.com/evandromr/176d6db06acfc76ab46261d035695113

To do:
-------
 

- [x] tests
- [x] rebin in frequency
- [x] rebin in time
- [ ] other tracing methods?



Patch:
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index bd8b7a61..cedf7b15 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -12,8 +12,7 @@ from stingray.gti import bin_intervals_from_gtis, check_gtis
 from stingray.utils import simon
 from stingray.crossspectrum import Crossspectrum, AveragedCrossspectrum
 
-
-__all__ = ["Powerspectrum", "AveragedPowerspectrum"]
+__all__ = ["Powerspectrum", "AveragedPowerspectrum", "DynamicalPowerspectrum"]
 
 
 def classical_pvalue(power, nspec):
@@ -79,7 +78,7 @@ def classical_pvalue(power, nspec):
 
     # If the power is really big, it's safe to say it's significant,
     # and the p-value will be nearly zero
-    if (power*nspec) > 30000:
+    if (power * nspec) > 30000:
         simon("Probability of no signal too miniscule to calculate.")
         return 0.0
 
@@ -100,10 +99,10 @@ def _pavnosigfun(power, nspec):
     while m >= 0:
 
         s = 0.0
-        for i in range(int(m)-1):
-            s += np.log(float(m-i))
+        for i in range(int(m) - 1):
+            s += np.log(float(m - i))
 
-        logterm = m*np.log(pn/2) - pn/2 - s
+        logterm = m * np.log(pn / 2) - pn / 2 - s
         term = np.exp(logterm)
         ratio = sum / term
 
@@ -117,7 +116,6 @@ def _pavnosigfun(power, nspec):
 
 
 class Powerspectrum(Crossspectrum):
-
     def __init__(self, lc=None, norm='frac', gti=None):
         """
         Make a Periodogram (power spectrum) from a (binned) light curve.
@@ -131,9 +129,9 @@ class Powerspectrum(Crossspectrum):
         lc: lightcurve.Lightcurve object, optional, default None
             The light curve data to be Fourier-transformed.
 
-        norm: {"leahy" | "rms"}, optional, default "rms"
+        norm: {"leahy" | "frac" | "abs" | "none" }, optional, default "frac"
             The normaliation of the periodogram to be used. Options are
-            "leahy" or "rms", default is "rms".
+            "leahy", "frac", "abs" and "none", default is "frac".
 
         Other Parameters
         ----------------
@@ -144,7 +142,7 @@ class Powerspectrum(Crossspectrum):
 
         Attributes
         ----------
-        norm: {"leahy" | "rms"}
+        norm: {"leahy" | "frac" | "abs" | "none"}
             the normalization of the periodogram
 
         freq: numpy.ndarray
@@ -209,10 +207,10 @@ class Powerspectrum(Crossspectrum):
         powers = self.power[minind:maxind]
 
         if self.norm.lower() == 'leahy':
-            rms = np.sqrt(np.sum(powers)/self.nphots)
+            rms = np.sqrt(np.sum(powers) / self.nphots)
 
         elif self.norm.lower() == "frac":
-            rms = np.sqrt(np.sum(powers*self.df))
+            rms = np.sqrt(np.sum(powers * self.df))
         else:
             raise TypeError("Normalization not recognized!")
 
@@ -238,9 +236,9 @@ class Powerspectrum(Crossspectrum):
         delta_rms: float
             the error on the fractional rms amplitude
         """
-        p_err = scipy.stats.chi2(2.0*self.m).var() * powers / self.m
-        drms_dp = 1 / (2*np.sqrt(np.sum(powers)*self.df))
-        delta_rms = np.sum(p_err*drms_dp*self.df)
+        p_err = scipy.stats.chi2(2.0 * self.m).var() * powers / self.m
+        drms_dp = 1 / (2 * np.sqrt(np.sum(powers) * self.df))
+        delta_rms = np.sum(p_err * drms_dp * self.df)
         return delta_rms
 
     def classical_significances(self, threshold=1, trial_correction=False):
@@ -298,7 +296,7 @@ class Powerspectrum(Crossspectrum):
         # calculate p-values for all powers
         # leave out zeroth power since it just encodes the number of photons!
         pv = np.array([classical_pvalue(power, self.m)
-                      for power in self.power])
+                       for power in self.power])
 
         # if trial correction is used, then correct the threshold for
         # the number of powers in the power spectrum
@@ -315,7 +313,6 @@ class Powerspectrum(Crossspectrum):
 
 
 class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
-
     def __init__(self, lc=None, segment_size=None, norm="frac", gti=None):
         """
         Make an averaged periodogram from a light curve by segmenting the light
@@ -334,9 +331,9 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
             of the segment_size, then any fraction left-over at the end of the
             time series will be lost.
 
-        norm: {"leahy" | "rms"}, optional, default "rms"
+        norm: {"leahy" | "frac" | "abs" | "none" }, optional, default "frac"
             The normaliation of the periodogram to be used. Options are
-            "leahy" or "rms", default is "rms".
+            "leahy", "frac", "abs" and "none", default is "frac".
 
 
         Other Parameters
@@ -348,7 +345,7 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
 
         Attributes
         ----------
-        norm: {"leahy" | "rms"}
+        norm: {"leahy" | "frac" | "abs" | "none"}
             the normalization of the periodogram
 
         freq: numpy.ndarray
@@ -417,3 +414,194 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
             nphots_all.append(np.sum(lc_seg.counts))
 
         return power_all, nphots_all
+
+
+class DynamicalPowerspectrum(AveragedPowerspectrum):
+    def __init__(self, lc, segment_size, norm="frac", gti=None):
+        """
+        Parameters
+        ----------
+        lc : lightcurve.Lightcurve object
+            The time series of which the Dynamical powerspectrum is
+            to be calculated.
+
+        segment_size : float, default 1
+             Length of the segment of light curve, default value is 1 second.
+
+        norm: {"leahy" | "frac" | "abs" | "none" }, optional, default "frac"
+            The normaliation of the periodogram to be used. Options are
+            "leahy", "frac", "abs" and "none", default is "frac".
+
+        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!
+
+        Attributes
+        ----------
+        segment_size: float
+            The size of each segment to average. Note that if the total
+            duration of each Lightcurve 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` attirbutes
+
+        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
+        """
+        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!")
+        self.segment_size = segment_size
+        self.norm = norm
+        self.gti = gti
+        self._make_matrix(lc)
+
+    def _make_matrix(self, lc):
+        """Create the matrix with freq and time columns."""
+        ps_all, _ = AveragedPowerspectrum._make_segment_spectrum(
+            self, lc, self.segment_size)
+        self.dyn_ps = np.array([ps.power for ps in ps_all]).T
+
+        self.freq = ps_all[0].freq
+        self.time = np.arange(lc.time[0] - 0.5 * lc.dt + 0.5 * self.segment_size,
+                              lc.time[-1] + 0.5 * lc.dt, self.segment_size)
+
+        # Assign lenght 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
+
+        if len(self.time) > self.dyn_ps.shape[0]:
+            self.time = self.time[:-1]
+
+    def rebin_frequency(self, df_new, method="sum"):
+        """
+        Rebin the Dynamic Power Spectrum to a new frequency resolution.
+        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 Dynamica 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.
+        """
+        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)
+
+        self.freq = freq_new
+        self.dyn_ps = np.array(dynspec_new).T
+        self.df = df_new
+
+    def trace_maximum(self, min_freq=None, max_freq=None, sigmaclip=False):
+        """
+        Return the indices of the maximum powers in each segment 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 Dynamica 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!")
+
+        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)
+
+        self.time = time_new
+        self.dyn_ps = np.array(dynspec_new)
+        self.dt = dt_new


Test patch:
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index f035266a..41fdbdc5 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -1,17 +1,16 @@
-
 import numpy as np
 
 from astropy.tests.helper import pytest
 
 from stingray import Lightcurve
-from stingray import Powerspectrum, AveragedPowerspectrum
+from stingray import Powerspectrum, AveragedPowerspectrum, \
+    DynamicalPowerspectrum
 from stingray.powerspectrum import classical_pvalue
 
 np.random.seed(20150907)
 
 
 class TestPowerspectrum(object):
-
     @classmethod
     def setup_class(cls):
         tstart = 0.0
@@ -88,9 +87,9 @@ class TestPowerspectrum(object):
         """
         ps = Powerspectrum(lc=self.lc)
         nn = ps.n
-        pp = ps.unnorm_power / np.float(nn)**2
-        p_int = np.sum(pp[:-1]*ps.df) + (pp[-1]*ps.df)/2
-        var_lc = np.var(self.lc.counts) / (2.*self.lc.tseg)
+        pp = ps.unnorm_power / np.float(nn) ** 2
+        p_int = np.sum(pp[:-1] * ps.df) + (pp[-1] * ps.df) / 2
+        var_lc = np.var(self.lc.counts) / (2. * self.lc.tseg)
         assert np.isclose(p_int, var_lc, atol=0.01, rtol=0.01)
 
     def test_frac_normalization_is_standard(self):
@@ -108,8 +107,8 @@ class TestPowerspectrum(object):
         of the light curve squared.
         """
         ps = Powerspectrum(lc=self.lc, norm="frac")
-        ps_int = np.sum(ps.power[:-1]*ps.df) + ps.power[-1]*ps.df/2
-        std_lc = np.var(self.lc.counts) / np.mean(self.lc.counts)**2
+        ps_int = np.sum(ps.power[:-1] * ps.df) + ps.power[-1] * ps.df / 2
+        std_lc = np.var(self.lc.counts) / np.mean(self.lc.counts) ** 2
         print(ps_int)
         print(std_lc)
         assert np.isclose(ps_int, std_lc, atol=0.01, rtol=0.01)
@@ -138,8 +137,8 @@ class TestPowerspectrum(object):
         square of the number of data points in the light curve
         """
         ps = Powerspectrum(lc=self.lc, norm="Leahy")
-        ps_var = (np.sum(self.lc.counts)/ps.n**2.) * \
-            (np.sum(ps.power[:-1]) + ps.power[-1]/2.)
+        ps_var = (np.sum(self.lc.counts) / ps.n ** 2.) * \
+                 (np.sum(ps.power[:-1]) + ps.power[-1] / 2.)
 
         assert np.isclose(ps_var, np.var(self.lc.counts), atol=0.01)
 
@@ -173,6 +172,7 @@ class TestPowerspectrum(object):
         ps = Powerspectrum(lc=self.lc, norm="Leahy")
         # replace powers
         ps.power = np.ones_like(ps.power) * 2.0
+
         rebin_factor = 2
         bin_ps = ps.rebin(rebin_factor*ps.df)
 
@@ -202,9 +202,9 @@ class TestPowerspectrum(object):
         """
         ps = Powerspectrum(lc=self.lc, norm="Leahy")
         bin_ps = ps.rebin(df)
-        assert np.isclose(bin_ps.freq[1]-bin_ps.freq[0], bin_ps.df,
+        assert np.isclose(bin_ps.freq[1] - bin_ps.freq[0], bin_ps.df,
                           atol=1e-4, rtol=1e-4)
-        assert np.isclose(bin_ps.freq[0], (ps.freq[0]-ps.df*0.5+bin_ps.df*0.5),
+        assert np.isclose(bin_ps.freq[0], (ps.freq[0] - ps.df * 0.5 + bin_ps.df * 0.5),
                           atol=1e-4, rtol=1e-4)
 
     def test_classical_significances_runs(self):
@@ -220,7 +220,7 @@ class TestPowerspectrum(object):
         ps = Powerspectrum(lc=self.lc, norm="leahy")
 
         # change the powers so that just one exceeds the threshold
-        ps.power = np.zeros_like(ps.power)+2.0
+        ps.power = np.zeros_like(ps.power) + 2.0
 
         index = 1
         ps.power[index] = 10.0
@@ -246,7 +246,7 @@ class TestPowerspectrum(object):
     def test_pvals_is_numpy_array(self):
         ps = Powerspectrum(lc=self.lc, norm="leahy")
         # change the powers so that just one exceeds the threshold
-        ps.power = np.zeros_like(ps.power)+2.0
+        ps.power = np.zeros_like(ps.power) + 2.0
 
         index = 1
         ps.power[index] = 10.0
@@ -261,7 +261,6 @@ class TestPowerspectrum(object):
 
 
 class TestAveragedPowerspectrum(object):
-
     @classmethod
     def setup_class(cls):
         tstart = 0.0
@@ -271,7 +270,7 @@ class TestAveragedPowerspectrum(object):
         time = np.arange(tstart + 0.5*dt, tend + 0.5*dt, dt)
 
         mean_count_rate = 1000.0
-        mean_counts = mean_count_rate*dt
+        mean_counts = mean_count_rate * dt
 
         poisson_counts = np.random.poisson(mean_counts,
                                            size=time.shape[0])
@@ -292,7 +291,7 @@ class TestAveragedPowerspectrum(object):
         assert ps.m == nseg
 
     def test_segments_with_leftover(self):
-        segment_size = self.lc.tseg/2. - 1.
+        segment_size = self.lc.tseg / 2. - 1.
         ps = AveragedPowerspectrum(self.lc, segment_size)
         assert np.isclose(ps.segment_size, segment_size)
         assert ps.m == 2
@@ -331,7 +330,7 @@ class TestAveragedPowerspectrum(object):
         time = np.arange(tstart + 0.5*dt, tend + 0.5*dt, dt)
 
         mean_count_rate = 1000.0
-        mean_counts = mean_count_rate*dt
+        mean_counts = mean_count_rate * dt
 
         lc_all = []
         for n in range(n_lcs):
@@ -353,10 +352,10 @@ class TestAveragedPowerspectrum(object):
             aps = AveragedPowerspectrum(lc=self.lc, segment_size=1,
                                         norm="Leahy")
             bin_aps = aps.rebin(df)
-            assert np.isclose(bin_aps.freq[1]-bin_aps.freq[0], bin_aps.df,
+            assert np.isclose(bin_aps.freq[1] - bin_aps.freq[0], bin_aps.df,
                               atol=1e-4, rtol=1e-4)
             assert np.isclose(bin_aps.freq[0],
-                              (aps.freq[0]-aps.df*0.5+bin_aps.df*0.5),
+                              (aps.freq[0] - aps.df * 0.5 + bin_aps.df * 0.5),
                               atol=1e-4, rtol=1e-4)
 
         def test_rebin_averagedps(self):
@@ -379,10 +378,10 @@ class TestAveragedPowerspectrum(object):
         tend = 1.0
         dt = 0.0001
 
-        time = np.linspace(tstart, tend, int((tend-tstart)/dt))
+        time = np.linspace(tstart, tend, int((tend - tstart) / dt))
 
         mean_count_rate = 1000.0
-        mean_counts = mean_count_rate*dt
+        mean_counts = mean_count_rate * dt
 
         lc_all = []
         for n in range(n_lcs):
@@ -403,7 +402,7 @@ class TestAveragedPowerspectrum(object):
         n = 100
         lc_all = []
         for i in range(n):
-            time = np.arange(0.0, 10.0, 10./100000)
+            time = np.arange(0.0, 10.0, 10. / 100000)
             counts = np.random.poisson(1000, size=time.shape[0])
             lc = Lightcurve(time, counts)
             lc_all.append(lc)
@@ -411,11 +410,10 @@ class TestAveragedPowerspectrum(object):
         ps = AveragedPowerspectrum(lc_all, 10.0, norm="leahy")
 
         assert np.isclose(np.mean(ps.power), 2.0, atol=1e-3, rtol=1e-3)
-        assert np.isclose(np.std(ps.power), 2.0/np.sqrt(n), atol=0.1, rtol=0.1)
+        assert np.isclose(np.std(ps.power), 2.0 / np.sqrt(n), atol=0.1, rtol=0.1)
 
 
 class TestClassicalSignificances(object):
-
     def test_function_runs(self):
         power = 2.0
         nspec = 1.0
@@ -493,10 +491,9 @@ class TestClassicalSignificances(object):
         pval1 = classical_pvalue(power1, nspec)
         pval2 = classical_pvalue(power2, nspec)
 
-        assert pval1-pval2 > 0.0
+        assert pval1 - pval2 > 0.0
 
     def test_pvalue_must_decrease_with_increasing_nspec(self):
-
         power = 3.0
         nspec1 = 1.0
         nspec2 = 10.0
@@ -504,10 +501,159 @@ class TestClassicalSignificances(object):
         pval1 = classical_pvalue(power, nspec1)
         pval2 = classical_pvalue(power, nspec2)
 
-        assert pval1-pval2 > 0.0
+        assert pval1 - pval2 > 0.0
 
     def test_very_large_powers_produce_zero_prob(self):
         power = 31000.0
         nspec = 1
         pval = classical_pvalue(power, nspec)
         assert np.isclose(pval, 0.0)
+
+
+class TestDynamicalPowerspectrum(object):
+    def setup_class(cls):
+        # generate timestamps
+        timestamps = np.linspace(1, 100, 10000)
+        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
+        lc = Lightcurve(timestamps, signal)
+        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 = DynamicalPowerspectrum(self.lc, segment_size=0)
+
+    def test_with_long_seg_size(self):
+        with pytest.raises(ValueError):
+            dps = DynamicalPowerspectrum(self.lc, segment_size=1000)
+
+    def test_matrix(self):
+        dps = DynamicalPowerspectrum(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):
+        dps = DynamicalPowerspectrum(self.lc, segment_size=3)
+        max_pos = dps.trace_maximum()
+
+        assert max(dps.freq[max_pos]) <= 1 / self.lc.dt
+        assert min(dps.freq[max_pos]) >= 1 / dps.segment_size
+
+    def test_trace_maximum_with_boundaries(self):
+        dps = DynamicalPowerspectrum(self.lc, segment_size=3)
+        minfreq = 21
+        maxfreq = 24
+        max_pos = dps.trace_maximum(min_freq=minfreq, max_freq=maxfreq)
+
+        assert max(dps.freq[max_pos]) <= maxfreq
+        assert min(dps.freq[max_pos]) >= minfreq
+
+    def test_size_of_trace_maximum(self):
+        dps = DynamicalPowerspectrum(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 = 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
+        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([1.5, 5.5, 9.5, 13.5])
+        rebin_dps = np.array([[0.7962963, 1.16402116, 0.28571429, 0.65625]])
+        dps = DynamicalPowerspectrum(self.lc_test, segment_size=segment_size)
+        dps.rebin_time(dt_new=dt_new)
+        assert np.allclose(dps.time, rebin_time)
+        assert np.allclose(dps.dyn_ps, rebin_dps)
+        assert np.isclose(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_dps = np.array([[5.76369293e-06],
+                              [7.07524761e-05],
+                              [6.24846189e+00],
+                              [5.77470465e-05],
+                              [1.76918128e-05]])
+        dps = DynamicalPowerspectrum(self.lc, segment_size=segment_size)
+        dps.rebin_frequency(df_new=df_new)
+        assert np.allclose(dps.freq, rebin_freq)
+        assert np.allclose(dps.dyn_ps, rebin_dps)
+        assert np.isclose(dps.df, df_new)
+
+    def test_rebin_time_mean_method(self):
+        segment_size = 3
+        dt_new = 4.0
+        rebin_time = np.array([1.5, 5.5, 9.5, 13.5])
+        rebin_dps = np.array([[0.59722222, 0.87301587, 0.21428571, 0.4921875]])
+        dps = DynamicalPowerspectrum(self.lc_test, segment_size=segment_size)
+        dps.rebin_time(dt_new=dt_new, method='mean')
+        assert np.allclose(dps.time, rebin_time)
+        assert np.allclose(dps.dyn_ps, rebin_dps)
+        assert np.isclose(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]])
+        dps = DynamicalPowerspectrum(self.lc, segment_size=segment_size)
+        dps.rebin_frequency(df_new=df_new, method="mean")
+        assert np.allclose(dps.freq, rebin_freq)
+        assert np.allclose(dps.dyn_ps, rebin_dps)
+        assert np.isclose(dps.df, df_new)
+
+    def test_rebin_time_average_method(self):
+        segment_size = 3
+        dt_new = 4.0
+        rebin_time = np.array([1.5, 5.5, 9.5, 13.5])
+        rebin_dps = np.array([[0.59722222, 0.87301587, 0.21428571, 0.4921875]])
+        dps = DynamicalPowerspectrum(self.lc_test, segment_size=segment_size)
+        dps.rebin_time(dt_new=dt_new, method='average')
+        assert np.allclose(dps.time, rebin_time)
+        assert np.allclose(dps.dyn_ps, rebin_dps)
+        assert np.isclose(dps.dt, dt_new)
+
+    def test_rebin_frequency_average_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]])
+        dps = DynamicalPowerspectrum(self.lc, segment_size=segment_size)
+        dps.rebin_frequency(df_new=df_new, method="average")
+        assert np.allclose(dps.freq, rebin_freq)
+        assert np.allclose(dps.dyn_ps, rebin_dps)
+        assert np.isclose(dps.df, df_new)