meg data analysis

py_neuromodulation/nm_bursts.py

@@ -1,6 +1,7 @@
 import enum
 import numpy as np
 from typing import Iterable
+from scipy import signal
 
 from py_neuromodulation import nm_features_abc, nm_filter
 
@@ -9,7 +10,6 @@ class Burst(nm_features_abc.Feature):
     def __init__(
         self, settings: dict, ch_names: Iterable[str], sfreq: float
     ) -> None:
-
         self.s = settings
         self.sfreq = sfreq
         self.ch_names = ch_names
@@ -34,7 +34,9 @@ def __init__(
             )
         ).astype(int)
 
-        self.num_max_samples_ring_buffer = int(self.sfreq * self.time_duration_s)
+        self.num_max_samples_ring_buffer = int(
+            self.sfreq * self.time_duration_s
+        )
 
         self.bandpass_filter = nm_filter.BandPassFilter(
             f_ranges=self.f_ranges,
@@ -62,35 +64,40 @@ def test_settings(
         ch_names: Iterable[str],
         sfreq: int | float,
     ):
-        assert (
-            isinstance(settings["burst_settings"]["threshold"], (float, int))
+        assert isinstance(
+            settings["burst_settings"]["threshold"], (float, int)
         ), f"burst settings threshold needs to be type int or float, got: {settings['burst_settings']['threshold']}"
         assert (
             0 < settings["burst_settings"]["threshold"] < 100
         ), f"burst setting threshold needs to be between 0 and 100, got: {settings['burst_settings']['threshold']}"
-        assert (
-            isinstance(settings["burst_settings"]["time_duration_s"], (float, int))
+        assert isinstance(
+            settings["burst_settings"]["time_duration_s"], (float, int)
         ), f"burst settings time_duration_s needs to be type int or float, got: {settings['burst_settings']['time_duration_s']}"
         assert (
             settings["burst_settings"]["time_duration_s"] > 0
         ), f"burst setting time_duration_s needs to be greater than 0, got: {settings['burst_settings']['time_duration_s']}"
 
         for fband_burst in settings["burst_settings"]["frequency_bands"]:
-            assert (
-                fband_burst in list(settings["frequency_ranges_hz"].keys())
+            assert fband_burst in list(
+                settings["frequency_ranges_hz"].keys()
             ), f"bursting {fband_burst} needs to be defined in settings['frequency_ranges_hz']"
 
-        for burst_feature in settings["burst_settings"]["burst_features"].keys():
+        for burst_feature in settings["burst_settings"][
+            "burst_features"
+        ].keys():
             assert isinstance(
-                settings["burst_settings"]["burst_features"][burst_feature], bool
+                settings["burst_settings"]["burst_features"][burst_feature],
+                bool,
             ), (
                 f"bursting feature {burst_feature} needs to be type bool, "
                 f"got: {settings['burst_settings']['burst_features'][burst_feature]}"
             )
-    def calc_feature(self, data: np.array, features_compute: dict) -> dict:
 
+    def calc_feature(self, data: np.array, features_compute: dict) -> dict:
         # filter_data returns (n_channels, n_fbands, n_samples)
-        filtered_data = self.bandpass_filter.filter_data(data)
+        filtered_data = np.abs(
+            signal.hilbert(self.bandpass_filter.filter_data(data), axis=2)
+        )
         for ch_idx, ch_name in enumerate(self.ch_names):
             for fband_idx, fband_name in enumerate(self.fband_names):
                 new_dat = filtered_data[ch_idx, fband_idx, :]
@@ -99,7 +106,7 @@ def calc_feature(self, data: np.array, features_compute: dict) -> dict:
                 else:
                     self.data_buffer[ch_name][fband_name] = np.concatenate(
                         (self.data_buffer[ch_name][fband_name], new_dat), axis=0
-                    )[-self.num_max_samples_ring_buffer:]
+                    )[-self.num_max_samples_ring_buffer :]
 
                 # calc features
                 burst_thr = np.percentile(

py_neuromodulation/nm_plots.py

@@ -449,7 +449,7 @@ def plot_all_features(
 
     cols_plt = [c for c in df.columns if c != "time"]
     if normalize is True:
-        data_plt = stats.zscore(df[cols_plt])
+        data_plt = stats.zscore(df[cols_plt], nan_policy="omit")
     else:
         data_plt = df[cols_plt]
 

pyproject.toml

@@ -4,7 +4,7 @@ requires = ["flit_core >=3.2,<4"]
 
 [project]
 name = "py_neuromodulation"
-version = "0.0.1"
+version = "0.0.2"
 authors = [{name = "Timon Merk", email="timon.merk@charite.de"}]
 classifiers = [
   "Development Status :: 2 - Pre-Alpha",

test_meg_data.py

@@ -0,0 +1,98 @@
+import numpy as np
+from matplotlib import pyplot as plt
+
+import py_neuromodulation as py_nm
+
+from py_neuromodulation import (
+    nm_analysis,
+    nm_define_nmchannels,
+    nm_plots,
+    nm_IO,
+)
+
+dat = nm_IO.loadmat(r"E:\Downloads\subj_1_off.mat")["data"]
+
+data = dat["trial"]
+sfreq = dat["hdr"]["Fs"]
+ch_names = dat["label"][[0, 1, 8, 11]]
+ch_types = dat["hdr"]["chantype"][[0, 1, 8, 11]]
+
+nm_channels = nm_define_nmchannels.set_channels(
+    ch_names=ch_names,
+    ch_types=ch_types,
+    reference="default",
+    bads=None,
+    new_names="default",
+    used_types=("lfp", "emg", "eog"),
+    target_keywords=None,
+)
+
+settings = py_nm.nm_settings.get_default_settings()
+settings = py_nm.nm_settings.reset_settings(settings)
+
+settings["features"]["fft"] = True
+settings["features"]["bursts"] = True
+settings["features"]["sharpwave_analysis"] = True
+settings["features"]["fooof"] = True
+
+settings["frequency_ranges_hz"] = {
+    "theta": [4, 8],
+    "alpha": [8, 12],
+    "low beta": [12, 20],
+    "high beta": [20, 30],
+    "low gamma": [30, 60],
+}
+
+settings["features"]["raw_hjorth"] = True
+settings["features"]["bispectrum"] = True
+settings["features"]["coherence"] = True
+settings["coherence"]["channels"] = [[ch_names[0], ch_names[2]]]
+settings["coherence"]["frequency_bands"] = ["high beta", "alpha"]
+settings["sharpwave_analysis_settings"]["estimator"]["mean"] = []
+
+for sw_feature in list(
+    settings["sharpwave_analysis_settings"]["sharpwave_features"].keys()
+):
+    settings["sharpwave_analysis_settings"]["sharpwave_features"][
+        sw_feature
+    ] = True
+    settings["sharpwave_analysis_settings"]["estimator"]["mean"].append(
+        sw_feature
+    )
+
+
+stream = py_nm.Stream(
+    settings=settings,
+    nm_channels=nm_channels,
+    verbose=True,
+    sfreq=sfreq,
+    line_noise=50,
+)
+
+features = stream.run(
+    data=data[[0, 1, 8, 11], : int(sfreq * 50)],
+    # out_path_root=PATH_OUT,
+    folder_name="meg",
+)
+
+analyzer = nm_analysis.Feature_Reader(
+    feature_dir="C:\\code\\py_neuromodulation",  # stream.PATH_OUT,  #
+    feature_file="meg",  # stream.PATH_OUT_folder_name,  #
+)
+
+analyzer.plot_all_features(ch_used="M1l")
+plt.clim(-3, 3)
+
+nm_plots.plot_corr_matrix(
+    feature=analyzer.feature_arr.filter(regex="M1l"),
+    ch_name="M1l",
+    feature_names=analyzer.feature_arr.filter(regex="M1l").columns,
+    feature_file=analyzer.feature_file,
+    show_plot=True,
+    figsize=(15, 15),
+)
+
+
+plt.plot(features["time"], features["M1r_fft_theta"])
+plt.plot(features["time"], features["M1r_fft_alpha"])
+plt.plot(features["time"], features["M1r_fft_low beta"])