---
title: "Spectral Analysis of Sleep EEG"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Spectral Analysis of Sleep EEG}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  eval     = FALSE
)
```

**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](https://mrpheus.circadia-lab.uk/articles/spectral-demo.html)

---

## Data

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

---

## Setup

```r
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.

```r
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.

```r
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:

```r
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:

```r
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:

```r
# 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.

```r
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

```r
# 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):

```r
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):

```r
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:**

```r
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.

```r
# 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.