Analysis

What filters cost, and what a detector finds

First, put a 20 µV shift lasting one second through a high-pass at 0.5 hertz, causal and zero-phase, and a 50 µV step through the 58 to 62 hertz notch. The device's high-pass order isn't stated, so a second-order Butterworth stands in for it.

import numpy as np
from scipy.signal import butter, iirnotch, lfilter, filtfilt, sosfilt, sosfiltfilt

fs = 500
t = np.arange(0, 6, 1 / fs)
slow = 20.0 * ((t >= 2) & (t < 3))            # a 20 uV shift lasting one second
blip = np.zeros_like(t); blip[t >= 4.5] = 50.0  # a 50 uV step, as when an electrode shifts

# A stand-in for the device's 0.5 Hz high-pass: second-order Butterworth.
hp = butter(2, 0.5, "highpass", fs=fs, output="sos")
causal, zero_phase = sosfilt(hp, slow), sosfiltfilt(hp, slow)
at = lambda x, s: x[int(s * fs)]
for when, s in (("0.5 s before it begins", 1.5), ("0.9 s into the shift", 2.9),
                ("just after it ends", 3.02), ("0.5 s after it ends", 3.5)):
    print(f"{when:24s} causal {at(causal, s):5.1f} uV   zero-phase {at(zero_phase, s):5.1f} uV")

# The 58 to 62 Hz band, nulled at 60 Hz, run as two sections: how long a step makes it ring.
b, a = iirnotch(60, 60 / 4, fs=fs)
rung = lfilter(b, a, lfilter(b, a, blip)) - blip
tail = np.abs(rung[t >= 4.5])
last = np.nonzero(tail > 0.01 * 50)[0].max() / fs
print(f"the notch rings for {1000 * last:.0f} ms after the step, above 1% of its size")

You should seeThe causal filter undoes the shift within a second, overshooting to −3.6 µV, and when it ends drops to −21.1 µV, a deflection that was never in the signal. The zero-phase filter keeps more of it, 7.3 µV, and shows −1.8 µV half a second before it begins. The notch rings for about 100 milliseconds after the step.

Then flag the Alpha Flower forehead recordings in 1-second windows, for blinks by their peak-to-peak size and for muscle by their power from 30 to 45 hertz, and check the flags against the blinks the analyst wrote down. Change the thresholds and watch both kinds of error.

import numpy as np
import mne
from scipy.signal import welch

def flags(path, ptp_uv=100, z=5):
    """Each 1-s window flagged for a blink (peak to peak) or for muscle (30 to 45 Hz power)."""
    raw = mne.io.read_raw_bdf(path, preload=True, verbose="error")
    fs = int(raw.info["sfreq"])
    x = 1e6 * raw.get_data(picks="eeg")
    w = x[:, : x.shape[1] // fs * fs].reshape(x.shape[0], -1, fs)   # channels, windows, samples
    ptp = np.ptp(w, axis=2).max(axis=0)
    f, p = welch(w, fs=fs, nperseg=fs // 2, axis=2)
    high = p[:, :, (f >= 30) & (f < 45)].mean(axis=2).max(axis=0)
    mad = np.median(np.abs(high - np.median(high)))
    return ptp > ptp_uv, high > np.median(high) + z * 1.4826 * mad, w, fs

def compare(flagged, marks, slack=1):
    """Marks with a flagged window within `slack` s, marks missed, flagged windows unmarked."""
    hit = [m for m in marks if flagged[max(0, int(m) - slack):int(m) + slack + 1].any()]
    near = np.zeros(len(flagged), bool)
    for m in marks:
        near[max(0, int(m) - slack):int(m) + slack + 1] = True
    return len(hit), len(marks) - len(hit), int((flagged & ~near).sum())

blink, muscle, w, fs = flags("w1-fpz.bdf")      # the path the program printed
marks = [3.5, 9, 15.5, 21, 27, 33.5, 40]        # the analyst's blink times from the sheet
found, missed, unmarked = compare(blink, marks)
print(f"{blink.sum()} windows flagged for blinks, {muscle.sum()} for muscle, of {len(blink)}")
print(f"marks found {found}, missed {missed}; flagged windows with no mark near {unmarked}")

def share(seg):
    f, p = welch(seg, fs=fs, nperseg=2 * fs)
    band = lambda lo, hi: p[:, (f >= lo) & (f < hi)].sum(axis=1)
    return float(np.mean(band(8, 13) / band(4, 30)))

keep = ~(blink | muscle)
for name, block in (("open", slice(3, 45)), ("closed", slice(48, 90))):
    every = np.concatenate(list(w[:, block].transpose(1, 0, 2)), axis=1)
    kept = np.concatenate([w[:, i] for i in range(len(keep))[block] if keep[i]], axis=1)
    print(f"{name}: share {share(every):.3f} with every window, "
          f"{share(kept):.3f} without the flagged")