
# Imports
from pathlib import Path

import numpy as np
import seaborn as sns
import matplotlib
import matplotlib.pyplot as plt

from eegclassify import main, load, clean, features, preprocess, transform
from eegwatch.devices.muse import CHANNELS_MUSE

%matplotlib inline
plt.rcParams['figure.figsize'] = [16, 4]
plt.rcParams['figure.dpi'] = 300

sfreq = 256
N_SAMPLES_PLOT = 400


%%javascript
document.title='erb-thesis/Signal - Jupyter'  // Set the document title to be able to track time spent working on the notebook with ActivityWatch


# Pick the last recording, or a test file
use_last_recording = False
if use_last_recording:
    recording_dir = Path("/home/erb/.eegnb/data/test/local/museS/subject0000/session001")
    assert recording_dir.exists()
    files = list(recording_dir.glob("*.csv"))
else:
    files = [load.TEST_EEG_FILES_MUSE[0]]

# Load the data
df = load.load_eeg(files)
print("Recording ended: ", df['timestamp'].max())
X = df.drop(columns=['timestamp']).to_numpy()
print(X.shape)

Recording ended:  2020-09-30 09:19:39.743000064+00:00
(15204, 4)


ax = sns.lineplot(data=X[-N_SAMPLES_PLOT:, :], dashes=False)
ax.set_title("Raw signal (pre-filtering)")
ax.legend(labels=CHANNELS_MUSE, loc='upper right');
ax.set_xlim(0, N_SAMPLES_PLOT)
ax.set_ylim(-125, 125)

(-125.0, 125.0)


X_f = clean.filter(X, sfreq)
ax = sns.lineplot(data=X_f[-N_SAMPLES_PLOT:, :], dashes=False)
ax.set_title("Filtered signal")
ax.legend(labels=CHANNELS_MUSE, loc='upper right');
ax.set_xlim(0, N_SAMPLES_PLOT)
ax.set_ylim(-125, 125)

Setting up band-pass filter from 3 - 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 3.00
- Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 2.00 Hz)
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 423 samples (1.652 sec)

(-125.0, 125.0)


# Check the stddev of each channel
# TODO: refactor signal check function to actually use this
X_check = X_f[-N_SAMPLES_PLOT:, :]
print("mean\t", np.mean(X_check, axis=0))
print("std\t", np.std(X_check, axis=0))
print("maxabs\t", np.max(np.abs(X_check), axis=0))

mean	 [ 0.46516693 -0.04659825  0.04369187  0.20437387]
std	 [21.51714701  6.95791439  8.81053102 19.64094595]
maxabs	 [74.13276658 23.78804101 42.90280673 67.95813093]


# std below this value is good
std_thres = 40

# std below this value is perfect
std_perfect = 10

def check(X: np.ndarray) -> float:
    std = np.std(X, axis=0)
    return std
    
def score(std: float) -> float:
    qual = std_thres - np.clip(std, std_perfect, std_thres)
    qual = qual / (std_thres - std_perfect)
    return np.clip(qual, 0, 1)

assert score(std_thres) == 0
assert score(std_perfect) == 1


# Plot signal quality as timeline plot
from eegclassify.plot import TimelineFigure

tl = TimelineFigure(title="Signal quality", figsize=(16, 4))

cmap = matplotlib.cm.get_cmap('RdYlGn')
chunk_len = 250

for i, channel in enumerate(CHANNELS_MUSE):
    events = []
    for si in range(chunk_len, X_f.shape[0], chunk_len):
        begin = si - chunk_len
        end = si
        segment = X_f[list(range(begin, end)), i]
        quality = score(check(segment))
        # Channels should be listed in top-down order, as in the muse-lsl viewer
        chidx = len(CHANNELS_MUSE) - 1 - i
        events.append((begin, end, cmap(quality), ''))
        
    tl.add_bar(events, channel)

tl.plot()

Signal¶