(1) function mean(  X::Matrix{T}, 
                    wl::S, 
                    mark::Vector{Vector{S}};
        overlapping :: Bool = false,
        offset :: S = 0,
        weights :: Union{Vector{Vector{R}}, Symbol}=:none) 
    where {T<:Real, S<:Int}

(2) function mean(  o::EEG; 
        overlapping :: Bool = false,
        offset :: S = 0,
        weights :: Union{Vector{Vector{T}}, Symbol} = :none,
        mark :: Union{Vector{Vector{S}}, Nothing} = nothing) 
    where {T<:Real, S<:Int}

Estimate the weighted mean ERPs (event-related potentials), as the standard arithmetic mean (default) or using the multivariate regression method, as detailed in (Congedo et al., 2016).

Tutorials

xxx

METHOD (1)

Arguments

• X: the whole EEG recording, a matrix of size $T×N$, where $T$ and $N$ denotes the number of samples and channels (sensors), respectively
• wl: the window (trial or ERP) length in samples
• mark: the marker vectors.

Optional Keyword Arguments

• overlapping: see overlapping
• offset: see offset
• weights: can be used to obtain weighted means. By default, equal weights are used. It can be:
  • a vector of vectors of non-negative real weights for the trials, with the same shape as mark, where the empty vectors of mark are ignored
  • :a : adaptive weights computed as the inverse of the squared Frobenius norm of the trials data, along the lines of (Congedo et al., 2016).

Return

A vector of mean ERPs, one for each non-empty vectors in mark. Each mean is a matrix of size $wl×n$, where $n$ is the number of electrodes.

METHOD (2)

The same as method (1), but taking as input an Eegle.InOut.EEG structure o, which has fields providing the recording (X), the ERP duration in samples (wl) and the markers (mark).

Different markers can be used instead by passing marker vectors with the mark kwarg.

See Eegle.ERPs.stim2mark, Eegle.ERPs.trialsWeights

Examples

using Eegle # or using Eegle.ERPs

# Method (1)

# Number of channels, sampling rate and window length
N, sr, wl = 19, 128, 128

# number of tags per class
nm = [30, 45, 28]

X = randn(sr*100, N)
mark = [[rand(1:sr*100-wl) for i=1:m] for m∈nm]

# compute the arithmetic means for all classes with adaptive weights
𝐌 = mean(X, wl, mark; weights=:a)

# compute the same means for class 1 and 3 only
𝐌 = mean(X, wl, [mark[1], mark[3]]; weights=:a

# compute the same mean only for the first class
# and return it as a matrix (not as a vector of matrices)
M = mean(X, wl, [mark[1]]; weights=:a)[1]

# compute the multivariate regression means with adaptive weighting
𝐌 = mean(X, wl, mark; overlapping=true, weights=:a)

# Method (2)

# read the example file for the P300 BCI paradigm
o = readNY(EXAMPLE_P300_1)

# compute means (adaptive weights and multivariate regression)
𝐌 = mean(o; overlapping=true, weights=:a)

# target average ERP
T_ERP = 𝐌[findfirst(isequal("target"), o.clabels)] 

# nontarget average ERP
NT_ERP = 𝐌[findfirst(isequal("nontarget"), o.clabels)] 

    function stim2mark( stim::Vector{S}, 
                        wl::S;
        offset::S=0, 
        code=nothing) 
    where S <: Int

Convert a stimulation vector into marker vectors.

Arguments

• stim: the stimulation vector to be converted
• wl: the window (trial or ERP) length in samples.

Optional Keyword Arguments

• code: by default, the output will hold as many marker vectors as the largest tag (integers) in stim, which may or may not hold instances of all integers up to the largest. If there are missing integers, the corresponding marker vector will be empty. Alternatively, a vector of tags coding the classes of stimulations in stim can be passed as kwarg code. In this case, arbitrary non-zero tags can be used (even negative) and the number of marker vectors will be equal to the number of unique integers in code. If code is provided, the marker vectors are arranged in the order given there, otherwise the first vector corresponds to the tag 1, the second to tag 2, etc. In any case, in each vector, the samples are sorted in ascending order.

Return

A vector of $z$ marker vectors, where $z$ is the number of classes, i.e., the highest integer in stim or the number of non-zero elements in code if it is provided.

See mark2stim

Examples

using Eegle # or using Eegle.ERPs

sr, wl = 128, 256 # sampling rate, window length of trials
ns = sr*100 # number of samples of the recording

# simulate a valid stimulations vector for three classes
stim = vcat([rand()<0.01 ? rand(1:3) : 0 for i = 1:ns-wl], zeros(Int, wl))

mark = stim2mark(stim, wl)

stim2 = mark2stim(mark, ns) # is identical to stim

stim == stim2

    function mark2stim( mark::Vector{Vector{S}}, 
                        ns::S;
        offset::S=0, code=nothing) 
    where S <: Int

Reverse transformation of stim2mark.

If code is provided, it must not contain 0.

Examples see stim2mark

    function merge( mark::Vector{Vector{S}}, 
                    mergeClasses::Vector{Vector{S}})
    where S <: Int

Merge the vectors of marker vectors mark and sort the markers within each class. Return another marker vectors. The merging pattern is determined by mergeClasses.

As an example, suppose mark holds 4 vectors of markers and mergeClasses=[[1, 2], [3, 4]], then the result will hold two markers vectors, vectors 1 and 2 of mark concatenated and sorted and vectors 3 and 4 of in mark concatenated and sorted. Empty mark vectors will be ignored.

This can be used to merge classes in ERP and BCI experiments.

Examples

using Eegle # or using Eegle.ERPs
mark =  [   [128, 367], 
            [245, 765, 986],
            [467, 880, 1025, 1456],
            [728, 1230, 1330, 1550, 1980],  
        ]

merged = merge(mark, [[1, 2], [3, 4]])

# return: 2-element Vector{Vector{Int64}}:
#           [128, 245, 367, 765, 986]
#           [467, 728, 880, 1025, 1230, 1330, 1456, 1550, 1980]

    function trials(X::Matrix{R}, 
                    mark::Vector{Vector{S}}, 
                    wl::S;
        shape::Symbol = :cat
        weights::Union{Vector{R}, Nothing} = nothing,
        linComb::Union{Vector{R}, S, Nothing} = nothing,
        offset::S = 0) 
    where {R<:Real, S<:Int}

Extract trials of duration wl from a tagged EEG recording X. Optionally, multiply them by weights and compute a linear combination across sensors thereof.

Arguments

• X: the whole EEG recording, a matrix of size $T×N$, where $T$ is the number of samples and $N$ the number of channels (sensors)
• mark: a marker vectors
• wl: the window (trial, e.g., ERP) length in samples.

Optional Keyword Arguments

• shape: see below
• weights: optional weights to be multiplied to the trials. Adaptive weights can be obtained passing the Eegle.ERPs.trialsWeights function.

• linComb: Optional linear combination to be applied to the trials, e.g., a spatial filter. It can be:
  • an integer: extract for each (weighted) trial only the data at the electrode indexed by linComb $∈[1,..,n]$ (linear combination by a one-hot vector)
  • a vector $f$ of $N$ real elements: extract for each (weighted) trial the linear combination $X_jf$.
• offset: see offset.

Return

• if shape ≠ :cat: a vector of vectors of trials or of linear combinations thereof
• if shape == :cat (default): all trials or the linear combinations thereof concatenated in a single vector.

Each extracted trial is a $wl×N$ matrix if linComb is nothing (default), otherwise it a vector of $wl$ elements holding the linear combination.

Examples

using Eegle # or using Eegle.ERPs

# Example P300 BCI session 
o = readNY(EXAMPLE_P300_1)

# Extract 1-s trials starting at samples specified in `mark`
mark = [ [245, 658, 987], [258, 758, 1987]  ]

# since `shape=:cat` (default),
# `𝐗` will hold the six trials concatenated 
𝐗 = trials(o.X, mark, o.sr)

# extract only the time-series at electrodes Cz
𝐗 = trials(o.X, mark, o.sr; 
            linComb = findfirst(isequal("Cz"), o.sensors))

# `𝐗[1]`, `𝐗[2]`, ... are in this case vectors, not matrices

# suppose `f` is a spatial filter
f = randn(o.ne)

# extract the filtered trials
𝐗 = trials(o.X, mark, o.sr; linComb = f)
 
# `𝐗[1]`, `𝐗[2]`, ... are vectors in this case too

# extract the trials in separated vectors for each class
𝐗 = trials(o.X, mark, o.sr; shape=:byClass)

    function trialsWeights( X::Matrix{R}, 
                            mark::Vector{Vector{S}},
                            wl::S;
        M::Union{Matrix{R}, Nothing} = nothing,
        offset::S = 0) 
    where {R<:Real, S<:Int}

Compute adaptive weights for trials as the inverse of their squared Frobenius norm, along the lines of (Congedo et al., 2016). The method is unsupervised, i.e., agnostic to class labels, but a supervised version is available using the M arguments.

Arguments

• X: the whole EEG recording, a matrix of size $T×N$, where $T$ is the number of samples and $N$ the number of channels (sensors), respectively
• mark: the marker vectors. For empty vectors, an empty vector is returned
• wl: the window (trial or ERP) length in samples.

Optional Keyword Arguments

• M: if it is a vector holding z matrices of size $T×N$, where z is the number of vectors in mark (classes of stimulations), then the weights are computed as the inverse of the squared norms of $X_{j(i)}-M_i$ for all trials $X_{j(i)}$, $j \in \{1, \ldots, k\}$, $i \in \{i, \ldots, z\}$ for each class $i$ separately.
• offset: see offset.

Examples

using Eegle # or using Eegle.ERPs

# Example P300 BCI session 
o = readNY(EXAMPLE_P300_1)

weights = trialsWeights(o.X, o.mark, o.wl) 

# weights[1] and weights[2] are the weights for class
# "non-target" and "target", respectively.

# check that the mean of weights within each class is 1
mean(weights[1])
mean(weights[2])

# Supervised weights using the mean ERPs
M = mean(o.X, o.wl, o.mark; overlapping=true, weights=:a)
weights = trialsWeights(o.X, o.mark, o.wl; M) 
# for the keyword argument, in Julia `M`, it is the same as `M=M`

    function reject(X::Matrix{R}, 
                    stim::Vector{Int}, 
                    wl::S;
        offset::S = 0,
        upperLimit::Union{R, S} = 1.2,
        returnDetails::Bool=false) 
    where {R<:Real, S<:Int}

Automatic rejection of artefacted trials in tagged EEG data by automatic amplitude thresholding.

Arguments

• X: the whole EEG recording, a matrix of size $T×N$, where $T$ is the number of samples and $N$ the number of electrodes
• stim: a stimulation vector
• wl: the trial length, i.e., the ERPs or trial duration, in samples.

Optioanl Keyword Arguments

• offset: see offset
• upperLimit: modulate the definition of the upper threshold (see below). a reasonable value ∈[1, 1.6]
• upperLimit: determine the output (see below).

Description

Let $v$ be the natural logarithm of the field root mean square (FRMS, see Eegle.Processing.globalFieldRMS) of X sorted in ascending order.

The lower threshold $l$ is defined as the tenth value of $v$ (robust minimum estimator).

The upper threshold $h$ is defined as

$h=m+((m-l)u)$,

where:

• $m$ is the mean of the $2wl$ central values of $v$, taken as a robust central tendency estimator
• $u$ is kwarg upperlimit (default=1.2).

All trials in which at least one sample of the log-FRMS exceeds $h$ (there are samples with over-threshold amplitude) or in which $l$ exceeds the log-FRMS all over the trial (almost zero everywhere) are rejected.

Return

• if returnDetails is false (default), a 5-tuple holding the following objects:
  • the stimulation vector stim with the tags corresponding to rejected trials set to zero (accepted trials),
  • the stimulation vector stim with the tags corresponding to accepted trials set to zero (rejected trials),
  • the first object as marker vectors,
  • the second object as marker vectors,
  • the number of rejected trials per class as a vector of integers.
• if returnDetails is true, a 9-tuple holding the above 5 objects and, in addition:
  • the log-FMRS (not sorted),
  • the mean $m$,
  • the lower threshold $l$,
  • the upper threshold $h$.

Examples

using Eegle # or using Eegle.ERPs

cleanstim, rejecstim, cleanmark, rejecmark, rejected = reject(X, stim, wl; upperLimit=1.5)

R = reject(X, stim, wl; upperLimit=1.5, returnDetails = true) # R is a tuple of 9 objects

norm((R[1].+R[2]).-stim)==0 # should be true

Tutorials

Tutorial ML 1, Tutorial ML 2