Repository: StingraySoftware/stingray

Issue:
Cross and Auto-correlation normalisation
Hi,

I'm not entirely sure the auto and cross-correlation functions returned by Stingray are normalised, as the output seems to depend on the amplitude of the input light-curves. I think it might be good to include an option, such that the returned correlation functions are =1 for a perfect correlation, =-1 for perfect anti-correlation, and -1<ccf/acf<1 otherwise. So essentially dividing through by the sqrt of the variance of each light-curve (see e.g Gaskell & Peterson 1987, Gardner & Done 2017)

The attached figure demonstrates the normalisation issue more clearly, where the cross-correlation has been calculated for two sin curves, offset slightly from each other (top panel). It can be seen that when changing the amplitude of the input light-curves (in this case by a factor of 0.25), the amplitude of the resulting cross-correlation also changes (blue and orange curves respectively). Also shown is the result from a test code (green), following Gardner & Done (2017), showing the desired normalisation. (I'd be happy to send you this code - although it's worth noting it still needs a fair bit of streamlining and optimising; so not currently very user friendly!)

![StingrayCrossCorrTest](https://user-images.githubusercontent.com/99818820/154451239-8f6a4024-c714-4373-88cf-d0dd0990da58.png)

Cheers,
Scott





Discussion:

Patch:
diff --git a/stingray/crosscorrelation.py b/stingray/crosscorrelation.py
index 36bb81ad..0b238162 100644
--- a/stingray/crosscorrelation.py
+++ b/stingray/crosscorrelation.py
@@ -39,6 +39,11 @@ class CrossCorrelation(object):
         See the relevant ``scipy`` documentation [scipy-docs]_
         for more details.
 
+    norm: {``none``, ``variance``}
+        if "variance", the cross correlation is normalized so that perfect
+        correlation gives 1, and perfect anticorrelation gives -1. See
+        Gaskell \& Peterson 1987, Gardner \& Done 2017
+
     Attributes
     ----------
     lc1: :class:`stingray.Lightcurve`
@@ -70,15 +75,19 @@ class CrossCorrelation(object):
     auto: bool
         An internal flag to indicate whether this is a cross-correlation or an auto-correlation.
 
+    norm: {``none``, ``variance``}
+        The normalization specified in input
+
     References
     ----------
     .. [scipy-docs] https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.signal.correlate.html
 
     """
 
-    def __init__(self, lc1=None, lc2=None, cross=None, mode='same'):
+    def __init__(self, lc1=None, lc2=None, cross=None, mode='same', norm="none"):
 
         self.auto = False
+        self.norm = norm
         if isinstance(mode, str) is False:
             raise TypeError("mode must be a string")
 
@@ -187,10 +196,22 @@ class CrossCorrelation(object):
         lc2_counts = self.lc2.counts - np.mean(self.lc2.counts)
 
         # Calculates cross-correlation of two lightcurves
-        self.corr = signal.correlate(lc1_counts, lc2_counts, self.mode)
-        self.n = len(self.corr)
+        self.corr = \
+            signal.correlate(lc1_counts, lc2_counts, self.mode)
+
+        self.n = np.size(self.corr)
         self.time_shift, self.time_lags, self.n = self.cal_timeshift(dt=self.dt)
 
+        # Normalization that makes the maximum correlation equal to 1, and
+        # maximum anticorrelation -1.
+        if self.norm == "variance":
+            # Note that self.corr is normalized so that the maximum is
+            # proportional to the number of bins in the first input
+            # light curve. Hence, the division by the lc size
+            variance1 = np.var(lc1.counts) - np.mean(lc1.counts_err)**2
+            variance2 = np.var(lc2.counts) - np.mean(lc2.counts_err)**2
+            self.corr = self.corr / np.sqrt(variance1 * variance2) / lc1_counts.size
+
     def cal_timeshift(self, dt=1.0):
         """
         Calculate the cross correlation against all possible time lags, both positive and negative.
@@ -212,7 +233,6 @@ class CrossCorrelation(object):
         """
         if self.dt is None:
             self.dt = dt
-
         if self.corr is None:
             if (self.lc1 is None or self.lc2 is None) and (self.cross is None):
                 raise StingrayError('Please provide either two lightcurve objects or \


Test patch:
diff --git a/stingray/tests/test_crosscorrelation.py b/stingray/tests/test_crosscorrelation.py
index a2333d74..07052fe5 100644
--- a/stingray/tests/test_crosscorrelation.py
+++ b/stingray/tests/test_crosscorrelation.py
@@ -14,6 +14,30 @@ import matplotlib.pyplot as plt
 
 
 
+class TestCrossCorrelationBase(object):
+    @classmethod
+    def setup_class(cls):
+        dt = 0.01
+        length = 10000
+        gti = [[0, length]]
+        times = np.arange(dt/2, length, dt)
+        freq = 1 / 50
+        flux1 = 0.5 + 0.5 * np.sin(2 * np.pi * freq * times)
+        flux2 = 0.5 + 0.5 * np.sin(2 * np.pi * freq * (times - 20))
+
+        cls.lc1 = Lightcurve(times, flux1, dt=dt, err_dist="gauss", gti=gti, skip_checks=True)
+        cls.lc2 = Lightcurve(times, flux2, dt=dt, err_dist="gauss", gti=gti, skip_checks=True)
+
+    def test_crosscorr(self):
+        cr = CrossCorrelation(self.lc1, self.lc2)
+        assert np.isclose(cr.time_shift, -20, atol=0.1)
+
+    def test_crosscorr_norm(self):
+        cr = CrossCorrelation(self.lc1, self.lc2, norm="variance")
+        assert np.isclose(cr.time_shift, -20, atol=0.1)
+        assert np.isclose(np.max(cr.corr), 1, atol=0.01)
+        assert np.isclose(np.min(cr.corr), -1, atol=0.01)
+
 class TestCrossCorrelation(object):
     @classmethod
     def setup_class(cls):