Load data

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

file = Path('../data/eeg/muse/subject0000/session001/recording_2021-07-15-11.32.37.PPG.csv')
df = pd.read_csv(file)
df['timestamps'] = pd.to_datetime(df['timestamps'],unit='s')
df.set_index('timestamps')

	PPG1	PPG2	PPG3
timestamps
2021-07-15 11:32:51.844000000	75768.0	173634.0	8.0
2021-07-15 11:32:51.860000000	78230.0	320708.0	0.0
2021-07-15 11:32:51.876000000	78175.0	321547.0	0.0
2021-07-15 11:32:51.891000064	78032.0	321479.0	0.0
2021-07-15 11:32:51.907000064	78092.0	321532.0	0.0
...	...	...	...
2021-07-15 12:25:40.984999936	81048.0	325373.0	50.0
2021-07-15 12:25:41.000999936	81059.0	325383.0	7.0
2021-07-15 12:25:41.016000000	81143.0	325321.0	30.0
2021-07-15 12:25:41.032000000	81181.0	325421.0	9.0
2021-07-15 12:25:41.048000000	81151.0	325296.0	0.0

202830 rows × 3 columns

Plot raw data

df.plot('timestamps', 'PPG1')

<AxesSubplot:xlabel='timestamps'>

fs = 1 / np.mean([(t2 - t1).total_seconds() for t1, t2 in zip(df['timestamps'].iloc[:-1], df['timestamps'].iloc[1:])])

from eegclassify.clean import filter
import scipy.signal

fs = 1 / np.mean([(t2 - t1).total_seconds() for t1, t2 in zip(df['timestamps'].iloc[:-1], df['timestamps'].iloc[1:])])
print(f"Sampling freq: {fs}")

def spectral(s: pd.Series):
    ppg = np.array([s.values.tolist()]).T
    ppg_filtered = filter(ppg, n_chans=1, low=30 / 60, high=200 / 60, sfreq=fs, verbose=False)
    f, Pxx = scipy.signal.welch(ppg_filtered.T, fs, nperseg=60*fs)
    f = 60 * f  # Convert to BPM
    return f, Pxx

Sampling freq: 64.00075032990458

f, Pxx = spectral(df['PPG1'].iloc[0:int(60 * fs)])
minidx = int(np.argmax(f > 40))  # plot only certain frequencies
maxidx = int(np.argmax(f > 120))  # plot only certain frequencies

plt.plot(f[minidx:maxidx], Pxx.T[minidx:maxidx])
plt.title("Heartrate")
plt.xlabel("Beats per minute")
plt.ylim(0);

Pulse

Compute the mean pulse

step = int(60 * fs)  # step has to be large enough to capture several heartbeats

pulse_coll = []
s = df['PPG1']
for t in range(0, len(df), step):
    if t + step > len(s):
        continue
    f, Pxx = spectral(s.iloc[t:t+step])
    pulse_coll += [f[np.argmax(Pxx)]]
    
pulse = np.array(pulse_coll)
plt.plot(pulse)

[<matplotlib.lines.Line2D at 0x7fa8b6e9a430>]