Analyses
For one pilot: the demo's detector rebuilt and run on each look, update by update. Then forced choices on windows of 0.5 to 4 seconds, by the demo's score and by canonical correlation with sinusoids at each rate and its harmonics (Lin and others, 2006), with Wolpaw's bits a minute and half a second added to each choice for moving the eyes. Then each rate's signal-to-noise ratio at O1 and at O2, with the eyes on the center. Change the w1 in the labels for each pilot.
import numpy as np
import mne
from scipy.signal import butter, sosfiltfilt
from scipy.signal.windows import hann
RATES = {"top": 8, "right": 10, "bottom": 12, "left": 15}
NEAR = (-2, -1.5, -1, 1, 1.5, 2)
def load(path):
raw = mne.io.read_raw_bdf(path, preload=True, verbose="error")
return raw.info["sfreq"], 1e6 * raw.get_data(picks="eeg")
def power(x, fs, hz):
"""Power at exactly hz, Hann-windowed, averaged over the channels, as the demo does."""
n = x.shape[1]
xs = (x - x.mean(axis=1, keepdims=True)) * hann(n, sym=False)
k = np.exp(-2j * np.pi * hz * np.arange(n) / fs)
return float(np.mean(np.abs(xs @ k) ** 2)) / n ** 2
def snr(x, fs, hz):
return power(x, fs, hz) / np.mean([power(x, fs, hz + d) for d in NEAR])
def score(x, fs, hz):
return snr(x, fs, hz) + 0.5 * max(0.0, snr(x, fs, 2 * hz) - 1)
def cca(x, fs, hz, harmonics=3):
"""The largest canonical correlation of the channels with sinusoids at hz's harmonics."""
t = np.arange(x.shape[1]) / fs
y = np.array([f(2 * np.pi * k * hz * t) for k in range(1, harmonics + 1)
for f in (np.sin, np.cos)])
qx = np.linalg.qr((x - x.mean(axis=1, keepdims=True)).T)[0]
qy = np.linalg.qr((y - y.mean(axis=1, keepdims=True)).T)[0]
return np.linalg.svd(qx.T @ qy, compute_uv=False)[0]
def windows(x, fs, seconds, step=0.25):
n = int(seconds * fs)
return [x[:, i:i + n] for i in range(0, x.shape[1] - n + 1, int(step * fs))]
def held(adjusted, threshold=2.5, margin=1.25, hold=2, release=2):
"""The demo's detector, update by update: the light it holds after each, or None."""
current = pending = None
count = empty = 0
out = []
for s in adjusted:
best, second = sorted(s.values())[::-1][:2]
top = max(s, key=s.get) if best >= threshold and best >= margin * second else None
if top is None:
pending, count, empty = None, 0, empty + 1
current = None if empty >= release else current
else:
empty = 0
if top == current:
pending, count = None, 0
elif top == pending:
count += 1
if count >= hold:
current, pending, count = top, None, 0
else:
pending, count = top, 1
out.append(current)
return out
def bits_a_minute(p, seconds, n=4):
"""Wolpaw's information transfer rate: n choices, each right with probability p."""
if p <= 1 / n:
return 0.0
bits = np.log2(n) + p * np.log2(p)
if p < 1:
bits += (1 - p) * np.log2((1 - p) / (n - 1))
return bits * 60 / seconds
PILOT = "w1" # run it in the folder the program wrote the files to
fs, steady = load(f"{PILOT}-steady.bdf") # the center of a still screen: the baseline
looks = {d: load(f"{PILOT}-look-{hz}.bdf")[1] for d, hz in RATES.items()}
# 1. The demo, exactly: 2-s windows four times a second, divided by the baseline, then held.
base = {d: max(1.0, np.mean([score(w, fs, hz) for w in windows(steady, fs, 2)]))
for d, hz in RATES.items()}
for d, x in looks.items():
h = held([{k: score(w, fs, hz) / base[k] for k, hz in RATES.items()}
for w in windows(x, fs, 2)])
right, nothing = np.mean([v == d for v in h]), np.mean([v is None for v in h])
print(f"looking {d:6s} held: right {right:.2f}, wrong {1 - right - nothing:.2f}, "
f"nothing {nothing:.2f}")
# 2. Forced choices on windows of each length: the demo's score against canonical correlation,
# which works on the signal band-passed from 6 to 50 Hz, as Lin and others (2006) did.
GAZE = 0.5 # seconds to move the eyes, added to each choice
bp = butter(4, [6, 50], "bandpass", fs=fs, output="sos")
passed = {d: sosfiltfilt(bp, x, axis=1) for d, x in looks.items()}
for seconds in (0.5, 1, 2, 4):
base = {d: max(1.0, np.mean([score(w, fs, hz) for w in windows(steady, fs, seconds)]))
for d, hz in RATES.items()}
right, n = {"score": 0, "cca": 0}, 0
for d, x in looks.items():
for w, v in zip(windows(x, fs, seconds), windows(passed[d], fs, seconds)):
s = {k: score(w, fs, hz) / base[k] for k, hz in RATES.items()}
c = {k: cca(v, fs, hz) for k, hz in RATES.items()}
right["score"] += max(s, key=s.get) == d
right["cca"] += max(c, key=c.get) == d
n += 1
print(f"{seconds} s: " + " ".join(
f"{k} {v / n:.2f} right, {bits_a_minute(v / n, seconds + GAZE):.0f} bits a minute"
for k, v in right.items()))
# 3. Eyes on the center, lights flickering: each rate's response at O1 and at O2.
for site in ("o1", "o2"):
fs, x = load(f"{PILOT}-center-{site}.bdf")
print(site, " ".join(f"{d} {np.mean([snr(w, fs, hz) for w in windows(x, fs, 4, 1)]):.2f}"
for d, hz in RATES.items()))
Then report, for your group's pilots and the class's:
- How often the rebuilt detector holds the right light, a wrong one, or none, for each light.
- Accuracy and bits a minute against the window, for both decoders, and the window where bits a minute peak.
- The 10 and 15 hertz responses at O1 against O2, with the 8 and 12 hertz lights, which sit on the vertical midline, as controls.
- Which light is hardest, and whether the pilot's alpha accounts for it.
Sources
Regan, D. (1966). Some characteristics of average steady-state and transient responses evoked by modulated light. Electroencephalography and Clinical Neurophysiology.
Barrett, G., Blumhardt, L., Halliday, A. M., Halliday, E., and Kriss, A. (1976). A paradox in the lateralisation of the visual evoked response. Nature.
Wolpaw, J. R., Birbaumer, N., McFarland, D. J., Pfurtscheller, G., and Vaughan, T. M. (2002). Brain-computer interfaces for communication and control. Clinical Neurophysiology.
Fisher, R. S., Harding, G., Erba, G., Barkley, G. L., and Wilkins, A. (2005). Photic- and pattern-induced seizures: a review for the Epilepsy Foundation of America Working Group. Epilepsia.
Kelly, S. P., Lalor, E. C., Finucane, C., McDarby, G., and Reilly, R. B. (2005). Visual spatial attention control in an independent brain-computer interface. IEEE Transactions on Biomedical Engineering.
Lin, Z., Zhang, C., Wu, W., and Gao, X. (2006). Frequency recognition based on canonical correlation analysis for SSVEP-based BCIs. IEEE Transactions on Biomedical Engineering.
Yuan, P., and others (2013). A study of the existing problems of estimating the information transfer rate in online brain-computer interfaces. Journal of Neural Engineering.
Chen, X., Wang, Y., Nakanishi, M., Gao, X., Jung, T.-P., and Gao, S. (2015). High-speed spelling with a noninvasive brain-computer interface. Proceedings of the National Academy of Sciences.
Norcia, A. M., Appelbaum, L. G., Ales, J. M., Cottereau, B. R., and Rossion, B. (2015). The steady-state visual evoked potential in vision research: a review. Journal of Vision.
Nakanishi, M., and others (2018). Enhancing detection of SSVEPs for a high-speed brain speller using task-related component analysis. IEEE Transactions on Biomedical Engineering.