Chapter 2 · Middle school
What two numbers describe every rhythm, from a shaken rope to your heartbeat to your brain?
A wave goes up and down again and again. Two numbers describe it.
Its frequency is how many times it repeats each second. It is measured in hertz: 10 hertz is ten waves every second.
Its amplitude is how big it is, from the middle to the top.
Slide the frequency, and more waves fit into the same second. Slide the amplitude, and each wave grows or shrinks, but repeats just as often.
Try it
Two people hold a long rope or a spring across the floor. One shakes it steadily while a third counts the waves for ten seconds. Shake it slowly, then faster, then bigger. Divide each count by ten for the frequency, in hertz.
| Shake | Waves in 10 seconds | Hertz |
|---|---|---|
| Slow | ||
| Faster | ||
| Bigger |
You should seeShaking faster raises the count, and the frequency with it. Shaking bigger makes taller waves, while the count stays about the same.
Your body is full of rhythms. Resting, your heart beats about once a second: about 1 hertz. You breathe more slowly, about 15 times a minute, which is a quarter of a hertz.
Your brain has rhythms too, and much faster ones. The rhythm you will measure in Alpha Flower is about 10 hertz, ten times faster than your heartbeat.
Try it
Find your pulse on the inside of your wrist with two fingers. Count the beats for 15 seconds, and multiply by four for beats a minute. Then count your breaths for a whole minute. Divide each by 60 for the frequency, in hertz.
| Rhythm | Count | A minute | Hertz |
|---|---|---|---|
| Heartbeat, 15 seconds | |||
| Breathing, 60 seconds |
You should seeA heartbeat of about 1 hertz and breathing well under 1 hertz, both far slower than alpha's 10.
A wave repeats 30 times in 3 seconds. What is its frequency?
Your heart beat 18 times in 15 seconds. How many beats is that a second?
Shaking the rope faster and shaking it bigger are different changes. Which one changes the frequency?
Chapter 2 · High school
What does it take to describe a rhythm, to measure it without fooling yourself, and to take a messy signal apart into the rhythms in it?
A sine wave is described by three numbers. Its frequency is how many times it repeats each second, in hertz. Its amplitude is how big it is, from the middle to the top. Its phase is where in its cycle it is at a given moment.
Two waves of the same frequency and amplitude can still cancel each other, if their phases are half a cycle apart. That is why neurons in step add up, and neurons out of step don't.
Try it
Two people hold a long rope or a spring across the floor. One shakes it steadily while a third counts the waves for ten seconds. Shake it slowly, then faster, then bigger.
| Shake | Waves in 10 seconds | Frequency, hertz | What changed |
|---|---|---|---|
| Slow | |||
| Faster | |||
| Bigger |
You should seeShaking faster raises the frequency. Shaking bigger raises the amplitude, and leaves the frequency about the same.
A digital device measures a signal at regular instants. How many times a second is its sample rate. To capture a frequency, it must sample more than twice in each cycle.
Sample a fast wave too slowly, and the samples trace a slower wave that isn't there. That is aliasing, the effect that makes a wheel in a film seem to spin backward.
The IntoMind One samples 500 times a second, so in principle it can capture rhythms below 250 hertz. Its filters remove most of what lies above 200 hertz.
Try it
On graph paper, draw a 10 hertz sine wave over one second. Mark its value at 50 evenly spaced instants, then at 12, then at 8. Join each set of marks with a smooth curve, and count the curve's waves in the second.
PredictSampled 8 times a second, the 10 hertz wave will look like hertz.
| Samples a second | The curve looks like, hertz | More than twice the frequency? |
|---|---|---|
| 50 | ||
| 12 | ||
| 8 |
You should seeAt 50 a second, the real wave. At 12 and at 8, a slow wave that isn't there.
Any signal, however messy, can be built by adding sine waves of different frequencies, amplitudes and phases. A spectrum shows how much of each frequency a signal contains.
In EEG from the back of the head with the eyes closed, the alpha rhythm stands up as a peak near 10 hertz, which you will measure in Alpha Flower.
Try it
In a spreadsheet, put the times 0 to 1 second down column A, in steps of 0.002. In column B, a 10 hertz wave: =SIN(2*PI()*10*A2). In column C, a 25 hertz wave half as big: =0.5*SIN(2*PI()*25*A2). Add them in column D, and plot it. Then sketch its spectrum by hand: which frequencies, and how much of each. Last, add noise to every sample with +0.5*(RAND()-0.5), and sketch it again.
You should seeA sum that looks messy, and a spectrum with just two peaks: 10 hertz, and 25 hertz half as tall. Noise adds a low floor at every frequency.
After the BCI Flight Simulator.
When you look at a light that flickers at a steady speed, the back of your brain follows it. In the spectrum, that shows as a peak at the light's speed, and smaller ones at two and three times it, its harmonics.
Try it
Before you look at the figure, sketch the spectrum you would expect from the back of the head while looking at the 12 hertz light, with the eyes relaxed. Mark the light's peak, its harmonics, and alpha. Then: where would the 10 hertz light's peak be hard to see, and why?
Shaking the rope faster and shaking it bigger are different changes. Which one is frequency, and which is amplitude?
The IntoMind One samples 500 times a second. What is the fastest rhythm it can capture, and why?
Looking at the 12 hertz light, where in the spectrum would you expect peaks? Why is the 10 hertz light the hardest to see?
Herrmann, C. S. (2001). Human EEG responses to 1–100 Hz flicker: resonance phenomena in visual cortex and their potential correlation to cognitive phenomena. Experimental Brain Research.
Chapter 2 · College
How do you take a sampled signal apart into its rhythms, and what does each choice in the estimate, the window, the segment and the band, do to the number you report?
A sinusoid is x(t) = A sin(2π f t + φ). Its frequency f is how many cycles it completes each second, its amplitude A how far it swings from its middle, and its phase φ where in its cycle it is at time zero.
EEG is not one sinusoid, but any stretch of it can be written as a sum of them. A spectrum says how much of each frequency the sum contains.
A device measures the voltage fs times a second, its sample rate. Only frequencies below fs/2, the Nyquist frequency, are represented faithfully. A frequency f above it folds back to |f − k fs| for the nearest whole k, a false slower wave: aliasing.
So a device filters before it samples. The demos run the IntoMind One at 500 samples a second, and as shipped it filters on the device: a high-pass at 0.5 hertz, notches at the mains frequencies and their harmonics, and a low-pass at 200 hertz, safely below the Nyquist frequency of 250.
The discrete Fourier transform of N samples gives the amplitude and phase of N/2 + 1 sinusoids at the frequencies k fs/N, from 0 to fs/2. Their spacing, the frequency resolution, is fs/N = 1/T, one over the duration: two seconds of signal resolve 0.5 hertz, whatever the sample rate.
The power at each frequency, divided by the spacing, estimates the power spectral density, in squared microvolts per hertz. Its area between two frequencies is the band power.
Cutting out N samples treats the signal as if it stopped sharply at both ends. A frequency that falls between two bins then spreads its power over many: leakage. A window tapers the segment to zero at its ends, which confines the spread at the cost of a wider peak. The Hann window, ½(1 − cos(2πn/N)), is the common choice.
Try it
Take the discrete Fourier transform of two seconds of a 10 hertz sine, then of a 10.25 hertz one, which falls halfway between bins, each with and without a Hann window. Scale by the window's sum, so a sine of amplitude 1 reads 1.
import numpy as np
fs = 500
t = np.arange(0, 2, 1 / fs) # 1,000 samples: bins every 0.5 Hz
for f in (10, 10.25):
x = np.sin(2 * np.pi * f * t)
for name, w in (("none", np.ones_like(t)), ("Hann", np.hanning(len(t)))):
amp = np.abs(np.fft.rfft(x * w)) * 2 / w.sum()
print(f, name, "peak", round(amp.max(), 2),
"bins 9 to 11.5 Hz", np.round(amp[18:24], 2),
"largest above 20 Hz", round(amp[40:].max(), 3))You should seeAt 10 hertz, one bin of 1 with no window, and a wider peak with Hann. At 10.25 hertz, a peak of 0.64 spread across many bins with no window, against 0.85 with Hann and almost nothing far away.
One periodogram is a noisy estimate: its value at each frequency scatters as widely as the value itself, however long the recording. Welch's method cuts the signal into overlapping segments, windows each, takes their periodograms and averages them. K independent segments shrink the scatter by about √K.
The segment length T sets the trade: resolution 1/T against the number of segments to average. Band power is the area under the estimate from f1 to f2. A share, such as alpha's power over the power from 4 to 30 hertz, cancels anything that scales every frequency alike.
Try it
Make a minute of noise with the 1/f spectrum of EEG, add a 10 hertz rhythm, and estimate its spectrum with segments of 1, 2 and 4 seconds. Plot the three, and compute alpha's share from each.
import numpy as np
from scipy.signal import welch
fs, n = 500, 60 * 500
rng = np.random.default_rng(1)
spec = rng.normal(size=n // 2 + 1) + 1j * rng.normal(size=n // 2 + 1)
f = np.fft.rfftfreq(n, 1 / fs)
spec[1:] /= np.sqrt(f[1:]) # power falling as 1/f
spec[0] = 0
x = np.fft.irfft(spec, n)
x = 20 * x / x.std() + 10 * np.sin(2 * np.pi * 10 * np.arange(n) / fs) # microvolts
def band(f, p, lo, hi):
pick = (f >= lo) & (f < hi)
return np.trapezoid(p[pick], f[pick])
for seconds in (1, 2, 4):
f, p = welch(x, fs=fs, nperseg=seconds * fs) # Hann windows, half overlap
print(seconds, "s:", round(band(f, p, 8, 13) / band(f, p, 4, 30), 3))You should seeShares that agree to within a few hundredths, and spectra that differ: smooth but blurred with one-second segments, sharp but ragged with four.
A 9 hertz wave is sampled 12 times a second. At what frequency does it appear, and why?
What frequency resolution do 2 seconds at 500 samples a second give, and which frequencies are there?
Why does a Hann window help a 10.25 hertz sine, and what does it cost?
Four seconds of signal, in two-second segments overlapping by half: how many segments does Welch's method average?
Why remove the mean before estimating a spectrum?
Welch, P. D. (1967). The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics.
Harris, F. J. (1978). On the use of windows for harmonic analysis with the discrete Fourier transform. Proceedings of the IEEE.
Cohen, M. X. (2014). Analyzing Neural Time Series Data: Theory and Practice. MIT Press.
Chapter 2 · Graduate
How do you estimate a spectrum from a few seconds of EEG, with what bias and what variance, and how much of a band's power is a rhythm rather than the background it sits on?
The IntoMind One samples each channel 250, 500 or 1,000 times a second, and every demo uses 500. A frequency above half the rate, the Nyquist frequency, folds back below it: at 500 a second, 300 hertz would appear at 200. The device's own filters, which can be read and set, by default remove drift below 0.5 hertz, notch out the mains frequencies and their harmonics, and low-pass at 200 hertz at this rate.
For N samples xn at rate fs, with a window wn, the periodogram estimates the power spectral density as
P(f) = |Σn wn xn e−2πifn/fs|² / (fs Σn wn²)
doubled for positive frequencies other than zero and the Nyquist. Its bins are fs/N apart: 0.49 hertz for 1,024 samples at 500 a second. Its expectation is the true spectrum smeared by the window's own spectrum, so power leaks from strong frequencies into weak ones. A rectangular window's first sidelobe sits 13 decibels below its main lobe, a Hann window's 31, at the price of a main lobe twice as wide (Harris, 1978). Its variance does not shrink with N: at each frequency it is about the true value times a chi-squared variable with two degrees of freedom, halved, so its standard deviation equals its mean however long the record.
Welch's method averages the periodograms of K windowed segments, usually overlapping by half (Welch, 1967). For independent segments the variance falls as 1/K, and the resolution is set by the segment's length, not the record's.
A multitaper estimate (Thomson, 1982) keeps one segment of length T and averages its periodograms under K orthogonal Slepian tapers (Slepian, 1978), each as concentrated as possible in a band of half-width W. With NW = TW, about K = 2NW − 1 tapers are well concentrated, the variance falls as 1/K, and the resolution is 2W. Two seconds at 2 hertz of resolution give NW = 2 and three tapers. The record's length, divided between resolution and variance, is the whole budget, and for EEG it is set by how long the signal stays the same (Babadi and Brown, 2014). Prerau and others (2017) work through the choice for sleep EEG.
An EEG spectrum is an aperiodic component with peaks above it. Donoghue and others (2020) model it in log power as
log10 P(f) = b − log10(k + fχ) + Σ peaks
with an offset b, a knee k, an exponent χ, and Gaussian peaks. The exponent and the offset differ between people, with age and between states. A band's power, or a ratio of two bands' powers, moves with either of them while no rhythm changes, so a ratio conflates periodic and aperiodic activity (Donoghue, Dominguez and Voytek, 2020).
Alpha's peak frequency differs between people, from about 8 to 12 hertz in healthy adults, and falls with age. The individual alpha frequency is the peak, or the center of gravity, of eyes-closed power in a search window such as 7 to 14 hertz (Klimesch, 1999), with smoothed, automated versions for noisy spectra (Corcoran and others, 2018). A fixed band of 8 to 13 hertz measures a different part of each person's alpha.
A spectrogram slides one window along the signal. A wavelet transform scales the window with frequency: a Morlet wavelet of n cycles at frequency f spreads over σt = n/(2πf) in time and σf = f/n in frequency, and their product is fixed at 1/(2π) (Cohen, 2014). Averaging power over time or trials can make brief bursts look like a sustained rhythm, as Jones (2016) and van Ede and others (2018) point out.
Alpha Flower reports alpha's share of the power from 4 to 30 hertz, from one Hann-windowed periodogram of 2.05 seconds every quarter second. Its Python program averages three Hann-windowed periodograms of 2 seconds over the last 4, once a second. The analysis compares the two with a multitaper estimate.
Analysis
Simulate twenty minutes of a signal whose spectrum you know exactly: an aperiodic component falling as f−χ above 1 hertz, plus an alpha rhythm, a resonance at 10 hertz. Estimate alpha's share in 300 windows of 4 seconds three ways: as the demo does, as its Python program does, and with three Slepian tapers. Then make the background steeper with the alpha unchanged, and measure alpha above a fitted background.
import numpy as np
from scipy.signal import welch, lfilter
from scipy.signal.windows import dpss, hann
fs, rng = 500, np.random.default_rng(7)
f = np.linspace(0.5, 60, 20000) # a fine grid for the true spectrum
r, theta = 0.99, 2 * np.pi * 10 / fs # alpha: an AR(2) resonance at 10 Hz
ar = [1, -2 * r * np.cos(theta), r * r]
def simulate(chi, minutes=20):
"""Aperiodic noise falling as f^-chi above 1 Hz, plus alpha, and the true spectra."""
n = minutes * 60 * fs
shape = np.maximum(np.fft.rfftfreq(n, 1 / fs), 1) ** (-chi / 2)
spec = np.fft.rfft(rng.normal(size=n)) * shape
x = 100 * np.fft.irfft(spec, n) + lfilter([1], ar, rng.normal(0, 0.1, n))
z = np.exp(-2j * np.pi * f / fs)
alpha = 0.1 ** 2 * 2 / fs / np.abs(1 + ar[1] * z + ar[2] * z ** 2) ** 2
return x, 100 ** 2 * 2 / fs * np.maximum(f, 1) ** -chi, alpha
def band(freqs, psd, lo, hi, rule=np.trapezoid):
pick = (freqs >= lo) & (freqs < hi)
return rule(psd[pick], freqs[pick])
def share(freqs, psd, rule=np.trapezoid):
return band(freqs, psd, 8, 13, rule) / band(freqs, psd, 4, 30, rule)
def demo(x): # one periodic-Hann periodogram of the first 1,024 samples, bins summed
seg = x[:1024] - x[:1024].mean()
p = np.abs(np.fft.rfft(seg * hann(1024, sym=False))) ** 2
return share(np.fft.rfftfreq(1024, 1 / fs), p, rule=lambda y, x: y.sum())
def program(x): # Welch over 4 s in 2-s segments, as eyes_closed.py does
freqs, p = welch(x - x.mean(), fs=fs, nperseg=2 * fs)
return share(freqs, p)
def multitaper(x): # the same 1,024 samples, 3 Slepian tapers, about 2 Hz of resolution
seg = x[:1024] - x[:1024].mean()
p = (np.abs(np.fft.rfft(dpss(1024, 2, 3) * seg, axis=1)) ** 2).mean(axis=0)
return share(np.fft.rfftfreq(1024, 1 / fs), p)
for chi in (1.5, 2.0):
x, background, alpha = simulate(chi)
true = background + alpha
print(f"exponent {chi}: true share {share(f, true):.3f}, "
f"true alpha power {band(f, alpha, 6, 14):.1f} uV^2")
windows = x.reshape(-1, 4 * fs) # 300 windows of 4 s
for name, estimate in (("demo", demo), ("program", program), ("multitaper", multitaper)):
v = np.array([estimate(w) for w in windows])
bias = v.mean() - share(f, true)
print(f" {name:10s} share {v.mean():.3f} bias {bias:+.3f} sd {v.std():.3f}")
# Fit the background in log-log, away from alpha, and measure alpha above it.
freqs, p = welch(x, fs=fs, nperseg=2 * fs)
freqs, p = freqs[1:], p[1:]
fit = ((freqs >= 2) & (freqs < 6)) | ((freqs >= 15) & (freqs < 40))
slope, offset = np.polyfit(np.log10(freqs[fit]), np.log10(p[fit]), 1)
above = p - 10 ** (offset + slope * np.log10(freqs))
print(f" fitted exponent {-slope:.2f}, "
f"alpha above the fit {band(freqs, above, 6, 14):.1f} uV^2")You should seeA true share of 0.464. The demo's estimate scatters most, with a standard deviation of 0.13 from one window to the next, the program's least, 0.08, and the multitaper's in between, 0.10, on the same two seconds as the demo. Every estimator is low by 0.01 to 0.03. With the exponent at 2.0, the true share rises to 0.613 though the alpha is the same, and alpha above the fitted background stays near its true 14.3 µV², at 13.4 and 13.7.
Then report: the program's estimator with segments of 1 and 4 seconds, and every estimator with the alpha at 8.5 hertz. Which bias grows, and why?
At 500 samples a second, what is the bin spacing of a 1,024-sample periodogram? With a Hann window, about how long must a segment be to resolve two equal peaks at 9.5 and 10.5 hertz?
One reading of the demo scatters by about 0.13 around its mean. Why doesn't the demo estimate the share over 10 seconds instead, and what do its two thresholds, 0.36 to count the eyes as closed and 0.27 to count them open again, do about the scatter?
For 2.048 seconds and 2 hertz of resolution, what NW, and how many tapers? What would 4.096 seconds at the same resolution give?
Two people have identical alpha rhythms, but the second's aperiodic exponent is 2.0 rather than 1.5, with the same offset. Whose share is larger? What does that do to a study that compares two groups by their share?
A Morlet wavelet of 7 cycles: what are its spreads in time and frequency at 10 hertz, and at 40?
Welch, P. D. (1967). The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics.
Harris, F. J. (1978). On the use of windows for harmonic analysis with the discrete Fourier transform. Proceedings of the IEEE.
Slepian, D. (1978). Prolate spheroidal wave functions, Fourier analysis, and uncertainty, V: the discrete case. Bell System Technical Journal.
Thomson, D. J. (1982). Spectrum estimation and harmonic analysis. Proceedings of the IEEE.
Klimesch, W. (1999). EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis. Brain Research Reviews.
Babadi, B., and Brown, E. N. (2014). A review of multitaper spectral analysis. IEEE Transactions on Biomedical Engineering.
Cohen, M. X. (2014). Analyzing Neural Time Series Data: Theory and Practice. MIT Press.
Jones, S. R. (2016). When brain rhythms aren't “rhythmic”: implication for their mechanisms and meaning. Current Opinion in Neurobiology.
Gao, R., Peterson, E. J., and Voytek, B. (2017). Inferring synaptic excitation/inhibition balance from field potentials. NeuroImage.
Prerau, M. J., and others (2017). Sleep neurophysiological dynamics through the lens of multitaper spectral analysis. Physiology.
Corcoran, A. W., Alday, P. M., Schlesewsky, M., and Bornkessel-Schlesewsky, I. (2018). Toward a reliable, automated method of individual alpha frequency (IAF) quantification. Psychophysiology.
van Ede, F., Quinn, A. J., Woolrich, M. W., and Nobre, A. C. (2018). Neural oscillations: sustained rhythms or transient burst-events? Trends in Neurosciences.
Donoghue, T., and others (2020). Parameterizing neural power spectra into periodic and aperiodic components. Nature Neuroscience.
Donoghue, T., Dominguez, J., and Voytek, B. (2020). Electrophysiological frequency band ratio measures conflate periodic and aperiodic neural activity. eNeuro.