Repository: StingraySoftware/stingray
Pull request #400 (2019-05-10 20:35:23)

[WIP] Fix of bug in definition of absolute rms normalization of the PSD
This is basically PR #388, except hopefully without merge conflicts?

@abigailStev: Can you try to check out this branch (not fork it), and then push to this PR directly?


Patch:
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index e751b2bf..a7ecfe82 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -190,8 +190,8 @@ class Crossspectrum(object):
             Two light curves used for computing the cross spectrum.
         """
         if lc1 is not lc2 and isinstance(lc1, Lightcurve):
-            self.pds1 = Crossspectrum(lc1, lc1, norm='none', power_type=self.power_type)
-            self.pds2 = Crossspectrum(lc2, lc2, norm='none', power_type=self.power_type)
+            self.pds1 = Crossspectrum(lc1, lc1, norm='none')
+            self.pds2 = Crossspectrum(lc2, lc2, norm='none')
 
     def _make_crossspectrum(self, lc1, lc2):
         """
@@ -320,7 +320,7 @@ class Crossspectrum(object):
         fourier_2 = scipy.fftpack.fft(lc2.counts)  # do Fourier transform 2
 
         freqs = scipy.fftpack.fftfreq(lc1.n, lc1.dt)
-        cross = fourier_1[freqs > 0] * np.conj(fourier_2[freqs > 0])
+        cross = np.multiply(fourier_1[freqs > 0], np.conj(fourier_2[freqs > 0]))
 
         return freqs[freqs > 0], cross
 
@@ -423,7 +423,7 @@ class Crossspectrum(object):
         """
 
         # The "effective" counts/bin is the geometrical mean of the counts/bin
-        # of the two light curves
+        # of the two light curves. Same goes for counts/second in meanrate.
 
         log_nphots1 = np.log(self.nphots1)
         log_nphots2 = np.log(self.nphots2)
@@ -431,6 +431,8 @@ class Crossspectrum(object):
         actual_nphots = np.float64(np.sqrt(np.exp(log_nphots1 + log_nphots2)))
         actual_mean = np.sqrt(self.meancounts1 * self.meancounts2)
 
+        meanrate = np.sqrt((self.nphots1 * self.nphots2) / tseg)
+
         assert actual_mean > 0.0, \
             "Mean count rate is <= 0. Something went wrong."
 
@@ -444,16 +446,14 @@ class Crossspectrum(object):
             raise ValueError("`power_type` not recognized!")
 
         if self.norm.lower() == 'leahy':
-            c = c_num
-            power = c * 2. / actual_nphots
+            power = c_num * 2. / actual_nphots
 
         elif self.norm.lower() == 'frac':
             c = c_num / np.float(self.n ** 2.)
             power = c * 2. * tseg / (actual_mean ** 2.0)
 
         elif self.norm.lower() == 'abs':
-            c = c_num / np.float(self.n ** 2.)
-            power = c * (2. * tseg)
+            power = c_num * 2. * meanrate / actual_nphots
 
         elif self.norm.lower() == 'none':
             power = unnorm_power


Test patch:
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index b92c7fd3..e0847bb5 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -102,10 +102,9 @@ class TestCrossspectrum(object):
 
         counts1 = np.random.poisson(0.01, size=time.shape[0])
         counts2 = np.random.negative_binomial(1, 0.09, size=time.shape[0])
-
         self.lc1 = Lightcurve(time, counts1, gti=[[tstart, tend]], dt=dt)
         self.lc2 = Lightcurve(time, counts2, gti=[[tstart, tend]], dt=dt)
-
+        self.rate1 = 100.  # mean count rate (counts/sec) of light curve 1
         self.cs = Crossspectrum(self.lc1, self.lc2)
 
     def test_make_empty_crossspectrum(self):
@@ -201,20 +200,32 @@ class TestCrossspectrum(object):
         assert type(new_cs) == type(self.cs)
         new_cs.time_lag()
 
+    def test_norm_abs(self):
+        # Testing for a power spectrum of lc1
+        cs = Crossspectrum(lc1=self.lc1, lc2=self.lc1, norm='abs')
+        assert len(cs.power) == 4999
+        assert cs.norm == 'abs'
+        abs_noise = 2. * self.rate1  # expected Poisson noise level
+        print(np.mean(cs.power), abs_noise)
+        assert np.isclose(np.mean(cs.power[1:]), abs_noise, atol=30)
+
     def test_norm_leahy(self):
-        cs = Crossspectrum(self.lc1, self.lc2, norm='leahy')
+        # Testing for a power spectrum of lc1
+        cs = Crossspectrum(lc1=self.lc1, lc2=self.lc1, norm='leahy')
         assert len(cs.power) == 4999
         assert cs.norm == 'leahy'
+        leahy_noise = 2.0  # expected Poisson noise level
+        print(np.mean(cs.power), leahy_noise)
+        assert np.isclose(np.mean(cs.power[1:]), leahy_noise, atol=0.2)
 
     def test_norm_frac(self):
-        cs = Crossspectrum(self.lc1, self.lc2, norm='frac')
+        # Testing for a power spectrum of lc1
+        cs = Crossspectrum(lc1=self.lc1, lc2=self.lc1, norm='frac')
         assert len(cs.power) == 4999
         assert cs.norm == 'frac'
-
-    def test_norm_abs(self):
-        cs = Crossspectrum(self.lc1, self.lc2, norm='abs')
-        assert len(cs.power) == 4999
-        assert cs.norm == 'abs'
+        frac_noise = 2. / self.rate1  # expected Poisson noise level
+        print(np.mean(cs.power), frac_noise)
+        assert np.isclose(np.mean(cs.power[1:]), frac_noise, atol=0.005)
 
     def test_failure_when_normalization_not_recognized(self):
         with pytest.raises(ValueError):
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index fad3137d..54d9e269 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -102,7 +102,7 @@ class TestPowerspectrum(object):
 
     def test_frac_normalization_correct(self):
         """
-        In rms normalization, the integral of the powers should be
+        In fractional rms normalization, the integral of the powers should be
         equal to the variance of the light curve divided by the mean
         of the light curve squared.
         """
@@ -162,7 +162,11 @@ class TestPowerspectrum(object):
         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):
+    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, 1e5)
         counts = np.random.poisson(1000, size=time.shape[0])
 
@@ -202,6 +206,23 @@ class TestPowerspectrum(object):
             rms_ps, rms_err = ps.compute_rms(min_freq=ps.freq[0],
                                              max_freq=ps.freq[-1])
 
+    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.
+        """
+        time = np.linspace(0, 1., 1e4)
+        counts = np.random.poisson(0.01, size=time.shape[0])
+
+        lc = Lightcurve(time, counts)
+        ps = Powerspectrum(lc, norm="abs")
+        print(lc.counts/lc.tseg)
+        abs_noise = 2. * 100  # expected Poisson noise level;
+                              # hardcoded value from above
+        print(np.mean(ps.power[1:]), abs_noise)
+        assert np.isclose(np.mean(ps.power[1:]), abs_noise, atol=30)
+
     def test_fractional_rms_error(self):
         """
         TODO: Need to write a test for the fractional rms error.