From c6e31d9556769575681033b9b94b2addabf1a03b Mon Sep 17 00:00:00 2001 From: sangfrois Date: Tue, 16 May 2023 13:39:18 -0400 Subject: [PATCH] pulse --- pipeline/data_in.py | 9 ++- pipeline/manager.py | 8 ++- pipeline/processors.py | 126 +++++++++++++++++++++++++++++++++++++++-- 3 files changed, 133 insertions(+), 10 deletions(-) diff --git a/pipeline/data_in.py b/pipeline/data_in.py index a837cf0..bbbf01a 100644 --- a/pipeline/data_in.py +++ b/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) diff --git a/pipeline/manager.py b/pipeline/manager.py index 3681e47..eaabc07 100644 --- a/pipeline/manager.py +++ b/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=[ diff --git a/pipeline/processors.py b/pipeline/processors.py index 0ca47ee..474f8d5 100644 --- a/pipeline/processors.py +++ b/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