add fooof documentation

neuromodulation · Nov 15, 2023 · df7558a · df7558a
1 parent e95465d
commit df7558a
Show file tree

Hide file tree

Showing 3 changed files with 46 additions and 14 deletions.
diff --git a/docs/source/usage.rst b/docs/source/usage.rst
@@ -427,6 +427,9 @@ aperiodic component. *fooof* specific parameters, e.g. *knee* or *max_n_peaks* a
        "knee": true
    }
 
+.. note::
+    When using the knee parameter, the *knee_frequency* is returned for every fit. See also the fooof `Aperiodic Component Fitting Notebook <https://fooof-tools.github.io/fooof/auto_tutorials/plot_05-AperiodicFitting.html#sphx-glr-auto-tutorials-plot-05-aperiodicfitting-py>`_. 
+
 Nonlinear measures for dynamical systems (nolds)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

diff --git a/py_neuromodulation/nm_fooof.py b/py_neuromodulation/nm_fooof.py
@@ -79,40 +79,67 @@ def calc_feature(
             except Exception as e:
                 print(e)
                 print(f"failing spectrum: {spectrum}")
+
+            if self.fm.fooofed_spectrum_ is None:
+                FIT_PASSED = False
+            else:
+                FIT_PASSED = True
+
             if self.settings_fooof["aperiodic"]["exponent"]:
                 features_compute[f"{ch_name}_fooof_a_exp"] = (
-                    np.nan_to_num(self.fm.get_params("aperiodic_params", "exponent"))
-                    if self.fm.fooofed_spectrum_ is not None
+                    np.nan_to_num(
+                        self.fm.get_params("aperiodic_params", "exponent")
+                    )
+                    if FIT_PASSED is True
                     else None
                 )
 
             if self.settings_fooof["aperiodic"]["offset"]:
                 features_compute[f"{ch_name}_fooof_a_offset"] = (
-                    np.nan_to_num(self.fm.get_params("aperiodic_params", "offset"))
-                    if self.fm.fooofed_spectrum_ is not None
+                    np.nan_to_num(
+                        self.fm.get_params("aperiodic_params", "offset")
+                    )
+                    if FIT_PASSED is True
                     else None
                 )
-            
+
             if self.settings_fooof["aperiodic"]["knee"]:
-                features_compute[f"{ch_name}_fooof_a_knee"] = (
-                    np.nan_to_num(self.fm.get_params("aperiodic_params", "knee"))
-                    if self.fm.fooofed_spectrum_ is not None
-                    else None
-                )
+                if FIT_PASSED is False:
+                    knee_freq = None
+                else:
+                    if self.fm.get_params("aperiodic_params", "exponent") != 0:
+                        knee_fooof = self.fm.get_params(
+                            "aperiodic_params", "knee"
+                        )
+                        knee_freq = np.nan_to_num(
+                            knee_fooof
+                            ** (
+                                1
+                                / self.fm.get_params(
+                                    "aperiodic_params", "exponent"
+                                )
+                            )
+                        )
+                    else:
+                        knee_freq = None
+
+                features_compute[
+                    f"{ch_name}_fooof_a_knee_frequency"
+                ] = knee_freq
 
             peaks_bw = (
                 self.fm.get_params("peak_params", "BW")
-                if self.fm.fooofed_spectrum_ is not None
+                if FIT_PASSED is True
                 else None
             )
             peaks_cf = (
                 self.fm.get_params("peak_params", "CF")
-                if self.fm.fooofed_spectrum_ is not None
+                if FIT_PASSED is True
                 else None
             )
             peaks_pw = (
                 self.fm.get_params("peak_params", "PW")
-                if self.fm.fooofed_spectrum_ is not None
+                if FIT_PASSED is True
                 else None
             )
 
@@ -122,7 +149,6 @@ def calc_feature(
                 peaks_pw = [peaks_pw]
 
             for peak_idx in range(self.max_n_peaks):
-
                 if self.settings_fooof["periodic"]["band_width"]:
                     features_compute[f"{ch_name}_fooof_p_{peak_idx}_bw"] = (
                         peaks_bw[peak_idx] if peak_idx < len(peaks_bw) else None

diff --git a/tests/test_bursts.py b/tests/test_bursts.py
@@ -87,6 +87,9 @@ def test_bursting_duration():
         f = bursts.calc_feature(data, {})
 
     np.random.seed(0)
+    # the percentile of the hilbert transform of a continuous oscillation will be high
+    # select better max amplitude
+    bursts = nm_bursts.Burst(settings, ch_names, sfreq)
     f_burst = bursts.calc_feature(
         beta_wave + np.random.random([NUM_CH, 1 * sfreq]), {}
     )