Episode 6: Averaging an Evoked Response

Measure how much averaging removes, with a Trigger, a Baseline and an EpochAverage on a signal that is nothing but noise.

An evoked response is small, a few microvolts against tens of microvolts of ongoing EEG, several times smaller than what it sits in. That is why nobody looks at a single trial: you average many, and the response emerges because it is the same every time while everything else is not. That is a claim you can measure, and this episode measures it. The generator produces only noise, with no evoked response in it at all, so whatever averaging removes is exactly the part that should shrink.

Noise, and markers from the script

A Generator with no signal, only noise, and a Marker fired by the script, standing in for a paradigm.

26    p = gp.Pipeline()
27
28    # Ongoing "EEG": noise and nothing else
29    source = gp.Generator(sampling_rate=SAMPLING_RATE, channel_count=1,
30                          signal_amplitude=0.0, noise_amplitude=10.0)
31
32    # Markers fired by this script, standing in for a paradigm
33    marker = gp.Marker()
Trigger, Baseline, EpochAverage

The three nodes an ERP pipeline always uses. The Trigger cuts an epoch around each marker. The Baseline subtracts each epoch’s pre-stimulus level, so every trial starts from a common zero; without it a slow drift between trials survives averaging as an offset. The EpochAverage averages the epochs as they arrive, cumulatively here, and says how many went in.

35    # An epoch around each marker, its pre-stimulus level subtracted,
36    # then averaged as they arrive
37    trigger = gp.Trigger(time_pre=TIME_PRE, time_post=TIME_POST, target=1)
38    baseline = gp.Baseline()
39    average = gp.EpochAverage(mode=gp.EpochAverage.CUMULATIVE)
Collect three things

An epoch arrives when a trigger fires, not on the sampling clock, and a Collector fed epochs stacks them into (time, channel, trial): one trial per epoch for the single epochs, and one trial per update for the running average. The count comes from the EpochAverage’s second output.

41    # Epochs stack into (time, channel, trial): one trial per epoch, and
42    # for the running average one trial per update
43    single = gp.Collector(name="single_epochs")
44    averaged = gp.Collector(name="average")
45    counted = gp.Collector(name="count")
46
47    p.connect(source, trigger)
48    p.connect(marker, trigger[gp.Trigger.PORT_TRIGGER])
49    p.connect(trigger, baseline)
50    p.connect(baseline, single)
51    p.connect(baseline, average)
52    p.connect(average, averaged)
53    p.connect(average[gp.EpochAverage.PORT_COUNT], counted)
Fire the markers

The script waits before the first marker: a marker arriving sooner than time_pre after the start has no pre-stimulus data to cut. Then one marker per epoch length, so no two epochs share a sample.

55    p.start()
56    # Let the buffer fill before the first marker, so the first epoch
57    # has its pre-stimulus samples
58    time.sleep(TIME_PRE + 0.3)
59    for _ in range(EPOCHS):
60        marker.emit(1)
61        time.sleep(INTERVAL)
62    time.sleep(TIME_POST + 0.3)
63    p.stop()
64    p.close()
Measure

The final average is the last trial of the running average, data[..., -1], not data[-1], which is the last sample of every update. Its root mean square against that of the single epochs is the ratio to compare with the square root of the count.

66    epochs = single.result
67    mean = averaged.result
68    counts = counted.result
69
70    n_used = int(np.max(counts.data)) if counts is not None else 0
71    single_rms = rms(epochs.data) if epochs is not None else float("nan")
72    # The last trial of the running average is the final average
73    average_rms = (rms(mean.data[..., -1]) if mean is not None
74                   else float("nan"))
75
76    print(f"{EPOCHS} markers fired, {n_used} epochs averaged")
77    print()
78    print(f"  rms across single epochs   {single_rms:8.3f}")
79    print(f"  rms of the final average   {average_rms:8.3f}")
80    if n_used > 0 and average_rms > 0:
81        print(f"  ratio                      {single_rms / average_rms:8.2f}"
82              f"   (sqrt({n_used}) = {np.sqrt(n_used):.2f})")
Result

Run the script. It takes about twenty seconds:

40 markers fired, 40 epochs averaged

  rms across single epochs   ...
  rms of the final average   ...
  ratio                      ...   (sqrt(40) = 6.32)

The ratio lands near the square root of the number of epochs, which is the whole reason the technique works. Read the other way round it becomes a planning tool: averaging buys signal-to-noise as the square root of N, for noise that is independent from trial to trial, which the generator’s is exactly and real EEG only roughly, so going from 40 trials to 160 is not four times better, it is two. That is the arithmetic behind every protocol that asks a subject to sit through hundreds of repetitions.

Three things to keep in mind

Keep the epochs apart. Fired every 0.1 s instead of every 0.5 s, the epochs overlap and share samples. On this generator’s white noise that costs little, since each sample is still independent of the next; on real EEG it does, because the ongoing activity is correlated from one moment to the next, so overlapping trials are no longer independent and the ratio stops following the square root of N.

Record the count. An average without its N is not a result anyone can interpret, and the count is also what tells you how many epochs the rejection thresholds, reject_level and reject_peak_to_peak, threw away.

The Baseline belongs between the Trigger and the average, and this episode puts it there; an ERP measured against a moving baseline is not measured against anything.

Summary

In this episode, you learned how to:

  • build the ERP chain: Trigger, Baseline, EpochAverage

  • collect epochs, a running average and its count as trials

  • measure what averaging removes, and plan a trial count from it

Proceed to Episode 7, where a decision, not a signal, is shown as a number.

Nodes in this episode: Generator, Marker, Trigger, Baseline, EpochAverage, Collector

File s5e6_averaging_an_evoked_response.py – View file on GitHub

 1"""Season 5, Episode 6: how much averaging removes, measured on noise."""
 2
 3import time
 4
 5import numpy as np
 6
 7import gpype as gp
 8
 9SAMPLING_RATE = 250
10EPOCHS = 40
11TIME_PRE = 0.2  # seconds of epoch before the marker
12TIME_POST = 0.3  # and after it
13# One epoch's length between markers, so no two epochs share a sample:
14# overlapping epochs are not independent trials
15INTERVAL = TIME_PRE + TIME_POST
16
17
18def rms(block) -> float:
19    """Root mean square over everything in the block."""
20    array = np.asarray(block, dtype=float)
21    return float(np.sqrt(np.mean(np.square(array))))
22
23
24if __name__ == "__main__":
25
26    p = gp.Pipeline()
27
28    # Ongoing "EEG": noise and nothing else
29    source = gp.Generator(sampling_rate=SAMPLING_RATE, channel_count=1,
30                          signal_amplitude=0.0, noise_amplitude=10.0)
31
32    # Markers fired by this script, standing in for a paradigm
33    marker = gp.Marker()
34
35    # An epoch around each marker, its pre-stimulus level subtracted,
36    # then averaged as they arrive
37    trigger = gp.Trigger(time_pre=TIME_PRE, time_post=TIME_POST, target=1)
38    baseline = gp.Baseline()
39    average = gp.EpochAverage(mode=gp.EpochAverage.CUMULATIVE)
40
41    # Epochs stack into (time, channel, trial): one trial per epoch, and
42    # for the running average one trial per update
43    single = gp.Collector(name="single_epochs")
44    averaged = gp.Collector(name="average")
45    counted = gp.Collector(name="count")
46
47    p.connect(source, trigger)
48    p.connect(marker, trigger[gp.Trigger.PORT_TRIGGER])
49    p.connect(trigger, baseline)
50    p.connect(baseline, single)
51    p.connect(baseline, average)
52    p.connect(average, averaged)
53    p.connect(average[gp.EpochAverage.PORT_COUNT], counted)
54
55    p.start()
56    # Let the buffer fill before the first marker, so the first epoch
57    # has its pre-stimulus samples
58    time.sleep(TIME_PRE + 0.3)
59    for _ in range(EPOCHS):
60        marker.emit(1)
61        time.sleep(INTERVAL)
62    time.sleep(TIME_POST + 0.3)
63    p.stop()
64    p.close()
65
66    epochs = single.result
67    mean = averaged.result
68    counts = counted.result
69
70    n_used = int(np.max(counts.data)) if counts is not None else 0
71    single_rms = rms(epochs.data) if epochs is not None else float("nan")
72    # The last trial of the running average is the final average
73    average_rms = (rms(mean.data[..., -1]) if mean is not None
74                   else float("nan"))
75
76    print(f"{EPOCHS} markers fired, {n_used} epochs averaged")
77    print()
78    print(f"  rms across single epochs   {single_rms:8.3f}")
79    print(f"  rms of the final average   {average_rms:8.3f}")
80    if n_used > 0 and average_rms > 0:
81        print(f"  ratio                      {single_rms / average_rms:8.2f}"
82              f"   (sqrt({n_used}) = {np.sqrt(n_used):.2f})")