diff --git a/docs/api.rst b/docs/api.rst
index 8a84549f0..cebc3df21 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -5,12 +5,57 @@ Stingray API
 
 Library of Time Series Methods For Astronomical X-ray Data.
 
+Base Class
+==========
+Most `stingray`` classes are subclasses of a single class, :class:`stingray.StingrayObject`, which
+implements most of the I/O functionality and common behavior (e.g. strategies to combine data and
+make operations such as the sum, difference, or negation). This class is not intended to be
+instantiated directly, but rather to be used as a base class for other classes. Any class wanting
+to inherit from :class:`stingray.StingrayObject` should define a ``main_array_attr`` attribute, which
+defines the name of the attribute that will be used to store the "independent variable" main data array.
+For example, for all time series-like objects, the main array is the time array, while for the
+periodograms the main array is the frequency array.
+All arrays sharing the length (not the shape: they might be multi-dimensional!) of the main array are called
+"array attributes" and are accessible through the ``array_attrs`` method.
+When applying a mask or any other selection to a :class:`stingray.StingrayObject`,
+all array attributes are filtered in the same way. Some array-like attributes might have the same length
+by chance, in this case the user or the developer should add these to the ``not_array_attr`` attribute.
+For example, :class:`stingray.StingrayTimeseries` has ``gti`` among the not_array_attrs, since it is an
+array but not related 1-to-1 to the main array, even if in some cases it might happen to have the same numbers
+of elements of the main array, which is ``time``.
+
+StingrayObject
+--------------
+
+.. autoclass:: stingray.StingrayObject
+   :members:
+
+----
+
 Data Classes
 ============
 
 These classes define basic functionality related to common data types and typical methods
 that apply to these data types, including basic read/write functionality. Currently
-implemented are :class:`stingray.Lightcurve` and :class:`stingray.events.EventList`.
+implemented are :class:`stingray.StingrayTimeseries`, :class:`stingray.Lightcurve` and
+:class:`stingray.events.EventList`.
+
+All time series-like data classes inherit from :class:`stingray.StingrayTimeseries`, which
+implements most of the common functionality. The main data array is stored in the ``time``
+attribute.
+Good Time Intervals (GTIs) are stored in the ``gti`` attribute, which is a list of 2-tuples or 2-lists
+containing the start and stop times of each GTI. The ``gti`` attribute is not an array attribute, since
+it is not related 1-to-1 to the main array, even if in some cases it might happen to have the same number
+of elements of the main array. It is by default added to the ``not_array_attr`` attribute.
+
+
+StingrayTimeseries
+------------------
+
+.. autoclass:: stingray.StingrayTimeseries
+   :members:
+
+----
 
 Lightcurve
 ----------
@@ -88,7 +133,7 @@ Dynamical Powerspectrum
 .. autoclass:: stingray.DynamicalPowerspectrum
    :members:
    :inherited-members:
-   
+
 ----
 
 CrossCorrelation
diff --git a/docs/changes/780.feature.rst b/docs/changes/780.feature.rst
new file mode 100644
index 000000000..29058279f
--- /dev/null
+++ b/docs/changes/780.feature.rst
@@ -0,0 +1 @@
+Introduce power colors à la Heil et al. 2015
diff --git a/docs/changes/782.feature.rst b/docs/changes/782.feature.rst
new file mode 100644
index 000000000..1ccafdc8b
--- /dev/null
+++ b/docs/changes/782.feature.rst
@@ -0,0 +1 @@
+Add function to randomize data in small bad time intervals
diff --git a/docs/changes/787.trivial.rst b/docs/changes/787.trivial.rst
new file mode 100644
index 000000000..a94fd889c
--- /dev/null
+++ b/docs/changes/787.trivial.rst
@@ -0,0 +1 @@
+More informative GTI messages
diff --git a/docs/index.rst b/docs/index.rst
index c6ae9caf1..a58f9c2df 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -24,14 +24,16 @@ Current Capabilities
 1. Data handling and simulation
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-* loading event lists from fits files of a few missions (RXTE/PCA, NuSTAR/FPM, XMM-Newton/EPIC, NICER/XTI)
-* constructing light curves from event data, various operations on light curves (e.g. addition, subtraction, joining, and truncation)
+* loading event lists from fits files (and generally good handling of OGIP-compliant missions, like RXTE/PCA, NuSTAR/FPM, XMM-Newton/EPIC, NICER/XTI)
+* constructing light curves and time series from event data
+* various operations on time series (e.g. addition, subtraction, joining, and truncation)
 * simulating a light curve with a given power spectrum
 * simulating a light curve from another light curve and a 1-d (time) or 2-d (time-energy) impulse response
 * simulating an event list from a given light curve _and_ with a given energy spectrum
 * Good Time Interval operations
+* Filling gaps in light curves with statistically sound fake data
 
-2. Fourier methods
+1. Fourier methods
 ~~~~~~~~~~~~~~~~~~
 * power spectra and cross spectra in Leahy, rms normalization, absolute rms and no normalization
 * averaged power spectra and cross spectra
@@ -44,6 +46,7 @@ Current Capabilities
 * bispectra; *needs testing*
 * (Bayesian) quasi-periodic oscillation searches
 * Lomb-Scargle periodograms and cross spectra
+* Power Colors
 
 3. Other time series methods
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -61,7 +64,6 @@ Other future additions we are currently implementing are:
 * bicoherence
 * phase-resolved spectroscopy of quasi-periodic oscillations
 * Fourier-frequency-resolved spectroscopy
-* power colours
 * full HEASARC-compatible mission support
 * pulsar searches with :math:`H`-test
 * binary pulsar searches
diff --git a/setup.cfg b/setup.cfg
index 256ed04c2..cfd51c42a 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -59,6 +59,7 @@ docs =
     docutils
     sphinx-astropy
     nbsphinx>=0.8.3,!=0.8.8
+    nbconvert<7.14
     pandoc
     ipython
     towncrier<22.12.0
diff --git a/stingray/base.py b/stingray/base.py
index 4e76cba6d..9b0cc62f5 100644
--- a/stingray/base.py
+++ b/stingray/base.py
@@ -15,8 +15,27 @@
 from astropy.time import Time, TimeDelta
 from astropy.units import Quantity
 
-from stingray.utils import sqsum
 from .io import _can_save_longdouble, _can_serialize_meta
+from .utils import (
+    sqsum,
+    assign_value_if_none,
+    make_nd_into_arrays,
+    make_1d_arrays_into_nd,
+    get_random_state,
+    find_nearest,
+    rebin_data,
+)
+from .gti import (
+    create_gti_mask,
+    check_gtis,
+    cross_two_gtis,
+    join_gtis,
+    gti_border_bins,
+    get_btis,
+    merge_gtis,
+    check_separate,
+    append_gtis,
+)
 
 from typing import TYPE_CHECKING, Type, TypeVar, Union
 
@@ -478,7 +497,6 @@ def to_pandas(self) -> DataFrame:
 
         """
         from pandas import DataFrame
-        from .utils import make_nd_into_arrays
 
         data = {}
         array_attrs = self.array_attrs() + [self.main_array_attr] + self.internal_array_attrs()
@@ -515,7 +533,6 @@ def from_pandas(cls: Type[Tso], ts: DataFrame) -> Tso:
 
         """
         import re
-        from .utils import make_1d_arrays_into_nd
 
         cls = cls()
 
@@ -1023,7 +1040,6 @@ def __getitem__(self, index):
         ts_new : :class:`StingrayObject` object
             The new :class:`StingrayObject` object with the set of selected data.
         """
-        from .utils import assign_value_if_none
 
         if isinstance(index, (int, np.integer)):
             start = index
@@ -1077,7 +1093,7 @@ class StingrayTimeseries(StingrayObject):
 
     dt: float
         The time resolution of the time series. Can be a scalar or an array attribute (useful
-        for non-uniformly sampled data or events from different instruments)
+        for non-evenly sampled data or events from different instruments)
 
     mjdref : float
         The MJD used as a reference for the time array.
@@ -1112,7 +1128,7 @@ class StingrayTimeseries(StingrayObject):
 
     dt: float
         The time resolution of the measurements. Can be a scalar or an array attribute (useful
-        for non-uniformly sampled data or events from different instruments)
+        for non-evenly sampled data or events from different instruments)
 
     mjdref : float
         The MJD used as a reference for the time array.
@@ -1200,8 +1216,6 @@ def gti(self, value):
 
     @property
     def mask(self):
-        from .gti import create_gti_mask
-
         if self._mask is None:
             self._mask = create_gti_mask(self.time, self.gti, dt=self.dt)
         return self._mask
@@ -1292,7 +1306,6 @@ def apply_gtis(self, new_gti=None, inplace: bool = True):
 
         """
         # I import here to avoid the risk of circular imports
-        from .gti import check_gtis, create_gti_mask
 
         if new_gti is None:
             new_gti = self.gti
@@ -1324,7 +1337,6 @@ def split_by_gti(self, gti=None, min_points=2):
         list_of_tss : list
             A list of :class:`StingrayTimeseries` objects, one for each GTI segment
         """
-        from .gti import gti_border_bins, create_gti_mask
 
         if gti is None:
             gti = self.gti
@@ -1530,8 +1542,6 @@ def _operation_with_other_obj(
             other = other.change_mjdref(self.mjdref)
 
         if not np.array_equal(self.gti, other.gti):
-            from .gti import cross_two_gtis
-
             warnings.warn(
                 "The good time intervals in the two time series are different. Data outside the "
                 "common GTIs will be discarded."
@@ -1634,8 +1644,6 @@ def __getitem__(self, index):
         >>> assert np.allclose(ts[2].counts, [33])
         >>> assert np.allclose(ts[:2].counts, [11, 22])
         """
-        from .utils import assign_value_if_none
-        from .gti import cross_two_gtis
 
         new_ts = super().__getitem__(index)
         step = 1
@@ -1721,7 +1729,6 @@ def truncate(self, start=0, stop=None, method="index"):
 
     def _truncate_by_index(self, start, stop):
         """Private method for truncation using index values."""
-        from .gti import cross_two_gtis
 
         new_ts = self.apply_mask(slice(start, stop))
 
@@ -1835,7 +1842,6 @@ def _join_timeseries(self, others, strategy="intersection", ignore_meta=[]):
         `ts_new` : :class:`StingrayTimeseries` object
             The resulting :class:`StingrayTimeseries` object.
         """
-        from .gti import check_separate, cross_gtis, append_gtis
 
         new_ts = type(self)()
 
@@ -1869,8 +1875,6 @@ def _join_timeseries(self, others, strategy="intersection", ignore_meta=[]):
 
         all_objs = [self] + others
 
-        from .gti import merge_gtis
-
         # Check if none of the GTIs was already initialized.
         all_gti = [obj._gti for obj in all_objs if obj._gti is not None]
 
@@ -2039,7 +2043,6 @@ def rebin(self, dt_new=None, f=None, method="sum"):
         ts_new: :class:`StingrayTimeseries` object
             The :class:`StingrayTimeseries` object with the new, binned time series.
         """
-        from .utils import rebin_data
 
         if f is None and dt_new is None:
             raise ValueError("You need to specify at least one between f and " "dt_new")
@@ -2134,6 +2137,162 @@ def sort(self, reverse=False, inplace=False):
             mask = mask[::-1]
         return self.apply_mask(mask, inplace=inplace)
 
+    def fill_bad_time_intervals(
+        self,
+        max_length=None,
+        attrs_to_randomize=None,
+        buffer_size=None,
+        even_sampling=None,
+        seed=None,
+    ):
+        """Fill short bad time intervals with random data.
+
+        .. warning::
+            This method is only appropriate for *very short* bad time intervals. The simulated data
+            are basically white noise, so they are able to alter the statistical properties of
+            variable data. For very short gaps in the data, the effect of these small
+            injections of white noise should be negligible. How short depends on the single case,
+            the user is urged not to use the method as a black box and make simulations to measure
+            its effect. If you have long bad time intervals, you should use more advanced
+            techniques, not currently available in Stingray for this use case, such as Gaussian
+            Processes. In particular, please verify that the values of ``max_length`` and
+            ``buffer_size`` are adequate to your case.
+
+        To fill the gaps in all but the time points (i.e., flux measures, energies), we take the
+        ``buffer_size`` (by default, the largest value between 100 and the estimated samples in
+        a ``max_length``-long gap) valid data points closest to the gap and repeat them randomly
+        with the same empirical statistical distribution. So, if the `my_fancy_attr` attribute, in
+        the 100 points of the buffer, has 30 times 10, 10 times 9, and 60 times 11, there will be
+        *on average* 30% of 10, 60% of 11, and 10% of 9 in the simulated data.
+
+        Times are treated differently depending on the fact that the time series is evenly
+        sampled or not. If it is not, the times are simulated from a uniform distribution with the
+        same count rate found in the buffer. Otherwise, times just follow the same grid used
+        inside GTIs. Using the evenly sampled or not is decided based on the ``even_sampling``
+        parameter. If left to ``None``, the time series is considered evenly sampled if
+        ``self.dt`` is greater than zero and the median separation between subsequent times is
+        within 1% of the time resolution.
+
+        Other Parameters
+        ----------------
+        max_length : float
+            Maximum length of a bad time interval to be filled. If None, the criterion is bad
+            time intervals shorter than 1/100th of the longest good time interval.
+        attrs_to_randomize : list of str, default None
+            List of array_attrs to randomize. ``If None``, all array_attrs are randomized.
+            It should not include ``time`` and ``_mask``, which are treated separately.
+        buffer_size : int, default 100
+            Number of good data points to use to calculate the means and variance the random data
+            on each side of the bad time interval
+        even_sampling : bool, default None
+            Force the treatment of the data as evenly sampled or not. If None, the data are
+            considered evenly sampled if ``self.dt`` is larger than zero and the median
+            separation between subsequent times is within 1% of ``self.dt``.
+        seed : int, default None
+            Random seed to use for the simulation. If None, a random seed is generated.
+
+        """
+
+        rs = get_random_state(seed)
+
+        if attrs_to_randomize is None:
+            attrs_to_randomize = self.array_attrs() + self.internal_array_attrs()
+            for attr in ["time", "_mask"]:
+                if attr in attrs_to_randomize:
+                    attrs_to_randomize.remove(attr)
+
+        attrs_to_leave_alone = [
+            a
+            for a in self.array_attrs() + self.internal_array_attrs()
+            if a not in attrs_to_randomize
+        ]
+
+        if max_length is None:
+            max_length = np.max(self.gti[:, 1] - self.gti[:, 0]) / 100
+
+        btis = get_btis(self.gti, self.time[0], self.time[-1])
+        if len(btis) == 0:
+            logging.info("No bad time intervals to fill")
+            return copy.deepcopy(self)
+        filtered_times = self.time[self.mask]
+
+        new_times = [filtered_times.copy()]
+        new_attrs = {}
+        mean_data_separation = np.median(np.diff(filtered_times))
+        if even_sampling is None:
+            # The time series is considered evenly sampled if the median separation between
+            # subsequent times is within 1% of the time resolution
+            even_sampling = False
+            if self.dt > 0 and np.isclose(mean_data_separation, self.dt, rtol=0.01):
+                even_sampling = True
+            logging.info(f"Data are {'not' if not even_sampling else ''} evenly sampled")
+
+        if even_sampling:
+            est_samples_in_gap = int(max_length / self.dt)
+        else:
+            est_samples_in_gap = int(max_length / mean_data_separation)
+
+        if buffer_size is None:
+            buffer_size = max(100, est_samples_in_gap)
+
+        added_gtis = []
+
+        total_filled_time = 0
+        for bti in btis:
+            length = bti[1] - bti[0]
+            if length > max_length:
+                continue
+            logging.info(f"Filling bad time interval {bti} ({length:.4f} s)")
+            epsilon = 1e-5 * length
+            added_gtis.append([bti[0] - epsilon, bti[1] + epsilon])
+            filt_low_t, filt_low_idx = find_nearest(filtered_times, bti[0])
+            filt_hig_t, filt_hig_idx = find_nearest(filtered_times, bti[1], side="right")
+            if even_sampling:
+                local_new_times = np.arange(bti[0] + self.dt / 2, bti[1], self.dt)
+                nevents = local_new_times.size
+            else:
+                low_time_arr = filtered_times[max(filt_low_idx - buffer_size, 0) : filt_low_idx]
+                high_time_arr = filtered_times[filt_hig_idx : buffer_size + filt_hig_idx]
+
+                ctrate = (
+                    np.count_nonzero(low_time_arr) / (filt_low_t - low_time_arr[0])
+                    + np.count_nonzero(high_time_arr) / (high_time_arr[-1] - filt_hig_t)
+                ) / 2
+                nevents = rs.poisson(ctrate * (bti[1] - bti[0]))
+                local_new_times = rs.uniform(bti[0], bti[1], nevents)
+            new_times.append(local_new_times)
+
+            for attr in attrs_to_randomize:
+                low_arr = getattr(self, attr)[max(buffer_size - filt_low_idx, 0) : filt_low_idx]
+                high_arr = getattr(self, attr)[filt_hig_idx : buffer_size + filt_hig_idx]
+                if attr not in new_attrs:
+                    new_attrs[attr] = [getattr(self, attr)[self.mask]]
+                new_attrs[attr].append(rs.choice(np.concatenate([low_arr, high_arr]), nevents))
+            for attr in attrs_to_leave_alone:
+                if attr not in new_attrs:
+                    new_attrs[attr] = [getattr(self, attr)[self.mask]]
+                if attr == "_mask":
+                    new_attrs[attr].append(np.ones(nevents, dtype=bool))
+                else:
+                    new_attrs[attr].append(np.zeros(nevents) + np.nan)
+            total_filled_time += length
+
+        logging.info(f"A total of {total_filled_time} s of data were simulated")
+
+        new_gtis = join_gtis(self.gti, added_gtis)
+        new_times = np.concatenate(new_times)
+        order = np.argsort(new_times)
+        new_obj = type(self)()
+        new_obj.time = new_times[order]
+
+        for attr in self.meta_attrs():
+            setattr(new_obj, attr, getattr(self, attr))
+
+        for attr, values in new_attrs.items():
+            setattr(new_obj, attr, np.concatenate(values)[order])
+        new_obj.gti = new_gtis
+        return new_obj
+
     def plot(
         self,
         attr,
@@ -2145,6 +2304,7 @@ def plot(
         save=False,
         filename=None,
         plot_btis=True,
+        axis_limits=None,
     ):
         """
         Plot the time series using ``matplotlib``.
@@ -2168,6 +2328,10 @@ def plot(
             could be ``['Time (s)', 'Counts (s^-1)']``
         ax : ``matplotlib.pyplot.axis`` object
             Axis to be used for plotting. Defaults to creating a new one.
+        axis_limits : list, tuple, string, default ``None``
+            Parameter to set axis properties of the ``matplotlib`` figure. For example
+            it can be a list like ``[xmin, xmax, ymin, ymax]`` or any other
+            acceptable argument for the``matplotlib.pyplot.axis()`` method.
         title : str, default ``None``
             The title of the plot.
         marker : str, default '-'
@@ -2182,17 +2346,23 @@ def plot(
             Plot the bad time intervals as red areas on the plot
         """
         import matplotlib.pyplot as plt
-        from .gti import get_btis
 
         if ax is None:
             plt.figure(attr)
             ax = plt.gca()
 
-        if labels is None:
+        valid_labels = (isinstance(labels, Iterable) and not isinstance(labels, str)) and len(
+            labels
+        ) == 2
+        if labels is not None and not valid_labels:
+            warnings.warn("``labels`` must be an iterable with two labels for x and y axes.")
+
+        if labels is None or not valid_labels:
             labels = ["Time (s)"] + [attr]
 
-        ylabel = labels[1]
         xlabel = labels[0]
+        ylabel = labels[1]
+        # Default values for labels
 
         ax.plot(self.time, getattr(self, attr), marker, ds="steps-mid", label=attr, zorder=10)
 
@@ -2202,11 +2372,15 @@ def plot(
                 getattr(self, attr),
                 yerr=getattr(self, attr + "_err"),
                 fmt="o",
+                zorder=10,
             )
 
         ax.set_ylabel(ylabel)
         ax.set_xlabel(xlabel)
 
+        if axis_limits is not None:
+            ax.set_xlim(axis_limits[0], axis_limits[1])
+            ax.set_ylim(axis_limits[2], axis_limits[3])
         if title is not None:
             ax.set_title(title)
 
@@ -2221,7 +2395,14 @@ def plot(
             tend = max(self.time[-1] + self.dt / 2, self.gti[-1, 1])
             btis = get_btis(self.gti, tstart, tend)
             for bti in btis:
-                plt.axvspan(bti[0], bti[1], alpha=0.5, color="r", zorder=10)
+                plt.axvspan(
+                    bti[0],
+                    bti[1],
+                    alpha=0.5,
+                    facecolor="r",
+                    zorder=1,
+                    edgecolor="none",
+                )
         return ax
 
 
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index d00cf83e6..e3e124a48 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -21,7 +21,7 @@
 from .fourier import avg_cs_from_iterables, error_on_averaged_cross_spectrum
 from .fourier import avg_cs_from_events, poisson_level
 from .fourier import fftfreq, fft, normalize_periodograms, raw_coherence
-from .fourier import get_flux_iterable_from_segments
+from .fourier import get_flux_iterable_from_segments, power_color
 
 from scipy.special import factorial
 
@@ -2028,6 +2028,9 @@ class DynamicalCrossspectrum(AveragedCrossspectrum):
         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``).
+
+    m: int
+        The number of averaged powers in each spectral bin (initially 1, it changes after rebinning).
     """
 
     def __init__(self, data1, data2, segment_size, norm="frac", gti=None, sample_time=None):
@@ -2075,16 +2078,21 @@ def _make_matrix(self, data1, data2):
             save_all=True,
         )
         self.dyn_ps = np.array(avg.cs_all).T
+        conv = avg.cs_all / avg.unnorm_cs_all
+        self.unnorm_conversion = np.nanmean(conv)
         self.freq = avg.freq
         current_gti = avg.gti
+        self.nphots1 = avg.nphots1
+        self.nphots2 = avg.nphots2
 
         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
+        self.m = 1
 
-    def rebin_frequency(self, df_new, method="sum"):
+    def rebin_frequency(self, df_new, method="average"):
         """
         Rebin the Dynamic Power Spectrum to a new frequency resolution.
         Rebinning is an in-place operation, i.e. will replace the existing
@@ -2101,14 +2109,14 @@ def rebin_frequency(self, df_new, method="sum"):
             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
+        method: {"sum" | "mean" | "average"}, optional, default "average"
+            This keyword argument sets whether the powers 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(
+            freq_new, bin_counts, bin_err, step = rebin_data(
                 self.freq, data, dx_new=df_new, method=method
             )
             dynspec_new.append(bin_counts)
@@ -2116,9 +2124,11 @@ def rebin_frequency(self, df_new, method="sum"):
         new_dynspec_object.freq = freq_new
         new_dynspec_object.dyn_ps = np.array(dynspec_new).T
         new_dynspec_object.df = df_new
+        new_dynspec_object.m = int(step) * self.m
+
         return new_dynspec_object
 
-    def rebin_time(self, dt_new, method="sum"):
+    def rebin_time(self, dt_new, method="average"):
         """
         Rebin the Dynamic Power Spectrum to a new time resolution.
 
@@ -2137,8 +2147,8 @@ def rebin_time(self, dt_new, method="sum"):
             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
+        method: {"sum" | "mean" | "average"}, optional, default "average"
+            This keyword argument sets whether the powers in the new bins
             should be summed or averaged.
 
         Returns
@@ -2156,7 +2166,7 @@ def rebin_time(self, dt_new, method="sum"):
 
         dynspec_new = []
         for data in self.dyn_ps:
-            time_new, bin_counts, _, _ = rebin_data(
+            time_new, bin_counts, _, step = rebin_data(
                 self.time, data, dt_new, method=method, dx=self.dt
             )
             dynspec_new.append(bin_counts)
@@ -2164,6 +2174,69 @@ def rebin_time(self, dt_new, method="sum"):
         new_dynspec_object.time = time_new
         new_dynspec_object.dyn_ps = np.array(dynspec_new)
         new_dynspec_object.dt = dt_new
+        new_dynspec_object.m = int(step) * self.m
+        return new_dynspec_object
+
+    def rebin_by_n_intervals(self, n, method="average"):
+        """
+        Rebin the Dynamic Power Spectrum to a new time resolution, by summing contiguous intervals.
+
+        This is different from meth:`DynamicalPowerspectrum.rebin_time` in that it averages ``n``
+        consecutive intervals regardless of their distance in time. ``rebin_time`` will instead
+        average intervals that are separated at most by a time ``dt_new``.
+
+        Parameters
+        ----------
+        n: int
+            The number of intervals to be combined into one.
+
+        method: {"sum" | "mean" | "average"}, optional, default "average"
+            This keyword argument sets whether the powers 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 not np.issubdtype(type(n), np.integer):
+            warnings.warn("n must be an integer. Casting to int")
+
+            n = int(n)
+
+        if n == 1:
+            return copy.deepcopy(self)
+        if n < 1:
+            raise ValueError("n must be >= 1")
+
+        new_dynspec_object = copy.deepcopy(self)
+
+        dynspec_new = []
+        time_new = []
+        for i, data in enumerate(self.dyn_ps.T):
+            if i % n == 0:
+                count = 1
+                bin_counts = data
+                bin_times = self.time[i]
+            else:
+                count += 1
+                bin_counts += data
+                bin_times += self.time[i]
+            if count == n:
+                if method in ["mean", "average"]:
+                    bin_counts /= n
+                dynspec_new.append(bin_counts)
+                bin_times /= n
+                time_new.append(bin_times)
+
+        new_dynspec_object.time = time_new
+        new_dynspec_object.dyn_ps = np.array(dynspec_new).T
+        new_dynspec_object.dt *= n
+        new_dynspec_object.m = n * self.m
+
         return new_dynspec_object
 
     def trace_maximum(self, min_freq=None, max_freq=None):
@@ -2198,6 +2271,125 @@ def trace_maximum(self, min_freq=None, max_freq=None):
 
         return np.array(max_positions)
 
+    def power_colors(
+        self,
+        freq_edges=[1 / 256, 1 / 32, 0.25, 2.0, 16.0],
+        freqs_to_exclude=None,
+        poisson_power=0,
+    ):
+        """
+        Compute the power colors of the dynamical power spectrum.
+
+        See Heil et al. 2015, MNRAS, 448, 3348
+
+        Parameters
+        ----------
+        freq_edges: iterable
+            The edges of the frequency bins to be used for the power colors.
+
+        freqs_to_exclude : 1-d or 2-d iterable, optional, default None
+            The ranges of frequencies to exclude from the calculation of the power color.
+            For example, the frequencies containing strong QPOs.
+            A 1-d iterable should contain two values for the edges of a single range. (E.g.
+            ``[0.1, 0.2]``). ``[[0.1, 0.2], [3, 4]]`` will exclude the ranges 0.1-0.2 Hz and
+            3-4 Hz.
+
+        poisson_power : float or iterable, optional, default 0
+            The Poisson noise level of the power spectrum. If iterable, it should have the same
+            length as ``freq``. (This might apply to the case of a power spectrum with a
+            strong dead time distortion)
+
+        Returns
+        -------
+        pc0: np.ndarray
+        pc0_err: np.ndarray
+        pc1: np.ndarray
+        pc1_err: np.ndarray
+            The power colors for each spectrum and their respective errors
+        """
+        power_colors = []
+        if np.iscomplexobj(self.dyn_ps):
+            warnings.warn("When using power_colors, complex powers will be cast to real.")
+
+        for ps in self.dyn_ps.T:
+            power_colors.append(
+                power_color(
+                    self.freq,
+                    ps.real,
+                    freq_edges=freq_edges,
+                    freqs_to_exclude=freqs_to_exclude,
+                    df=self.df,
+                    poisson_power=poisson_power,
+                    m=self.m,
+                )
+            )
+
+        pc0, pc0_err, pc1, pc1_err = np.array(power_colors).T
+        return pc0, pc0_err, pc1, pc1_err
+
+    def compute_rms(self, min_freq, max_freq, poisson_noise_level=0):
+        """
+        Compute the fractional rms amplitude in the power spectrum
+        between two frequencies.
+
+        Parameters
+        ----------
+        min_freq: float
+            The lower frequency bound for the calculation.
+
+        max_freq: float
+            The upper frequency bound for the calculation.
+
+        Other parameters
+        ----------------
+        poisson_noise_level : float, default is None
+            This is the Poisson noise level of the PDS with same
+            normalization as the PDS.
+
+        Returns
+        -------
+        rms: float
+            The fractional rms amplitude contained between ``min_freq`` and
+            ``max_freq``.
+
+        rms_err: float
+            The error on the fractional rms amplitude.
+
+        """
+        from .fourier import rms_calculation
+
+        dyn_ps_unnorm = self.dyn_ps / self.unnorm_conversion
+        poisson_noise_unnrom = poisson_noise_level / self.unnorm_conversion
+        if not hasattr(self, "nphots"):
+            nphots = (self.nphots1 * self.nphots2) ** 0.5
+        else:
+            nphots = self.nphots
+        T = self.dt * self.m
+        M_freqs = self.m
+        K_freqs = 1
+        minind = self.freq.searchsorted(min_freq)
+        maxind = self.freq.searchsorted(max_freq)
+        min_freq = self.freq[minind]
+        freq_bins = maxind - minind
+        rmss = []
+        rms_errs = []
+        for unnorm_powers in dyn_ps_unnorm.T:
+            r, re = rms_calculation(
+                unnorm_powers[minind:maxind],
+                min_freq,
+                max_freq,
+                nphots,
+                T,
+                M_freqs,
+                K_freqs,
+                freq_bins,
+                poisson_noise_unnrom,
+                deadtime=0.0,
+            )
+            rmss.append(r)
+            rms_errs.append(re)
+        return rmss, rms_errs
+
 
 def crossspectrum_from_time_array(
     times1,
diff --git a/stingray/fourier.py b/stingray/fourier.py
index 608fb9cb9..a39d1da5e 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -15,6 +15,309 @@
 from .utils import fft, fftfreq, histogram, show_progress, sum_if_not_none_or_initialize
 
 
+def integrate_power_in_frequency_range(
+    frequency,
+    power,
+    frequency_range,
+    power_err=None,
+    df=None,
+    m=1,
+    poisson_power=0,
+):
+    """
+    Integrate the power in a given frequency range.
+
+    Parameters
+    ----------
+    frequency : iterable
+        The frequencies of the power spectrum
+    power : iterable
+        The power at each frequency
+    frequency_range : iterable of length 2
+        The frequency range to integrate
+    power_err : iterable, optional, default None
+        The power error bar at each frequency
+    df : float or float iterable, optional, default None
+        The frequency resolution of the input data. If None, it is calculated
+        from the median difference of input frequencies.
+    m : int, optional, default 1
+        The number of segments and/or contiguous frequency bins averaged to obtain power.
+        Only needed if ``power_err`` is None
+    poisson_power : float, optional, default 0
+        The Poisson noise level of the power spectrum.
+
+    Returns
+    -------
+    power_integrated : float
+        The integrated power
+    power_integrated_err : float
+        The error on the integrated power
+
+    """
+    frequency = np.array(frequency, dtype=float)
+    power = np.array(power)
+    if not np.iscomplexobj(power):
+        power = power.astype(float)
+    if not isinstance(poisson_power, Iterable):
+        poisson_power = np.ones_like(frequency) * poisson_power
+    if df is None:
+        df = np.median(np.diff(frequency))
+    if not isinstance(df, Iterable):
+        df = np.ones_like(frequency) * df
+
+    # The frequency_range is in the middle or at the edge of each bin. When only a part of the
+    # bin is included, we need to add the fraction of the bin that is included.
+    frequency_mask = (frequency + df / 2 > frequency_range[0]) & (
+        frequency - df / 2 < frequency_range[1]
+    )
+
+    freqs_to_integrate = frequency[frequency_mask]
+    poisson_power = poisson_power[frequency_mask]
+    correction_ratios = np.ones_like(freqs_to_integrate)
+    dfs_to_integrate = df[frequency_mask]
+
+    # The first and last bins are only partially included. We need to correct for the
+    # fraction of the bin actually included.
+    correction_ratios[0] = (
+        freqs_to_integrate[0] + 0.5 * dfs_to_integrate[0] - frequency_range[0]
+    ) / dfs_to_integrate[0]
+    correction_ratios[-1] = (
+        frequency_range[-1] - freqs_to_integrate[-1] + 0.5 * dfs_to_integrate[-1]
+    ) / dfs_to_integrate[-1]
+    dfs_to_integrate = dfs_to_integrate * correction_ratios
+
+    powers_to_integrate = power[frequency_mask]
+
+    if power_err is None:
+        power_err_to_integrate = powers_to_integrate / np.sqrt(m)
+    else:
+        power_err_to_integrate = np.asarray(power_err)[frequency_mask]
+
+    power_integrated = np.sum((powers_to_integrate - poisson_power) * dfs_to_integrate)
+    power_err_integrated = np.sqrt(np.sum((power_err_to_integrate * dfs_to_integrate) ** 2))
+    return power_integrated, power_err_integrated
+
+
+def power_color(
+    frequency,
+    power,
+    power_err=None,
+    freq_edges=[1 / 256, 1 / 32, 0.25, 2.0, 16.0],
+    df=None,
+    m=1,
+    freqs_to_exclude=None,
+    poisson_power=0,
+    return_log=False,
+):
+    """
+    Calculate two power colors from a power spectrum.
+
+    Power colors are an alternative to spectral colors to understand the spectral state of an
+    accreting source. They are defined as the ratio of the power in two frequency ranges,
+    analogously to the colors calculated from electromagnetic spectra.
+
+    This function calculates two power colors, using the four frequency ranges contained
+    between the five frequency edges in ``freq_edges``. Given [f0, f1, f2, f3, f4], the
+    two power colors are calculated as the following ratios of the integrated power
+    (which are variances):
+
+    + PC0 = Var([f0, f1]) / Var([f2, f3])
+
+    + PC1 = Var([f1, f2]) / Var([f3, f4])
+
+    Errors are calculated using simple error propagation from the integrated power errors.
+
+    See Heil et al. 2015, MNRAS, 448, 3348
+
+    Parameters
+    ----------
+    frequency : iterable
+        The frequencies of the power spectrum
+    power : iterable
+        The power at each frequency
+
+    Other Parameters
+    ----------------
+    power_err : iterable
+        The power error bar at each frequency
+    freq_edges : iterable, optional, default ``[0.0039, 0.031, 0.25, 2.0, 16.0]``
+        The five edges defining the four frequency intervals to use to calculate the power color.
+        If empty, the power color is calculated using the frequencies from Heil et al. 2015.
+    df : float or float iterable, optional, default None
+        The frequency resolution of the input data. If None, it is calculated
+        from the median difference of input frequencies.
+    m : int, optional, default 1
+        The number of segments and/or contiguous frequency bins averaged to obtain power
+    freqs_to_exclude : 1-d or 2-d iterable, optional, default None
+        The ranges of frequencies to exclude from the calculation of the power color.
+        For example, the frequencies containing strong QPOs.
+        A 1-d iterable should contain two values for the edges of a single range. (E.g.
+        ``[0.1, 0.2]``). ``[[0.1, 0.2], [3, 4]]`` will exclude the ranges 0.1-0.2 Hz and 3-4 Hz.
+    poisson_power : float or iterable, optional, default 0
+        The Poisson noise level of the power spectrum. If iterable, it should have the same
+        length as ``frequency``. (This might apply to the case of a power spectrum with a
+        strong dead time distortion
+    return_log : bool, optional, default False
+        Return the base-10 logarithm of the power color and the errors
+
+    Returns
+    -------
+    PC0 : float
+        The first power color
+    PC0_err : float
+        The error on the first power color
+    PC1 : float
+        The second power color
+    PC1_err : float
+        The error on the second power color
+    """
+    freq_edges = np.asarray(freq_edges)
+    if len(freq_edges) != 5:
+        raise ValueError("freq_edges must have 5 elements")
+
+    frequency = np.asarray(frequency)
+    power = np.asarray(power)
+
+    if df is None:
+        df = np.median(np.diff(frequency))
+    input_frequency_low_edges = frequency - df / 2
+    input_frequency_high_edges = frequency + df / 2
+
+    if freq_edges.min() < input_frequency_low_edges[0]:
+        raise ValueError("The minimum frequency is larger than the first frequency edge")
+    if freq_edges.max() > input_frequency_high_edges[-1]:
+        raise ValueError("The maximum frequency is lower than the last frequency edge")
+
+    if power_err is None:
+        power_err = power / np.sqrt(m)
+    else:
+        power_err = np.asarray(power_err)
+
+    if freqs_to_exclude is not None:
+        if len(np.shape(freqs_to_exclude)) == 1:
+            freqs_to_exclude = [freqs_to_exclude]
+
+        if (
+            not isinstance(freqs_to_exclude, Iterable)
+            or len(np.shape(freqs_to_exclude)) != 2
+            or np.shape(freqs_to_exclude)[1] != 2
+        ):
+            raise ValueError("freqs_to_exclude must be of format [[f0, f1], [f2, f3], ...]")
+        for f0, f1 in freqs_to_exclude:
+            frequency_mask = (input_frequency_low_edges > f0) & (input_frequency_high_edges < f1)
+            idx0, idx1 = np.searchsorted(frequency, [f0, f1])
+            power[frequency_mask] = np.mean([power[idx0], power[idx1]])
+
+    var00, var00_err = integrate_power_in_frequency_range(
+        frequency,
+        power,
+        freq_edges[:2],
+        power_err=power_err,
+        df=df,
+        m=m,
+        poisson_power=poisson_power,
+    )
+    var01, var01_err = integrate_power_in_frequency_range(
+        frequency,
+        power,
+        freq_edges[2:4],
+        power_err=power_err,
+        df=df,
+        m=m,
+        poisson_power=poisson_power,
+    )
+    var10, var10_err = integrate_power_in_frequency_range(
+        frequency,
+        power,
+        freq_edges[1:3],
+        power_err=power_err,
+        df=df,
+        m=m,
+        poisson_power=poisson_power,
+    )
+    var11, var11_err = integrate_power_in_frequency_range(
+        frequency,
+        power,
+        freq_edges[3:5],
+        power_err=power_err,
+        df=df,
+        m=m,
+        poisson_power=poisson_power,
+    )
+    pc0 = var00 / var01
+    pc1 = var10 / var11
+    pc0_err = pc0 * (var00_err / var00 + var01_err / var01)
+    pc1_err = pc1 * (var10_err / var10 + var11_err / var11)
+    if return_log:
+        pc0_err = 1 / pc0 * pc0_err
+        pc1_err = 1 / pc1 * pc1_err
+        pc0 = np.log10(pc0)
+        pc1 = np.log10(pc1)
+    return pc0, pc0_err, pc1, pc1_err
+
+
+def hue_from_power_color(pc0, pc1, center=[4.51920, 0.453724]):
+    """Measure the angle of a point in the log-power color diagram wrt the center.
+
+    Angles are measured in radians, **in the clockwise direction**, with respect to a line oriented
+    at -45 degrees wrt the horizontal axis.
+
+    See Heil et al. 2015, MNRAS, 448, 3348
+
+    Parameters
+    ----------
+    pc0 : float
+        The (linear, not log!) power color in the first frequency range
+    pc1 : float
+        The (linear, not log!) power color in the second frequency range
+
+    Other Parameters
+    ----------------
+    center : iterable, optional, default [4.51920, 0.453724]
+        The coordinates of the center of the power color diagram
+
+    Returns
+    -------
+    hue : float
+        The angle of the point wrt the center, in radians
+    """
+    pc0 = np.log10(pc0)
+    pc1 = np.log10(pc1)
+
+    center = np.log10(np.asarray(center))
+
+    return hue_from_logpower_color(pc0, pc1, center=center)
+
+
+def hue_from_logpower_color(log10pc0, log10pc1, center=(np.log10(4.51920), np.log10(0.453724))):
+    """Measure the angle of a point in the log-power color diagram wrt the center.
+
+    Angles are measured in radians, **in the clockwise direction**, with respect to a line oriented
+    at -45 degrees wrt the horizontal axis.
+
+    See Heil et al. 2015, MNRAS, 448, 3348
+
+    Parameters
+    ----------
+    log10pc0 : float
+        The log10 power color in the first frequency range
+    log10pc1 : float
+        The log10 power color in the second frequency range
+
+    Other Parameters
+    ----------------
+    center : iterable, optional, default ``log10([4.51920, 0.453724])``
+        The coordinates of the center of the power color diagram
+
+    Returns
+    -------
+    hue : float
+        The angle of the point wrt the center, in radians
+    """
+    hue = 3 / 4 * np.pi - np.arctan2(log10pc1 - center[1], log10pc0 - center[0])
+    return hue
+
+
 def positive_fft_bins(n_bin, include_zero=False):
     """
     Give the range of positive frequencies of a complex FFT.
@@ -1207,7 +1510,7 @@ def local_show_progress(a):
 
         # If the user wants to normalize using the mean of the total
         # lightcurve, normalize it here
-        cs_seg = unnorm_power
+        cs_seg = copy.deepcopy(unnorm_power)
         if not use_common_mean:
             mean = n_ph / n_bin
 
@@ -1593,7 +1896,7 @@ def local_show_progress(a):
                 unnorm_pd1 = unnorm_pd1[fgt0]
                 unnorm_pd2 = unnorm_pd2[fgt0]
 
-        cs_seg = unnorm_power
+        cs_seg = copy.deepcopy(unnorm_power)
         p1_seg = unnorm_pd1
         p2_seg = unnorm_pd2
 
diff --git a/stingray/io.py b/stingray/io.py
index 1cdd7d83a..465394a72 100644
--- a/stingray/io.py
+++ b/stingray/io.py
@@ -2,7 +2,8 @@
 import math
 import copy
 import os
-import pickle
+import sys
+import traceback
 import warnings
 from collections.abc import Iterable
 
@@ -596,7 +597,8 @@ def load_events_and_gtis(
     if modekey is not None and modekey in probe_header:
         mode = probe_header[modekey].strip()
 
-    gtistring = get_key_from_mission_info(db, "gti", "GTI,STDGTI", instr, mode)
+    if gtistring is None:
+        gtistring = get_key_from_mission_info(db, "gti", "GTI,STDGTI", instr, mode)
     if hduname is None:
         hduname = get_key_from_mission_info(db, "events", "EVENTS", instr, mode)
 
@@ -655,10 +657,12 @@ def load_events_and_gtis(
                 accepted_gtistrings=accepted_gtistrings,
                 det_numbers=det_number,
             )
-        except Exception:  # pragma: no cover
+        except Exception as e:  # pragma: no cover
             warnings.warn(
-                "No extensions found with a valid name. "
-                "Please check the `accepted_gtistrings` values.",
+                (
+                    f"No valid GTI extensions found. \nError: {str(e)}\n"
+                    "GTIs will be set to the entire time series."
+                ),
                 AstropyUserWarning,
             )
             gti_list = np.array([[t_start, t_stop]], dtype=np.longdouble)
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index 8602482ef..f53933387 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -1608,11 +1608,14 @@ def plot(
         self,
         witherrors=False,
         labels=None,
-        axis=None,
+        ax=None,
         title=None,
         marker="-",
         save=False,
         filename=None,
+        axis_limits=None,
+        axis=None,
+        plot_btis=True,
     ):
         """
         Plot the light curve using ``matplotlib``.
@@ -1629,11 +1632,14 @@ def plot(
         labels : iterable, default ``None``
             A list of tuple with ``xlabel`` and ``ylabel`` as strings.
 
-        axis : list, tuple, string, default ``None``
+        axis_limits : list, tuple, string, default ``None``
             Parameter to set axis properties of the ``matplotlib`` figure. For example
             it can be a list like ``[xmin, xmax, ymin, ymax]`` or any other
             acceptable argument for the``matplotlib.pyplot.axis()`` method.
 
+        axis : list, tuple, string, default ``None``
+            Deprecated in favor of ``axis_limits``, same functionality.
+
         title : str, default ``None``
             The title of the plot.
 
@@ -1647,39 +1653,37 @@ def plot(
 
         filename : str
             File name of the image to save. Depends on the boolean ``save``.
-        """
 
-        fig = plt.figure()
-        if witherrors:
-            fig = plt.errorbar(self.time, self.counts, yerr=self.counts_err, fmt=marker)
-        else:
-            fig = plt.plot(self.time, self.counts, marker)
-
-        if labels is not None:
-            try:
-                plt.xlabel(labels[0])
-                plt.ylabel(labels[1])
-            except TypeError:
-                utils.simon("``labels`` must be either a list or tuple with " "x and y labels.")
-                raise
-            except IndexError:
-                utils.simon("``labels`` must have two labels for x and y " "axes.")
-                # Not raising here because in case of len(labels)==1, only
-                # x-axis will be labelled.
+        ax : ``matplotlib.pyplot.axis`` object
+            Axis to be used for plotting. Defaults to creating a new one.
 
+        plot_btis : bool
+            Plot the bad time intervals as red areas on the plot
+        """
         if axis is not None:
-            plt.axis(axis)
-
-        if title is not None:
-            plt.title(title)
-
-        if save:
-            if filename is None:
-                plt.savefig("out.png")
-            else:
-                plt.savefig(filename)
-        else:
-            plt.show(block=False)
+            warnings.warn(
+                "The ``axis`` argument is deprecated in favor of ``axis_limits``. "
+                "Please use that instead.",
+                DeprecationWarning,
+            )
+            axis_limits = axis
+
+        flux_attr = "counts"
+        if not self.input_counts:
+            flux_attr = "countrate"
+
+        return super().plot(
+            flux_attr,
+            witherrors=witherrors,
+            labels=labels,
+            ax=ax,
+            title=title,
+            marker=marker,
+            save=save,
+            filename=filename,
+            plot_btis=plot_btis,
+            axis_limits=axis_limits,
+        )
 
     @classmethod
     def read(
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 82b92faa7..53495b239 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -210,7 +210,8 @@ def compute_rms(
             M_freq = self.m
             K_freq = self.k
             freq_bins = maxind - minind
-        T = self.dt * self.n
+
+        T = self.dt * self.n * 2
 
         if white_noise_offset is not None:
             powers = self.power[minind:maxind]
@@ -237,6 +238,7 @@ def compute_rms(
                 poisson_noise_unnorm = unnormalize_periodograms(
                     poisson_noise_level, self.dt, self.n, self.nphots, norm=self.norm
                 )
+
             return rms_calculation(
                 self.unnorm_power[minind:maxind],
                 min_freq,
@@ -1013,6 +1015,9 @@ class DynamicalPowerspectrum(DynamicalCrossspectrum):
         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``).
+
+    m: int
+        The number of averaged cross spectra.
     """
 
     def __init__(self, lc, segment_size, norm="frac", gti=None, sample_time=None):
@@ -1055,6 +1060,8 @@ def _make_matrix(self, lc):
             gti=self.gti,
             save_all=True,
         )
+        conv = avg.cs_all / avg.unnorm_cs_all
+        self.unnorm_conversion = np.nanmean(conv)
         self.dyn_ps = np.array(avg.cs_all).T
         self.freq = avg.freq
         current_gti = avg.gti
@@ -1064,6 +1071,57 @@ def _make_matrix(self, lc):
         self.time = tstart + 0.5 * self.segment_size
         self.df = avg.df
         self.dt = self.segment_size
+        self.meanrate = avg.nphots / avg.n / avg.dt
+        self.nphots = avg.nphots
+        self.m = 1
+
+    def power_colors(
+        self,
+        freq_edges=[1 / 256, 1 / 32, 0.25, 2.0, 16.0],
+        freqs_to_exclude=None,
+        poisson_power=None,
+    ):
+        """
+        Return the power colors of the dynamical power spectrum.
+
+        Parameters
+        ----------
+        freq_edges: iterable
+            The edges of the frequency bins to be used for the power colors.
+
+        freqs_to_exclude : 1-d or 2-d iterable, optional, default None
+            The ranges of frequencies to exclude from the calculation of the power color.
+            For example, the frequencies containing strong QPOs.
+            A 1-d iterable should contain two values for the edges of a single range. (E.g.
+            ``[0.1, 0.2]``). ``[[0.1, 0.2], [3, 4]]`` will exclude the ranges 0.1-0.2 Hz and
+            3-4 Hz.
+
+        poisson_level : float or iterable, optional
+            Defaults to the theoretical Poisson noise level (e.g. 2 for Leahy normalization).
+            The Poisson noise level of the power spectrum. If iterable, it should have the same
+            length as ``frequency``. (This might apply to the case of a power spectrum with a
+            strong dead time distortion)
+
+        Returns
+        -------
+        pc0: np.ndarray
+        pc0_err: np.ndarray
+        pc1: np.ndarray
+        pc1_err: np.ndarray
+            The power colors for each spectrum and their respective errors
+        """
+        if poisson_power is None:
+            poisson_power = poisson_level(
+                norm=self.norm,
+                meanrate=self.meanrate,
+                n_ph=self.nphots,
+            )
+
+        return super().power_colors(
+            freq_edges=freq_edges,
+            freqs_to_exclude=freqs_to_exclude,
+            poisson_power=poisson_power,
+        )
 
 
 def powerspectrum_from_time_array(
diff --git a/stingray/pulse/overlapandsave/ols.py b/stingray/pulse/overlapandsave/ols.py
index 831372bd0..d6ac909b2 100644
--- a/stingray/pulse/overlapandsave/ols.py
+++ b/stingray/pulse/overlapandsave/ols.py
@@ -196,7 +196,7 @@ def prepareh(h, nfft: List[int], rfftn=None):
         The FFT-transformed, conjugate filter array
     """
     rfftn = rfftn or np.fft.rfftn
-    return np.conj(rfftn(flip(np.conj(h)), nfft))
+    return np.conj(rfftn(flip(np.conj(h)), nfft, axes=np.arange(len(nfft))))
 
 
 def slice2range(s: slice):
@@ -365,7 +365,9 @@ def olsStep(
         for (start, length, nh, border) in zip(starts, lengths, nh, border)
     )
     xpart = padEdges(x, slices, mode=mode, **kwargs)
-    output = irfftn(rfftn(xpart, nfft) * hfftconj, nfft)
+    output = irfftn(
+        rfftn(xpart, nfft, axes=np.arange(len(nfft))) * hfftconj, nfft, axes=np.arange(len(nfft))
+    )
     return output[tuple(slice(0, s) for s in lengths)]
 
 
diff --git a/stingray/tests/test_base.py b/stingray/tests/test_base.py
index bbeb1ae04..56624dc5a 100644
--- a/stingray/tests/test_base.py
+++ b/stingray/tests/test_base.py
@@ -1291,3 +1291,133 @@ def test_join_ignore_attr(self):
         assert np.allclose(ts_new.time, [1, 2, 3, 4, 5, 7])
         assert not hasattr(ts_new, "instr")
         assert ts_new.mission == (1, 2)
+
+
+class TestFillBTI(object):
+    @classmethod
+    def setup_class(cls):
+        cls.rand_time = np.sort(np.random.uniform(0, 1000, 100000))
+        cls.rand_ener = np.random.uniform(0, 100, 100000)
+        cls.gti = [[0, 900], [950, 1000]]
+        blablas = np.random.normal(0, 1, cls.rand_ener.size)
+        cls.ev_like = StingrayTimeseries(
+            cls.rand_time, energy=cls.rand_ener, blablas=blablas, gti=cls.gti
+        )
+        time_edges = np.linspace(0, 1000, 1001)
+        counts = np.histogram(cls.rand_time, bins=time_edges)[0]
+        blablas = np.random.normal(0, 1, 1000)
+        cls.lc_like = StingrayTimeseries(
+            time=time_edges[:-1] + 0.5, counts=counts, blablas=blablas, gti=cls.gti, dt=1
+        )
+
+    def test_no_btis_returns_copy(self):
+        ts = StingrayTimeseries([1, 2, 3], energy=[4, 6, 8], gti=[[0.5, 3.5]])
+        ts_new = ts.fill_bad_time_intervals()
+        assert ts == ts_new
+
+    def test_event_like(self):
+        ev_like_filt = copy.deepcopy(self.ev_like)
+        # I introduce a small gap in the GTIs
+        ev_like_filt.gti = np.asarray([[0, 498], [500, 520], [522, 700], [702, 900], [950, 1000]])
+        ev_new = ev_like_filt.fill_bad_time_intervals()
+
+        assert np.allclose(ev_new.gti, self.gti)
+
+        # Now, I set the same GTIs as the original event list, and the data
+        # should be the same
+        ev_new.gti = ev_like_filt.gti
+
+        new_masked, filt_masked = ev_new.apply_gtis(), ev_like_filt.apply_gtis()
+        for attr in ["time", "energy", "blablas"]:
+            assert np.allclose(getattr(new_masked, attr), getattr(filt_masked, attr))
+
+    def test_lc_like(self):
+        lc_like_filt = copy.deepcopy(self.lc_like)
+        # I introduce a small gap in the GTIs
+        lc_like_filt.gti = np.asarray([[0, 498], [500, 520], [522, 700], [702, 900], [950, 1000]])
+        lc_new = lc_like_filt.fill_bad_time_intervals()
+        assert np.allclose(lc_new.gti, self.gti)
+
+        lc_like_gtifilt = self.lc_like.apply_gtis(inplace=False)
+        # In this case, the time array should also be the same as the original
+        assert np.allclose(lc_new.time, lc_like_gtifilt.time)
+
+        # Now, I set the same GTIs as the original event list, and the data
+        # should be the same
+        lc_new.gti = lc_like_filt.gti
+
+        new_masked, filt_masked = lc_new.apply_gtis(), lc_like_filt.apply_gtis()
+        for attr in ["time", "counts", "blablas"]:
+            assert np.allclose(getattr(new_masked, attr), getattr(filt_masked, attr))
+
+    def test_ignore_attrs_ev_like(self):
+        ev_like_filt = copy.deepcopy(self.ev_like)
+        # I introduce a small gap in the GTIs
+        ev_like_filt.gti = np.asarray([[0, 498], [500, 900], [950, 1000]])
+        ev_new0 = ev_like_filt.fill_bad_time_intervals(seed=1234)
+        ev_new1 = ev_like_filt.fill_bad_time_intervals(seed=1234, attrs_to_randomize=["energy"])
+        assert np.allclose(ev_new0.gti, ev_new1.gti)
+        assert np.allclose(ev_new0.time, ev_new1.time)
+
+        assert np.count_nonzero(np.isnan(ev_new0.blablas)) == 0
+        assert np.count_nonzero(np.isnan(ev_new1.blablas)) > 0
+        assert np.count_nonzero(np.isnan(ev_new1.energy)) == 0
+
+    def test_ignore_attrs_lc_like(self):
+        lc_like_filt = copy.deepcopy(self.lc_like)
+        # I introduce a small gap in the GTIs
+        lc_like_filt.gti = np.asarray([[0, 498], [500, 900], [950, 1000]])
+        lc_new0 = lc_like_filt.fill_bad_time_intervals(seed=1234)
+        lc_new1 = lc_like_filt.fill_bad_time_intervals(seed=1234, attrs_to_randomize=["counts"])
+        assert np.allclose(lc_new0.gti, lc_new1.gti)
+        assert np.allclose(lc_new0.time, lc_new1.time)
+
+        assert np.count_nonzero(np.isnan(lc_new0.blablas)) == 0
+        assert np.count_nonzero(np.isnan(lc_new1.blablas)) > 0
+        assert np.count_nonzero(np.isnan(lc_new1.counts)) == 0
+
+    def test_forcing_non_uniform(self):
+        ev_like_filt = copy.deepcopy(self.ev_like)
+        # I introduce a small gap in the GTIs
+        ev_like_filt.gti = np.asarray([[0, 498], [500, 900], [950, 1000]])
+        # Results should be exactly the same
+        ev_new0 = ev_like_filt.fill_bad_time_intervals(even_sampling=False, seed=201903)
+        ev_new1 = ev_like_filt.fill_bad_time_intervals(even_sampling=None, seed=201903)
+        for attr in ["time", "energy"]:
+            assert np.allclose(getattr(ev_new0, attr), getattr(ev_new1, attr))
+
+    def test_forcing_uniform(self):
+        lc_like_filt = copy.deepcopy(self.lc_like)
+        # I introduce a small gap in the GTIs
+        lc_like_filt.gti = np.asarray([[0, 498], [500, 900], [950, 1000]])
+        # Results should be exactly the same
+        lc_new0 = lc_like_filt.fill_bad_time_intervals(even_sampling=True, seed=201903)
+        lc_new1 = lc_like_filt.fill_bad_time_intervals(even_sampling=None, seed=201903)
+        for attr in ["time", "counts", "blablas"]:
+            assert np.allclose(getattr(lc_new0, attr), getattr(lc_new1, attr))
+
+    def test_bti_close_to_edge_event_like(self):
+        ev_like_filt = copy.deepcopy(self.ev_like)
+        # I introduce a small gap in the GTIs
+        ev_like_filt.gti = np.asarray([[0, 0.5], [1, 900], [950, 1000]])
+        ev_new = ev_like_filt.fill_bad_time_intervals()
+        assert np.allclose(ev_new.gti, self.gti)
+
+        ev_like_filt = copy.deepcopy(self.ev_like)
+        # I introduce a small gap in the GTIs
+        ev_like_filt.gti = np.asarray([[0, 900], [950, 999], [999.5, 1000]])
+        ev_new = ev_like_filt.fill_bad_time_intervals()
+        assert np.allclose(ev_new.gti, self.gti)
+
+    def test_bti_close_to_edge_lc_like(self):
+        lc_like_filt = copy.deepcopy(self.lc_like)
+        # I introduce a small gap in the GTIs
+        lc_like_filt.gti = np.asarray([[0, 0.5], [1, 900], [950, 1000]])
+        lc_new = lc_like_filt.fill_bad_time_intervals()
+        assert np.allclose(lc_new.gti, self.gti)
+
+        lc_like_filt = copy.deepcopy(self.lc_like)
+        # I introduce a small gap in the GTIs
+        lc_like_filt.gti = np.asarray([[0, 900], [950, 999], [999.5, 1000]])
+        lc_new = lc_like_filt.fill_bad_time_intervals()
+        assert np.allclose(lc_new.gti, self.gti)
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index 08ec13b86..5647a3d97 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -1408,6 +1408,30 @@ 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)
+        with pytest.warns(UserWarning, match="When using power_colors, complex "):
+            dps_ev.power_colors(freq_edges=[1 / 5, 1 / 2, 1, 2.0, 16.0])
+
+    def test_rms_is_correct(self):
+        lc = copy.deepcopy(self.lc)
+        lc.counts = np.random.poisson(lc.counts)
+        dps = DynamicalCrossspectrum(lc, lc, segment_size=10, norm="leahy")
+        rms, rmse = dps.compute_rms(1 / 5, 16.0, poisson_noise_level=2)
+        from stingray.powerspectrum import AveragedPowerspectrum
+
+        ps = AveragedPowerspectrum()
+        ps.freq = dps.freq
+        ps.power = dps.dyn_ps.T[0]
+        ps.unnorm_power = ps.power / dps.unnorm_conversion
+        ps.df = dps.df
+        ps.m = dps.m
+        ps.n = dps.freq.size
+        ps.dt = lc.dt
+        ps.norm = dps.norm
+        ps.k = 1
+        ps.nphots = (dps.nphots1 * dps.nphots2) ** 0.5
+        rms2, rmse2 = ps.compute_rms(1 / 5, 16.0, poisson_noise_level=2)
+        assert np.isclose(rms[0], rms2)
+        assert np.isclose(rmse[0], rmse2, rtol=0.01)
 
     def test_works_with_events_and_its_complex(self):
         lc = copy.deepcopy(self.lc)
@@ -1473,18 +1497,68 @@ def test_rebin_small_df(self):
         with pytest.raises(ValueError):
             dps.rebin_frequency(df_new=dps.df / 2.0)
 
-    def test_rebin_time_default_method(self):
+    def test_rebin_time_sum_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, method="sum")
+        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_n(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_by_n_intervals(n=2, method="sum")
+        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_n_average(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_by_n_intervals(n=2, method="average")
+        assert np.allclose(new_dps.time, rebin_time)
+        assert np.allclose(new_dps.dyn_ps, rebin_dps / 2)
+        assert np.isclose(new_dps.dt, dt_new)
+
+    def test_rebin_n_warns_for_non_integer(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)
+        with pytest.warns(UserWarning, match="n must be an integer. Casting to int"):
+            new_dps = dps.rebin_by_n_intervals(n=2.1, method="sum")
         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):
+    def test_rebin_n_fails_for_n_lt_1(self):
+        segment_size = 3
+
+        dps = DynamicalCrossspectrum(self.lc_test, self.lc_test, segment_size=segment_size)
+        with pytest.raises(ValueError, match="n must be >= 1"):
+            _ = dps.rebin_by_n_intervals(n=0)
+        dps2 = dps.rebin_by_n_intervals(n=1)
+        assert np.allclose(dps2.dyn_ps, dps.dyn_ps)
+
+    def test_rebin_n_copies_for_n_1(self):
+        segment_size = 3
+
+        dps = DynamicalCrossspectrum(self.lc_test, self.lc_test, segment_size=segment_size)
+        dps2 = dps.rebin_by_n_intervals(n=1)
+        assert np.allclose(dps2.dyn_ps, dps.dyn_ps)
+
+    def test_rebin_frequency_sum_method(self):
         segment_size = 50
         df_new = 10.0
         rebin_freq = np.array([5.01, 15.01, 25.01, 35.01])
@@ -1501,7 +1575,7 @@ def test_rebin_frequency_default_method(self):
         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)
+        new_dps = dps.rebin_frequency(df_new=df_new, method="sum")
         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)
diff --git a/stingray/tests/test_fourier.py b/stingray/tests/test_fourier.py
index 216ee49b0..91b8df6dc 100644
--- a/stingray/tests/test_fourier.py
+++ b/stingray/tests/test_fourier.py
@@ -293,7 +293,19 @@ def test_avg_cs_use_common_mean_similar_stats(self, norm):
             silent=True,
             return_subcs=True,
         )
+
         assert np.isclose(out_comm["power"].std(), out["power"].std(), rtol=0.1)
+        assert np.isclose(out_comm["unnorm_power"].std(), out["unnorm_power"].std(), rtol=0.1)
+        # Run the same check on the subcs
+        assert np.isclose(out_comm.meta["subcs"].std(), out.meta["subcs"].std(), rtol=0.1)
+        assert np.isclose(
+            out_comm.meta["unnorm_subcs"].std(), out.meta["unnorm_subcs"].std(), rtol=0.1
+        )
+        # Now verify that the normalizations are consistent between single power and subcs
+        assert np.isclose(
+            out_comm["unnorm_power"].std() / out_comm["power"].std(),
+            out_comm.meta["unnorm_subcs"].std() / out_comm.meta["subcs"].std(),
+        )
 
     @pytest.mark.parametrize("use_common_mean", [True, False])
     @pytest.mark.parametrize("norm", ["frac", "abs", "none", "leahy"])
@@ -663,3 +675,123 @@ def test_impose_symmetry_lsft():
     assert np.all((imp_sym_fast.real) == np.flip(imp_sym_fast.real))
     assert np.all((imp_sym_fast.imag) == (-np.flip(imp_sym_fast.imag)))
     assert np.all(freqs_new_slow == freqs_new_fast)
+
+
+class TestIntegration(object):
+    @classmethod
+    def setup_class(cls):
+        cls.freq = [0, 1, 2, 3]
+        cls.power = [2, 2, 2, 2]
+        cls.power_err = [1, 1, 1, 1]
+
+    def test_power_integration_middle_bin(self):
+        freq_range = [1, 2]
+        pow, powe = integrate_power_in_frequency_range(self.freq, self.power, freq_range)
+        assert np.isclose(pow, 2)
+        assert np.isclose(powe, np.sqrt(2))
+
+    def test_power_integration_precise(self):
+        freq_range = [0.5, 2.5]
+        df = 1
+        pow, powe = integrate_power_in_frequency_range(self.freq, self.power, freq_range, df=df)
+        assert np.allclose(pow, 4)
+        assert np.allclose(powe, 2 * np.sqrt(2))
+
+    def test_power_integration_poisson(self):
+        freq_range = [0.5, 2.5]
+        for poisson_power in (1, np.ones_like(self.power)):
+            pow, powe = integrate_power_in_frequency_range(
+                self.freq, self.power, freq_range, poisson_power=poisson_power
+            )
+            assert np.allclose(pow, 2)
+            assert np.allclose(powe, 2 * np.sqrt(2))
+
+    def test_power_integration_err(self):
+        freq_range = [0.5, 2.5]
+        pow, powe = integrate_power_in_frequency_range(
+            self.freq, self.power, freq_range, power_err=self.power_err
+        )
+        assert np.allclose(pow, 4)
+        assert np.allclose(powe, np.sqrt(2))
+
+    def test_power_integration_m(self):
+        freq_range = [0.5, 2.5]
+        pow, powe = integrate_power_in_frequency_range(self.freq, self.power, freq_range, m=4)
+        assert np.allclose(pow, 4)
+        assert np.allclose(powe, np.sqrt(2))
+
+
+class TestPowerColor(object):
+    @classmethod
+    def setup_class(cls):
+        cls.freq = np.arange(0.0001, 17, 0.00001)
+        cls.power = 1 / cls.freq
+
+    def test_power_color(self):
+        pc0, _, pc1, _ = power_color(self.freq, self.power)
+        # The colors calculated with these frequency edges on a 1/f spectrum should be 1
+        assert np.isclose(pc0, 1)
+        assert np.isclose(pc1, 1)
+
+    def test_return_log(self):
+        pc0, _, pc1, _ = power_color(self.freq, self.power, return_log=True)
+        # The colors calculated with these frequency edges on a 1/f spectrum should be 1
+        assert np.isclose(pc0, 0, atol=0.001)
+        assert np.isclose(pc1, 0, atol=0.001)
+
+    def test_bad_edges(self):
+        good = self.freq > 1 / 255  # the smallest frequency is 1/256
+        with pytest.raises(ValueError, match="The minimum frequency is larger "):
+            power_color(self.freq[good], self.power[good])
+
+        good = self.freq < 15  # the smallest frequency is 1/256
+        with pytest.raises(ValueError, match="The maximum frequency is lower "):
+            power_color(self.freq[good], self.power[good])
+
+        with pytest.raises(ValueError, match="freq_edges must have 5 elements"):
+            power_color(self.freq, self.power, freq_edges=[1])
+        with pytest.raises(ValueError, match="freq_edges must have 5 elements"):
+            power_color(self.freq, self.power, freq_edges=[1, 2, 3, 4, 5, 6])
+
+    def test_bad_excluded_interval(self):
+        for fte in ([1, 1.1, 3.0], [4], [[1, 1.1, 3.0]], 0, [[[1, 3]]]):
+            with pytest.raises(ValueError, match="freqs_to_exclude must be of "):
+                power_color(self.freq, self.power, freqs_to_exclude=fte)
+
+    def test_excluded_frequencies(self):
+        pc0, _, pc1, _ = power_color(self.freq, self.power, freqs_to_exclude=[1, 1.1])
+        # The colors calculated with these frequency edges on a 1/f spectrum should be 1
+        # The excluded frequency interval is small enough that the approximation should work
+        assert np.isclose(pc0, 1, atol=0.001)
+        assert np.isclose(pc1, 1, atol=0.001)
+
+    def test_with_power_err(self):
+        pc0, pc0_err, pc1, pc1_err = power_color(
+            self.freq,
+            self.power,
+            power_err=self.power / 2,
+        )
+        pc0e, pc0e_err, pc1e, pc1e_err = power_color(
+            self.freq,
+            self.power,
+            power_err=self.power,
+        )
+        assert np.isclose(pc0, 1, atol=0.001)
+        assert np.isclose(pc1, 1, atol=0.001)
+        assert np.isclose(pc0e, 1, atol=0.001)
+        assert np.isclose(pc1e, 1, atol=0.001)
+        assert np.isclose(pc0e_err / pc0_err, 2, atol=0.001)
+        assert np.isclose(pc1e_err / pc1_err, 2, atol=0.001)
+
+    def test_hue(self):
+        center = (4.51920, 0.453724)
+        log_center = np.log10(np.asarray(center))
+        for angle in np.radians(np.arange(0, 380, 20)):
+            factor = np.random.uniform(0.1, 10)
+            x = factor * np.cos(3 / 4 * np.pi - angle) + log_center[0]
+            y = factor * np.sin(3 / 4 * np.pi - angle) + log_center[1]
+            hue = hue_from_power_color(10**x, 10**y, center)
+            # Compare the angles in a safe way
+            c2 = (np.sin(hue) - np.sin(angle)) ** 2 + (np.cos(hue) - np.cos(angle)) ** 2
+            angle_diff = np.arccos((2.0 - c2) / 2.0)
+            assert np.isclose(angle_diff, 0, atol=0.001)
diff --git a/stingray/tests/test_io.py b/stingray/tests/test_io.py
index b67ae2754..86ef6c6e7 100644
--- a/stingray/tests/test_io.py
+++ b/stingray/tests/test_io.py
@@ -71,7 +71,7 @@ def test_read_whole_mission_info(self):
     def test_event_file_read_and_automatic_sort(self):
         """Test event file reading."""
         fname = os.path.join(datadir, "monol_testA_calib.evt")
-        with pytest.warns(AstropyUserWarning, match="No extensions found with a"):
+        with pytest.warns(AstropyUserWarning, match="No valid GTI extensions"):
             evdata = load_events_and_gtis(fname)
         fname_unsrt = os.path.join(datadir, "monol_testA_calib_unsrt.evt")
         with pytest.warns(UserWarning, match="not sorted. Sorting them for you"):
@@ -93,7 +93,7 @@ def test_event_file_read_additional_warns_uncal(self):
     def test_event_file_read_additional_energy_cal(self):
         """Test event file reading."""
         fname = os.path.join(datadir, "monol_testA_calib.evt")
-        with pytest.warns(AstropyUserWarning, match="No extensions found with a"):
+        with pytest.warns(AstropyUserWarning, match="No valid GTI extensions"):
             vals = load_events_and_gtis(fname, additional_columns=["energy"])
         # These energies were calibrated with a different calibration than
         # returned from rough_calibration, on purpose! (notice the +1.)
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index 68e03d71d..6b770bbd3 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -1148,9 +1148,10 @@ def test_plot_simple(self):
     def test_plot_wrong_label_type(self):
         lc = Lightcurve(self.times, self.counts)
 
-        with pytest.raises(TypeError):
-            with pytest.warns(UserWarning, match="must be either a list or tuple") as w:
-                lc.plot(labels=123)
+        with pytest.warns(
+            UserWarning, match="``labels`` must be an iterable with two labels "
+        ) as w:
+            lc.plot(labels=123)
         plt.close("all")
 
     def test_plot_labels_index_error(self):
@@ -1158,7 +1159,9 @@ def test_plot_labels_index_error(self):
         with pytest.warns(UserWarning) as w:
             lc.plot(labels=("x"))
 
-            assert np.any(["must have two labels" in str(wi.message) for wi in w])
+            assert np.any(
+                ["``labels`` must be an iterable with two labels " in str(wi.message) for wi in w]
+            )
         plt.close("all")
 
     def test_plot_default_filename(self):
@@ -1175,14 +1178,19 @@ def test_plot_custom_filename(self):
         os.unlink("lc.png")
         plt.close("all")
 
-    def test_plot_axis(self):
+    def test_plot_axis_arg(self):
         lc = Lightcurve(self.times, self.counts)
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore", category=UserWarning)
+        with pytest.warns(DeprecationWarning, match="argument is deprecated in favor"):
             lc.plot(axis=[0, 1, 0, 100])
         assert plt.fignum_exists(1)
         plt.close("all")
 
+    def test_plot_axis_limits_arg(self):
+        lc = Lightcurve(self.times, self.counts)
+        lc.plot(axis_limits=[0, 1, 0, 100])
+        assert plt.fignum_exists(1)
+        plt.close("all")
+
     def test_plot_title(self):
         lc = Lightcurve(self.times, self.counts)
         with warnings.catch_warnings():
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index 45561fa3f..a43d4c20d 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -1034,6 +1034,29 @@ def test_works_with_events(self):
 
         dps_ev = DynamicalPowerspectrum(ev, segment_size=10, sample_time=self.lc.dt)
         assert np.allclose(dps.dyn_ps, dps_ev.dyn_ps)
+        dps_ev.power_colors(freq_edges=[1 / 5, 1 / 2, 1, 2.0, 16.0])
+
+    def test_rms_is_correct(self):
+        lc = copy.deepcopy(self.lc)
+        lc.counts = np.random.poisson(lc.counts)
+        dps = DynamicalPowerspectrum(lc, segment_size=10, norm="leahy")
+        rms, rmse = dps.compute_rms(1 / 5, 16.0, poisson_noise_level=2)
+        from stingray.powerspectrum import AveragedPowerspectrum
+
+        ps = AveragedPowerspectrum()
+        ps.freq = dps.freq
+        ps.power = dps.dyn_ps.T[0]
+        ps.unnorm_power = ps.power / dps.unnorm_conversion
+        ps.df = dps.df
+        ps.m = dps.m
+        ps.n = dps.freq.size
+        ps.dt = lc.dt
+        ps.norm = dps.norm
+        ps.k = 1
+        ps.nphots = dps.nphots
+        rms2, rmse2 = ps.compute_rms(1 / 5, 16.0, poisson_noise_level=2)
+        assert np.isclose(rms[0], rms2)
+        assert np.isclose(rmse[0], rmse2, rtol=0.01)
 
     def test_with_long_seg_size(self):
         with pytest.raises(ValueError):
@@ -1087,13 +1110,13 @@ def test_rebin_small_df(self):
         with pytest.raises(ValueError):
             dps.rebin_frequency(df_new=dps.df / 2.0)
 
-    def test_rebin_time_default_method(self):
+    def test_rebin_time_sum_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 = DynamicalPowerspectrum(self.lc_test, segment_size=segment_size)
-        new_dps = dps.rebin_time(dt_new=dt_new)
+        new_dps = dps.rebin_time(dt_new=dt_new, method="sum")
         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)
@@ -1115,7 +1138,7 @@ def test_rebin_frequency_default_method(self):
         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)
+        new_dps = dps.rebin_frequency(df_new=df_new, method="sum")
         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)
diff --git a/stingray/utils.py b/stingray/utils.py
index ce0adb298..7f0b16312 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -1299,7 +1299,7 @@ def nearest_power_of_two(x):
     return x_nearest
 
 
-def find_nearest(array, value):
+def find_nearest(array, value, side="left"):
     """
     Return the array value that is closest to the input value (Abigail Stevens:
     Thanks StackOverflow!)
@@ -1313,6 +1313,11 @@ def find_nearest(array, value):
     value : int or float
         The value you want to find the closest to in the array.
 
+    Other Parameters
+    ----------------
+    side : str
+        Look at the ``numpy.searchsorted`` documentation for more information.
+
     Returns
     -------
     array[idx] : int or float
@@ -1322,7 +1327,7 @@ def find_nearest(array, value):
         The index of the array of the closest value.
 
     """
-    idx = np.searchsorted(array, value, side="left")
+    idx = np.searchsorted(array, value, side=side)
     if idx == len(array) or np.fabs(value - array[idx - 1]) < np.fabs(value - array[idx]):
         return array[idx - 1], idx - 1
     else: