[BACKPORT] [BUG] Fix normalisation procedure for antisymmetrised PAC (#115) (#139)

Author: tsbinns

changelog.md

@@ -4,7 +4,7 @@
 - No changes.
 
 
-## [Version 1.2](https://pybispectra.readthedocs.io/en/1.2.2/)
+## [Version 1.2](https://pybispectra.readthedocs.io/en/1.2.3/)
 
 ##### Enhancements
 - Added general `Bispectrum` and `Threenorm` classes for computing with flexible kmn channel combinations.
@@ -14,6 +14,7 @@
 - Fixed error where bandpass filter settings for the SSD method in `SpatioSpectralFilter` were not being applied correctly.
 - Fixed error where `indices` in `ResultsCFC`, `ResultsTDE`, and `ResultsGeneral` classes were not being mapped to results correctly.
 - Fixed error where NumPy integers and floats were not being recognised as valid types.
+- Fixed error where univariate normalisation of antisymmetrised PAC was not being applied correctly.
 
 ##### API
 - Changed the default value of `min_ratio` in `SpatioSpectralFilter.get_transformed_data()` from `1.0` to `-inf`.

docs/source/refs.bib

@@ -44,6 +44,18 @@ @article{Chella2014
   doi={10.1016/j.neuroimage.2013.12.064}
 }
 
+@article{Chella2016,
+  title={Bispectral pairwise interacting source analysis for identifying systems of cross-frequency interacting brain sources from electroencephalographic or magnetoencephalographic signals},
+  author={Chella, Federico and Pizzella, Vittorio and Zappasodi, Filippo and Nolte, Guido and Marzetti, Laura},
+  journal={Physical Review E},
+  volume={93},
+  number={5},
+  pages={052420},
+  year={2016},
+  publisher={APS},
+  doi={10.1103/PhysRevE.93.052420}
+}
+
 @inbook{Chen2004,
   title={Time {D}elay {E}stimation},
   author={Chen, Jingdong and Huang, Yiteng and Benesty, Jacob},

examples/plot_compute_pac.py

@@ -54,16 +54,16 @@
 #
 # where the resulting values lie in the range :math:`[0, 1]`. Furthermore, PAC can be
 # antisymmetrised by subtracting the results from those found using the transposed
-# bispectrum, :math:`\textbf{B}_{xyx}`, :footcite:`Chella2014`
+# bispectrum, :math:`\textbf{B}_{yxy}`, :footcite:`Chella2014`
 #
 # :math:`\textrm{PAC}_{\textrm{antisym}}(\textbf{x}_{f_1},\textbf{y}_{f_2})=|
-# \textbf{B}_{xyy}-\textbf{B}_{xyx}|` .
+# \textbf{B}_{xyy}(f_1,f_2)-\textbf{B}_{yxy}(f_1,f_2)|` .
 #
 # In the context of analysing PAC between two signals, antisymmetrisation allows you to
 # correct for spurious estimates of coupling arising from interactions within the
 # signals themselves in instances of source mixing, providing a more robust connectivity
 # metric :footcite:`PellegriniPreprint`. The same principle applies for the
-# antisymmetrisation of the bicoherence.
+# antisymmetrisation of the bicoherence :footcite:`Chella2016`.
 
 ########################################################################################
 # Loading data and computing Fourier coefficients

src/pybispectra/cfc/pac.py

@@ -140,10 +140,17 @@ def compute(
 
         where the resulting values lie in the range :math:`[0, 1]`. Furthermore, PAC can
         be antisymmetrised by subtracting the results from those found using the
-        transposed bispectrum, :math:`\textbf{B}_{xyx}`, :footcite:`Chella2014`
+        transposed bispectrum, :math:`\textbf{B}_{yxy}`, :footcite:`Chella2014`
 
         :math:`\textrm{PAC}_{\textrm{antisym}}(\textbf{x}_{f_1},\textbf{y}_{f_2})=
-        |\textbf{B}_{xyy}-\textbf{B}_{xyx}|` .
+        |\textbf{B}_{xyy}(f_1,f_2)-\textbf{B}_{yxy}(f_1,f_2)|` .
+
+        A modified approach is used for the normalisation of antisymmetrised PAC
+        :footcite:`Chella2016`
+
+        :math:`\textrm{PAC}_{\textrm{norm,antisym}}(\textbf{x}_{f_1},\textbf{y}_{f_2})=
+        \Large|\frac{\textbf{B}_{xyy}(f_1,f_2)-\textbf{B}_{yxy}(f_1,f_2)}
+        {\textbf{N}_{xyy}(f_1,f_2)+\textbf{N}_{yxy}(f_1,f_2)}|` .
 
         If the seed and target for a given connection is the same channel and
         antisymmetrisation is being performed, :obj:`numpy.nan` values are returned.
@@ -171,7 +178,6 @@ def compute(
         self._compute_bispectrum()
         if self._return_threenorm:
             self._compute_threenorm()
-            self._bicoherence = self._bispectrum / self._threenorm
         self._compute_pac()
         self._store_results()
 
@@ -189,7 +195,6 @@ def _reset_attrs(self) -> None:
 
         self._bispectrum = None
         self._threenorm = None
-        self._bicoherence = None
 
         self._pac_nosym_nonorm = None
         self._pac_nosym_threenorm = None
@@ -340,10 +345,13 @@ def _compute_pac(self) -> None:
 
         if self._return_threenorm:
             if self._return_nosym:
-                self._pac_nosym_threenorm = np.abs(self._bicoherence[0])
+                self._pac_nosym_threenorm = np.abs(
+                    self._bispectrum[0] / self._threenorm[0]
+                )
             if self._return_antisym:
                 self._pac_antisym_threenorm = np.abs(
-                    self._bicoherence[0] - self._bicoherence[1]
+                    (self._bispectrum[0] - self._bispectrum[1])
+                    / (self._threenorm[0] + self._threenorm[1])
                 )
                 for con_i, seed, target in zip(
                     range(self._n_cons), self._seeds, self._targets