Adding a processor for Pulse signal using serial

pipeline/data_in.py

@@ -155,9 +155,10 @@ class EEGRecording(EEGStream):
 
     Parameters:
         raw (str, BaseRaw): file-name of a raw EEG file or an instance of mne.io.BaseRaw
+        buffer_seconds (int): the number of seconds to buffer incoming data
     """
 
-    def __init__(self, raw: Union[str, BaseRaw]):
+    def __init__(self, raw: Union[str, BaseRaw], buffer_seconds: int = 5):
         # load raw EEG data
         if not isinstance(raw, BaseRaw):
             raw = read_raw(raw)
@@ -169,12 +170,13 @@ def __init__(self, raw: Union[str, BaseRaw]):
         self.mock_stream.start()
 
         # start the LSL client
-        super(EEGRecording, self).__init__(host=host)
+        super(EEGRecording, self).__init__(host=host, buffer_seconds=buffer_seconds)
 
     @staticmethod
     def make_eegbci(
         subjects: Union[int, List[int]] = 1,
         runs: Union[int, List[int]] = [1, 2],
+        buffer_seconds: int = 5,
     ):
         """
         Static utility function to instantiate an EEGRecording instance using
@@ -186,7 +188,8 @@ def make_eegbci(
         Parameters:
             subjects (int, List[int]): which subject(s) to load data from
             runs (int, List[int]): which run(s) to load from the corresponding subject
+            buffer_seconds (int): the number of seconds to buffer incoming data
         """
         raw = concatenate_raws([read_raw(p) for p in eegbci.load_data(subjects, runs)])
         eegbci.standardize(raw)
-        return EEGRecording(raw)
+        return EEGRecording(raw, buffer_seconds=buffer_seconds)

pipeline/manager.py

@@ -106,6 +106,9 @@ def run(self):
     mngr = Manager(
         data_in={
             "file": data_in.EEGRecording.make_eegbci(),
+            "plant": data_in.SerialStream(sfreq=100, buffer_seconds=5),
+            "pulse": data_in.SerialStream(sfreq=100, buffer_seconds=60),
+            "muse": data_in.EEGStream(host=""),
         },
         processors=[
             processors.PSD(label="delta"),
@@ -115,8 +118,9 @@ def run(self):
             processors.PSD(label="gamma"),
             processors.LempelZiv(),
             processors.Ratio("/file/alpha", "/file/theta", "alpha/theta"),
-            processors.Biocolor(channels={"file": ["C3"]}),
-            processors.Biotuner(channels={"file": ["C3", "C4", "O1", "O2"]}),
+            processors.Biocolor(channels={"muse": ["AF7"]}),
+            processors.Biotuner(channels={"plant": ["serial"]}),
+            processors.Pulse(channels={"pulse": ["serial"]})
         ],
         normalization=normalization.StaticBaselineNormal(duration=30),
         data_out=[

pipeline/processors.py

@@ -3,7 +3,8 @@
 import threading
 import time
 from typing import Callable, Dict, List, Optional, Tuple
-
+import neurokit2 as nk
+import pandas as pd
 import mne
 import numpy as np
 from antropy import lziv_complexity, spectral_entropy
@@ -501,9 +502,14 @@ def process(
         result = {}
         for i, hsvs in enumerate(latest_hsvs):
             for j, hsv in enumerate(hsvs):
-                result[f"{self.label}/ch{i}_peak{j}_hue"] = hsv[0]
-                result[f"{self.label}/ch{i}_peak{j}_sat"] = hsv[1]
-                result[f"{self.label}/ch{i}_peak{j}_val"] = hsv[2]
+                if info["nchan"] > 1:
+                    result[f"{self.label}/ch{i}_peak{j}_hue"] = hsv[0]
+                    result[f"{self.label}/ch{i}_peak{j}_sat"] = hsv[1]
+                    result[f"{self.label}/ch{i}_peak{j}_val"] = hsv[2]
+                else:
+                    result[f"{self.label}/peak{j}_hue"] = hsv[0]
+                    result[f"{self.label}/peak{j}_sat"] = hsv[1]
+                    result[f"{self.label}/peak{j}_val"] = hsv[2]
         return result
 
 
@@ -637,7 +643,7 @@ def process(
         result = {}
         normalization_mask = {}
         for i in range(len(metrics)):
-            ch_prefix = f"ch{i}_"
+            ch_prefix = f"ch{i}_" if info["nchan"] > 1 else ""
             result[f"{self.label}/{ch_prefix}harmsim"] = metrics[i]["harmsim"]
             normalization_mask[f"{self.label}/{ch_prefix}harmsim"] = True
             result[f"{self.label}/{ch_prefix}cons"] = metrics[i]["cons"]
@@ -667,3 +673,113 @@ def process(
                 result[f"{self.label}/harm_conn/{i}/{j}"] = harm_conn[i][j]
                 normalization_mask[f"{self.label}/harm_conn/{i}/{j}"] = False
         return result
+
+
+class Pulse(Processor):
+    """
+    Feature extractor for Pulse metrics (HRV).
+
+    Parameters:
+        label (str): label under which to save the extracted feature
+        channels (Dict[str, List[str]]): channel list for each input stream
+        extraction_frequency (float, optional): the frequency in Hz at which to run the peak extraction loop
+    """
+
+    def __init__(
+        self,
+        label: str = "pulse",
+        channels: Dict[str, List[str]] = None,
+        extraction_frequency: float = 1 / 5,
+    ):
+        super(Pulse, self).__init__(label, channels)
+        self.sfreq = None
+        self.latest_raw = None
+        self.latest_hrv = None
+        self.raw_lock = threading.Lock()
+        self.features_lock = threading.Lock()
+        self.extraction_frequency = extraction_frequency
+        self.extraction_thread = threading.Thread(
+            target=self.extraction_loop, daemon=True
+        )
+        self.extraction_thread.start()
+
+    def extraction_loop(self):
+        """
+        This function runs the biotuner_realtime function in a loop in a separate thread.
+        It continuously grabs the latest raw data, processes it, and updates the latest_hsvs.
+        """
+        import warnings
+        warnings.filterwarnings("ignore")
+        while True:
+            with self.raw_lock:
+                raw = self.latest_raw
+
+            if raw is None:
+                time.sleep(0.05)
+                continue
+            try:
+                if raw.shape[0] > 1:
+                    print("got more than one channel")
+                ppg, info_ppg = nk.ppg_process(raw[0], sampling_rate=self.sfreq)
+                hrv_df = nk.hrv(info_ppg, sampling_rate=self.sfreq)
+            except:
+                print("neurokit failed.")
+                continue
+
+            with self.features_lock:
+                self.latest_ppg = ppg
+                self.latest_info = info_ppg
+                self.latest_hrv = hrv_df
+
+            if self.extraction_frequency is not None:
+                sleep_time = 1 / self.extraction_frequency
+                time.sleep(sleep_time)
+
+    def process(
+        self,
+        raw: np.ndarray,
+        info: mne.Info,
+        processed: Dict[str, float],
+        intermediates: Dict[str, np.ndarray],
+    ):
+        """
+        This function computes the HRV from the pulse signal.
+
+        Parameters:
+            raw (np.ndarray): the raw PPG buffer with shape (Channels, Time)
+            info (mne.Info): info object containing e.g. channel names, sampling frequency, etc.
+            processed (Dict[str, float]): dictionary collecting extracted features
+            intermediates (Dict[str, np.ndarray]): dictionary containing intermediate representations
+        """
+        self.sfreq = info["sfreq"]
+
+        with self.raw_lock:
+            self.latest_raw = raw
+            # self.latest_info = info_ppg
+            # self.latest_ppg = ppg
+
+        with self.features_lock:
+            hrv_df = self.latest_hrv
+        if hrv_df is None:
+            hrv_df = pd.DataFrame({"HRV_SDNN": [0], "HRV_RMSSD": [0], "HRV_MeanNN": [0], "HRV_pNN50": [0],
+                                   "HRV_LF": [0], "HRV_HF": [0], "HRV_LFHF": [0],
+                                   "HRV_SD1": [0], "HRV_SD2": [0], "HRV_SD1SD2": [0], "HRV_ApEn":[0], "HRV_SampEn":[0],
+                                   "HRV_DFA_alpha1": [0], "HRV_DFA_alpha2": [0], })
+
+        result = {}
+        result[f"{self.label}/sdnn"] = hrv_df["HRV_SDNN"].values[0]
+        result[f"{self.label}/rmssd"] = hrv_df["HRV_RMSSD"].values[0]
+        result[f"{self.label}/meanHR"] = hrv_df["HRV_MeanNN"].values[0]
+        result[f"{self.label}/pnn50"] = hrv_df["HRV_pNN50"].values[0]
+        result[f"{self.label}/lf"] = hrv_df["HRV_LF"].values[0]
+        result[f"{self.label}/hf"] = hrv_df["HRV_HF"].values[0]
+        result[f"{self.label}/lfhf"] = hrv_df["HRV_LFHF"].values[0]
+        result[f"{self.label}/sd1"] = hrv_df["HRV_SD1"].values[0]
+        result[f"{self.label}/sd2"] = hrv_df["HRV_SD2"].values[0]
+        result[f"{self.label}/sd1sd2"] = hrv_df["HRV_SD1SD2"].values[0]
+        result[f"{self.label}/apen"] = hrv_df["HRV_ApEn"].values[0]
+        result[f"{self.label}/sampen"] = hrv_df["HRV_SampEn"].values[0]
+        result[f"{self.label}/dfa1"] = hrv_df["HRV_DFA_alpha1"].values[0]
+        result[f"{self.label}/dfa2"] = hrv_df["HRV_DFA_alpha2"].values[0]
+
+        return result