Spectral Analysis of Sleep EEG

mrpheus provides two complementary spectral analysis functions: compute_band_power() for per-epoch power in standard EEG frequency bands, and compute_temporal_bandpower() for tracking how band power evolves across the full recording night via a sliding window. Both use a Welch PSD implementation validated against YASA’s Python bandpower() function (absolute power MARE < 0.08 %, relative power MARE < 0.8 %).

Note — All code blocks use eval = FALSE because they require a local EDF file. For a fully-executed walkthrough including figures see the companion article: Spectral Analysis: Live Walkthrough

Data

The examples use SC4001E0-PSG.edf from the Sleep-EDF Cassette study (Kemp et al., 2000). See the PSG Preprocessing Pipeline vignette for download instructions.

Setup

library(mrpheus)
library(dplyr)
library(ggplot2)

rec    <- read_edf("SC4001E0-PSG.edf")
psg    <- prepare_psg(rec) |> preprocess_psg()
stages <- stage_epochs(psg)   # data frame with columns: epoch, stage, ...

Band power per epoch

compute_band_power() estimates the Welch PSD for each epoch and integrates power within six standard EEG bands.

bp <- compute_band_power(psg)
bp
#> # A tibble: 2,650 × 9
#>   epoch channel      delta theta alpha sigma   beta  gamma total_power
#>   <int> <chr>        <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>       <dbl>
#> 1     1 EEG Fpz-Cz  342.  28.4  12.1  4.81   3.21   1.05        391.
#> ...

Default band definitions match YASA’s bandpower():

Band	Range
delta	0.5–4 Hz
theta	4–8 Hz
alpha	8–12 Hz
sigma	12–16 Hz
beta	16–30 Hz
gamma	30–45 Hz

Relative power

Set relative = TRUE to express each band as a fraction of the total band-summed power. This is the standard for between-subject comparisons as it removes the effect of overall signal amplitude.

bp_rel <- compute_band_power(psg, relative = TRUE)

Band power by sleep stage

Join the band power table with the scored hypnogram to compute mean power per sleep stage:

stage_bp <- bp_rel |>
  left_join(stages |> select(epoch, stage), by = "epoch") |>
  group_by(stage) |>
  summarise(across(delta:gamma, mean, na.rm = TRUE))

stage_bp

Custom bands

The bands argument accepts any named list of frequency ranges:

bp_custom <- compute_band_power(
  psg,
  bands = list(
    slow_osc  = c(0.5, 1),
    delta     = c(1,   4),
    theta     = c(4,   8),
    spindle   = c(11, 16)   # narrower sigma for spindle-specific power
  )
)

Welch parameters

The defaults (win_sec = 4, noverlap = 200, nfft = 1024) match YASA. Override them for higher frequency resolution or different smoothing:

# Higher frequency resolution (longer window)
bp_hires <- compute_band_power(psg, win_sec = 8, noverlap = 400L, nfft = 2048L)

Temporal band power

compute_temporal_bandpower() applies a sliding window over the full recording, concatenating epochs within each window before computing the PSD. This reveals how spectral composition changes across the night — the classic homeostatic decline of slow-wave (delta) power and the cyclic appearance of REM-associated theta and alpha.

tbp <- compute_temporal_bandpower(psg, hypno = stages$stage)
tbp
#> # A tibble: 220 × 7
#>   time_hours epoch_start epoch_end dominant_stage band    power relative_power
#>        <dbl>       <int>     <int>          <int> <chr>   <dbl>          <dbl>
#> 1      0           1        10              2 delta 312.           0.791
#> ...

Each row represents one (window, band) combination. By default windows are 10 epochs (5 minutes) wide with 50 % overlap.

Customising the window

# Finer temporal resolution: 6-epoch (3 min) windows, 50 % overlap
tbp <- compute_temporal_bandpower(
  psg,
  hypno         = stages$stage,
  window_epochs = 6L,
  step_epochs   = 3L
)

# Non-overlapping windows
tbp <- compute_temporal_bandpower(
  psg,
  hypno         = stages$stage,
  window_epochs = 10L,
  step_epochs   = 10L
)

Visualising temporal dynamics

The long-format output is designed for direct use with ggplot2.

Delta power across the night (slow-wave homeostasis):

tbp |>
  filter(band == "delta") |>
  ggplot(aes(time_hours, relative_power)) +
  geom_line(colour = "#014370") +
  geom_smooth(method = "lm", se = FALSE,
              colour = "#B83E2C", linewidth = 0.6, linetype = "dashed") +
  labs(
    x = "Time (hours from recording onset)",
    y = "Relative delta power",
    title = "Slow-wave homeostasis: delta power across the night"
  ) +
  theme_minimal()

All bands as a faceted panel (matches the notebook’s temporal plot):

band_colours <- c(
  delta = "#014370", theta = "#1B6799", alpha = "#9BDFE2",
  sigma = "#FFA75D", beta  = "#FC544A", gamma = "#FFECD4"
)

tbp |>
  ggplot(aes(time_hours, relative_power, colour = band)) +
  geom_line(alpha = 0.8) +
  facet_wrap(~band, scales = "free_y", ncol = 2) +
  scale_colour_manual(values = band_colours, guide = "none") +
  labs(
    x = "Time (hours)",
    y = "Relative power",
    title = "EEG band power across the night"
  ) +
  theme_minimal()

Colouring by dominant sleep stage:

stage_labels <- c("0" = "Wake", "1" = "N1", "2" = "N2",
                  "3" = "N3",   "4" = "REM")

tbp |>
  filter(band == "delta") |>
  mutate(stage_label = stage_labels[as.character(dominant_stage)]) |>
  ggplot(aes(time_hours, relative_power, colour = stage_label)) +
  geom_point(size = 1.5, alpha = 0.8) +
  scale_colour_brewer(palette = "Set2", name = "Stage") +
  labs(x = "Time (hours)", y = "Relative delta power") +
  theme_minimal()

Combining band power and staging

A common analysis pattern: compute per-stage mean band power, then fit a linear model to test how it relates to clinical variables.

# Per-stage mean relative band power
stage_summary <- bp_rel |>
  left_join(stages |> select(epoch, stage), by = "epoch") |>
  group_by(stage) |>
  summarise(across(delta:gamma, \(x) mean(x, na.rm = TRUE)), .groups = "drop") |>
  mutate(stage = factor(stage, levels = 0:4,
                        labels = c("Wake", "N1", "N2", "N3", "REM")))

# Bar chart of mean relative power per stage
stage_summary |>
  tidyr::pivot_longer(delta:gamma, names_to = "band", values_to = "power") |>
  ggplot(aes(stage, power, fill = band)) +
  geom_col(position = "dodge") +
  scale_fill_brewer(palette = "Blues") +
  labs(x = "Sleep stage", y = "Mean relative power", fill = "Band") +
  theme_minimal()

Parity with YASA

compute_band_power() was validated against YASA’s bandpower() on SC4001E0-PSG.edf (2,650 epochs, channel EEG Fpz-Cz):

Band	Abs. MARE	Rel. MARE
delta	0.079 %	0.771 %
theta	0.059 %	0.750 %
alpha	0.048 %	0.519 %
sigma	0.054 %	0.747 %
beta	0.032 %	0.748 %
gamma	0.034 %	0.748 %

The residual ~0.75 % relative error is due to floating-point accumulation across six bands in two different language runtimes and is not fixable without artificial precision alignment. It is scientifically negligible for any EEG analysis.

References

Kemp B, Zwinderman A, Tuk B, Kamphuisen H, Oberyé J (2000). Analysis of a sleep-dependent neuronal feedback loop: the slow-wave microcontinuity of the EEG. IEEE Transactions on Biomedical Engineering, 47(9), 1185–1194.

Vallat R, Walker MP (2021). An open-source, high-performance tool for automated sleep staging. eLife, 10, e70092.