EEGToolkit.jl

A struct representing time series data.

Fields

• x::Vector{<:AbstractFloat}: Time series data.
• fs::Integer: Sampling rate.

segment(v::Vector{T}, L::Int; overlap::Union{<:AbstractFloat,Integer}=0, symmetric=false) where {T}

Splits a vector v into segments of length L with an overlap overlap expressed as a fraction of L. The overlap defaults to 0 (no overlap). Returns a vector $v$ of vectors - i.e. Vector{Vector{T}} - with $\vec{v_i}$ the $i$th segment in the split.

The function always attempts to capture the whole vector, even if the final split is not of length L. For example,

> x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]
> segment(x, 5)
2-element Vector{Vector{Int64}}:
[1, 2, 3, 4, 5]
[6, 7, 8, 9, 0]

> segment(x, 7)
2-element Vector{Vector{Int64}}:
[1, 2, 3, 4, 5, 6, 7]
[8, 9, 0]

Set symmetric=true to ensure that, if this occurs, the last split is dropped.

> segment(x, 3; symmetric=true)
3-element Vector{Vector{Int64}}:
[1, 2, 3]
[4, 5, 6]
[7, 8, 9]

If L is equal to the segment length, segment raises a warning and returns a vector with only the original vector: [v]. The return value ensures type-safety but the warning is raised because splitting a vector over its length is potentially a programming mistake.

segment(ts::TimeSeries, L::Int; kargs...)

Wrapper to segment the vector ts.x in the time series ts.

epoch(signal::TimeSeries, n::Integer; epoch_length::Integer=30)

Returns a vector [x₁, …, xₖ] with all values xᵢ corresponding to the nth epoch in the signal.

epoch(signal::TimeSeries, n::Integer, m::Integer)

Returns a vector [x₁, …, xₖ] with all indexes corresponding to epochs n, n+1, …, m of the EEG. The default sampling rate is used to compute the indexes.

plot_ts(ts::TimeSeries, s::Integer, e::Integer; norm=false, ylab="Amplitude (uV)")

Plots TimeSeries from epoch s to epoch e. The series many be normalized.

plot_ts(ts::TimeSeries, s::Integer; kargs...)

Plots TimeSeries at epoch s.

plot_ts(ts::TimeSeries; norm=false, ylab="Amplitude (uV)")

Plots TimeSeries. The series may be normalized.

seconds_to_time(seconds::Union{AbstractFloat, Integer})

Helper function: maps a time in seconds to a Time object.

A struct for the EEG data type. An EEG is simply conceived as a collection of labeled time series.

Fields

• signals::Dict{String, TimeSeries}: A dictionary mapping signal labels (strings) to arrays of floating-point values.

Constructors

EEG(file::String; id::String=""): Instantiates an EEG from an EDF file (file).

Example

eeg_data = EEG("path/to/edf_data/data.edf")

remove_channel!(eeg::EEG, channel::String)

Removes a channel from the EEG.

remove_channel!(eeg::EEG, channels::Vector{String})

Removes a list of channels from the EEG.

plot_eeg(eeg::EEG, s::Integer, e::Integer; channels::Vector{String}=[""], spacing::AbstractFloat=1.5)

Plots EEG channels from epoch s to epoch e. Specific channels may be selected with the channels karg. The spacing argument is an added factor in the normalization of the EEG signals - the vertical distance between each signal in the plot grows proportionally to spacing.

function artifactreject(signal::TimeSeries, anomdict::Dict{Int, Vector{Int}}; epochlength::Integer=30, subepochlength::Integer=5)::Vector{Vector{<:AbstractFloat}};

This function removes from a signal the sub-epochs which contain artifacts. It requires a TimeSeries and a Dict{Int, Vector{Int}}, hereby termed anom_dict (for anomaly dictionary).

anom_dict is understood to be such that anom_dict[i] = [n₁, …, nₖ] means the ith epoch has artifacts at sub-epochs n₁, ..., nₖ.

The return value is a segmented signal (Vector{Vector<:AbstractFloat}}), each of whose segments corresponds to an epoch with its artifcat-contaminated sub-epochs removed. In other words, if the result holds the return value of this function, result[i] contains what is left from the ith epoch after removing its contaminated sub-epochs. It is possible that result[i] is empty, if all sub-epochs of epoch i contained artifacts.

artifact_reject(signal::TimeSeries, anoms::Vector{Integer}; epoch_length::Integer=30)

An anomaly vector $\vec{x} \in \mathbb{N}^{n}$ is a sorted vector whose values are those epochs in an TimeSeries that contain anomalies or artifacts. This function segments the TimeSeries and filters out all epochs containing artifacts.

function nrem(staging::Staging, n::Integer=30, m::Integer=10)

Finds the k underlying NREM periods in a staging vector. Returns a vector of vectors V s.t. the ith vector in V contains the epochs which comprise the ith NREM period. Thus, the length of V is k the number of NREM periods.

The n parameter is the number of (not necessarily consecutive) epochs which comprise a NREM period. The m parameter is the number of REM or wakefulness epochs required to mark the end of a NREM period (if n NREM epochs were parsed before) or to disregard the current sequence and begin parsing from the next NREM epoch.

The staging parameter is a vector containing only the symbols 1, …, 6, ? where 5 denotes REM, 6 denotes wakefulness, and ? denotes unscored/unstaged.

sigma_index(x::Vector{<:AbstractFloat}, fs::Integer)

The $\sigma$-index algorithm (Huupponen et al., 2007) find abnormally high amplitude values in the spindle frequency band. Per each 1 second window of the EEG, it computes

• the maximum amplitude in the spindle frequency band, which we call $S_{max}$
• the average amplitude in the low alpha and theta frequencies, which we call

$\alpha_{mean}, \theta_{mean}$

• the maximum alpha amplitude $\alpha_{max}$

The $\sigma$-index of each window is defined to be zero if $\alpha_max > S_{max}$, and otherwise

\[f(S_{max}, \alpha_{mean}, \phi_{mean}) = \frac{2S_{max}}{ \alpha_{mean} + \theta_{mean} } \]

Higher values are indicative of a higher spindle probability. The rejection threshold recommended in the original paper is $\lambda = 4.5$.

relative_spindle_power(x::Vector{<:AbstractFloat}, fs::Integer)

The Relative Spindle Power (RSP) algorithm (Devuyst et al., 2011) also detects abnormal values along the spindle frequency band. For every 1 second window, the amplitude spectrum $S(t)$ is computed, and the RSP is defined as

\[RSP(t) = \frac{\int_{11}^{16} S(t, f) df}{\int_{0.5}^{40} S(t, f) df}\]

This definition is more intelligible than the that of the $\sigma$-index, insofar as it represents the ratio of the total power in the spindle band with respect to the total power in the $[0.5, 40]$ frequency range. It is evident by definition that $0 \leq RSP \leq 1$. Higher values are indicative of a higher spindle probability (it should be clear that $RSP$ is not a probability itself). The rejection threshold recommended in the original paper is $\lambda = 0.22$.

Structure for amplitude spectrum estimations. Estimations are by default one sided, with frequencies ranging from [0, fₛ/2]. The formula used is

\[\frac{2|H(f)|}{\sum_i w_i}\]

with $w_i$ a Hanning window.

Fields

• freq::Vector{<:AbstractFloat}: Frequency range of the spectrum
• spectrum::Vector{<:AbstractFloat}: Estimated spectral amplitude

Constructors

AmplitudeSpectrum(x::Vector{<:AbstractFloat}, sampling_rate::Integer, pad::Integer) : Computes a direct PSD over a signal x with a given sampling_rate.

Structure for PSD estimations. Estimations are by default one sided, with frequencies ranging from [0, fₛ/2].

The default formula is

\[\frac{2|H(f)|^2}{\zeta \sum_i w_i^2}\]

with $w_1, \ldots, w_n$ a Hanning window and $\zeta$ a normalization factor which defaults to $1$.

Barlett or Welch's mehtod can be used, where the formula becomes

\[\frac{1}{M K \varphi} \sum_i^M \left[ \frac{2|H_i(f)|^2}{ \sum_i w_i^2} \right]\]

where $w_1, \ldots, w_n$ a Hanning window, $M$ the number of segments, $K$ the number of samples per segment, $H_i(f)$ the FFT of the $i$th segment of the signal, and $\varphi$ a normalization factor defaulting to 1.

Fields

• freq::Vector{<:AbstractFloat}: Frequency range of the spectrum
• spectrum::Vector{<:AbstractFloat} : Estimated spectral density in dB.

Constructors

• PSD(x::Vector{<:AbstractFloat}, fs::Integer; window_function::Function = hanning, pad::Integer=0, normalization::Real=1): Computes PSD estimation of a signal x with sampling rate fs. A window_function is applied to the signal, defaulting to a Hanning window. The signal may be padded to an optional length pad (defaults to zero, i.e. no padding). A Real normalization is added to the denominator, which defaults to 1.
• PSD(segs::Vector{Vector{T}}, fs::Integer; window_function=hanning, normalization::Real=1) where {T<:AbstractFloat}: Computes the average spectrum of the segment vectors segs. The estimation is normalized with a normalization that defaults to $1$. The window_function is applied to all segments prior to their PSD estimation.

PSD(x::Vector{<:AbstractFloat}, fs::Integer, seg_length::Integer; overlap::Union{<:AbstractFloat,Integer} = 0.5, window_function::Function=hanning, normalization::Real=1): Segments the signal x into segments of length seg_length, with an overlap of overlap which defaults to 0.5 (half the segment length). After signal segmentation, the average spectrum of the resulting segments is returned. A window_function is applied to all segments prior to their PSD estimation, defaulting to a Hanning window. The final estimation is normalized with a normalization factor that defaults to 1.

• PSD(ts::TimeSeries; kargs...) : Wrapper to apply the first or second constructor to a TimeSeries signal.
• PSD(ts::TimeSeries, seg_length::Integer; kargs...): Wrapper to apply the third constructor to a TimeSeries signal.
• PSD(freq::Vector{<:AbstractFloat}, spectrum::Vector{<:AbstractFloat}): Direct constructor.

plot_psd(psd::PSD; freq_lim=30.0)

Plot a PSD with x-axis being frequency and y-axis being estimated power spectrum.

A spectrogram is a matrix $S^{M \times F}$ where $M$ is the number of windows in the windowing of a signal and $F$ is the length of the spectrum vector in any given window (i.e. the frequency resolution). It is useful to observe spectral changes in time or to compute the spectrum of time-regions of interest (e.g. only NREM periods in a sleep EEG). The information available in direct PSD can be inferred from the spectrogram with ease.

For instance, let $f_1, f_2, \ldots, f_k$ be a strictly increasing sequence of frequencies. Assume these frequencies correspond to the column indexes $c_1, c_2, \ldots, c_k$ of $S$. Then the mean power in the frequency range $[f_1, f_k]$ is

\[\frac{1}{M} \sum_{i=1}^{M}\left[\frac{1}{c_k - c_1}\sum_{j=c_1}^{c_k} S_{ij}\right] = \frac{1}{M\big(c_k - c_1\big)}\sum_{i=1}^{M}\sum_{j=c_1}^{c_k} S_{ij}\]

In this package, mean power in a frequency range is computed with the mean_band_power function.

Fields

• time::Vector : Time domain
• freq::Vector{<:AbstractFloat}: Frequency domain
• spectrums::Matrix{<:AbstractFloat}: Power spectrum. Rows are time and columns are frequency; the value in spectrums[row, freq] is the power at time window row for frequency freq.
• segment_length::Integer : Length of each segment in time.
• aggregated_spectra::Vector{<:AbstractFloat} : Spectral average (mean of rows)

Constructors

• Spectrogram(segs::Vector{Vector{T}}, psd_function::Function; dB = false) where {T<:AbstractFloat}: Given a sequence of windows $w_1, \ldots, w_k$ contained in the segs argument, computes the PSD within each window using a custom psd_function.
• Spectrogram(signal::Vector{<:AbstractFloat}, window_length::Integer, psd_function::Function; overlap::Union{AbstractFloat, Integer}=0, dB=false): Splits a signal into (potentially overlapping) segments of length window_length and computes the Spectrogram over this windowing using the first constructor. A custom psd_function is used within each window. Symmetry is enforced over the split signal, meaning that if the last segment is of length not equal to the rest, it is dropped. Thus, all windows are of equal length.
• function Spectrogram(ts::TimeSeries, window_length::Integer, psd_function::Function; kargs...): Wrapper constructor for a TimeSeries object.

plot_spectrogram(spec::Spectrogram; freq_lim::AbstractFloat=30.0, type::Int=1, color=:nipy_spectral)

Plots a spectogram spec either in 2d (type = 1) or 3d (type = 2). An optional frequency limit (freq_lim) may be set (defaults to 30Hz). The color palette color may be set; defaults to nipy_spectral.

freq_band(spec::Union{PSD}, lower::AbstractFloat, upper::AbstractFloat)

Given a PSD, returns a Vector{<:AbstractFloat} with the powers within the frequency band [lower, upper].

freq_band(spec::Spectrogram, lower::AbstractFloat, upper::AbstractFloat, window::Integer)

Given a Spectrogram, returns a Vector{<:AbstractFloat} with the powers within a frequency band [lower, upper] of a specific window (row of the spectrogram).

freq_band(spec::Spectrogram, lower::AbstractFloat, upper::AbstractFloat)

Given a Spectrogram, returns a Matrix{<:AbstractFloat} with the powers within a frequency band [lower, upper] across all time windows.

mean_band_power(spec::Spectrogram, lower::AbstractFloat, upper::AbstractFloat)

Given a Spectrogram, returns the mean power in a given frequency band [lower, upper]. This function effectively computes

\[\frac{1}{M\big(c_k - c_1\big)}\sum_{i=1}^{M}\sum_{j=c_1}^{c_k} S_{ij}\]

mean_band_power(spec::PSD, lower::AbstractFloat, upper::AbstractFloat)

Given a PSD, returns the mean power in a given frequency band [lower, upper].

total_band_power(psd::PSD, lower::AbstractFloat, upper::AbstractFloat)

Given a PSD, computes the total power in the frequency band [lower, upper].

analyze_eeg(signal::Vector{<:AbstractFloat}, fs::Integer)::Spectrogram

Perform a standardized analysis of an EEG signal. This analysis procedure succesfully replicated results from Washington State University in collaboration with the developer's laboratory at UPenn.

The standardized procedure is as follows: split the signal into 30-sec epochs, each of which is split into 5-sec sub-epochs. Each epoch's spectrum is the aggregated spectra from its sub-epochs; the signal's spectrum is the aggregated spectra from its epochs. Rectangular window is used. Overlap is set to zero. No normalization is performed.

analyze_eeg(signal::Vector{<:AbstractFloat}, fs::Integer, epoch_indexes::Vector{<:Integer})::Spectrogram

Perform a standardized analysis of the specified epochs of an EEG signal. This analysis procedure succesfully replicated results from Washington State University in collaboration with the developer's laboratory at UPenn.

The standardized procedure is as follows: split the signal into 30-sec epochs, each of which is split into 5-sec sub-epochs. Each epoch's spectrum is the aggregated spectra from its sub-epochs; the signal's spectrum is the aggregated spectra from its epochs.

function detect_artifacts(eeg::EEG, channel_name::String, seg_length::Int)::Nothing

Performs artifact detection in the eeg channel named channel_name. Stores resulting vector of Artifact objects in the _artifacts dictionary of the eeg with key channel_name.

Artifact detection is performed on each segment of the signal segmented with seg_length. On each segment, the CAPA algorithm (Fisch et. al 2022) is used to detect epidemic distributional changes in the mean value of the segment. The penalty is an integer value β such that β ln(n) (with n the segment's length) is the penalty used by the CAPA algorithm to penalize the introduction of artifacts. The β penalty defaults to 24, which is much higher than the value recommended in Fisch et. al but matched human supervision on 78Hz sleep EEGs at the developer's laboratory.

function plot_artifacts_in_epochs(signal::TimeSeries, artifacts::ArtifactData, from::Integer, to::Integer; annotate::Bool=false)

Given a signal and a non-empty vector of artifacts, plots the signal and highlights the existing artifacts from epoch from to epoch to. If annotate is set to true, artifacts are annotated with their mean change and test statistic.

function plotartifactsinepochs(eeg::EEG, channelname::String,

Given an eeg and an artifact-detected channel_name, plots the existing artifacts in the channel from epoch from to epoch to. If annotate is set to true, artifacts are annotated with their mean change and test statistic.

resample(signal::TimeSeries, new_fs::Float64)::TimeSeries

Resamples the input signal to the given new_fs, using rational factor resampling (L/M). Returns a new TimeSeries with the resampled signal and updated sampling rate.

Given an integer $n$, finds the least $m = 2^k$ s.t. $m \geq n$.

zero_pad(v::Vector{T}, desired_length::Integer) where {T<:AbstractFloat}

Zero-pads a numeric vector v to a desired_length