diff --git a/docs/changes/865.feature.rst b/docs/changes/865.feature.rst
new file mode 100644
index 000000000..79c571338
--- /dev/null
+++ b/docs/changes/865.feature.rst
@@ -0,0 +1 @@
+Use the formulas from Ingram+2019 consistently when calculating cross spectral errors, including lags
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 5673fbcd8..74a6cf705 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -503,6 +503,11 @@ class Crossspectrum(StingrayObject):
         Skip initial checks, for speed or other reasons (you need to trust your
         inputs!)
 
+    channels_overlap: bool, default False
+        If True, the reference band contains all the photons of the subject band.
+        This happens, for example, when calculating covariance spectra (see, e.g.,
+        the docs for ``CovarianceSpectrum``). This will generally be false in a
+        single cross spectrum.
 
     Attributes
     ----------
@@ -553,6 +558,7 @@ def __init__(
         fullspec=False,
         skip_checks=False,
         save_all=False,
+        channels_overlap=False,
     ):
         self._type = None
         # for backwards compatibility
@@ -560,7 +566,7 @@ def __init__(
             data1 = lc1
         if data2 is None:
             data2 = lc2
-
+        self.channels_overlap = channels_overlap
         empty = data1 is None and data2 is None
         good_input = not empty
 
@@ -1219,6 +1225,7 @@ def from_time_array(
         silent=False,
         fullspec=False,
         use_common_mean=True,
+        channels_overlap=False,
     ):
         """Calculate AveragedCrossspectrum from two arrays of event times.
 
@@ -1260,6 +1267,11 @@ def from_time_array(
         power_type : str, default 'all'
             If 'all', give complex powers. If 'abs', the absolute value; if 'real',
             the real part
+        channels_overlap: bool, default False
+            If True, the reference band contains all the photons of the subject band.
+            This happens, for example, when calculating covariance spectra (see, e.g.,
+            the docs for ``CovarianceSpectrum``). This will generally be false in a
+            single cross spectrum.
         """
 
         return crossspectrum_from_time_array(
@@ -1273,6 +1285,7 @@ def from_time_array(
             silent=silent,
             fullspec=fullspec,
             use_common_mean=use_common_mean,
+            channels_overlap=channels_overlap,
         )
 
     @staticmethod
@@ -1287,6 +1300,7 @@ def from_events(
         fullspec=False,
         use_common_mean=True,
         gti=None,
+        channels_overlap=False,
     ):
         """Calculate AveragedCrossspectrum from two event lists
 
@@ -1328,6 +1342,11 @@ def from_events(
             input object GTIs! If you're getting errors regarding your GTIs,
             don't use this and only give GTIs to the input objects before
             making the cross spectrum.
+        channels_overlap: bool, default False
+            If True, the reference band contains all the photons of the subject band.
+            This happens, for example, when calculating covariance spectra (see, e.g.,
+            the docs for ``CovarianceSpectrum``). This will generally be false in a
+            single cross spectrum.
         """
 
         return crossspectrum_from_events(
@@ -1341,6 +1360,7 @@ def from_events(
             fullspec=fullspec,
             use_common_mean=use_common_mean,
             gti=gti,
+            channels_overlap=channels_overlap,
         )
 
     @staticmethod
@@ -1354,6 +1374,7 @@ def from_lightcurve(
         fullspec=False,
         use_common_mean=True,
         gti=None,
+        channels_overlap=False,
     ):
         """Calculate AveragedCrossspectrum from two light curves
 
@@ -1392,6 +1413,11 @@ def from_lightcurve(
             input object GTIs! If you're getting errors regarding your GTIs,
             don't  use this and only give GTIs to the input objects before
             making the cross spectrum.
+        channels_overlap: bool, default False
+            If True, the reference band contains all the photons of the subject band.
+            This happens, for example, when calculating covariance spectra (see, e.g.,
+            the docs for ``CovarianceSpectrum``). This will generally be false in a
+            single cross spectrum.
         """
         return crossspectrum_from_lightcurve(
             lc1,
@@ -1403,6 +1429,7 @@ def from_lightcurve(
             fullspec=fullspec,
             use_common_mean=use_common_mean,
             gti=gti,
+            channels_overlap=channels_overlap,
         )
 
     @staticmethod
@@ -1418,6 +1445,7 @@ def from_stingray_timeseries(
         fullspec=False,
         use_common_mean=True,
         gti=None,
+        channels_overlap=False,
     ):
         """Calculate AveragedCrossspectrum from two light curves
 
@@ -1460,6 +1488,11 @@ def from_stingray_timeseries(
             input object GTIs! If you're getting errors regarding your GTIs,
             don't  use this and only give GTIs to the input objects before
             making the cross spectrum.
+        channels_overlap: bool, default False
+            If True, the reference band contains all the photons of the subject band.
+            This happens, for example, when calculating covariance spectra (see, e.g.,
+            the docs for ``CovarianceSpectrum``). This will generally be false in a
+            single cross spectrum.
         """
         return crossspectrum_from_timeseries(
             ts1,
@@ -1473,6 +1506,7 @@ def from_stingray_timeseries(
             fullspec=fullspec,
             use_common_mean=use_common_mean,
             gti=gti,
+            channels_overlap=channels_overlap,
         )
 
     @staticmethod
@@ -1487,6 +1521,7 @@ def from_lc_iterable(
         fullspec=False,
         use_common_mean=True,
         gti=None,
+        channels_overlap=False,
     ):
         """Calculate AveragedCrossspectrum from two light curves
 
@@ -1531,6 +1566,11 @@ def from_lc_iterable(
         save_all : bool, default False
             If True, save the cross spectrum of each segment in the ``cs_all``
             attribute of the output :class:`Crossspectrum` object.
+        channels_overlap: bool, default False
+            If True, the reference band contains all the photons of the subject band.
+            This happens, for example, when calculating covariance spectra (see, e.g.,
+            the docs for ``CovarianceSpectrum``). This will generally be false in a
+            single cross spectrum.
         """
 
         return crossspectrum_from_lc_iterable(
@@ -1544,6 +1584,7 @@ def from_lc_iterable(
             fullspec=fullspec,
             use_common_mean=use_common_mean,
             gti=gti,
+            channels_overlap=channels_overlap,
         )
 
     def _initialize_from_any_input(
@@ -1601,6 +1642,7 @@ def _initialize_from_any_input(
                 use_common_mean=use_common_mean,
                 gti=gti,
                 save_all=save_all,
+                channels_overlap=self.channels_overlap,
             )
         elif isinstance(data1, Lightcurve):
             spec = crossspectrum_from_lightcurve(
@@ -1614,6 +1656,7 @@ def _initialize_from_any_input(
                 use_common_mean=use_common_mean,
                 gti=gti,
                 save_all=save_all,
+                channels_overlap=self.channels_overlap,
             )
             spec.lc1 = data1
             spec.lc2 = data2
@@ -1632,6 +1675,7 @@ def _initialize_from_any_input(
                 gti=gti,
                 use_common_mean=use_common_mean,
                 save_all=save_all,
+                channels_overlap=self.channels_overlap,
             )
         else:  # pragma: no cover
             raise TypeError(f"Bad inputs to Crosssspectrum: {type(data1)}")
@@ -1736,6 +1780,67 @@ def deadtime_correct(
 
         return new_spec
 
+    def _calculate_errors(self):
+        """Calculate the errors on cross powers and lags.
+
+        Uses different formulas if the reference band contains the photons of the subject band.
+        This might happen, for example, when calculating covariance spectra using a large
+        reference band.
+        See the details in the documentation of
+        :func:`stingray.fourier.error_on_averaged_cross_spectrum`.
+
+        Please note that we have dedicated methods for covariance spectra and other variability
+        versus energy spectra in `stingray.varenergyspectrum`, even though they only work for
+        input event lists at the moment.
+        """
+        P1noise = poisson_level(norm="none", meanrate=self.countrate1, n_ph=self.nphots1)
+        P2noise = poisson_level(norm="none", meanrate=self.countrate2, n_ph=self.nphots2)
+
+        dRe, dIm, dphi, _ = error_on_averaged_cross_spectrum(
+            self.unnorm_power,
+            self.pds1.unnorm_power,
+            self.pds2.unnorm_power,
+            self.m,
+            P1noise,
+            P2noise,
+            common_ref=self.channels_overlap,
+        )
+        bad = np.isnan(dRe) | np.isnan(dIm)
+
+        if np.any(bad):
+            warnings.warn(
+                "Some error bars in the Averaged Crossspectrum are invalid."
+                "Defaulting to sqrt(2 / M) in Leahy norm, rescaled to the appropriate norm."
+            )
+
+        Nph = np.sqrt(self.nphots1 * self.nphots2)
+
+        assert self.nph == Nph
+
+        default_err = np.sqrt(2 / self.m) * Nph / 2
+        if isinstance(default_err, Iterable):
+            default_err = np.array(default_err)[bad]
+
+        dRe[bad] = default_err
+        dIm[bad] = default_err
+
+        power_err = dRe + 1.0j * dIm
+
+        self.unnorm_power_err = power_err
+
+        self.power_err = normalize_periodograms(
+            power_err, self.dt, self.n, self.mean, n_ph=Nph, variance=self.variance, norm=self.norm
+        )
+
+        self.pds1.power_err = self.pds1.power / np.sqrt(self.pds1.m)
+        self.pds2.power_err = self.pds2.power / np.sqrt(self.pds2.m)
+        self.pds1.unnorm_power_err = self.pds1.unnorm_power / np.sqrt(self.pds1.m)
+        self.pds2.unnorm_power_err = self.pds2.unnorm_power / np.sqrt(self.pds2.m)
+        self.lag_err = dphi
+
+        assert hasattr(self, "df")
+        assert hasattr(self, "dt")
+
 
 class AveragedCrossspectrum(Crossspectrum):
     type = "crossspectrum"
@@ -1821,6 +1926,12 @@ class AveragedCrossspectrum(Crossspectrum):
         don't  use this and only give GTIs to the input objects before
         making the cross spectrum.
 
+    channels_overlap: bool, default False
+        If True, the reference band contains all the photons of the subject band.
+        This happens, for example, when calculating covariance spectra (see, e.g.,
+        the docs for ``CovarianceSpectrum``). This will generally be false in a
+        single cross spectrum.
+
     Attributes
     ----------
     freq: numpy.ndarray
@@ -1872,6 +1983,7 @@ def __init__(
         save_all=False,
         use_common_mean=True,
         skip_checks=False,
+        channels_overlap=False,
     ):
         self._type = None
         # for backwards compatibility
@@ -1902,6 +2014,7 @@ def __init__(
         self.save_all = save_all
         self.segment_size = segment_size
         self.show_progress = not silent
+        self.channels_overlap = channels_overlap
 
         if empty or not good_input:
             return self._initialize_empty()
@@ -2021,14 +2134,9 @@ def coherence(self):
     def phase_lag(self):
         """Return the fourier phase lag of the cross spectrum."""
         lag = np.angle(self.unnorm_power)
-        coh, uncert = self.coherence()
-
-        dum = (1.0 - coh) / (2.0 * coh)
-
-        dum[coh == 0] = 0.0
-
-        lag_err = np.sqrt(dum / self.m)
-        return lag, lag_err
+        if not hasattr(self, "lag_err") or self.lag_err.size != lag.size:
+            self._calculate_errors()
+        return lag, self.lag_err
 
     def time_lag(self):
         """Calculate time lag and uncertainty.
@@ -2577,6 +2685,7 @@ def crossspectrum_from_time_array(
     fullspec=False,
     use_common_mean=True,
     save_all=False,
+    channels_overlap=False,
 ):
     """Calculate AveragedCrossspectrum from two arrays of event times.
 
@@ -2614,6 +2723,11 @@ def crossspectrum_from_time_array(
     power_type : str, default 'all'
         If 'all', give complex powers. If 'abs', the absolute value; if 'real',
         the real part
+    channels_overlap: bool, default False
+        If True, the reference band contains all the photons of the subject band.
+        This happens, for example, when calculating covariance spectra (see, e.g.,
+        the docs for ``CovarianceSpectrum``). This will generally be false in a
+        single cross spectrum.
 
     Returns
     -------
@@ -2638,7 +2752,9 @@ def crossspectrum_from_time_array(
         return_subcs=save_all,
     )
 
-    return _create_crossspectrum_from_result_table(results, force_averaged=force_averaged)
+    return _create_crossspectrum_from_result_table(
+        results, force_averaged=force_averaged, channels_overlap=channels_overlap
+    )
 
 
 def crossspectrum_from_events(
@@ -2653,6 +2769,7 @@ def crossspectrum_from_events(
     use_common_mean=True,
     gti=None,
     save_all=False,
+    channels_overlap=False,
 ):
     """Calculate AveragedCrossspectrum from two event lists
 
@@ -2690,6 +2807,11 @@ def crossspectrum_from_events(
     gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
         Good Time intervals. Defaults to the common GTIs from the two input
         objects
+    channels_overlap: bool, default False
+        If True, the reference band contains all the photons of the subject band.
+        This happens, for example, when calculating covariance spectra (see, e.g.,
+        the docs for ``CovarianceSpectrum``). This will generally be false in a
+        single cross spectrum.
 
     Returns
     -------
@@ -2712,6 +2834,7 @@ def crossspectrum_from_events(
         fullspec=fullspec,
         use_common_mean=use_common_mean,
         save_all=save_all,
+        channels_overlap=channels_overlap,
     )
 
 
@@ -2726,6 +2849,7 @@ def crossspectrum_from_lightcurve(
     use_common_mean=True,
     gti=None,
     save_all=False,
+    channels_overlap=False,
 ):
     """Calculate AveragedCrossspectrum from two light curves
 
@@ -2764,6 +2888,11 @@ def crossspectrum_from_lightcurve(
         Save all intermediate spectra used for the final average. Use with care.
         This is likely to fill up your RAM on medium-sized datasets, and to
         slow down the computation when rebinning.
+    channels_overlap: bool, default False
+        If True, the reference band contains all the photons of the subject band.
+        This happens, for example, when calculating covariance spectra (see, e.g.,
+        the docs for ``CovarianceSpectrum``). This will generally be false in a
+        single cross spectrum.
 
     Returns
     -------
@@ -2788,6 +2917,7 @@ def crossspectrum_from_lightcurve(
         use_common_mean=use_common_mean,
         gti=gti,
         save_all=save_all,
+        channels_overlap=channels_overlap,
     )
 
 
@@ -2804,6 +2934,7 @@ def crossspectrum_from_timeseries(
     use_common_mean=True,
     gti=None,
     save_all=False,
+    channels_overlap=False,
 ):
     """Calculate AveragedCrossspectrum from two time series
 
@@ -2846,6 +2977,11 @@ def crossspectrum_from_timeseries(
         Save all intermediate spectra used for the final average. Use with care.
         This is likely to fill up your RAM on medium-sized datasets, and to
         slow down the computation when rebinning.
+    channels_overlap: bool, default False
+        If True, the reference band contains all the photons of the subject band.
+        This happens, for example, when calculating covariance spectra (see, e.g.,
+        the docs for ``CovarianceSpectrum``). This will generally be false in a
+        single cross spectrum.
 
     Returns
     -------
@@ -2882,7 +3018,9 @@ def crossspectrum_from_timeseries(
         return_subcs=save_all,
     )
 
-    return _create_crossspectrum_from_result_table(results, force_averaged=force_averaged)
+    return _create_crossspectrum_from_result_table(
+        results, force_averaged=force_averaged, channels_overlap=channels_overlap
+    )
 
 
 def crossspectrum_from_lc_iterable(
@@ -2897,6 +3035,7 @@ def crossspectrum_from_lc_iterable(
     use_common_mean=True,
     gti=None,
     save_all=False,
+    channels_overlap=False,
 ):
     """Calculate AveragedCrossspectrum from two light curves
 
@@ -2934,6 +3073,11 @@ def crossspectrum_from_lc_iterable(
     gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
         Good Time intervals. The GTIs of the input light curves are
         intersected with these.
+    channels_overlap: bool, default False
+        If True, the reference band contains all the photons of the subject band.
+        This happens, for example, when calculating covariance spectra (see, e.g.,
+        the docs for ``CovarianceSpectrum``). This will generally be false in a
+        single cross spectrum.
 
     Returns
     -------
@@ -2980,10 +3124,12 @@ def iterate_lc_counts(iter_lc):
         return_auxil=True,
         return_subcs=save_all,
     )
-    return _create_crossspectrum_from_result_table(results, force_averaged=force_averaged)
+    return _create_crossspectrum_from_result_table(
+        results, force_averaged=force_averaged, channels_overlap=channels_overlap
+    )
 
 
-def _create_crossspectrum_from_result_table(table, force_averaged=False):
+def _create_crossspectrum_from_result_table(table, force_averaged=False, channels_overlap=False):
     """Copy the columns and metadata from the results of
     ``stingray.fourier.avg_cs_from_XX`` functions into
     `AveragedCrossspectrum` or `Crossspectrum` objects.
@@ -3002,6 +3148,12 @@ def _create_crossspectrum_from_result_table(table, force_averaged=False):
     ----------------
     force_averaged : bool, default False
 
+    channels_overlap: bool, default False
+        If True, the reference band contains all the photons of the subject band.
+        This happens, for example, when calculating covariance spectra (see, e.g.,
+        the docs for ``CovarianceSpectrum``). This will generally be false in a
+        single cross spectrum.
+
     Returns
     -------
     spec : `AveragedCrossspectrum` or `Crossspectrum`
@@ -3049,49 +3201,11 @@ def _create_crossspectrum_from_result_table(table, force_averaged=False):
         if attr.endswith("2"):
             setattr(cs.pds2, attr[:-1], val)
 
-    # I start from unnormalized, and I normalize after correcting for bad error values
-    P1noise = poisson_level(norm="none", meanrate=cs.countrate1, n_ph=cs.nphots1)
-    P2noise = poisson_level(norm="none", meanrate=cs.countrate2, n_ph=cs.nphots2)
-
-    dRe, dIm, _, _ = error_on_averaged_cross_spectrum(
-        cs.unnorm_power,
-        cs.pds1.unnorm_power,
-        cs.pds2.unnorm_power,
-        cs.m,
-        P1noise,
-        P2noise,
-        common_ref=False,
-    )
-
-    bad = np.isnan(dRe) | np.isnan(dIm)
-
-    if np.any(bad):
-        warnings.warn(
-            "Some error bars in the Averaged Crossspectrum are invalid."
-            "Defaulting to sqrt(2 / M) in Leahy norm, rescaled to the appropriate norm."
-        )
-
-    Nph = np.sqrt(cs.nphots1 * cs.nphots2)
-    default_err = np.sqrt(2 / cs.m) * Nph / 2
-
-    dRe[bad] = default_err
-    dIm[bad] = default_err
-
-    power_err = dRe + 1.0j * dIm
-
-    cs.unnorm_power_err = power_err
-
     mean = table.meta["mean"]
     nph = table.meta["nphots"]
-    cs.power_err = normalize_periodograms(
-        power_err, cs.dt, cs.n, mean, n_ph=nph, variance=cs.variance, norm=cs.norm
-    )
-
-    cs.pds1.power_err = cs.pds1.power / np.sqrt(cs.pds1.m)
-    cs.pds2.power_err = cs.pds2.power / np.sqrt(cs.pds2.m)
-    cs.pds1.unnorm_power_err = cs.pds1.unnorm_power / np.sqrt(cs.pds1.m)
-    cs.pds2.unnorm_power_err = cs.pds2.unnorm_power / np.sqrt(cs.pds2.m)
+    cs.mean = mean
+    cs.channels_overlap = channels_overlap
 
-    assert hasattr(cs, "df")
-    assert hasattr(cs, "dt")
+    cs.nph = nph
+    cs._calculate_errors()
     return cs
diff --git a/stingray/fourier.py b/stingray/fourier.py
index 3e88527c9..8e14dc00f 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -1084,7 +1084,14 @@ def error_on_averaged_cross_spectrum(
     cross_power, seg_power, ref_power, n_ave, seg_power_noise, ref_power_noise, common_ref=False
 ):
     """
-    Error on cross spectral quantities, From Ingram 2019.
+    Error on cross spectral quantities.
+
+    `Ingram 2019<https://ui.adsabs.harvard.edu/abs/2019MNRAS.489.3927I/abstract>`__ details
+    the derivation of the corrected formulas when the subject band is included in the reference
+    band (and photons are presents in both bands, so this does _not_ involve overlapping bands from
+    different detectors).
+    The keyword argument ``common_ref`` indicates if the reference band also contains the photons
+    of the subject band.
 
     Note: this is only valid for a very large number of averaged powers.
     Beware if n_ave < 50 or so.
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index b12caeadf..4b5b9dbe8 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -1676,3 +1676,186 @@ def test_shift_and_add(self):
         assert np.array_equal(output.m, [2, 3, 3, 3, 2])
         assert np.array_equal(output.power, [2.0, 2.0, 5.0, 2.0, 1.5])
         assert np.allclose(output.freq, [0.05, 0.15, 0.25, 0.35, 0.45])
+
+
+class TestAveragedCrossspectrumOverlap(object):
+    def setup_class(self):
+        tstart = 0.0
+        tend = 1.0
+        dt = np.longdouble(0.0001)
+
+        time = np.arange(tstart + 0.5 * dt, tend + 0.5 * dt, dt)
+
+        counts1 = np.random.poisson(1, size=time.shape[0])
+        counts2 = np.random.poisson(1, size=time.shape[0]) + counts1
+
+        self.lc1 = Lightcurve(time, counts1, gti=[[tstart, tend]], dt=dt)
+        self.lc2 = Lightcurve(time, counts2, gti=[[tstart, tend]], dt=dt)
+
+        self.cs = AveragedCrossspectrum(
+            self.lc1, self.lc2, segment_size=1, save_all=True, channels_overlap=True
+        )
+
+    def test_save_all(self):
+        cs = AveragedCrossspectrum(
+            self.lc1, self.lc2, segment_size=1, save_all=True, channels_overlap=True
+        )
+        assert hasattr(self.cs, "cs_all")
+
+    def test_rebin_with_valid_type_attribute(self):
+        new_df = 2
+        aps = AveragedCrossspectrum(
+            self.lc1, self.lc2, segment_size=1, norm="leahy", channels_overlap=True
+        )
+
+        assert aps.rebin(df=new_df)
+
+    @pytest.mark.parametrize("err_dist", ["poisson", "gauss"])
+    def test_with_iterable_of_lightcurves(self, err_dist):
+        def iter_lc(lc, n):
+            "Generator of n parts of lc."
+            t0 = int(len(lc) / n)
+            t = t0
+            i = 0
+            while True:
+                lc_seg = lc[i:t]
+                yield lc_seg
+                if t + t0 > len(lc):
+                    break
+                else:
+                    i, t = t, t + t0
+
+        lc1 = copy.deepcopy(self.lc1)
+        lc2 = copy.deepcopy(self.lc2)
+        lc1.err_dist = lc2.err_dist = err_dist
+        with pytest.warns(UserWarning) as record:
+            cs = AveragedCrossspectrum(
+                iter_lc(self.lc1, 1), iter_lc(self.lc2, 1), segment_size=1, channels_overlap=True
+            )
+        message = "The averaged Cross spectrum from a generator "
+
+        assert np.any([message in r.message.args[0] for r in record])
+
+    def test_with_multiple_lightcurves_variable_length(self):
+        gti = [[0, 0.05], [0.05, 0.5], [0.555, 1.0]]
+        lc1 = copy.deepcopy(self.lc1)
+        lc1.gti = gti
+        lc2 = copy.deepcopy(self.lc2)
+        lc2.gti = gti
+
+        lc1_split = lc1.split_by_gti()
+        lc2_split = lc2.split_by_gti()
+
+        cs = AveragedCrossspectrum(
+            lc1_split,
+            lc2_split,
+            segment_size=0.05,
+            norm="leahy",
+            silent=True,
+            channels_overlap=True,
+        )
+
+    def test_coherence(self):
+        with pytest.warns(UserWarning) as w:
+            coh = self.cs.coherence()
+
+            assert len(coh[0]) == 4999
+            assert len(coh[1]) == 4999
+            assert issubclass(w[-1].category, UserWarning)
+
+    def test_normalize_crossspectrum(self):
+        cs1 = Crossspectrum(self.lc1, self.lc2, norm="leahy", channels_overlap=True)
+        cs2 = Crossspectrum(
+            self.lc1, self.lc2, norm="leahy", power_type="all", channels_overlap=True
+        )
+        cs3 = Crossspectrum(
+            self.lc1, self.lc2, norm="leahy", power_type="real", channels_overlap=True
+        )
+        cs4 = Crossspectrum(
+            self.lc1, self.lc2, norm="leahy", power_type="absolute", channels_overlap=True
+        )
+        assert np.allclose(cs1.power.real, cs3.power)
+        assert np.all(np.isclose(np.abs(cs2.power), cs4.power, atol=0.0001))
+
+    def test_normalize_crossspectrum_with_method_inplace(self):
+        cs1 = AveragedCrossspectrum.from_lightcurve(
+            self.lc1, self.lc2, segment_size=1, norm="abs", channels_overlap=True
+        )
+        cs2 = cs1.to_norm("leahy", inplace=True)
+        cs3 = cs1.to_norm("leahy", inplace=False)
+        assert cs3 is not cs1
+        assert cs2 is cs1
+
+    @pytest.mark.parametrize("norm1", ["leahy", "abs", "frac", "none"])
+    @pytest.mark.parametrize("norm2", ["leahy", "abs", "frac", "none"])
+    def test_normalize_crossspectrum_with_method(self, norm1, norm2):
+        cs1 = AveragedCrossspectrum.from_lightcurve(
+            self.lc1, self.lc2, segment_size=1, norm=norm1, channels_overlap=True
+        )
+        cs2 = AveragedCrossspectrum.from_lightcurve(
+            self.lc1, self.lc2, segment_size=1, norm=norm2, channels_overlap=True
+        )
+        cs3 = cs2.to_norm(norm1)
+        for attr in ["power", "power_err", "unnorm_power", "unnorm_power_err"]:
+            assert np.allclose(getattr(cs1, attr), getattr(cs3, attr))
+            assert np.allclose(getattr(cs1.pds1, attr), getattr(cs3.pds1, attr))
+            assert np.allclose(getattr(cs1.pds2, attr), getattr(cs3.pds2, attr))
+
+    @pytest.mark.parametrize("f", [None, 1.5])
+    @pytest.mark.parametrize("norm", ["leahy", "abs", "frac", "none"])
+    def test_rebin_factor_rebins_all_attrs(self, f, norm):
+        cs1 = AveragedCrossspectrum.from_lightcurve(
+            self.lc1, self.lc2, segment_size=1, norm=norm, channels_overlap=True
+        )
+        # N.B.: if f is not None, df gets ignored.
+        new_cs = cs1.rebin(df=1.5, f=f)
+        N = new_cs.freq.size
+        for attr in ["power", "power_err", "unnorm_power", "unnorm_power_err"]:
+            assert hasattr(new_cs, attr) and getattr(new_cs, attr).size == N
+            assert hasattr(new_cs.pds1, attr) and getattr(new_cs.pds1, attr).size == N
+            assert hasattr(new_cs.pds2, attr) and getattr(new_cs.pds2, attr).size == N
+
+        for attr in cs1.meta_attrs():
+            if attr not in ["df", "gti", "m"]:
+                assert getattr(cs1, attr) == getattr(new_cs, attr)
+
+    @pytest.mark.parametrize("norm", ["leahy", "abs", "frac", "none"])
+    def test_rebin_factor_log_rebins_all_attrs(self, norm):
+        cs1 = AveragedCrossspectrum.from_lightcurve(
+            self.lc1, self.lc2, segment_size=1, norm=norm, channels_overlap=True
+        )
+        new_cs = cs1.rebin_log(0.03)
+        N = new_cs.freq.size
+        for attr in ["power", "power_err", "unnorm_power", "unnorm_power_err"]:
+            assert hasattr(new_cs, attr) and getattr(new_cs, attr).size == N
+            assert hasattr(new_cs.pds1, attr) and getattr(new_cs.pds1, attr).size == N
+            assert hasattr(new_cs.pds2, attr) and getattr(new_cs.pds2, attr).size == N
+
+        for attr in cs1.meta_attrs():
+            if attr not in ["df", "gti", "m", "k"]:
+                assert np.all(getattr(cs1, attr) == getattr(new_cs, attr))
+
+    def test_rebin(self):
+        new_cs = self.cs.rebin(df=1.5)
+        assert hasattr(new_cs, "dt") and new_cs.dt is not None
+        assert new_cs.df == 1.5
+        new_cs.time_lag()
+
+    def test_rebin_factor(self):
+        new_cs = self.cs.rebin(f=1.5)
+        assert hasattr(new_cs, "dt") and new_cs.dt is not None
+        assert new_cs.df == self.cs.df * 1.5
+        new_cs.time_lag()
+
+    def test_rebin_log(self):
+        # For now, just verify that it doesn't crash
+        new_cs = self.cs.rebin_log(f=0.1)
+        assert hasattr(new_cs, "dt") and new_cs.dt is not None
+        assert isinstance(new_cs, type(self.cs))
+        new_cs.time_lag()
+
+    def test_rebin_log_returns_complex_values_and_errors(self):
+        # For now, just verify that it doesn't crash
+        new_cs = self.cs.rebin_log(f=0.1)
+        assert np.iscomplexobj(new_cs.power[0])
+        assert np.iscomplexobj(new_cs.power_err[0])