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 [37]:
# Import packages
import neurokit as nk
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib
# 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("hls")) # Better colours
In [38]:
# Download data
df = pd.read_csv("https://raw.githubusercontent.com/neuropsychology/NeuroKit.py/master/examples/Bio/bio_100Hz.csv")
# Plot it
df.plot()
Out[38]:
df contains 2.5 minutes of data recorded at 100Hz (2.5 x 60 x 100 = 15000 data points). There are 4 channels, EDA, ECG, RSP and the Photosensor used to localize events. In the present case, there are four events, corresponding to emotionally negative and neutral pictures presented for 3 seconds.
Biosignals processing can be done quite easily using NeuroKit with the bio_process() function. Simply provide the appropriate biosignal channels and additional channels that you want to keep (for example, the photosensor), and bio_process() will take care of the rest. It will returns a dict containing a dataframe df, including the raw as well as processed signals, and features relevant to each provided signal.
In [39]:
# Process the signals
bio = nk.bio_process(ecg=df["ECG"], rsp=df["RSP"], eda=df["EDA"], add=df["Photosensor"], sampling_rate=100)
# Plot the processed dataframe, normalizing all variables for viewing purpose
nk.z_score(bio["df"]).plot()
Out[39]:
Woah, there's a lot going on there! From 3 variables of interest (ECG, RSP and EDA), bio_process() produced 18 signals. Moreover, the bio dict contains three other dicts, ECG, RSP and EDA containing other features that cannot be simply added in a dataframe. Let's see what we can do with that.
In [40]:
bio["ECG"]["Average_Signal_Quality"] # Get average quality
Out[40]:
As you can see, the average quality of the ECG signal is 99%. See this TO BE DONE tutorial for how to record a good signal.
Let's take a look at each individual heart beat, synchronized by their R peak. You can plot all of them by doing the following:
In [41]:
pd.DataFrame(bio["ECG"]["Cardiac_Cycles"]).plot(legend=False) # Plot all the heart beats
Out[41]:
A large number of HRV indices can be found by checking out bio["ECG"]["HRV"].
One of the most popular RSA algorithm (P2T) implementation can be found in the main data frame.
In [42]:
nk.z_score(bio["df"][["ECG_Filtered", "RSP_Filtered", "RSA"]])[1000:2500].plot()
Out[42]:
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.
This experiment consisted of 4 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 [43]:
condition_list = ["Negative", "Neutral", "Neutral", "Negative"]
First, we must find events onset within our photosensor's signal using the find_events() function. Specify a cut direction (should it select events that are higher or lower than the treshold).
In [44]:
events = nk.find_events(df["Photosensor"], cut="lower")
events
Out[44]:
As we can see, find_events() returns a dict containing onsets and durations for each event. Here, each event lasts for approximately 300 data points (= 3 seconds sampled at 100Hz).
Then, we have to split our dataframe in epochs, i.e. segments of data around the event. We set our epochs to start one second before the event start (onset=-100) and to last for 700 data points, in our case equal to 7 s (since the signal is sampled at 100Hz).
In [45]:
epochs = nk.create_epochs(bio["df"], events["onsets"], duration=700, onset=-100)
Let's plot the first epoch.
In [46]:
nk.z_score(epochs[0][["ECG_Filtered", "EDA_Filtered", "Photosensor"]]).plot()
Out[46]:
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 [47]:
data = {} # Initialize an empty dict
for epoch_index in epochs:
data[epoch_index] = {} # Initialize an empty dict for the current epoch
epoch = epochs[epoch_index]
# ECG
baseline = epoch["ECG_RR_Interval"].ix[-100:0].mean() # Baseline
rr_max = epoch["ECG_RR_Interval"].ix[0:400].max() # Maximum RR interval
data[epoch_index]["HRV_MaxRR"] = rr_max - baseline # Corrected for baseline
# EDA - SCR
scr_max = epoch["SCR_Peaks"].ix[0:600].max() # Maximum SCR peak
if np.isnan(scr_max):
scr_max = 0 # If no SCR, consider the magnitude, i.e. that the value is 0
data[epoch_index]["SCR_Magnitude"] = scr_max
data = pd.DataFrame.from_dict(data, orient="index") # Convert to a dataframe
data["Condition"] = condition_list # Add the conditions
data # Print
Out[47]:
In [48]:
sns.boxplot(x="Condition", y="HRV_MaxRR", data=data)
Out[48]:
In accord with the litterature, we observe that the RR interval is higher in the negative than in the neutral condition.
In [49]:
sns.boxplot(x="Condition", y="SCR_Magnitude", data=data)
Out[49]:
In the same line, the skin conductance response (SCR) is higher in the negative condition compared to the neutral one. Overall, these results suggest that the acquired biosignals are sensitive to the cognitive processing of emotional stimuli.