# Apply AASM-recommended bandpass filter
eeg_filt <- eegFilter(
  eeg,
  lowcut = 0.3,
  highcut = 35,
  order = 4,
  assay_name = "raw",
  output_assay = "filtered"
)

# Re-referencing to linked mastoids (A1+A2)/2 is common for sleep
# This would require mastoid channels in real data
# For this example, we'll use the average reference
eeg_filt <- eegRereference(
  eeg_filt,
  ref_type = "average",
  assay_name = "filtered",
  output_assay = "reref"
)

# Flag 30-second epochs whose peak amplitude exceeds +/-150 uV as artifacts
reref <- assay(eeg_filt, "reref")
epoch_len <- 30 * samplingRate(eeg_filt)
n_epochs <- nrow(reref) %/% epoch_len
artifact_epochs <- Filter(function(e) {
  seg <- reref[((e - 1) * epoch_len + 1):(e * epoch_len), , drop = FALSE]
  max(abs(seg)) > 150
}, seq_len(n_epochs))
artifact_epochs <- unlist(artifact_epochs)

# Mark artifact epochs for exclusion from analysis
metadata(eeg_filt)$artifact_epochs <- artifact_epochs

# Perform automatic sleep staging in 30-second epochs
stages <- eegSleepStage(
  eeg_filt,
  epoch_sec = 30,  # AASM standard epoch length
  assay_name = "reref"
)

# eegSleepStage() returns a data.frame with one row per epoch
head(stages)

# Sequence of stage labels and their distribution
head(stages$stage)
table(stages$stage)

# Spectral power features used for classification
grep("_power$", names(stages), value = TRUE)

# Store staging results in metadata (used by eegSleepMetrics)
metadata(eeg_filt)$sleep_stages <- stages

# Detect sleep spindles across all channels
spindles <- eegSpindleDetect(
  eeg_filt,
  freq_range = c(11, 16),      # Sigma band
  min_duration_ms = 500,       # Minimum duration
  max_duration_ms = 2000,      # Maximum duration
  threshold_sd = 2.0,          # Amplitude threshold in SD of the background
  assay_name = "reref"
)

# eegSpindleDetect() returns a data.frame of spindle events
# (channel, start/end sample, duration_ms, peak_sample, peak_amplitude,
#  frequency_hz)
head(spindles)

# Spindle count per channel (channel index)
table(spindles$channel)

# Overall spindle density (spindles per minute of recording)
rec_minutes <- nrow(assay(eeg_filt, "reref")) / samplingRate(eeg_filt) / 60
nrow(spindles) / rec_minutes

# Visualize spindle onsets over time (hours)
if (nrow(spindles) > 0) {
  onset_hours <- spindles$start_sample / samplingRate(eeg_filt) / 3600
  hist(onset_hours,
       xlab = "Time (hours)", ylab = "Spindle Count",
       main = "Sleep Spindle Distribution")
}

# Detect K-complexes (most prominent in central leads C3/C4)
kcomplexes <- eegKcomplexDetect(
  eeg_filt,
  min_neg_amplitude = 75,   # μV negative deflection
  min_duration_ms = 500,
  max_duration_ms = 1500,
  assay_name = "reref"
)

# eegKcomplexDetect() returns a data.frame of K-complex events
# (channel, negative/positive peak samples and amplitudes, duration_ms)
head(kcomplexes)

# Count per channel (channel index; C3/C4 are indices 1/2)
table(kcomplexes$channel)

# K-complex density (events per minute of recording)
nrow(kcomplexes) / rec_minutes

# Detect slow waves (maximum amplitude is over frontal regions, F3/F4)
slow_waves <- eegSlowWaveDetect(
  eeg_filt,
  freq_range = c(0.5, 2.0),   # Delta band slow waves
  min_amplitude = 75,         # μV (AASM criterion)
  assay_name = "reref"
)

# eegSlowWaveDetect() returns a data.frame of slow-wave events
# (channel, start/end sample, negative_peak, positive_peak, duration_ms, slope)
head(slow_waves)

# Slow wave density (waves per minute of recording)
slow_waves_density <- nrow(slow_waves) / rec_minutes
slow_waves_density

# SWA: mean negative-peak amplitude
mean_amplitude <- mean(abs(slow_waves$negative_peak))
mean_amplitude

# Slow wave slope (marker of sleep depth)
mean(slow_waves$slope)

# Full-night SWA dynamics (an exponential decline across sleep cycles) require
# a complete overnight recording. Here we summarise slow-wave amplitude over
# time by binning slow-wave events into successive segments of the recording.
if (nrow(slow_waves) > 0) {
  sw_time_min <- slow_waves$start_sample / samplingRate(eeg_filt) / 60
  n_bins <- 5
  bin <- cut(sw_time_min, breaks = n_bins, labels = FALSE)
  swa_by_bin <- tapply(abs(slow_waves$negative_peak), bin, mean)

  plot(seq_along(swa_by_bin), swa_by_bin,
       type = "b", pch = 19,
       xlab = "Recording segment", ylab = "Mean slow-wave amplitude (uV)",
       main = "Slow Wave Amplitude Across the Recording")
}

# 1. Load and preprocess data (short excerpt for a fast demo)
eeg <- make_eeg_sleep(n_time = 20 * 60 * 256, n_channels = 6, sr = 256)
eeg <- eegFilter(eeg, lowcut = 0.3, highcut = 35)
eeg <- eegRereference(eeg, ref_type = "average")

# 2. Automatic sleep staging
stages <- eegSleepStage(eeg, epoch_sec = 30)
metadata(eeg)$sleep_stages <- stages

# 3. Detect sleep-specific waveforms (each returns an event data.frame)
spindles <- eegSpindleDetect(eeg, freq_range = c(11, 16))
kcomplexes <- eegKcomplexDetect(eeg, min_neg_amplitude = 75)
slow_waves <- eegSlowWaveDetect(eeg, freq_range = c(0.5, 2.0), min_amplitude = 75)

# 4. Compute sleep architecture metrics (reads metadata$sleep_stages)
metrics <- eegSleepMetrics(eeg)

# 5. Generate visualization
eegPlotHypnogram(eeg, stages = stages)

# 6. Export comprehensive report
report <- list(
  demographics = list(age = 35, sex = "M"),
  sleep_metrics = metrics,
  spindle_count = nrow(spindles),
  kcomplex_count = nrow(kcomplexes),
  swa_mean = mean(abs(slow_waves$negative_peak)),
  n_epochs_scored = nrow(stages)
)

# Save report for clinical review
saveRDS(report, file.path(tempdir(), "sleep_study_report.rds"))

sessionInfo()