Repository: StingraySoftware/stingray

Issue:
Errors in rms_calculation are incorrect
When simulating many times a PDS with a given rms and measuring it with our `rms_calculation` in `fourier.py`, the error bars are underestimated (red lines indicate the calculated error, the blue ones the actual standard deviation of the measurements in the thousands of realization)
![image](https://github.com/StingraySoftware/stingray/assets/7190189/c27a9288-cf33-4505-a6ba-0faea0e18d1f)

I re-implemented the method from scratch, using clearer passages to normalize the spectrum correctly, and I obtain something like this: 
![image](https://github.com/StingraySoftware/stingray/assets/7190189/8a2daf5e-7ed5-4f0c-a503-2a5c3bc423c7)

So, there must be an error somewhere. Since the current code is rather clunky, I'm inclined to deprecate the old function and create a new one with clearer parameters (e.g. making it clear whether a parameter refers to the full observation or each segment of the averaged spectrum, like `nphots_per_segment`)

Discussion:

Patch:
diff --git a/docs/changes/792.bugfix.rst b/docs/changes/792.bugfix.rst
new file mode 100644
index 00000000..76409e8c
--- /dev/null
+++ b/docs/changes/792.bugfix.rst
@@ -0,0 +1,1 @@
+Fix rms computation and error bars
diff --git a/docs/notebooks b/docs/notebooks
index 1b09cfca..44828775 160000
--- a/docs/notebooks
+++ b/docs/notebooks
@@ -1,1 +1,1 @@
-Subproject commit 1b09cfcad4684ac6f55227e2e92a343b183ef435
+Subproject commit 448287755a373496cd5f81ae6ee904d3de0f4561
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index e3e124a4..75d40950 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -6,22 +6,20 @@ import numpy as np
 import scipy
 import scipy.optimize
 import scipy.stats
-from astropy import log
 import matplotlib.pyplot as plt
 
 from stingray.exceptions import StingrayError
-from stingray.gti import bin_intervals_from_gtis, check_gtis, cross_two_gtis
 from stingray.utils import rebin_data, rebin_data_log, simon
 
 from .base import StingrayObject
 from .events import EventList
 from .gti import cross_two_gtis, time_intervals_from_gtis
 from .lightcurve import Lightcurve
-from .utils import show_progress
 from .fourier import avg_cs_from_iterables, error_on_averaged_cross_spectrum
 from .fourier import avg_cs_from_events, poisson_level
-from .fourier import fftfreq, fft, normalize_periodograms, raw_coherence
+from .fourier import normalize_periodograms, raw_coherence
 from .fourier import get_flux_iterable_from_segments, power_color
+from .fourier import get_rms_from_unnorm_periodogram
 
 from scipy.special import factorial
 
@@ -853,7 +851,7 @@ class Crossspectrum(StingrayObject):
 
         for attr in ["power", "power_err"]:
             unnorm_attr = "unnorm_" + attr
-            if not hasattr(self, unnorm_attr):
+            if not hasattr(self, unnorm_attr) or getattr(self, unnorm_attr) is None:
                 continue
             power = normalize_periodograms(
                 getattr(self, unnorm_attr),
@@ -2356,36 +2354,33 @@ class DynamicalCrossspectrum(AveragedCrossspectrum):
             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
+
+        good = (self.freq >= min_freq) & (self.freq <= max_freq)
+
+        M_freq = self.m
+        if isinstance(self.m, Iterable):
+            M_freq = self.m[good]
+
         rmss = []
         rms_errs = []
+
         for unnorm_powers in dyn_ps_unnorm.T:
-            r, re = rms_calculation(
-                unnorm_powers[minind:maxind],
-                min_freq,
-                max_freq,
+            r, re = get_rms_from_unnorm_periodogram(
+                unnorm_powers[good],
                 nphots,
-                T,
-                M_freqs,
-                K_freqs,
-                freq_bins,
-                poisson_noise_unnrom,
-                deadtime=0.0,
+                self.df,
+                M=M_freq,
+                poisson_noise_unnorm=poisson_noise_unnrom,
+                segment_size=self.segment_size,
+                kind="frac",
             )
+
             rmss.append(r)
             rms_errs.append(re)
         return rmss, rms_errs
diff --git a/stingray/deadtime/fad.py b/stingray/deadtime/fad.py
index 6ea87679..01a45398 100644
--- a/stingray/deadtime/fad.py
+++ b/stingray/deadtime/fad.py
@@ -9,11 +9,11 @@ from astropy import log
 from astropy.table import Table
 
 from stingray.lightcurve import Lightcurve
-from ..crossspectrum import AveragedCrossspectrum, show_progress, get_flux_generator
+from ..crossspectrum import AveragedCrossspectrum, get_flux_generator
 from ..powerspectrum import AveragedPowerspectrum
 from ..fourier import normalize_periodograms, fft, fftfreq, positive_fft_bins
-
-from ..gti import cross_two_gtis, bin_intervals_from_gtis
+from ..utils import show_progress
+from ..gti import cross_two_gtis
 
 
 __all__ = ["calculate_FAD_correction", "get_periodograms_from_FAD_results", "FAD"]
diff --git a/stingray/fourier.py b/stingray/fourier.py
index a39d1da5..52d3e046 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -1050,6 +1050,163 @@ def estimate_intrinsic_coherence(cross_power, power1, power2, power1_noise, powe
     return new_coherence
 
 
+def get_rms_from_rms_norm_periodogram(power_sqrms, poisson_noise_sqrms, df, M, low_M_buffer_size=4):
+    r"""Calculate integrated rms spectrum (frac or abs).
+
+    If M=1, it starts by rebinning the powers slightly in order to get a slightly
+    better approximation for the error bars.
+
+    Parameters
+    ----------
+    power_sqrms: array-like
+        Powers, in units of fractional rms ($(rms/mean)^2 Hz{-1}$)
+    poisson_noise_sqrms: float
+        Poisson noise level, in units of fractional rms ($(rms/mean)^2 Hz{-1}$
+    df: float or ``np.array``, same dimension of ``power_sqrms``
+        The frequency resolution of each power
+    M: int or ``np.array``, same dimension of ``power_sqrms``
+        The number of powers averaged to obtain each value of power.
+
+    Other Parameters
+    ----------------
+    low_M_buffer_size : int, default 4
+        If M=1, the powers are rebinned to have a minimum of ``low_M_buffer_size`` powers
+        in each bin. This is done to get a better estimate of the error bars.
+    """
+    from stingray.utils import rebin_data
+
+    m_is_iterable = isinstance(M, Iterable)
+    df_is_iterable = isinstance(df, Iterable)
+
+    # if M is an iterable but all values are the same, let's simplify
+    if m_is_iterable and len(list(set(M))) == 1:
+        M = M[0]
+        m_is_iterable = False
+    # the same with df
+    if df_is_iterable and len(list(set(df))) == 1:
+        df = df[0]
+        df_is_iterable = False
+
+    low_M_values = M < 30
+
+    if np.any(M < 30):
+        quantity = "Some"
+        if not m_is_iterable or np.count_nonzero(low_M_values) == M.size:
+            quantity = "All"
+        warnings.warn(
+            f"{quantity} power spectral bins have M<30. The error bars on the rms might be wrong. "
+            "In some cases one might try to increase the number of segments, for example by "
+            "reducing the segment size, in order to obtain at least 30 segments."
+        )
+    # But they cannot be of different kind. There would be something wrong with the data
+    if m_is_iterable != df_is_iterable:
+        raise ValueError("M and df must be either both constant, or none of them.")
+
+    # If powers are not rebinned and M=1, we rebin them slightly. The error on the power is tricky,
+    # because powers follow a non-central chi squared distribution. If we combine powers with
+    # very different underlying signal level, the error bars will be completely wrong. But
+    # nearby powers have a higher chance of having similar values, and so, by combining them,
+    # we have a higher chance of obtaining sensible quasi-Gaussian error bars, easier to
+    # propagate through standard quadrature summation
+    if not m_is_iterable and M == 1 and power_sqrms.size > low_M_buffer_size:
+        _, local_power_sqrms, _, local_M = rebin_data(
+            np.arange(power_sqrms.size) * df,
+            power_sqrms,
+            df * low_M_buffer_size,
+            yerr=None,
+            method="average",
+            dx=df,
+        )
+        local_df = df * local_M
+        total_local_powers = local_M * local_power_sqrms.size
+        total_power_err = np.sqrt(np.sum(local_power_sqrms**2 * local_df**2 / local_M))
+        total_power_err *= power_sqrms.size / total_local_powers
+    else:
+        total_power_err = np.sqrt(np.sum(power_sqrms**2 * df**2 / M))
+
+    powers_sub = power_sqrms - poisson_noise_sqrms
+
+    total_power_sub = np.sum(powers_sub * df)
+
+    if total_power_sub < 0:
+        # By the definition, it makes no sense to define an error bar here.
+        warnings.warn("Poisson-subtracted power is below 0")
+        return 0.0, 0.0
+
+    rms = np.sqrt(total_power_sub)
+
+    high_snr_err = 0.5 / rms * total_power_err
+
+    return rms, high_snr_err
+
+
+def get_rms_from_unnorm_periodogram(
+    unnorm_powers,
+    nphots_per_segment,
+    df,
+    M=1,
+    poisson_noise_unnorm=None,
+    segment_size=None,
+    kind="frac",
+):
+    """Calculate the fractional rms amplitude from unnormalized powers.
+
+    We assume the powers come from an unnormalized Bartlett periodogram.
+    If so, the Poisson noise level is ``nphots_per_segment``, but the user
+    can specify otherwise (e.g. if the Poisson noise level is altered by dead time).
+    The ``segment_size`` and ``nphots_per_segment`` parameters refer to the length
+    and averaged counts of each segment of data used for the Bartlett periodogram.
+
+    Parameters
+    ----------
+    unnorm_powers : np.ndarray
+        The unnormalized power spectrum
+    nphots_per_segment : float
+        The averaged number of photons per segment of the data used for the Bartlett periodogram
+    df : float
+        The frequency resolution of the periodogram
+
+    Other parameters
+    ----------------
+    poisson_noise_unnorm : float
+        The unnormalized Poisson noise level
+    segment_size : float
+        The size of the segment, in seconds
+    M : int
+        The number of segments averaged to obtain the periodogram
+    kind : str
+        One of "frac" or "abs"
+    """
+    if segment_size is None:
+        segment_size = 1 / np.min(df)
+
+    if poisson_noise_unnorm is None:
+        poisson_noise_unnorm = nphots_per_segment
+
+    meanrate = nphots_per_segment / segment_size
+
+    def to_leahy(powers):
+        return powers * 2.0 / nphots_per_segment
+
+    def to_frac(powers):
+        return to_leahy(powers) / meanrate
+
+    def to_abs(powers):
+        return to_leahy(powers) * meanrate
+
+    if kind.startswith("frac"):
+        to_norm = to_frac
+    elif kind.startswith("abs"):
+        to_norm = to_abs
+    else:
+        raise ValueError("Only 'frac' or 'abs' rms are supported.")
+
+    poisson = to_norm(poisson_noise_unnorm)
+    powers = to_norm(unnorm_powers)
+
+    return get_rms_from_rms_norm_periodogram(powers, poisson, df, M)
+
+
 def rms_calculation(
     unnorm_powers,
     min_freq,
@@ -1059,7 +1216,7 @@ def rms_calculation(
     M_freqs,
     K_freqs,
     freq_bins,
-    poisson_noise_unnrom,
+    poisson_noise_unnorm,
     deadtime=0.0,
 ):
     """
@@ -1069,7 +1226,7 @@ def rms_calculation(
 
     Parameters
     ----------
-    unnrom_powers: array of float
+    unnorm_powers: array of float
         unnormalised power or cross spectrum, the array has already been
         filtered for the given frequency range
 
@@ -1100,7 +1257,7 @@ def rms_calculation(
         if it NOT rebinned freq_bins is the number of frequency bins
         in the given frequency range.
 
-    poisson_noise_unnrom : float
+    poisson_noise_unnorm : float
         This is the Poisson noise level unnormalised.
 
     Other parameters
@@ -1118,25 +1275,18 @@ def rms_calculation(
         The error on the fractional rms amplitude.
 
     """
-    rms_squared = (
-        np.sum((unnorm_powers - poisson_noise_unnrom) * 1 / T * K_freqs) * 2 * T / nphots**2
+    warnings.warn(
+        "The rms_calculation function is deprecated. Use get_rms_from_unnorm_periodogram instead.",
+        DeprecationWarning,
     )
+    rms_norm_powers = unnorm_powers * 2 * T / nphots**2
+    rms_poisson_noise = poisson_noise_unnorm * 2 * T / nphots**2
 
-    if rms_squared < 0.0:
-        rms_err = np.sqrt(
-            np.var((unnorm_powers - poisson_noise_unnrom) * 1 / T * K_freqs) * 2 * T / nphots**2
-        )
-        return 0.0, rms_err
-    rms = np.sqrt(rms_squared)
-
-    rms_noise_squared = (
-        poisson_noise_unnrom * (max_freq - min_freq) * 2 * T / nphots**2
-    )  # rms of the noise
-    rms_err_squared = (2 * rms_squared * rms_noise_squared + rms_noise_squared**2) / (
-        2 * np.sum(M_freqs) * freq_bins * rms_squared
-    )
-    rms_err = np.sqrt(rms_err_squared)
+    df = 1 / T * K_freqs
 
+    rms, rms_err = get_rms_from_rms_norm_periodogram(
+        rms_norm_powers, rms_poisson_noise, df, M_freqs
+    )
     return rms, rms_err
 
 
@@ -1641,7 +1791,7 @@ def avg_cs_from_iterables_quick(flux_iterable1, flux_iterable2, dt, norm="frac")
 
     """
     # Initialize stuff
-    unnorm_cross = unnorm_pds1 = unnorm_pds2 = None
+    unnorm_cross = None
     n_ave = 0
 
     sum_of_photons1 = sum_of_photons2 = 0
diff --git a/stingray/multitaper.py b/stingray/multitaper.py
index 6585ccdc..6de5bfab 100644
--- a/stingray/multitaper.py
+++ b/stingray/multitaper.py
@@ -171,6 +171,7 @@ class Multitaper(Powerspectrum):
             self.nphots = None
             self.k = 1
             self.jk_var_deg_freedom = None
+            self.segment_size = None
             return
         elif not isinstance(data, EventList):
             lc = data
@@ -182,6 +183,7 @@ class Multitaper(Powerspectrum):
         self.power_type = "real"
         self.fullspec = False
         self.k = 1
+        self.segment_size = None
 
         self._make_multitaper_periodogram(
             lc,
@@ -684,23 +686,6 @@ class Multitaper(Powerspectrum):
 
         return bin_mtp
 
-    def compute_rms(
-        self, min_freq, max_freq, poisson_noise_level=None, white_noise_offset=None, deadtime=0.0
-    ):
-        return Powerspectrum.compute_rms(
-            self,
-            min_freq,
-            max_freq,
-            poisson_noise_level=poisson_noise_level,
-            white_noise_offset=white_noise_offset,
-            deadtime=deadtime,
-        )
-
-    def classical_significances(self, threshold, trial_correction):
-        return Powerspectrum.classical_significances(
-            self, threshold=threshold, trial_correction=trial_correction
-        )
-
     def _fourier_multitaper_lomb_scargle(self, lc, NW=4, low_bias=True):
         """Compute the multitaper lomb-scargle spectral estimate.
 
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 53495b23..fe2fb2b2 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -1,4 +1,3 @@
-import copy
 import warnings
 from collections.abc import Generator, Iterable
 
@@ -17,14 +16,8 @@ from .lightcurve import Lightcurve
 from .fourier import avg_pds_from_iterable, unnormalize_periodograms
 from .fourier import avg_pds_from_events
 from .fourier import get_flux_iterable_from_segments
-from .fourier import rms_calculation, poisson_level
-
-try:
-    from tqdm import tqdm as show_progress
-except ImportError:
-
-    def show_progress(a, **kwargs):
-        return a
+from .fourier import poisson_level
+from .fourier import get_rms_from_unnorm_periodogram
 
 
 __all__ = ["Powerspectrum", "AveragedPowerspectrum", "DynamicalPowerspectrum"]
@@ -147,9 +140,7 @@ class Powerspectrum(Crossspectrum):
 
         return bin_ps
 
-    def compute_rms(
-        self, min_freq, max_freq, poisson_noise_level=None, white_noise_offset=None, deadtime=0.0
-    ):
+    def compute_rms(self, min_freq, max_freq, poisson_noise_level=None):
         """
         Compute the fractional rms amplitude in the power spectrum
         between two frequencies.
@@ -174,12 +165,6 @@ class Powerspectrum(Crossspectrum):
             this function using the same normalisation of the PDS
             and it will get subtracted from powers here.
 
-        white_noise_offset : float, default None
-            This is the white noise level, in Leahy normalization. In the ideal
-            case, this is 2. Dead time and other instrumental effects can alter
-            it. The user can fit the white noise level outside this function
-            and it will get subtracted from powers here.
-
         Returns
         -------
         rms: float
@@ -190,68 +175,37 @@ class Powerspectrum(Crossspectrum):
             The error on the fractional rms amplitude.
 
         """
-        minind = self.freq.searchsorted(min_freq)
-        maxind = self.freq.searchsorted(max_freq)
-        min_freq = self.freq[minind]
-
-        # To avoid corner case of searchsorted, where maxind goes out of the array
-        if maxind >= len(self.freq) - 1:
-            max_freq = self.freq[maxind - 1]
-        else:
-            max_freq = self.freq[maxind]
-
-        nphots = self.nphots
-        # distinguish the rebinned and non-rebinned case
-        if isinstance(self.m, Iterable):
-            M_freq = self.m[minind:maxind]
-            K_freq = self.k[minind:maxind]
-            freq_bins = 1
-        else:
-            M_freq = self.m
-            K_freq = self.k
-            freq_bins = maxind - minind
 
-        T = self.dt * self.n * 2
+        good = (self.freq >= min_freq) & (self.freq <= max_freq)
 
-        if white_noise_offset is not None:
-            powers = self.power[minind:maxind]
-            warnings.warn(
-                "the option white_noise_offset now deprecated and will be "
-                "removed in the next major release. The routine"
-                "is correct only with non-rebinned power-spectra.",
-                DeprecationWarning,
-            )
+        M_freq = self.m
+        K_freq = self.k
 
-            if self.norm.lower() == "leahy":
-                powers_leahy = powers.copy()
-            else:
-                powers_leahy = self.unnorm_power[minind:maxind].real * 2 / nphots
+        if isinstance(self.k, Iterable):
+            K_freq = self.k[good]
 
-            rms = np.sqrt(np.sum(powers_leahy - white_noise_offset) / nphots)
-            rms_err = self._rms_error(powers_leahy)
-            return rms, rms_err
+        if isinstance(self.m, Iterable):
+            M_freq = self.m[good]
 
+        if poisson_noise_level is None:
+            poisson_noise_unnorm = poisson_level("none", n_ph=self.nphots)
         else:
-            if poisson_noise_level is None:
-                poisson_noise_unnorm = poisson_level("none", n_ph=self.nphots)
-            else:
-                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,
-                max_freq,
-                self.nphots,
-                T,
-                M_freq,
-                K_freq,
-                freq_bins,
-                poisson_noise_unnorm,
-                deadtime,
+            poisson_noise_unnorm = unnormalize_periodograms(
+                poisson_noise_level, self.dt, self.n, self.nphots, norm=self.norm
             )
 
+        rms, rmse = get_rms_from_unnorm_periodogram(
+            self.unnorm_power[good],
+            self.nphots,
+            self.df * K_freq,
+            M=M_freq,
+            poisson_noise_unnorm=poisson_noise_unnorm,
+            segment_size=self.segment_size,
+            kind="frac",
+        )
+
+        return rms, rmse
+
     def _rms_error(self, powers):
         r"""
         Compute the error on the fractional rms amplitude using error


Test patch:
diff --git a/stingray/deadtime/tests/test_fad.py b/stingray/deadtime/tests/test_fad.py
index b17daa96..3fc5fd81 100644
--- a/stingray/deadtime/tests/test_fad.py
+++ b/stingray/deadtime/tests/test_fad.py
@@ -124,9 +124,6 @@ def test_fad_power_spectrum_compliant_objects(ctrate):
     cs_f = results["cs"].power
     ptot_f = results["ptot"].power
 
-    # Verify that compute_rms works
-    results["ptot"].compute_rms(0.01, 1)
-
     n = length / segment_size
     ncounts_per_intv1 = ncounts1 * segment_size / length
     ncounts_per_intv2 = ncounts2 * segment_size / length
@@ -277,9 +274,6 @@ def test_fad_power_spectrum_compliant_leahy_objects(ctrate):
     cs_f = results_out["cs"].power
     ptot_f = results_out["ptot"].power
 
-    # Verify that compute_rms works
-    results_out["ptot"].compute_rms(0.01, 1)
-
     n = length / segment_size
 
     pds_std_theor = 2 / np.sqrt(n)
diff --git a/stingray/tests/test_covariancespectrum.py b/stingray/tests/test_covariancespectrum.py
index a4d47149..fb8b6fb8 100644
--- a/stingray/tests/test_covariancespectrum.py
+++ b/stingray/tests/test_covariancespectrum.py
@@ -96,8 +96,9 @@ class TestCovariancespectrumwithEvents(object):
 
     def test_with_unsorted_event_list(self):
         event_list = np.array([[2, 1], [2, 3], [1, 2], [5, 2], [4, 1]])
-        with pytest.warns(UserWarning, match="must be sorted"):
+        with pytest.warns(UserWarning) as record:
             _ = Covariancespectrum(event_list, dt=1)
+        assert np.any(["must be sorted" in str(r.message) for r in record])
 
     def test_with_std_as_iterable(self):
         c = Covariancespectrum(self.event_list, dt=1, std=[1, 2])
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index 5647a3d9..ea5022fa 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -1,4 +1,5 @@
 import os
+import importlib
 import numpy as np
 import pytest
 import warnings
@@ -7,34 +8,20 @@ import scipy.special
 from astropy.io import fits
 from stingray import Lightcurve
 from stingray import Crossspectrum, AveragedCrossspectrum, DynamicalCrossspectrum
+from stingray import AveragedPowerspectrum
 from stingray.crossspectrum import cospectra_pvalue
 from stingray.crossspectrum import normalize_crossspectrum, normalize_crossspectrum_gauss
 from stingray.crossspectrum import coherence, time_lag
 from stingray import StingrayError
 from stingray.simulator import Simulator
-from stingray.fourier import poisson_level, raw_coherence
+from stingray.fourier import poisson_level
 
 from stingray.events import EventList
 import copy
 
-_HAS_XARRAY = _HAS_PANDAS = _HAS_H5PY = True
-
-try:
-    import xarray
-    from xarray import Dataset
-except ImportError:
-    _HAS_XARRAY = False
-
-try:
-    import pandas
-    from pandas import DataFrame
-except ImportError:
-    _HAS_PANDAS = False
-
-try:
-    import h5py
-except ImportError:
-    _HAS_H5PY = False
+_HAS_XARRAY = importlib.util.find_spec("xarray") is not None
+_HAS_PANDAS = importlib.util.find_spec("pandas") is not None
+_HAS_H5PY = importlib.util.find_spec("h5py") is not None
 
 np.random.seed(20160528)
 curdir = os.path.abspath(os.path.dirname(__file__))
@@ -450,7 +437,7 @@ class TestAveragedCrossspectrumEvents(object):
     def test_rebin_log(self):
         # For now, just verify that it doesn't crash
         new_cs = self.acs.rebin_log(f=0.1)
-        assert type(new_cs) == type(self.acs)
+        assert isinstance(new_cs, type(self.acs))
         new_cs.time_lag()
 
     def test_rebin_log_returns_complex_values(self):
@@ -758,18 +745,18 @@ class TestCrossspectrum(object):
     def test_rebin_log(self):
         # For now, just verify that it doesn't crash
         new_cs = self.cs.rebin_log(f=0.1)
-        assert type(new_cs) == type(self.cs)
+        assert isinstance(new_cs, type(self.cs))
         new_cs.time_lag()
 
-    def test_norm_abs(self):
+    def test_norm_abs_same_lc(self):
         # Testing for a power spectrum of lc1
         cs = Crossspectrum(self.lc1, self.lc1, norm="abs")
         assert len(cs.power) == 4999
         assert cs.norm == "abs"
         abs_noise = 2.0 * self.rate1  # expected Poisson noise level
-        assert np.isclose(np.mean(cs.power[1:]), abs_noise)
+        assert np.isclose(np.mean(cs.power[1:]), abs_noise, rtol=0.2)
 
-    def test_norm_leahy(self):
+    def test_norm_leahy_same_lc(self):
         # with pytest.warns(UserWarning) as record:
         cs = Crossspectrum(self.lc1, self.lc1, norm="leahy")
         assert len(cs.power) == 4999
@@ -777,7 +764,7 @@ class TestCrossspectrum(object):
         leahy_noise = 2.0  # expected Poisson noise level
         assert np.isclose(np.mean(cs.power[1:]), leahy_noise, rtol=0.02)
 
-    def test_norm_frac(self):
+    def test_norm_frac_same_lc(self):
         with pytest.warns(UserWarning) as record:
             cs = Crossspectrum(self.lc1, self.lc1, norm="frac")
         assert len(cs.power) == 4999
@@ -963,8 +950,8 @@ class TestCrossspectrum(object):
 
     def test_fullspec(self):
         csT = Crossspectrum(self.lc1, self.lc2, fullspec=True)
-        assert csT.fullspec == True
-        assert self.cs.fullspec == False
+        assert csT.fullspec is True
+        assert self.cs.fullspec is False
         assert csT.n == self.cs.n
         assert csT.n == len(csT.power)
         assert self.cs.n != len(self.cs.power)
@@ -1195,7 +1182,7 @@ class TestAveragedCrossspectrum(object):
         # 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 type(new_cs) == type(self.cs)
+        assert isinstance(new_cs, type(self.cs))
         new_cs.time_lag()
 
     def test_rebin_log_returns_complex_values_and_errors(self):
@@ -1413,10 +1400,12 @@ class TestDynamicalCrossspectrum(object):
 
     def test_rms_is_correct(self):
         lc = copy.deepcopy(self.lc)
-        lc.counts = np.random.poisson(lc.counts)
+        # Create a clear variable signal with an exponential decay
+        lc.counts = np.random.poisson(100000 * np.exp(-(lc.time) / 100))
+
         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
+        with pytest.warns(UserWarning, match="All power spectral bins have M<30"):
+            rms, rmse = dps.compute_rms(1 / 5, 16.0, poisson_noise_level=2)
 
         ps = AveragedPowerspectrum()
         ps.freq = dps.freq
@@ -1429,7 +1418,8 @@ class TestDynamicalCrossspectrum(object):
         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)
+        with pytest.warns(UserWarning, match="All power spectral bins have M<30"):
+            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)
 
diff --git a/stingray/tests/test_fourier.py b/stingray/tests/test_fourier.py
index 91b8df6d..51b9b560 100644
--- a/stingray/tests/test_fourier.py
+++ b/stingray/tests/test_fourier.py
@@ -1,15 +1,33 @@
 import os
-from pickle import FALSE
 
+from collections.abc import Iterable
 import pytest
+import numpy as np
+from astropy.table import Table
+
+from stingray.fourier import fft, fftfreq, normalize_abs, normalize_frac, poisson_level
+from stingray.fourier import (
+    get_flux_iterable_from_segments,
+    avg_pds_from_events,
+    avg_cs_from_events,
+)
+from stingray.fourier import normalize_periodograms, raw_coherence, estimate_intrinsic_coherence
+from stingray.fourier import bias_term, error_on_averaged_cross_spectrum, unnormalize_periodograms
+from stingray.fourier import impose_symmetry_lsft, lsft_slow, lsft_fast, rms_calculation
+from stingray.fourier import get_average_ctrate, normalize_leahy_from_variance
+from stingray.fourier import integrate_power_in_frequency_range, power_color, hue_from_power_color
+from stingray.fourier import get_rms_from_rms_norm_periodogram, get_rms_from_unnorm_periodogram
 
-from stingray.fourier import *
 from stingray.utils import check_allclose_and_print
+from astropy.modeling.models import Lorentz1D
 
 curdir = os.path.abspath(os.path.dirname(__file__))
 datadir = os.path.join(curdir, "data")
 
 
+rng = np.random.RandomState(137259723)
+
+
 def compare_tables(table1, table2, rtol=0.001, discard=[]):
     for key in table1.meta.keys():
         if key in discard:
@@ -36,7 +54,7 @@ def test_norm():
     N = 1000000
     dt = 0.2
     meanrate = mean / dt
-    lc = np.random.poisson(mean, N)
+    lc = rng.poisson(mean, N)
     pds = np.abs(fft(lc)) ** 2
     freq = fftfreq(N, dt)
     good = slice(1, N // 2)
@@ -99,8 +117,8 @@ class TestCoherence(object):
             Table.read(os.path.join(datadir, "sample_variable_series.fits"))["data"][:10000] * 1000
         )
         print(data.max(), data.min())
-        cls.data1 = np.random.poisson(data)
-        cls.data2 = np.random.poisson(data)
+        cls.data1 = rng.poisson(data)
+        cls.data2 = rng.poisson(data)
         ft1 = np.fft.fft(cls.data1)
         ft2 = np.fft.fft(cls.data2)
         dt = 0.01
@@ -138,7 +156,7 @@ class TestCoherence(object):
         C, P1, P2 = self.cross[:nbins], self.pds1[:nbins], self.pds2[:nbins]
         bsq = bias_term(P1, P2, self.p1noise, self.p2noise, self.N)
         # must be lower than bsq!
-        low_coh_cross = np.random.normal(bsq**0.5 / 10, bsq**0.5 / 100) + 0.0j
+        low_coh_cross = rng.normal(bsq**0.5 / 10, bsq**0.5 / 100) + 0.0j
         coh = raw_coherence(low_coh_cross, P1, P2, self.p1noise, self.p2noise, self.N)
         assert np.allclose(coh, 0)
         # Do it with a single number
@@ -176,13 +194,13 @@ class TestFourier(object):
         cls.ctrate = 10000
         cls.N = np.rint(cls.length / cls.dt).astype(int)
         cls.dt = cls.length / cls.N
-        cls.times = np.sort(np.random.uniform(0, cls.length, int(cls.length * cls.ctrate)))
+        cls.times = np.sort(rng.uniform(0, cls.length, int(cls.length * cls.ctrate)))
         cls.gti = np.asarray([[0, cls.length]])
         cls.counts, bins = np.histogram(cls.times, bins=np.linspace(0, cls.length, cls.N + 1))
         cls.errs = np.ones_like(cls.counts) * np.sqrt(cls.ctrate)
         cls.bin_times = (bins[:-1] + bins[1:]) / 2
         cls.segment_size = 5.0
-        cls.times2 = np.sort(np.random.uniform(0, cls.length, int(cls.length * cls.ctrate)))
+        cls.times2 = np.sort(rng.uniform(0, cls.length, int(cls.length * cls.ctrate)))
         cls.counts2, _ = np.histogram(cls.times2, bins=np.linspace(0, cls.length, cls.N + 1))
         cls.errs2 = np.ones_like(cls.counts2) * np.sqrt(cls.ctrate)
 
@@ -225,15 +243,15 @@ class TestFourier(object):
             assert np.allclose(fe, fc)
 
     def test_avg_pds_bad_input(self):
-        times = np.sort(np.random.uniform(0, 1000, 1))
+        times = np.sort(rng.uniform(0, 1000, 1))
         out_ev = avg_pds_from_events(times, self.gti, self.segment_size, self.dt)
         assert out_ev is None
 
     @pytest.mark.parametrize("return_subcs", [True, False])
     @pytest.mark.parametrize("return_auxil", [True, False])
     def test_avg_cs_bad_input(self, return_auxil, return_subcs):
-        times1 = np.sort(np.random.uniform(0, 1000, 1))
-        times2 = np.sort(np.random.uniform(0, 1000, 1))
+        times1 = np.sort(rng.uniform(0, 1000, 1))
+        times2 = np.sort(rng.uniform(0, 1000, 1))
         out_ev = avg_cs_from_events(
             times1,
             times2,
@@ -443,7 +461,7 @@ class TestNorms(object):
         good = freq > 0
         cls.good = good
         cls.meanrate = cls.mean / cls.dt
-        cls.lc = np.random.poisson(cls.mean, cls.N).astype(float)
+        cls.lc = rng.poisson(cls.mean, cls.N).astype(float)
         cls.nph = np.sum(cls.lc)
         cls.pds = (np.abs(np.fft.fft(cls.lc)) ** 2)[good]
         cls.cross = ((np.fft.fft(cls.lc)) ** 2)[good]
@@ -624,7 +642,7 @@ def test_lags(phlag):
     def func(time, phase=0):
         return 2 + np.sin(2 * np.pi * (time * freq - phase))
 
-    time = np.sort(np.random.uniform(0, 100, 3000))
+    time = np.sort(rng.uniform(0, 100, 3000))
     ft0 = lsft_slow(func(time, 0), time, np.array([freq]))
     ft1 = lsft_slow(func(time, phlag), time, np.array([freq]))
     measured_lag = (np.angle(ft1) - np.angle(ft0)) / 2 / np.pi
@@ -787,7 +805,7 @@ class TestPowerColor(object):
         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)
+            factor = rng.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)
@@ -795,3 +813,175 @@ class TestPowerColor(object):
             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)
+
+
+class TestRMS(object):
+    @classmethod
+    def setup_class(cls):
+        fwhm = 0.23456
+        cls.segment_size = 256
+        cls.df = 1 / cls.segment_size
+
+        cls.freqs = np.arange(cls.df, 1.54232, cls.df)
+        pds_shape_func = Lorentz1D(x_0=0, fwhm=fwhm)
+        cls.pds_shape_raw = pds_shape_func(cls.freqs)
+
+        pds_shape_func_qpo = Lorentz1D(x_0=0, fwhm=0.312567) + Lorentz1D(x_0=0.5, fwhm=0.1)
+        cls.pds_shape_qpo_raw = pds_shape_func_qpo(cls.freqs)
+
+    def _prepare_pds_for_rms_tests(self, rms, nphots, M, distort_poisson_by=1, with_qpo=False):
+        meanrate = nphots / self.segment_size
+        poisson_noise_rms = 2 / meanrate
+        pds_shape = self.pds_shape_raw if not with_qpo else self.pds_shape_qpo_raw
+
+        pds_shape_rms = pds_shape / np.sum(pds_shape * self.df) * rms**2
+        pds_shape_rms += poisson_noise_rms * distort_poisson_by
+
+        random_part = rng.chisquare(2 * M, size=pds_shape.size) / 2 / M
+        pds_rms_noisy = random_part * pds_shape_rms
+
+        pds_unnorm = pds_rms_noisy * meanrate / 2 * nphots
+        return pds_rms_noisy, pds_unnorm
+
+    @pytest.mark.parametrize("M", [100, 10000])
+    @pytest.mark.parametrize("nphots", [100_000, 1_000_000])
+    @pytest.mark.parametrize("rms", [0.05, 0.1, 0.32, 0.5])
+    @pytest.mark.parametrize("with_qpo", [False, True])
+    def test_rms(self, M, nphots, rms, with_qpo):
+        meanrate = nphots / self.segment_size
+        poisson_noise_rms = 2 / meanrate
+        pds_rms_noisy, pds_unnorm = self._prepare_pds_for_rms_tests(
+            rms, nphots, M, with_qpo=with_qpo
+        )
+
+        rms_from_unnorm, rmse_from_unnorm = get_rms_from_unnorm_periodogram(
+            pds_unnorm,
+            nphots,
+            self.df,
+            M=M,
+        )
+        rms_from_rms, rmse_from_rms = get_rms_from_rms_norm_periodogram(
+            pds_rms_noisy, poisson_noise_rms, self.df, M
+        )
+
+        assert np.isclose(rms_from_rms, rms, atol=3 * rmse_from_rms)
+        assert np.isclose(rms_from_unnorm, rms, atol=3 * rmse_from_unnorm)
+
+    @pytest.mark.parametrize("M", [100, 10000])
+    @pytest.mark.parametrize("nphots", [100_000, 1_000_000])
+    @pytest.mark.parametrize("rms", [0.05, 0.1, 0.32, 0.5])
+    def test_rms_abs(self, M, nphots, rms):
+        meanrate = nphots / self.segment_size
+        _, pds_unnorm = self._prepare_pds_for_rms_tests(rms, nphots, M)
+
+        rms_from_unnorm, rmse_from_unnorm = get_rms_from_unnorm_periodogram(
+            pds_unnorm, nphots, self.df, M=M, kind="abs"
+        )
+        assert np.isclose(rms_from_unnorm, rms * meanrate, atol=3 * rmse_from_unnorm * meanrate)
+
+    @pytest.mark.parametrize("M", [1, 10])
+    @pytest.mark.parametrize("nphots", [100_000, 1_000_000])
+    @pytest.mark.parametrize("rms", [0.05, 0.1, 0.32, 0.5])
+    def test_rms_M_low(self, M, nphots, rms):
+        """Test that the warning is raised when M is low."""
+        meanrate = nphots / self.segment_size
+        poisson_noise_rms = 2 / meanrate
+
+        pds_rms_noisy, pds_unnorm = self._prepare_pds_for_rms_tests(rms, nphots, M)
+
+        with pytest.warns(UserWarning, match="All power spectral bins have M<30."):
+            rms_from_unnorm, rmse_from_unnorm = get_rms_from_unnorm_periodogram(
+                pds_unnorm,
+                nphots,
+                self.df,
+                M=M,
+            )
+        with pytest.warns(UserWarning, match="All power spectral bins have M<30."):
+            rms_from_rms, rmse_from_rms = get_rms_from_rms_norm_periodogram(
+                pds_rms_noisy, poisson_noise_rms, self.df, M
+            )
+
+        assert np.isclose(rms_from_rms, rms, atol=3 * rmse_from_rms)
+        assert np.isclose(rms_from_unnorm, rms, atol=3 * rmse_from_unnorm)
+
+    @pytest.mark.parametrize("nphots", [100_000, 1_000_000])
+    def test_rms_low(self, nphots):
+        meanrate = nphots / self.segment_size
+        poisson_noise_rms = 2 / meanrate
+        M = 100
+
+        pds_rms_noisy, pds_unnorm = self._prepare_pds_for_rms_tests(
+            0, nphots, M, distort_poisson_by=0.9
+        )
+
+        with pytest.warns(UserWarning, match="Poisson-subtracted power is below 0"):
+            get_rms_from_unnorm_periodogram(
+                pds_unnorm,
+                nphots,
+                self.df,
+                M=M,
+                kind="frac",
+            )
+        with pytest.warns(UserWarning, match="Poisson-subtracted power is below 0"):
+            get_rms_from_rms_norm_periodogram(pds_rms_noisy, poisson_noise_rms, self.df, M)
+
+    def test_array_m_and_df(self):
+        # Very safe, high-rms dataset
+        nphots = 1_000_000
+        rms = 0.5
+        M = 1000
+
+        meanrate = nphots / self.segment_size
+        poisson_noise_rms = 2 / meanrate
+
+        pds_rms_noisy, _ = self._prepare_pds_for_rms_tests(rms, nphots, M)
+
+        M = np.zeros_like(pds_rms_noisy) + 100
+        df = np.zeros_like(pds_rms_noisy) + self.df
+
+        rms_from_rms, rmse_from_rms = get_rms_from_rms_norm_periodogram(
+            pds_rms_noisy, poisson_noise_rms, df, M
+        )
+
+        assert np.isclose(rms_from_rms, rms, atol=3 * rmse_from_rms)
+
+    def test_incompatible_m_and_df(self):
+        # Make df non constant
+        df = np.zeros_like(self.pds_shape_raw) + self.df
+        df[-1] = 2 * self.df
+
+        with pytest.raises(
+            ValueError, match="M and df must be either both constant, or none of them."
+        ):
+            get_rms_from_rms_norm_periodogram(self.pds_shape_raw, 2, df, M=100)
+
+    def test_invalid_kind(self):
+        # Make df non constant
+
+        with pytest.raises(ValueError, match="Only 'frac' or 'abs' rms are supported."):
+            get_rms_from_unnorm_periodogram(self.pds_shape_raw, 2, 0.1, M=100, kind="asdfkhf")
+
+    def test_deprecation_rms_calculation(self):
+        nphots = 1_000_000
+        rms = 0.5
+        M = 1000
+        _, pds_unnorm = self._prepare_pds_for_rms_tests(rms, nphots, M)
+        with pytest.warns(DeprecationWarning, match="The rms_calculation function is deprecated"):
+            rms, _ = rms_calculation(
+                pds_unnorm,
+                self.freqs.min(),
+                self.freqs.max(),
+                nphots,
+                self.segment_size,
+                M,
+                1,
+                len(self.freqs),
+                poisson_noise_unnorm=nphots,
+            )
+        rms_from_unnorm, rmse_from_unnorm = get_rms_from_unnorm_periodogram(
+            pds_unnorm,
+            nphots,
+            self.df,
+            M=M,
+        )
+        assert np.isclose(rms, rms_from_unnorm, atol=3 * rmse_from_unnorm)
diff --git a/stingray/tests/test_io.py b/stingray/tests/test_io.py
index 86ef6c6e..8c30e558 100644
--- a/stingray/tests/test_io.py
+++ b/stingray/tests/test_io.py
@@ -75,7 +75,8 @@ class TestIO(object):
             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"):
-            evdata_unsrt = load_events_and_gtis(fname_unsrt)
+            with pytest.warns(AstropyUserWarning, match="No valid GTI extensions"):
+                evdata_unsrt = load_events_and_gtis(fname_unsrt)
 
         for attr in "ev_list", "energy_list", "pi_list":
             assert np.allclose(getattr(evdata, attr), getattr(evdata_unsrt, attr))
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index 78d82e1f..b247b4bc 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -758,8 +758,9 @@ class TestLightcurve(object):
         lc1 = Lightcurve(self.times, self.counts)
         lc2 = Lightcurve(_times, _counts)
 
-        with pytest.warns(UserWarning, match="different bin widths") as w:
+        with pytest.warns(UserWarning) as record:
             lc1.join(lc2)
+        assert np.any(["different bin widths" in str(r.message) for r in record])
 
     def test_join_with_different_mjdref(self):
         shift = 86400.0  # day
@@ -809,8 +810,9 @@ class TestLightcurve(object):
         gti1 = [[4.5, 9.5]]
         lc0 = Lightcurve(time0, counts=count0, err=np.asarray(count0) / 2, dt=1, gti=gti0)
         lc1 = Lightcurve(time1, counts=count1, dt=1, gti=gti1)
-        with pytest.warns(UserWarning, match="The _counts_err array is empty in one of the"):
+        with pytest.warns(UserWarning) as record:
             lc = lc0.concatenate(lc1)
+        assert np.any(["The _counts_err array" in str(r.message) for r in record])
         assert np.allclose(lc.counts, count0 + count1)
         # Errors have been defined inside
         assert len(lc.counts_err) == len(lc.counts)
@@ -867,8 +869,9 @@ class TestLightcurve(object):
             warnings.simplefilter("ignore", category=UserWarning)
             lc2 = Lightcurve(_times, _counts, err_dist="gauss")
 
-        with pytest.warns(UserWarning, match="We are setting the errors to zero."):
+        with pytest.warns(UserWarning) as record:
             lc3 = lc1.join(lc2)
+            assert np.any(["We are setting the errors to zero." in str(r.message) for r in record])
         assert np.allclose(lc3.counts_err, np.zeros_like(lc3.time))
 
     def test_truncate_by_index(self):
@@ -1640,10 +1643,7 @@ class TestBexvar(object):
     @pytest.mark.skipif("not _HAS_ULTRANEST")
     def test_bexvar_with_dt_as_array(self):
         # create lightcurve with ``dt`` as an array
-        with pytest.warns(
-            UserWarning,
-            match="Some functionalities of Stingray Lightcurve will not work when `dt` is Iterable",
-        ):
+        with pytest.warns(UserWarning) as record:
             lc = Lightcurve(
                 time=self.time,
                 counts=self.src_counts,
@@ -1652,8 +1652,15 @@ class TestBexvar(object):
                 bg_counts=self.bg_counts,
                 bg_ratio=self.bg_ratio,
                 frac_exp=self.frac_exp,
+                skip_checks=True,
             )
-
+        assert np.any(
+            [
+                "Some functionalities of Stingray Lightcurve will not work when `dt` is Iterable"
+                in str(r.message)
+                for r in record
+            ]
+        )
         # provide time intervals externally to find bexvar
         log_cr_sigma_from_method = lc.bexvar()
         log_cr_sigma_result = np.load(self.fname_result, allow_pickle=True)[1]
@@ -1690,11 +1697,8 @@ class TestArraydt(object):
         bg_ratio = np.array([1, 1, 0.5, 1])
         frac_exp = np.array([1, 1, 1, 1])
         gti = np.array([[0.5, 4.5]])
-        with pytest.warns(
-            UserWarning,
-            match="Some functionalities of Stingray Lightcurve will not work when `dt` is Iterable",
-        ):
-            lc = Lightcurve(
+        with pytest.warns(UserWarning) as record:
+            Lightcurve(
                 time=times,
                 counts=counts,
                 dt=dt,
@@ -1704,15 +1708,19 @@ class TestArraydt(object):
                 bg_ratio=bg_ratio,
                 frac_exp=frac_exp,
             )
+            assert np.any(
+                [
+                    "Some functionalities of Stingray Lightcurve will not work when `dt` is Iterable"
+                    in str(r.message)
+                    for r in record
+                ]
+            )
 
         # demonstrate that we can create a Lightcurve object with dt being an array of floats
         # and without explicitly providing gtis.
 
-        with pytest.warns(
-            UserWarning,
-            match="Some functionalities of Stingray Lightcurve will not work when `dt` is Iterable",
-        ):
-            lc = Lightcurve(
+        with pytest.warns(UserWarning) as record:
+            Lightcurve(
                 time=times,
                 counts=counts,
                 dt=dt,
@@ -1721,13 +1729,24 @@ class TestArraydt(object):
                 bg_ratio=bg_ratio,
                 frac_exp=frac_exp,
             )
+            assert np.any(
+                [
+                    "Some functionalities of Stingray Lightcurve will not work when `dt` is Iterable"
+                    in str(r.message)
+                    for r in record
+                ]
+            )
 
     def test_warning_when_dt_is_array(self):
-        with pytest.warns(
-            UserWarning,
-            match="Some functionalities of Stingray Lightcurve will not work when `dt` is Iterable",
-        ):
+        with pytest.warns(UserWarning) as record:
             _ = Lightcurve(time=self.times, counts=self.counts, dt=self.dt)
+        assert np.any(
+            [
+                "Some functionalities of Stingray Lightcurve will not work when `dt` is Iterable"
+                in str(r.message)
+                for r in record
+            ]
+        )
 
     def test_truncate_by_index_when_dt_is_array(self):
         """
diff --git a/stingray/tests/test_multitaper.py b/stingray/tests/test_multitaper.py
index 073b9ff0..0f4d0c81 100644
--- a/stingray/tests/test_multitaper.py
+++ b/stingray/tests/test_multitaper.py
@@ -149,14 +149,16 @@ class TestMultitaper(object):
 
         lc = Lightcurve(time, counts=poisson_counts, dt=1, gti=[[0, 100]])
         mtp = Multitaper(lc, norm="leahy")
-        rms_mtp_l, rms_err_l = mtp.compute_rms(
-            min_freq=mtp.freq[1], max_freq=mtp.freq[-1], poisson_noise_level=0
-        )
+        with pytest.warns(UserWarning, match="M<30"):
+            rms_mtp_l, rms_err_l = mtp.compute_rms(
+                min_freq=mtp.freq[1], max_freq=mtp.freq[-1], poisson_noise_level=0
+            )
 
         mtp = Multitaper(lc, norm="frac")
-        rms_mtp, rms_err = mtp.compute_rms(
-            min_freq=mtp.freq[1], max_freq=mtp.freq[-1], poisson_noise_level=0
-        )
+        with pytest.warns(UserWarning, match="M<30"):
+            rms_mtp, rms_err = mtp.compute_rms(
+                min_freq=mtp.freq[1], max_freq=mtp.freq[-1], poisson_noise_level=0
+            )
         assert np.allclose(rms_mtp, rms_mtp_l, atol=0.01)
         assert np.allclose(rms_err, rms_err_l, atol=0.01)
 
@@ -236,14 +238,14 @@ class TestMultitaper(object):
 
     def test_get_adaptive_psd_with_less_tapers(self):
         with pytest.warns(UserWarning) as record:
-            mtp = Multitaper(lc=self.lc, NW=1.5, adaptive=True)
+            mtp = Multitaper(data=self.lc, NW=1.5, adaptive=True)
         assert np.any(["Not adaptively" in r.message.args[0] for r in record])
         assert mtp.multitaper_norm_power is not None
 
     @pytest.mark.parametrize("lombscargle", [False, True])
     def test_max_eigval_less_than_threshold(self, lombscargle):
         with pytest.warns(UserWarning) as record:
-            mtp = Multitaper(lc=self.lc, NW=0.5, low_bias=True, lombscargle=lombscargle)
+            mtp = Multitaper(data=self.lc, NW=0.5, low_bias=True, lombscargle=lombscargle)
         assert np.any(["not properly use low_bias" in r.message.args[0] for r in record])
         assert len(mtp.eigvals) > 0
 
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index a43d4c20..cd1f6c04 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -2,33 +2,19 @@ import os
 import numpy as np
 import copy
 import warnings
+import importlib
 
 import pytest
 import matplotlib.pyplot as plt
 from astropy.io import fits
-from astropy.table import Table
 from stingray import Lightcurve
 from stingray.events import EventList
 from stingray import Powerspectrum, AveragedPowerspectrum, DynamicalPowerspectrum
+from astropy.modeling.models import Lorentz1D
 
-_HAS_XARRAY = _HAS_PANDAS = _HAS_H5PY = True
-
-try:
-    import xarray
-    from xarray import Dataset
-except ImportError:
-    _HAS_XARRAY = False
-
-try:
-    import pandas
-    from pandas import DataFrame
-except ImportError:
-    _HAS_PANDAS = False
-
-try:
-    import h5py
-except ImportError:
-    _HAS_H5PY = False
+_HAS_XARRAY = importlib.util.find_spec("xarray") is not None
+_HAS_PANDAS = importlib.util.find_spec("pandas") is not None
+_HAS_H5PY = importlib.util.find_spec("h5py") is not None
 
 
 def clear_all_figs():
@@ -37,7 +23,8 @@ def clear_all_figs():
         plt.close(fig)
 
 
-np.random.seed(20150907)
+rng = np.random.RandomState(20150907)
+
 curdir = os.path.abspath(os.path.dirname(__file__))
 datadir = os.path.join(curdir, "data")
 
@@ -50,7 +37,7 @@ class TestAveragedPowerspectrumEvents(object):
         cls.dt = 0.0001
         cls.segment_size = tend - tstart
 
-        times = np.sort(np.random.uniform(tstart, tend, 1000))
+        times = np.sort(rng.uniform(tstart, tend, 1000))
         gti = np.array([[tstart, tend]])
 
         cls.events = EventList(times, gti=gti)
@@ -343,7 +330,7 @@ class TestAveragedPowerspectrumEvents(object):
         lc_all = []
         for i in range(n):
             time = np.arange(0.0, 10.0, 10.0 / 10000)
-            counts = np.random.poisson(1000, size=time.shape[0])
+            counts = rng.poisson(1000, size=time.shape[0])
             lc = Lightcurve(time, counts)
             lc_all.append(lc)
 
@@ -366,7 +353,7 @@ class TestPowerspectrum(object):
         mean_count_rate = 100.0
         mean_counts = mean_count_rate * dt
 
-        poisson_counts = np.random.poisson(mean_counts, size=time.shape[0])
+        poisson_counts = rng.poisson(mean_counts, size=time.shape[0])
 
         cls.lc = Lightcurve(time, counts=poisson_counts, dt=dt, gti=[[tstart, tend]])
 
@@ -463,156 +450,13 @@ class TestPowerspectrum(object):
         std_lc = np.var(self.lc.counts) / np.mean(self.lc.counts) ** 2
         assert np.isclose(ps_int, std_lc, atol=0.01, rtol=0.01)
 
-    def test_compute_rms_wrong_norm(self):
-        ps = Powerspectrum(self.lc)
-        ps.norm = "gibberish"
-        # This will now pass, due to changes on 2022-10-10
-        ps.compute_rms(0, 10)
-
-    def test_compute_rms_rebinning_is_consistent(self):
-        time = np.arange(0, 100, 1) + 0.5
-
-        poisson_counts = np.random.poisson(100.0, size=time.shape[0])
-
-        lc = Lightcurve(time, counts=poisson_counts, dt=1, gti=[[0, 100]])
-        ps = Powerspectrum(lc, norm="leahy")
-        ps_rebinned = ps.rebin_log()
-        rms, err_rms = ps.compute_rms(
-            min_freq=ps.freq[1], max_freq=ps.freq[-2], poisson_noise_level=0
-        )
-        rms_reb, err_rms_reb = ps_rebinned.compute_rms(
-            min_freq=ps.freq[1], max_freq=ps.freq[-2], poisson_noise_level=0
-        )
-        assert np.isclose(rms, rms_reb, atol=0.01, rtol=0.01)
-
-    def test_fractional_rms_in_frac_norm_is_consistent(self):
-        time = np.arange(0, 100, 1) + 0.5
-
-        poisson_counts = np.random.poisson(100.0, size=time.shape[0])
-
-        lc = Lightcurve(time, counts=poisson_counts, dt=1, gti=[[0, 100]])
-        ps = Powerspectrum(lc, norm="leahy")
-        rms_ps_l, rms_err_l = ps.compute_rms(
-            min_freq=ps.freq[1], max_freq=ps.freq[-1], poisson_noise_level=0
-        )
-
-        ps = Powerspectrum(lc, norm="frac")
-        rms_ps, rms_err = ps.compute_rms(
-            min_freq=ps.freq[1], max_freq=ps.freq[-1], poisson_noise_level=0
-        )
-        assert np.allclose(rms_ps, rms_ps_l, atol=0.01)
-        assert np.allclose(rms_err, rms_err_l, atol=0.01)
-
-    def test_fractional_rms_in_frac_norm_is_consistent_old(self):
-        with pytest.warns(DeprecationWarning, match="the option white_noise_offset"):
-            time = np.arange(0, 100, 1) + 0.5
-
-            poisson_counts = np.random.poisson(100.0, size=time.shape[0])
-
-            lc = Lightcurve(time, counts=poisson_counts, dt=1, gti=[[0, 100]])
-            ps = Powerspectrum(lc, norm="leahy")
-            rms_ps_l, rms_err_l = ps.compute_rms(
-                min_freq=ps.freq[1], max_freq=ps.freq[-1], white_noise_offset=0
-            )
-
-            ps = Powerspectrum(lc, norm="frac")
-            rms_ps, rms_err = ps.compute_rms(
-                min_freq=ps.freq[1], max_freq=ps.freq[-1], white_noise_offset=0
-            )
-            assert np.allclose(rms_ps, rms_ps_l, atol=0.01)
-            assert np.allclose(rms_err, rms_err_l, atol=0.01)
-
-    @pytest.mark.parametrize("norm", ["frac", "abs", "none"])
-    def test_fractional_rms_in_frac_norm_is_consistent_averaged_noPnoise(self, norm):
-        time = np.arange(0, 400, 1) + 0.5
-
-        poisson_counts = np.random.poisson(100.0, size=time.shape[0])
-
-        lc = Lightcurve(time, counts=poisson_counts, dt=1, gti=[[0, 400]])
-        ps = AveragedPowerspectrum(lc, norm="leahy", segment_size=100, silent=True)
-        rms_ps_l, rms_err_l = ps.compute_rms(
-            min_freq=ps.freq[1], max_freq=ps.freq[-1], poisson_noise_level=0
-        )
-
-        ps = AveragedPowerspectrum(lc, norm=norm, segment_size=100)
-        rms_ps, rms_err = ps.compute_rms(
-            min_freq=ps.freq[1], max_freq=ps.freq[-1], poisson_noise_level=0
-        )
-        assert np.allclose(rms_ps, rms_ps_l, atol=0.01)
-        assert np.allclose(rms_err, rms_err_l, atol=0.01)
-
-    @pytest.mark.parametrize("norm", ["frac", "abs", "none"])
-    def test_fractional_rms_in_frac_norm_is_consistent_averaged(self, norm):
-        time = np.arange(0, 400, 1) + 0.5
-
-        data = Table.read(os.path.join(datadir, "sample_variable_series.fits"))["data"][:400] * 1000
-        poisson_counts = np.random.poisson(data)
-
-        lc = Lightcurve(time, counts=poisson_counts, dt=1, gti=[[0, 400]])
-        ps = AveragedPowerspectrum(lc, norm="leahy", segment_size=100, silent=True)
-        rms_ps_l, rms_err_l = ps.compute_rms(min_freq=ps.freq[1], max_freq=ps.freq[-1])
-
-        ps = AveragedPowerspectrum(lc, norm=norm, segment_size=100)
-        rms_ps, rms_err = ps.compute_rms(min_freq=ps.freq[1], max_freq=ps.freq[-1])
-        assert np.allclose(rms_ps, rms_ps_l, atol=0.01)
-        assert np.allclose(rms_err, rms_err_l, atol=0.01)
-
-    @pytest.mark.parametrize("norm", ["frac", "abs", "none"])
-    def test_fractional_rms_in_frac_norm_is_consistent_averaged_freq_range(self, norm):
-        time = np.arange(0, 400, 1) + 0.5
-
-        data = Table.read(os.path.join(datadir, "sample_variable_series.fits"))["data"][:400] * 1000
-        poisson_counts = np.random.poisson(data)
-
-        lc = Lightcurve(time, counts=poisson_counts, dt=1, gti=[[0, 400]])
-        ps = AveragedPowerspectrum(lc, norm="leahy", segment_size=100, silent=True)
-        rms_ps_l, rms_err_l = ps.compute_rms(min_freq=ps.freq[5], max_freq=ps.freq[-5])
-
-        ps = AveragedPowerspectrum(lc, norm=norm, segment_size=100)
-        rms_ps, rms_err = ps.compute_rms(min_freq=ps.freq[5], max_freq=ps.freq[-5])
-        assert np.allclose(rms_ps, rms_ps_l, atol=0.01)
-        assert np.allclose(rms_err, rms_err_l, atol=0.01)
-
-    def test_fractional_rms_in_frac_norm_is_consistent_averaged_old(self):
-        with pytest.warns(DeprecationWarning, match="the option white_noise_offset"):
-            time = np.arange(0, 400, 1) + 0.5
-
-            poisson_counts = np.random.poisson(100.0, size=time.shape[0])
-
-            lc = Lightcurve(time, counts=poisson_counts, dt=1, gti=[[0, 400]])
-            ps = AveragedPowerspectrum(lc, norm="leahy", segment_size=100, silent=True)
-            rms_ps_l, rms_err_l = ps.compute_rms(
-                min_freq=ps.freq[1], max_freq=ps.freq[-1], white_noise_offset=0
-            )
-
-            ps = AveragedPowerspectrum(lc, norm="frac", segment_size=100)
-            rms_ps, rms_err = ps.compute_rms(
-                min_freq=ps.freq[1], max_freq=ps.freq[-1], white_noise_offset=0
-            )
-            assert np.allclose(rms_ps, rms_ps_l, atol=0.01)
-            assert np.allclose(rms_err, rms_err_l, atol=0.01)
-
-    def test_fractional_rms_in_frac_norm(self):
-        with pytest.warns(DeprecationWarning, match="the option white_noise_offset"):
-            time = np.arange(0, 400, 1) + 0.5
-
-            poisson_counts = np.random.poisson(100.0, size=time.shape[0])
-
-            lc = Lightcurve(time, counts=poisson_counts, dt=1, gti=[[0, 400]])
-            ps = AveragedPowerspectrum(lc, norm="frac", segment_size=100)
-            rms_ps, rms_err = ps.compute_rms(
-                min_freq=ps.freq[1], max_freq=ps.freq[-1], white_noise_offset=0
-            )
-            rms_lc = np.std(lc.counts) / np.mean(lc.counts)
-            assert np.isclose(rms_ps, rms_lc, atol=0.01)
-
     def test_leahy_norm_Poisson_noise(self):
         """
         In Leahy normalization, the poisson noise level (so, in the absence of
         a signal, the average power) should be equal to 2.
         """
         time = np.linspace(0, 10.0, 10**5)
-        counts = np.random.poisson(1000, size=time.shape[0])
+        counts = rng.poisson(1000, size=time.shape[0])
 
         lc = Lightcurve(time, counts)
         ps = Powerspectrum(lc, norm="leahy")
@@ -632,31 +476,15 @@ class TestPowerspectrum(object):
 
         assert np.isclose(ps_var, np.var(self.lc.counts), atol=0.01)
 
-    def test_fractional_rms_in_leahy_norm(self):
-        """
-        fractional rms should only be *approximately* equal the standard
-        deviation divided by the mean of the light curve. Therefore, we allow
-        for a larger tolerance in np.isclose()
-        """
-        with pytest.warns(DeprecationWarning):
-            ps = Powerspectrum(self.lc, norm="Leahy")
-            rms_ps, rms_err = ps.compute_rms(
-                min_freq=ps.freq[0], max_freq=ps.freq[-1], white_noise_offset=0
-            )
-
-            rms_lc = np.std(self.lc.counts) / np.mean(self.lc.counts)
-            assert np.isclose(rms_ps, rms_lc, atol=0.01)
-
     def test_abs_norm_Poisson_noise(self):
         """
         Poisson noise level for a light curve with absolute rms-squared
         normalization should be approximately 2 * the mean count rate of the
         light curve.
         """
-        np.random.seed(101)
 
         time = np.linspace(0, 1.0, 10**4)
-        counts = np.random.poisson(0.01, size=time.shape[0])
+        counts = rng.poisson(0.01, size=time.shape[0])
 
         lc = Lightcurve(time, counts)
         ps = Powerspectrum(lc, norm="abs")
@@ -787,7 +615,7 @@ class TestAveragedPowerspectrum(object):
         mean_count_rate = 1000.0
         mean_counts = mean_count_rate * dt
 
-        poisson_counts = np.random.poisson(mean_counts, size=time.shape[0])
+        poisson_counts = rng.poisson(mean_counts, size=time.shape[0])
 
         cls.lc = Lightcurve(time, counts=poisson_counts, gti=[[tstart, tend]], dt=dt)
 
@@ -870,7 +698,7 @@ class TestAveragedPowerspectrum(object):
 
         lc_all = []
         for n in range(n_lcs):
-            poisson_counts = np.random.poisson(mean_counts, size=len(time))
+            poisson_counts = rng.poisson(mean_counts, size=len(time))
 
             lc = Lightcurve(time, counts=poisson_counts, gti=[[tstart, tend]], dt=dt)
             lc_all.append(lc)
@@ -881,7 +709,7 @@ class TestAveragedPowerspectrum(object):
     def test_with_zero_counts(self):
         nbins = 100
         x = np.linspace(0, 10, nbins)
-        y0 = np.random.normal(loc=10, scale=0.5, size=int(0.4 * nbins))
+        y0 = rng.normal(loc=10, scale=0.5, size=int(0.4 * nbins))
         y1 = np.zeros(int(0.6 * nbins))
         y = np.hstack([y0, y1])
 
@@ -969,7 +797,7 @@ class TestAveragedPowerspectrum(object):
 
         lc_all = []
         for n in range(n_lcs):
-            poisson_counts = np.random.poisson(mean_counts, size=len(time))
+            poisson_counts = rng.poisson(mean_counts, size=len(time))
 
             lc = Lightcurve(time, counts=poisson_counts)
             lc_all.append(lc)
@@ -986,7 +814,7 @@ class TestAveragedPowerspectrum(object):
         lc_all = []
         for i in range(n):
             time = np.arange(0.0, 10.0, 10.0 / 10000)
-            counts = np.random.poisson(1000, size=time.shape[0])
+            counts = rng.poisson(1000, size=time.shape[0])
             lc = Lightcurve(time, counts)
             lc_all.append(lc)
 
@@ -1036,28 +864,6 @@ class TestDynamicalPowerspectrum(object):
         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):
             dps = DynamicalPowerspectrum(self.lc, segment_size=1000)
@@ -1190,7 +996,7 @@ class TestRoundTrip:
     def setup_class(cls):
         cls.cs = AveragedPowerspectrum()
         cls.cs.freq = np.arange(10)
-        cls.cs.power = np.random.uniform(0, 10, 10)
+        cls.cs.power = rng.uniform(0, 10, 10)
         cls.cs.m = 2
         cls.cs.nphots1 = 34
 
@@ -1252,3 +1058,64 @@ class TestRoundTrip:
         os.unlink(fname)
 
         self._check_equal(so, new_so)
+
+
+class TestRMS(object):
+    @classmethod
+    def setup_class(cls):
+        fwhm = 0.23456
+        cls.segment_size = 256
+        cls.df = 1 / cls.segment_size
+
+        cls.freqs = np.arange(cls.df, 1, cls.df)
+        dt = 0.5 / cls.freqs.max()
+
+        pds_shape_func = Lorentz1D(x_0=0, fwhm=fwhm)
+        cls.pds_shape_raw = pds_shape_func(cls.freqs)
+        cls.M = 1000
+        cls.nphots = 1_000_000
+        cls.rms = 0.5
+        meanrate = cls.nphots / cls.segment_size
+        cls.poisson_noise_rms = 2 / meanrate
+        pds_shape_rms = cls.pds_shape_raw / np.sum(cls.pds_shape_raw * cls.df) * cls.rms**2
+        pds_shape_rms += cls.poisson_noise_rms
+
+        random_part = rng.chisquare(2 * cls.M, size=cls.pds_shape_raw.size) / 2 / cls.M
+        pds_rms_noisy = random_part * pds_shape_rms
+
+        pds_unnorm = pds_rms_noisy * meanrate / 2 * cls.nphots
+        cls.pds = AveragedPowerspectrum()
+        cls.pds.freq = cls.freqs
+        cls.pds.unnorm_power = pds_unnorm
+        cls.pds.power = pds_rms_noisy
+        cls.pds.df = cls.df
+        cls.pds.m = cls.M
+        cls.pds.nphots = cls.nphots
+        cls.pds.norm = "frac"
+        cls.pds.dt = dt
+        cls.pds.n = cls.pds.freq.size
+
+    @pytest.mark.parametrize("norm", ["none", "frac", "leahy", "abs"])
+    def test_rms(self, norm):
+        pds = self.pds.to_norm(norm)
+        rms_from_ps, rmse_from_ps = pds.compute_rms(self.freqs.min(), self.freqs.max())
+        assert np.isclose(rms_from_ps, self.rms, atol=3 * rmse_from_ps)
+
+    def test_rms_M_low(self):
+        """Test that the warning is raised when M is low."""
+
+        pds = copy.deepcopy(self.pds)
+        pds.m = 23
+
+        with pytest.warns(UserWarning, match="All power spectral bins have M<30."):
+            pds.compute_rms(
+                self.freqs.min(), self.freqs.max(), poisson_noise_level=self.poisson_noise_rms
+            )
+
+    @pytest.mark.parametrize("norm", ["none", "frac", "leahy", "abs"])
+    def test_rms_rebinning(self, norm):
+        pds = self.pds.to_norm(norm)
+        pds = pds.rebin_log(0.04)
+        rms_from_ps, rmse_from_ps = pds.compute_rms(self.freqs.min(), self.freqs.max())
+
+        assert np.isclose(rms_from_ps, self.rms, atol=3 * rmse_from_ps)
diff --git a/stingray/tests/test_stats.py b/stingray/tests/test_stats.py
index d7306b63..d859ce97 100644
--- a/stingray/tests/test_stats.py
+++ b/stingray/tests/test_stats.py
@@ -186,7 +186,7 @@ class TestClassicalSignificances(object):
     def test_very_large_powers_produce_zero_prob(self):
         power = 31000.0
         nspec = 1
-        with pytest.warns(DeprecationWarning):
+        with pytest.warns((DeprecationWarning, UserWarning)):
             pval = classical_pvalue(power, nspec)
         assert np.isclose(pval, 0.0)