Repository: StingraySoftware/stingray (python)

Issue:
In `Powerspectrum`, the calculation of rms and rms error are wrong
I could see two problems with the method `compute_rms` and the related `_rms_error`:

1. It just sums over the powers, without eliminating the white noise level

2. The error is always assumed in Leahy normalization!


Pull request description:
Fix rms error
Resolves #285
Now it is possible to subtract white noise when calculating the rms.
Fixes (I think) an inconsistency in the normalization factors used when calculating error bars on the rms.
Need to check that the values coming out of this make sense in different conditions

Patch:
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 5366648f..4385b305 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -695,6 +695,7 @@ class AveragedCrossspectrum(Crossspectrum):
 
         power_avg /= np.float(m)
         power_err_avg = np.sqrt(power_err_avg) / m
+        unnorm_power_avg /= np.float(m)
 
         self.freq = self.cs_all[0].freq
         self.power = power_avg
@@ -737,7 +738,7 @@ class AveragedCrossspectrum(Crossspectrum):
                   "expected statistical distributions.")
 
         # Calculate average coherence
-        unnorm_power_avg = self.unnorm_power / self.m
+        unnorm_power_avg = self.unnorm_power
 
         num = np.absolute(unnorm_power_avg) ** 2
 
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index a0c20950..e49eaf74 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -182,7 +182,7 @@ class Powerspectrum(Crossspectrum):
 
         return bin_ps
 
-    def compute_rms(self, min_freq, max_freq):
+    def compute_rms(self, min_freq, max_freq, white_noise_offset=0.):
         """
         Compute the fractional rms amplitude in the periodgram
         between two frequencies.
@@ -195,6 +195,14 @@ class Powerspectrum(Crossspectrum):
         max_freq: float
             The upper frequency bound for the calculation
 
+        Other parameters
+        ----------------
+        white_noise_offset : float, default 0
+            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
@@ -205,16 +213,18 @@ class Powerspectrum(Crossspectrum):
         minind = self.freq.searchsorted(min_freq)
         maxind = self.freq.searchsorted(max_freq)
         powers = self.power[minind:maxind]
+        nphots = self.nphots
 
         if self.norm.lower() == 'leahy':
-            rms = np.sqrt(np.sum(powers) / self.nphots)
-
+            powers_leahy = powers.copy()
         elif self.norm.lower() == "frac":
-            rms = np.sqrt(np.sum(powers * self.df))
+            powers_leahy = \
+                self.unnorm_power[minind:maxind].real * 2 / nphots
         else:
             raise TypeError("Normalization not recognized!")
 
-        rms_err = self._rms_error(powers)
+        rms = np.sqrt(np.sum(powers_leahy - white_noise_offset) / nphots)
+        rms_err = self._rms_error(powers_leahy)
 
         return rms, rms_err
 
@@ -226,6 +236,10 @@ class Powerspectrum(Crossspectrum):
         strictly correct. We should be using the underlying power spectrum,
         but in the absence of an estimate of that, this will have to do.
 
+        $r = \sqrt{P}$
+
+        $\delta r = \frac{1}{2 * \sqrt{P}} \delta P$
+
         Parameters
         ----------
         powers: iterable
@@ -236,9 +250,16 @@ class Powerspectrum(Crossspectrum):
         delta_rms: float
             the error on the fractional rms amplitude
         """
-        p_err = scipy.stats.chi2(2.0 * self.m).var() * powers / self.m
-        drms_dp = 1 / (2 * np.sqrt(np.sum(powers) * self.df))
-        delta_rms = np.sum(p_err * drms_dp * self.df)
+        nphots = self.nphots
+        p_err = scipy.stats.chi2(2.0 * self.m).var() * powers / self.m / nphots
+
+        rms = np.sum(powers) / nphots
+        pow = np.sqrt(rms)
+
+        drms_dp = 1 / (2 * pow)
+
+        sq_sum_err = np.sqrt(np.sum(p_err**2))
+        delta_rms = sq_sum_err * drms_dp
         return delta_rms
 
     def classical_significances(self, threshold=1, trial_correction=False):


Test patch:
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index e8dfd46d..5c4803f6 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -113,11 +113,55 @@ class TestPowerspectrum(object):
         print(std_lc)
         assert np.isclose(ps_int, std_lc, 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=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=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)
+
+    def test_fractional_rms_in_frac_norm_is_consistent_averaged(self):
+        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=lc, norm="leahy", segment_size=100)
+        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=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):
-        ps = Powerspectrum(lc=self.lc, norm="frac")
+        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=lc, norm="frac", segment_size=100)
         rms_ps, rms_err = ps.compute_rms(min_freq=ps.freq[1],
-                                         max_freq=ps.freq[-1])
-        rms_lc = np.std(self.lc.counts) / np.mean(self.lc.counts)
+                                         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_correct(self):