Welcome to the course for biosignals processing using NeuroKit and python. You'll find the necessary files to run this example in the examples section.
In [1]:
# Import packages
import neurokit as nk
import pandas as pd
import numpy as np
import matplotlib
import seaborn as sns
# Plotting preferences
%matplotlib inline
matplotlib.rcParams['figure.figsize'] = [14.0, 10.0] # Bigger figures
sns.set_style("whitegrid") # White background
sns.set_palette(sns.color_palette("colorblind")) # Better colours
In [2]:
# Download resting-state data
df = pd.read_csv("https://raw.githubusercontent.com/neuropsychology/NeuroKit.py/master/examples/Bio/data/bio_rest.csv", index_col=0)
# Plot it
df.plot()
Out[2]:
df contains about 5 minutes of data recorded at 1000Hz. There are 4 channels, EDA, ECG, RSP and the Photosensor used to localize events. In the present case, there is only one event, one sequence of 5 min during which the participant was instructed to to nothing.
First thing that we're gonna do is crop that data according to the photosensor channel, to keep only the sequence of interest.
In [3]:
# We want to find events on the Photosensor channel, when it goes down (hence, cut is set to lower).
events = nk.find_events(df["Photosensor"], cut="lower")
print(events)
find_events returns a dict containing onsets and durations of each event. Here, it correctly detected only one event. Then, we're gonna crop our data according to that event. The create_epochs function returns a list containing epochs of data corresponding to each event. As we have only one event, we're gonna select the 0th element of that list.
In [4]:
df = nk.create_epochs(df, events["onsets"], duration=events["durations"], onset=0)
df = df[0] # Select the first (0th) element of that list.
Biosignals processing can be done quite easily using NeuroKit with the bio_process function. Simply provide the biosignal channels and additional channels that you want to keep (for example, the photosensor). bio_process returns a dict containing a dataframe df, including raw and processed signals, as well as features relevant to each provided signal.
In [5]:
bio = nk.bio_process(ecg=df["ECG"], rsp=df["RSP"], eda=df["EDA"], add=df["Photosensor"])
# Plot the processed dataframe
bio["df"].plot()
Out[5]:
Aside from this dataframe, bio contains also several features computed signal wise.
Many indices of HRV, a finely tuned measure of heart-brain communication, are computed.
In [6]:
bio["ECG"]["HRV"]
Out[6]:
TO BE DONE.
TO BE DONE.
The processing functions automatically extracts each individual heartbeat, synchronized by their R peak. You can plot all of them.
In [7]:
bio["ECG"]["Heart_Beats"]
Out[7]:
In [8]:
pd.DataFrame(bio["ECG"]["Heart_Beats"]).T.plot(legend=False) # Plot all the heart beats
Out[8]:
In [7]:
# Print all the HRV indices
bio["ECG_Features"]["ECG_HRV"]
Out[7]:
This experiment consisted of 8 events (when the photosensor signal goes down), which were 2 types of images that were shown to the participant: "Negative" vs "Neutral". The following list is the condition order.
In [12]:
condition_list = ["Negative", "Negative", "Neutral", "Neutral", "Neutral", "Negative", "Negative", "Neutral"]
First, we must find events onset within our photosensor's signal using the find_events() function. This function requires a treshold and a cut direction (should it select events that are higher or lower than the treshold).
In [8]:
events = nk.find_events(df["Photosensor"], treshold = 3, cut="lower")
events
Out[8]:
Then, we divise our dataframe in epochs, i.e. segments of data around the event. We set our epochs to start at the event start (onset=0) and to last for 5000 data points, in our case equal to 5 s (since the signal is sampled at 1000Hz).
In [10]:
epochs = nk.create_epochs(bio["Bio"], events["onsets"], duration=5000, onset=0)
We can then itereate through the epochs and store the interesting results in a new dict that will be, at the end, converted to a dataframe.
In [13]:
evoked = {} # Initialize an empty dict
for epoch in epochs:
evoked[epoch] = {} # Initialize an empty dict for the current epoch
evoked[epoch]["Heart_Rate"] = epochs[epoch]["Heart_Rate"].mean() # Heart Rate mean
evoked[epoch]["RSP_Rate"] = epochs[epoch]["RSP_Rate"].mean() # Respiration Rate mean
evoked[epoch]["EDA_Filtered"] = epochs[epoch]["EDA_Filtered"].mean() # EDA mean
evoked[epoch]["EDA_Max"] = max(epochs[epoch]["EDA_Filtered"]) # Max EDA value
# SRC_Peaks are scored np.nan (NaN values) in the absence of peak. We want to change it to 0
if np.isnan(epochs[epoch]["SCR_Peaks"].mean()):
evoked[epoch]["SCR_Peaks"] = 0
else:
evoked[epoch]["SCR_Peaks"] = epochs[epoch]["SCR_Peaks"].mean()
evoked = pd.DataFrame.from_dict(evoked, orient="index") # Convert to a dataframe
evoked["Condition"] = condition_list # Add the conditions
evoked # Print
Out[13]:
In [14]:
sns.boxplot(x="Condition", y="Heart_Rate", data=evoked)
Out[14]:
In [15]:
sns.boxplot(x="Condition", y="RSP_Rate", data=evoked)
Out[15]:
In [16]:
sns.boxplot(x="Condition", y="EDA_Filtered", data=evoked)
Out[16]:
In [17]:
sns.boxplot(x="Condition", y="EDA_Max", data=evoked)
Out[17]:
In [18]:
sns.boxplot(x="Condition", y="SCR_Peaks", data=evoked)
Out[18]:
In [ ]: