function filtfilt(  X::Matrix, 
                        sr::Int, 
                        responseType::DSP.FilterType; 
        designMethod::DSP.ZeroPoleGain=Butterworth(2))

Apply a digital filter in a forward-backward manner to obtain a linear phase response.

Arguments

• X: the $T×N$ EEG recording, where $T$ and $N$ denotes the number of samples and channels (sensors), respectively
• sr: the sampling rate of X
• responseType: a filter response type of the DSP.jl package: Lowpass, Highpass, Bandpass, or Bandstop
• designMethod: The filter design method of the DSP.jl package: Butterworth, Chebyshev1, Chebyshev2, Elliptic, or FIRWindow.

By default, designMethod is Butterworth(2), that is, a second-order Butterworth filter.

For the analysis in the frequency domain, see the FourierAnalysis package.

Examples

using Eegle

sr = 128 # sampling rate
X = randn(sr*8, 19)

filteredX = filtfilt(X, sr, Bandpass(1, 24))

filteredX = filtfilt(X, sr, Bandstop(49, 51))

filteredX = filtfilt(X, sr, Highpass(10))

filteredX = filtfilt(X, sr, Lowpass(10))

filteredX = filtfilt(X, sr, Bandstop(49, 51); 
                    designMethod = Chebyshev1(4, 0.5))

# Apply a filter bank
responses = [Bandpass(1, 4), Bandpass(4, 8), Bandpass(8, 12), Bandpass(12, 16)]
filterBank = [filtfilt(X, 128, r; designMethod = Butterworth(8)) for r ∈ responses]

    function centeringMatrix(N::Int)

The common average reference (CAR) operator for referencing EEG data potentials so that their mean across sensors (space) is zero at all samples.

Let $X$ be the $T×N$ EEG recording, where $T$ and $N$ denotes the number of samples and channels (sensors), respectively, and let $H_N$ be the $N×N$ recentering matrix, then

$Y=XH$

is the CAR (or centered) data.

$H_N$ is named the common average reference operator. It is given at p.67 by (Searle, 1982), as

$H_N = I_N - \frac{1}{N} \left( \mathbf{1}_N \mathbf{1}_N^\top \right)$

where $I_N$ is the N-dimensional identity matrix and $\mathbf{1}_N$ is the $N$-dimensional vector of ones.

Alias ℌ (U+0210C, with escape sequence "frakH")

Return the $N×N$ centering matrix.

Examples

using Eegle

X = randn(128, 19)

# CAR
X_car = X * centeringMatrix(size(X, 2))
# or
X_car = X * ℌ(size(X, 2))

# double centered data: zero mean across time and space
X_dc = ℌ(size(X, 1)) * X * ℌ(size(X, 2))

    function globalFieldPower(X::AbstractMatrix{T}; 
        func::Function = identity) 
    where T<:Real 

The global field power (GFP) is the sample-by-sample total EEG power.

Let $X$ be the $T×N$ EEG recording, where $T$ and $N$ denotes the number of samples and channels, respectively, and $x_t$ be the vector of $N$ potentials at sample $t∈{1,...,T}$, then the GFP at each sample is given by

$x_t^Tx_t$.

Function func can be applied element-wise to the output (none by default). Anonymous functions can be used.

Usually the GFP is computed on common average reference data — see centeringMatrix.

Return the vector comprising the $T$ GFP values.

See also globalFieldRMS

Examples

using Eegle

X=randn(128, 19)
# ℌ is an alias for centeringMatrix

# using an anonymous function
g = globalFieldPower(X * ℌ(size(X, 2)); func=x->sqrt(x/size(X, 2)))

    function globalFieldRMS(X::AbstractMatrix{T}; 
        func=identity) 
    where T<:Real

The global field root mean square (GFRMS) is the square root of the globalFieldPower once this has been divided by the number of electrodes.

Let $X$ be the $T×N$ EEG recording, where $T$ and $N$ denotes the number of samples and channels, respectively, and let $x_t$ be the vector of $N$ potentials at sample $t∈{1,...T}$, then the GFRMS at each sample is given by

$\sqrt{\frac{1}{N} (x_t^\top x_t)}$.

Function func can be applied element-wise to the output (none by default). Anonymous functions can be used.

Usually the GFRMS is computed on common average reference data — see centeringMatrix.

Return the vector comprising the $T$ GFRMS values.

See also globalFieldPower

Examples

using Eegle

X=randn(128, 19)

g = globalFieldRMS(X * ℌ(size(X, 2)))
# ℌ is an alias for centeringMatrix

# using an anonymous function
g = globalFieldRMSX * ℌ(size(X, 2)); func=x->x^2)

    function epoching(X::AbstractMatrix{T}, sr;
        wl::Int = round(Int, sr*1.5),
        slide::Int = 0,
        minSize::S = round(Int, sr*1.5),
        lowPass::Union{T, S} = 14,
        richReturn::Bool = false) 
    where {T<:Real, S<:Int}

Segment an EEG file in successive epochs and compute a vector of unit ranges delimiting the epochs. This is used to extract epochs from spontaneous EEG recording. For tagged data (e.g., ERPs and BCI data), use Eegle.ERPs.trials instead.

Two segmentation methods are possible, the standard fixed-length epoching and the adaptive epoching based on the local minima of the low-pass filtered global field root mean square (GFRMS).

• Standard (default): wl is set to a positive integer (by default is sr*1.5), which determines the length in samples of the epochs. A positive value of slide (default=0) determines the number of overlapping samples. By default there will be no overlapping.
• Adaptive: if wl= 0, the GFRMS is computed, low-pass filtered using lowPass (in Hz) as the cut-off (default = 14 Hz) and segmented ensuring that the minimum epoch size (in samples) is minSize, which default is the nuber of samples covering 1.5s. Set LowPass to sr/2 (Nyquist frequency) for no low-pass filtering of the GFRMS.

Return

• Standard: if richReturn=false (default) $r$, else the 3-tuple ($r$, 0 and wl)
• Adaptive: if richReturn=false (default) $r$, else the 3-tuple ($r$, $m$, $l$),

where $r$ is the computed vector of unit ranges (a Vector{UnitRange{Int64}} type), $m$ the vector with the low-pass filtered GFMRS and $l$ the vector of epoch lengths.

Examples

using Eegle

X=randn(6144, 19)
sr = 128

# standard 1s epoching with 50% overlap
ranges = epoching(X, sr;
        wl = sr,
        slide = sr ÷ 2)
# return (1:64, 65:128, ...)

# standard 4s epoching with no overlap
ranges = epoching(X, sr;
        wl = sr * 4)
# return (1:512, 513:1024, ...)

# adaptive epoching of θ (4Hz-7.5Hz) oscillations
Xθ = filtfilt(X, sr, Bandpass(4, 7.5))
ranges = epoching(Xθ, sr;
        minSize = round(Int, sr ÷ 4), # at least one θ cycle
        lowPass = 7.5)  # ignore minima due to higher frequencies

# Get the epochs from any of the above:
𝐗 = [X[r, :] for r ∈ ranges] # or 𝐗 = [Xθ[r, :] for r ∈ ranges]

# Get the covariance matrices of the epochs from any of the above
𝐂 = covmat(𝐗) # See CovarianceMatrices.jl

# If only the covariance matrices are needed,
# a more memory-efficient way avoiding the extraction of 𝐗 is
𝐂 = ℍVector(covmat([view(X, r, :) for r ∈ ranges]))
𝐂θ = ℍVector(covmat([view(Xθ, r, :) for r ∈ ranges]))

See Eegle.BCI.covmat