EEG and LFP spectra in Python: the Fourier transform from scratch

In lesson 5 the first thing we did to the raw voltage was throw away everything slower than 300 Hz, because the spikes lived above it. This lesson is about what we threw away. Below 300 Hz, the voltage near an electrode is the summed activity of thousands of neurons at once, the local field potential, and recorded from the scalp through skull and skin it is the EEG. It is where the brain’s rhythms live, and the tool that reveals them is the Fourier transform, which we will build from a dot product before using NumPy’s. By the end you will have reproduced, on real recordings, the first discovery ever made with EEG and the best-known rhythm in the rodent brain.

flowchart LR

A[Frequency: match a signal against sines] --> B[Power spectrum: average over segments]

B --> C[EEG: eyes open vs closed]

C --> D[Spectrogram: power over time]

B --> E[Mouse LFP: 16 depths]

E --> F[Band-pass: pull one rhythm out]

The lesson. The first box is the only new idea; the rest is applying it.

Step 1: what a frequency is, and how to measure it

Any signal can be written as a sum of sine waves at different frequencies, each with its own size and its own shift in time. That is a theorem, not a metaphor, and it means a recording can be described two ways: as a voltage at each moment, or as an amount at each frequency. The Fourier transform converts one into the other. Here is the whole idea in one function: to ask how much of a given frequency is in a signal, multiply the signal by a sine wave at that frequency and add up. Where the two line up, the products are positive; where they do not, they cancel. Do it with a cosine as well, because the rhythm might be shifted in time, and combine the two.

import numpy as np
import matplotlib.pyplot as plt

fs = 200                                                   # samples per second
t = np.arange(0, 4, 1 / fs)                                # four seconds
signal = 2 * np.sin(2 * np.pi * 3 * t) + 1 * np.sin(2 * np.pi * 10 * t + 1) + 0.5 * np.sin(2 * np.pi * 27 * t)

def amplitude_at(x, freq, fs):
    """How much of frequency `freq` is in x: match it against a cosine and a sine of that frequency."""
    tt = np.arange(len(x)) / fs
    c = np.sum(x * np.cos(2 * np.pi * freq * tt))            # how well x lines up with a cosine
    s = np.sum(x * np.sin(2 * np.pi * freq * tt))            # and with a sine (same frequency, shifted)
    return 2 * np.sqrt(c ** 2 + s ** 2) / len(x)             # combine them: the amplitude, whatever the phase

freqs = np.arange(0, 50, 0.25)
amps = np.array([amplitude_at(signal, f, fs) for f in freqs])
for f in freqs[amps > 0.3]:
    print(f"{f:5.2f} Hz: amplitude {amps[freqs == f][0]:.2f}")

spectrum = np.abs(np.fft.rfft(signal)) * 2 / len(signal)   # the same thing, every frequency at once, in a millisecond
f_fft = np.fft.rfftfreq(len(signal), 1 / fs)
print("FFT agrees at 10 Hz:", spectrum[f_fft == 10].round(2))
#  3.00 Hz: amplitude 2.00
# 10.00 Hz: amplitude 1.00
# 27.00 Hz: amplitude 0.50
# FFT agrees at 10 Hz: [1.]

A made-up signal built from three sine waves, and its spectrum, read back from the signal by the matching function. The three ingredients come out at the right frequencies with the right amplitudes. np.fft.rfft gives the same answer for every frequency at once, which is the only reason to use it instead of the loop.

The loop is slow and the fast Fourier transform, FFT, is the same calculation organised cleverly, which is why everyone uses it. But it is worth having written the loop once: a Fourier transform is nothing more than a list of dot products with sine waves, one per frequency, and everything from here on is that.

Step 2: the first rhythm

Hans Berger recorded the first human EEG in Jena in 1924 and published in 1929. His first observation is still the clearest in the field: when a resting subject closes their eyes, a 10-cycle-per-second wave appears over the back of the head and vanishes when they open them. He called it the alpha rhythm because it was the first thing he saw. Here it is, from a subject in the PhysioNet motor-imagery dataset, one minute with eyes open and one with eyes closed, at electrode Oz over visual cortex.

eeg = np.load("eeg_s001_open_closed.zip")
fs = float(eeg["fs"])                                      # 160 samples per second
names = list(eeg["channels"])
oz = names.index("Oz")                                     # the electrode at the back of the head, over visual cortex
open_eyes = eeg["eyes_open"][oz] - eeg["eyes_open"][oz].mean()
closed_eyes = eeg["eyes_closed"][oz] - eeg["eyes_closed"][oz].mean()
te = np.arange(len(open_eyes)) / fs

fig, axes = plt.subplots(2, 1, figsize=(10, 4.6), sharex=True, sharey=True)
axes[0].plot(te, open_eyes, lw=0.7); axes[0].set(ylabel="uV", title="Oz, eyes open")
axes[1].plot(te, closed_eyes, lw=0.7); axes[1].set(ylabel="uV", title="Oz, eyes closed", xlabel="time (s)", xlim=(10, 15))
plt.show()

Five seconds at Oz, eyes open above and eyes closed below. Even by eye the lower trace is more regular: count the waves in any one second of it and you get about ten.

Step 3: the power spectrum

To put a number on it, measure how much of each frequency is present. There is one practical wrinkle. A Fourier transform of a whole minute gives a very finely spaced and very noisy spectrum, because the brain does not hold a frequency steady for a minute. The fix, due to Peter Welch in 1967, is to cut the recording into short segments, transform each one, and average the power. Each segment is also tapered at the edges with a window so that the cut itself does not look like a sudden event. Ten lines.

def power_spectrum(x, fs, seconds=2.0):
    """Average power at each frequency, estimated from many short segments (Welch's method)."""
    n = int(seconds * fs)
    window = np.hanning(n)                                 # taper each segment so its edges do not look like jumps
    freqs = np.fft.rfftfreq(n, 1 / fs)
    power = np.zeros(len(freqs))
    segments = 0
    for start in range(0, len(x) - n + 1, n // 2):         # segments overlap by half
        seg = (x[start:start + n] - x[start:start + n].mean()) * window
        power += np.abs(np.fft.rfft(seg)) ** 2
        segments += 1
    return freqs, power / segments / (fs * np.sum(window ** 2))   # scaled to microvolts squared per Hz

f, p_open = power_spectrum(open_eyes, fs)
f, p_closed = power_spectrum(closed_eyes, fs)
alpha = (f >= 8) & (f <= 12)
print(f"alpha band power at Oz: eyes open {p_open[alpha].mean():5.1f}, eyes closed {p_closed[alpha].mean():5.1f}, ratio {p_closed[alpha].mean() / p_open[alpha].mean():.1f}")
print(f"peak of the eyes-closed spectrum between 4 and 30 Hz: {f[(f > 4) & (f < 30)][np.argmax(p_closed[(f > 4) & (f < 30)])]:.1f} Hz")
# alpha band power at Oz: eyes open  24.3, eyes closed 327.4, ratio 13.5
# peak of the eyes-closed spectrum between 4 and 30 Hz: 10.0 Hz

fig, ax = plt.subplots(figsize=(10, 4))
ax.semilogy(f, p_open, label="eyes open")
ax.semilogy(f, p_closed, label="eyes closed")
ax.axvspan(8, 12, color="orange", alpha=0.15)
ax.set(xlabel="frequency (Hz)", ylabel="power (uV$^2$ / Hz)", xlim=(0, 60))
ax.legend()
plt.show()

Power at each frequency for Oz, eyes open in grey and eyes closed in blue, on a log scale with the conventional bands shaded. Eyes closed, a peak appears at exactly 10 Hz with 13 times the alpha power, and a smaller bump near 20 Hz, its harmonic. Everything else is similar. The downward slope on both curves is the 1/f background that every brain recording has.

Thirteen times more alpha power with the eyes closed, with the peak at 10.0 Hz. That is Berger’s result, ninety-seven years on, in a few lines. Two more things in the figure deserve a look. Both spectra fall steadily with frequency: this 1/f shape is universal in brain signals and means that comparing power across bands is only meaningful relative to that slope, which is why the ratio between conditions is a safer number than the raw power. And the bands themselves are conventions, not physics: useful labels for regions of the spectrum where particular rhythms tend to appear.

Band	Frequency	Where you see it
delta	1 to 4 Hz	Deep sleep; large and slow, everywhere
theta	4 to 8 Hz (6 to 10 in rodents)	Hippocampus during movement and memory; frontal midline during effort
alpha	8 to 12 Hz	Visual cortex at rest with the eyes closed; the first rhythm ever recorded
beta	13 to 30 Hz	Motor cortex holding still; drops before a movement
gamma	30 to 100 Hz and up	Local processing; weak at the scalp, strong on the probe

Step 4: where it is

Alpha is a rhythm of visual cortex at rest, so it should be strongest at the back of the head and weaker at the front. The same function on four electrodes:

for ch in ["Oz", "O1", "Cz", "Fp1"]:
    i = names.index(ch)
    _, po = power_spectrum(eeg["eyes_open"][i], fs)
    _, pc = power_spectrum(eeg["eyes_closed"][i], fs)
    print(f"{ch:4s} alpha power, closed / open: {pc[alpha].mean() / po[alpha].mean():4.1f}")
# Oz   alpha power, closed / open: 13.5
# O1   alpha power, closed / open: 15.2
# Cz   alpha power, closed / open:  4.2
# Fp1  alpha power, closed / open:  2.6

Fifteen-fold at O1 and thirteen-fold at Oz, both at the back; four-fold at the vertex; under three-fold at the forehead. The rhythm is generated in and around visual cortex, and what Fp1 sees is the tail of it conducted through the head. The ICA primer showed the same conduction working the other way, carrying eye blinks backwards.

Step 5: power over time

A spectrum averages over the whole recording, which hides when things happen. A spectrogram keeps the time axis: compute a spectrum in a short window, slide the window along, and stack the results into an image with time across and frequency up. The function is the previous one without the averaging.

def spectrogram(x, fs, window_s=1.0, step_s=0.25):
    """Power spectrum in a sliding window: frequency x time."""
    n, step = int(window_s * fs), int(step_s * fs)
    taper = np.hanning(n)
    starts = np.arange(0, len(x) - n, step)
    freqs = np.fft.rfftfreq(n, 1 / fs)
    S = np.array([np.abs(np.fft.rfft((x[s:s + n] - x[s:s + n].mean()) * taper)) ** 2 for s in starts]).T
    return freqs, (starts + n / 2) / fs, S

both = np.concatenate([open_eyes, closed_eyes])            # one minute eyes open, then one minute eyes closed
f_s, t_s, S = spectrogram(both, fs)

fig, ax = plt.subplots(figsize=(10, 4))
ax.pcolormesh(t_s, f_s, np.log10(S), shading="auto")
ax.axvline(len(open_eyes) / fs, color="white", ls="--")
ax.set(ylim=(0, 40), xlabel="time (s)", ylabel="frequency (Hz)")
plt.show()

Oz over two minutes, the eyes-open minute followed by the eyes-closed minute, as a spectrogram: power in one-second windows, brighter is more. The 10 Hz band switches on at the dashed line and stays on, with brief gaps. The low-frequency glow at the bottom is the 1/f background.

This is how most EEG is actually looked at, because the interesting questions are about change: what happens when the eyes close, when a stimulus appears, when a seizure starts. The trade-off is built into the method: a one-second window cannot tell 10 Hz from 10.5 Hz, and a ten-second window cannot tell you which second the rhythm began in. You choose the window for the question.

Step 6: inside the brain

The scalp sees the whole cortex blurred together. A probe sees one place at a time, and the Neuropixels probe from lesson 5 has 384 sites along its length, so it sees many places at once. This file holds 40 seconds of the slow band, sampled at 2,500 Hz, from 16 sites spaced along a probe in the same mouse, the other probe from the lesson 5 session. The IBL’s histology puts the sites in visual cortex at the top, then white matter, then the hippocampus, then the thalamus at the tip. Compute a spectrum at every depth and ask where the theta band, 6 to 10 Hz in rodents, stands out from the background.

lf = np.load("ibl_ks023_lfp.zip")
lfp = lf["lfp"].astype(float) * float(lf["uv_per_bit"])    # samples x 16 channels, microvolts
fs_lfp = float(lf["fs"])                                   # 2500 samples per second
depth, region = lf["depth_um"], lf["area"]
print(len(lfp) / fs_lfp, "seconds;", "depths (um from tip):", depth.astype(int))
print("regions:", list(region))

theta_ratio = []
for k in range(lfp.shape[1]):
    f_l, p_l = power_spectrum(lfp[:, k], fs_lfp, seconds=2.0)
    theta = (f_l >= 6) & (f_l <= 10)
    delta = (f_l >= 1) & (f_l <= 4)
    theta_ratio.append(p_l[theta].mean() / p_l[delta].mean())
for d, r, ratio in zip(depth, region, theta_ratio):
    print(f"{d:5.0f} um  {r:8s} theta / delta = {ratio:.2f}")
# 40.0 seconds; depths (um from tip): [  40  300  540  800 1040 1300 1540 1800 2040 2300 2540 2800 3040 3300
#  3540 3800]
# regions: ['ml', 'VPM', 'PO', 'PO', 'PO', 'LP', 'LP', 'DG-mo', 'DG-po', 'DG-mo', 'CA1', 'CA1', 'scwm', 'VISa5', 'VISa4', 'VISa2/3']
#    40 um  ml       theta / delta = 0.09
#   300 um  VPM      theta / delta = 0.10
#   540 um  PO       theta / delta = 0.10
#   800 um  PO       theta / delta = 0.10
#  1040 um  PO       theta / delta = 0.11
#  1300 um  LP       theta / delta = 0.13
#  1540 um  LP       theta / delta = 0.15
#  1800 um  DG-mo    theta / delta = 0.20
#  2040 um  DG-po    theta / delta = 0.29
#  2300 um  DG-mo    theta / delta = 0.57
#  2540 um  CA1      theta / delta = 0.55
#  2800 um  CA1      theta / delta = 0.11
#  3040 um  scwm     theta / delta = 0.10
#  3300 um  VISa5    theta / delta = 0.08
#  3540 um  VISa4    theta / delta = 0.07
#  3800 um  VISa2/3  theta / delta = 0.06

Left: three seconds of LFP at four of the sixteen depths. Right: theta power relative to the slow background at every depth, with the brain region from the IBL's histology. Theta is strongest in the dentate gyrus and CA1 of the hippocampus, five times the ratio seen in cortex or thalamus, and in the traces the hippocampal sites are visibly swinging at about seven cycles per second while cortex and thalamus barely move. — Left: three seconds of LFP at four of the sixteen depths. Right: theta power relative to the slow background at every depth, with the brain region from the IBL’s histology. Theta is strongest in the dentate gyrus and CA1 of the hippocampus, five times the ratio seen in cortex or thalamus, and in the traces the hippocampal sites are visibly swinging at about seven cycles per second while cortex and thalamus barely move.

The hippocampus is the one part of this probe that oscillates, and it does so at about 7 Hz, which is the rodent theta rhythm: the signature of a mouse that is moving or attending, generated by a pacemaker in the septum and organising the hippocampus into cycles within which place cells fire in sequence. The peak of the ratio sits at 2,300 micrometres, in the dentate gyrus molecular layer just below CA1, which is where the textbooks put the largest theta amplitude, at the hippocampal fissure. Four hundred micrometres further up, in the CA1 cell layer itself, the ratio has already fallen to 0.11, because the field reverses polarity across the layer and nearly cancels there. The LFP is a spatial signal as much as a temporal one.

Step 7: keep what we threw away

Lesson 5 used a filter to remove the slow activity. Use the same tool to keep only a band of it: a band-pass filter from 6 to 10 Hz, applied to the theta-rich site and to a site at the probe tip, pulls the rhythm out of the raw trace so you can see it against the signal it came from.

from scipy.signal import butter, filtfilt
b, a = butter(3, [6 / (fs_lfp / 2), 10 / (fs_lfp / 2)], "band")   # keep only 6 to 10 Hz
k = int(np.argmax(theta_ratio))                                      # the channel where theta is strongest
theta_wave = filtfilt(b, a, lfp[:, k])
tl = np.arange(len(lfp)) / fs_lfp

fig, axes = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
axes[0].plot(tl, lfp[:, k], lw=0.7)
axes[0].plot(tl, theta_wave, lw=1.4)
axes[0].set(ylabel="uV", title=f"{region[k]}, {depth[k]:.0f} um: raw LFP and its 6 to 10 Hz part")
axes[1].plot(tl, lfp[:, 0], lw=0.7)
axes[1].plot(tl, filtfilt(b, a, lfp[:, 0]), lw=1.4)
axes[1].set(ylabel="uV", xlabel="time (s)", xlim=(10, 13), title=f"{region[0]}, {depth[0]:.0f} um: the same")
plt.show()

Top: the dentate gyrus site, raw LFP in grey and its 6 to 10 Hz component in blue: a near-sinusoidal theta wave of several hundred microvolts, riding under the sharper fluctuations. Bottom: the same filter at the probe tip, in a fibre tract below the thalamus, finds a wave fifty times smaller. A filter always returns something in its band; whether it means anything is what the spectrum tells you.

That last point is the warning for this lesson. A band-pass filter will produce a 6 to 10 Hz wave from any signal, including pure noise, because that is what it does. The lower trace shows it: a neat small sinusoid from a site with no theta rhythm at all. Filtered traces are for seeing; the spectrum is for deciding.

What you just did

Measured the amount of a frequency in a signal with a dot product, and understood the FFT as that, done for every frequency at once.
Written Welch’s method for a power spectrum and a spectrogram, each in about ten lines.
Reproduced Berger’s 1929 discovery: a 10 Hz alpha rhythm over visual cortex that appears when the eyes close, thirteen times stronger, and fades toward the front of the head.
Mapped a rhythm to a place in the brain: hippocampal theta at 7 Hz, strongest at the dentate gyrus, absent in cortex and thalamus on the same probe.
Used a band-pass filter to isolate a rhythm, and learned why a filtered trace is not evidence that the rhythm exists.

Exercises

Change the segment length in power_spectrum to 0.5 and then 8 seconds. What happens to the sharpness of the alpha peak, and to the roughness of the curve? This is the same trade-off as the spectrogram window.
Repeat step 3 for every one of the 64 channels and plot the closed/open alpha ratio on the schematic head from the ICA primer. Where exactly is the rhythm?
In the LFP file, compute the spectrum at the hippocampal sites in two-second windows across the 40 seconds. Does theta come and go? Mice run in bouts; theta should follow.
Harder: at the CA1 site at 2,800 um the theta ratio was low, and the text claims the field reverses polarity across the cell layer. Test it: band-pass both CA1 sites and plot them together. Are the theta waves in phase or opposite?

Data: EEG from the PhysioNet EEG Motor Movement/Imagery Dataset (Schalk et al. 2004; Goldberger et al. 2000), subject 1, runs 1 and 2, Open Data Commons Attribution licence. LFP from the International Brain Laboratory brain-wide map, CC-BY 4.0, session KS023 2019-12-10, probe 01, LF band, seconds 1500 to 1540, with brain regions from the IBL’s histology alignment. The code runs top to bottom in a few seconds.

Programming In Neuroscience

What the filter threw away: EEG, LFP and the Fourier transform

Step 1: what a frequency is, and how to measure it

Step 2: the first rhythm

Step 3: the power spectrum

Step 4: where it is

Step 5: power over time

Step 6: inside the brain

Step 7: keep what we threw away

What you just did

Exercises

Like this:

Related

Step 1: what a frequency is, and how to measure it

Step 2: the first rhythm

Step 3: the power spectrum

Step 4: where it is

Step 5: power over time

Step 6: inside the brain

Step 7: keep what we threw away

What you just did

Exercises

Share this:

How a neural network learns: a playground

How a brain works, how a language model works, and where they meet

Detecting a seizure: from scalp EEG to a Raspberry Pi monitor

Like this:

Related

Discover more from Programming In Neuroscience