Essential Motor Cortex Signal Processing: an ERP and functional connectivity MATLAB toolbox -- user guide version 2.0

TL;DR

This paper introduces version 2.0 of a MATLAB toolbox designed for comprehensive analysis of motor cortex signals, including ERP estimation, functional connectivity, and EMG quantification, aiding neuroscience research.

Contribution

It provides a free, open-source MATLAB toolbox implementing methods for ERP, cortical connectivity, and EMG analysis tailored for motor cortex studies.

Findings

Includes tools for ERP extraction and analysis.

Provides functions for cortical connectivity measurement.

Supports EMG quantification during motor tasks.

Abstract

The purpose of this document is to help individuals use the "Essential Motor Cortex Signal Processing MATLAB Toolbox". The toolbox implements various methods for three major aspects of investigating human motor cortex from Neuroscience view point: (1) ERP estimation and quantification, (2) Cortical Functional Connectivity analysis and (3) EMG quantification. The toolbox -- which is distributed under the terms of the GNU GENERAL PUBLIC LICENSE as a set of MATLAB R routines -- can be downloaded directly at the address: http://oset.ir/category.php?dir=Tools or from the public repository on GitHub, at address below: https://github.com/EsiSeraj/ERP Connectivity EMG Analysis The purpose of this toolbox is threefold: 1. Extract the event-related-potential (ERP) from preprocessed cerebral signals (i.e. EEG, MEG, etc.), identify and then quantify the event-related synchronization/desynchronization (ERS/ERD) events. Both time-course dynamics and time-frequency (TF) analyzes are included. 2. Measure, quantify and demonstrate the cortical functional connectivity (CFC) across scalp electrodes. These set of functions can also be applied to various types of cerebral signals (i.e. electric and magnetic). 3. Quantify electromyogram (EMG) recorded from active muscles during performing motor tasks.

Equations

$$PLV(f)=\left|\frac{1}{T}\sum _{t=1}^T\exp \left(j[{\varphi }_y(t,f)-{\varphi }_x(t,f)]\right)\right|$$

$$|MSC(f){|}^2=\frac{|PS{D}_{xy}(f){|}^2}{PS{D}_{xx}(f)PS{D}_{yy}(f)}=\frac{|{\sum }_{i=1}^N{X}_i(f)Y_i^{\ast }(f){|}^2}{{\sum }_{i=1}^N|X_i(f){|}^2{\sum }_{i=1}^N|Y_i(f){|}^2}$$

Full text

The toolbox – which is distributed under the terms of the GNU GENERAL PUBLIC LICENSE as a set of MATLAB® routines – can be downloaded directly at the address:

http://oset.ir/category.php?dir=Tools

.

or from the public repository on GitHub, at address below:

https://github.com/EsiSeraj/ERP_Connectivity_EMG_Analysis

.

The purpose of this toolbox is threefold:

Extract the event-related-potential (ERP) from preprocessed cerebral signals (i.e. EEG, MEG, etc.), identify and then quantify the event-related synchronization/desynchronization (ERS/ERD) events. Both time-course dynamics and time-frequency (TF) analyzes are included. 2. 2.

Measure, quantify and demonstrate the cortical functional connectivity (CFC) across scalp electrodes. These set of functions can also be applied to various types of cerebral signals (i.e. electric and magnetic). 3. 3.

Quantify electromyogram (EMG) recorded from active muscles during performing motor tasks.

Key-words: Event-Related Potential, ERP, Event-Related Synchronization, ERS, Event-Related Desynchronization, ERD, ERP Time Dynamics, Time-Frequency Analysis, Cortical Functional Connectivity, CFC, Electroencephalogram, EEG, Electromyogram, EMG, EMG Quantification, MATLAB Function, Free Toolbox, User Guide, Manual

1. 1 Introduction

2. 1.1 License - No Warranty

3. 1.1.1 The Open-Source Electrophysiological Toolbox (OSET)

4. 1.2 Citation

5. 1.3 Download and Utilization

6. 1.4 Getting Help

7. 2 User Guide

8. 2.1 Overview

9. 2.2 Fundamentals

10. 2.2.1 Data: Setup, Equipments and Format

11. 2.2.2 Preconditioning

12. 2.2.3 Time-course ERP Dynamics

13. 2.2.4 Time-frequency ERP Analysis

14. 2.2.5 Time-course Cortical Functional Connectivity

15. 2.2.6 Pair-wise Cortical Functional Connectivity

16. 2.2.7 EMG Quantification

17. 3 Reference Manual

18. 3.1 bdf2mat_main.m

19. 3.2 trigger_synch.m

20. 3.3 trigger_avg_erp.m

21. 3.4 trigger_avg_TF_erp.m

22. 3.5 erp_quantification.m

23. 3.6 TCPLV.m

24. 3.7 PWPLV.m

25. 3.8 PWCoherence.m

26. 3.9 emg_quantification.m

27. 3.10 drift_reject.m

28. 3.11 emg_onset.m

29. 3.12 phase_est.m

30. 3.13 PLV_PhaseSeq.m

31. 3.14 sig_trend.m

32. 3.15 task_separator.m

33. 4 Acknowledgment

34. 5 References

1 Introduction

This document is meant to help individuals use the ”Essential Motor Cortex Signal Processing MATLAB Toolbox”. Using the most popular reference articles in literature, the toolbox implements various methods for three major aspects of neuro-physiological investigation of human motor cortex: (1) ERP estimation and quantification (e.g. based on [1, 2, 3, 4]), (2) Cortical Functional Connectivity analysis (e.g. based on [5, 6, 7, 8, 9]) and (3) EMG quantification (e.g. based on [10, 11]).

The purpose of this toolbox is threefold:

Extract the event-related potential (ERP) from preprocessed cerebral signals (i.e. EEG, MEG, etc.), identify and then quantify the event-related synchronization/desynchronization (ERS/ ERD) events. Both time-course dynamics and time-frequency (TF) analysis are included. 2. 2.

Measure, quantify and demonstrate the cortical functional connectivity (CFC) across scalp electrodes. These set of functions can also be applied to various types of cerebral signals (i.e. electric and magnetic). 3. 3.

Quantify electromyogram (EMG) recorded from active muscles during performing motor tasks.

A primary goal of this toolset is to ease understanding the routines and help individuals alter our codes according to their study and/or implement their own techniques. For this purpose, herein, we first present a detailed tutorial (e.g. Section 2.2.) on signal preconditioning and interpretable-implementation of introduced methods, from a practical view-point.

1.1 License - No Warranty

This program is free software; you can redistribute it and/or modify it under the terms of the GNU GENERAL PUBLIC LICENSE as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU GENERAL PUBLIC LICENSE for more details. You should have received a copy of the GNU GENERAL PUBLIC LICENSE along with this program; if not, see $\langle$ http://www.gnu.org/licenses/ $\rangle$ or write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

1.1.1 The Open-Source Electrophysiological Toolbox (OSET)

Open Source Electrophysiological Toolbox (OSET) is a collection of electrophysiological data and open source codes for biosignal generation, modeling, processing, and filtering. OSET, version 3.1, 2014 Released under the GNU GENERAL PUBLIC LICENSE [12]. Copyright© 2012.

As a progressive general-purpose open-source toolset, OSET is one of the main sources to access the functions and codes of current motor cortex signal processing toolbox and also its documentation and dependencies.

1.2 Citation

Within the limits of the GNU GENERAL PUBLIC LICENSE, you can use the toolbox as you please; however, if you use the toolbox in a work of your own that you wish to publish, you need to make sure to cite this user manual and the original studies properly, as shown below. This way you will contribute to helping other scholars find these items.

• •

Esmaeil Seraj and Karthiga Mahalingam “Essential Motor Cortex Signal Processing: an ERP and functional connectivity MATLAB toolbox - User Guide Version 1.0,” arXiv Preprint, June 2019 [Online].

1.3 Download and Utilization

The latest version of the toolbox can be downloaded directly from OSET at the address:

http://oset.ir/category.php?dir=Tools

or from the public repository on GitHub, at:

https://github.com/EsiSeraj/ERP_Connectivity_EMG_Analysis

The functions and m-files can be downloaded separately as you need or all together in a compressed file named ERP_Connectivity_EMG_Analysis_Toolbox. Once you have downloaded and uncompressed the toolbox, 35 files represented in following Tables shall be appeared in your chosen directory. Table (1) represents all m-files (i.e. core functions, internal computational functions and demo test files) included in this toolset with their short descriptions and references. Additionally, a copy of both GNU general public license, this user manual and ten recorded sample data are also included (see Section 2.2.1 for details). Moreover, you might add the path of the directory in which you stored the toolbox to your MATLAB in order to easily use and apply the functions for your own dataset.

[FIGURE:]

1.4 Getting Help

If you have added the toolbox directory to the MATLAB® path you can simply type:

<doc   function_name>           or           <help   function_name>

in command window to get online help for the function you are using. Furthermore, you can also contact any of the authors111E. Seraj, K. Mahalingam: {eseraj3, kmahalingam}@gatech.edu directly to ask any related questions or discuss possible difficulties or errors you might encounter. Please feel free to contact in either case.

2 User Guide

2.1 Overview

This document is meant to help individuals use the ”Essential Motor Cortex Signal Processing MATLAB Toolbox” which implements various methods for three major aspects of investigating human motor cortex from Neuroscience view point: (1) ERP estimation and quantification (e.g. based on [1, 2, 3, 4]), (2) Cortical Functional Connectivity analysis (e.g. based on [5, 6, 7, 8, 9]) and (3) EMG quantification (e.g. based on [10, 11]). Measurements and methodologies are all derived based on the procedures suggested by the most popular reference articles in that category, as aforementioned.

The current version of toolbox covers six well-known and widely used approaches within Neuroscience community for analyzing motor cortex potentials and EMG signals, as follows:

• •

Preconditioning: Preprocessing the signals, including baseline drift removal, artifact rejection, movement onset detection (i.e. for non-cued movements), trigger synchronization (i.e. essential step for ERP extraction through non–synchronized trials of EEG data) and etc.

• •

Time-course ERP Dynamics: Estimation and quantification of ERP (i.e. ERS and ERD events).

• •

Time-frequency ERP Analysis: Estimation and representation of time-frequency (TF) maps through various widely-accepted TF analysis methods such as Short-time Fourier Transform (STFT), Continuous Wavelet Transform (CWT), Narrow-band Channelization (NBCH) and etc.

• •

Time-course Cortical Functional Connectivity: Estimation and quantification of local/large -scale functional connectivity between arbitrary pairs of electrodes through Phase Locking Value (PLV) [7, 9] and Magnitude Squared Coherence (MSC) [6].

• •

Pair-wise Cortical Functional Connectivity: Estimation and quantification of local/large -scale functional connectivity maps across scalp electrodes through PLV and MSC.

• •

EMG Quantification: Estimation and quantification of electromyogram (EMG) records of active muscles during performing motor tasks.

2.2 Fundamentals

2.2.1 Data: Setup, Equipments and Format

ActiveTwo is a very well-known commonly used high-resolution biosignal acquisition system that comes with advantageous capabilities such as active electrodes [15]. Active electrodes are those that does not require extensive skin (scalp) preparation due to the presence of pre-amplifiers in them. The pin electrodes must be connected to the scalp on one end and to the acquisition system on the other. The acquisition system is then connected to the PC to which it relays electrode data serially. Electrode gel has to be used on the scalp with electrodes to improve conductivity and to reduce noise. Care should be taken to reduce noise as much as possible by properly grounding all electrical appliances connected with the electrode system, making sure there are minimum head movements while recording etc.

Before starting data collection, all electrodes need to be checked for their proper functionality. This can be performed through the single bucket test [16]. Electrodes can be placed according to the 10-20 system [17], which is a very common approach for most studies. EEG cap has to be adjusted based on head circumference, nasion-inion, and ear-to-ear distance. Reference electrode areas are normally chosen to be the average between two ears [15] and thus CMS/DRL flat type electrodes should be placed behind earlobes. Finally, EMG electrodes should be placed directly on the active muscle the activity of which is to be recorded.

The EEG signals recorded by ActiveTwo are stored in *.bdf format files. We provided a very useful MATLAB m-file by Gleb Tcheslavski [14] that converts *.bdf format files to *.mat MATLAB data files.

2.2.2 Preconditioning

In ERP and functional connectivity analysis problems, preconditioning mostly refers to preprocessing the signals. This includes baseline drift removal, artifact rejection, movement onset detection from EMG signals, trigger synchronization (i.e. as described above) and etc.

Quick note on artifact rejection:

The concept of noise in this problem could be regarded to two different parts of recorded signals. First group are the common sources of noise observed in EEG recordings such as muscle activity, eye movement, electrocardiographic activity, instrumentation and equipment related artifacts and interference, slow drifts and amplifier saturation and etc. [2]. Moreover, background EEG activity is considered as the other source of noise [1]. In this case, ERP is considered as signal and the background EEG activity as noise [2, 1].

Artifact rejection in later case is normally easy, since as widely accepted, background ongoing EEG is considered as uncorrelated to signal (ERP) and thus, can be rejected by an ensemble averaging over trials [1, 2, 3]. The best way of dealing with the first group of noises however, is to not have any in first place. The easiest sources of noise to deal with are AC power lines, lighting and electronic equipment such as computers, displays and TVs, wireless routers, notebooks and also mobile phones. The basic step here is to simply remove any unnecessary sources of electro-magnetic (EM) noise from the recording room and, if possible, replace equipment using alternate current with direct current [18, 19]. Further, a significant number of noise sources can be rejected through different steps of low-pass, high-pass and sometimes notch filters.

Drift cancellation:

Despite having various sources of additional noise on EEG, one of the most important ones to deal with, in this application, is the baseline drift noise (a.k.a base line wander). The drift noise normally happens due to sweating, drifts of cap and/or electrodes and similar reasons, which could normally lead to amplifier saturation [2, 18] and thus incorrect and unreliable measurements. The effect of drifts on ERP estimation application is shown in Fig. 1 which is directly employed from [2]. As shown, although the effect of drift removal from 182 trials is relatively small, still is considerable for analyzing the relative decreases and increases in ERP power (ERD/ERS) quantification.

As suggested in [2], a number of issues are to be taken into consideration, such as detecting negative drifts as well as positive ones, using similar specs for the drift rejection method across trials and also choosing a proper-length (long enough) temporal window for drift calculation (preferably between 1 or two seconds) [2]. Nevertheless, the filtering specs in EEG artifact removal are highly dependent on the data and recording conditions and thus, should be chosen empirically [2, 20]. The strongest recommendation throughout the literature is to evaluate the filter before becoming committed to [2, 20]. An important consideration to be taken into account here is that the filtering procedure itself, depending on the characteristics and nature of the utilized filter and filtering method, can also significantly affect the data in the non-filtered frequency ranges and thus can affect the onset detection and/or amplitude of estimated ERP. Therefore, to prevent these problems, it is highly recommended that filtering should be limited to what is necessary and unavoidable [18].

In this toolbox, the function drift_reject.m is provided for this purpose which uses two stages of median or moving average temporal filtering to extract the base-line wander and reject it. To operate this function you can simply call function below with proper input parameters:

sig = drift_reject(raw_sig, L1, L2, approach)

where the inputs and outputs are described later in Section 3. Fig. 2 is an elaboration of how how this function works.

2.2.3 Time-course ERP Dynamics

For this part, the conventional approach, as proposed in [1], is implemented. Accordingly, the well-known five-step method presented below is employed for ERP extraction from EEG signals:

• •

Step #1 – Band-pass frequency filtering the EEG data in each trial

• •

Step #2 – Squaring the samples of filtered trials to find the power

• •

Step #3 – Ensemble averaging over synchronized trials according to trigger time

• •

Step #4 – Calculating the trend of output sequence in order to obtain ERP signal

• •

Step #5 – Quantifying the estimated ERP signal according to a reference period

The estimated signal in Step #4 is the desired ERP signal. Each step is elaborated below.

Step #1:

For this step, zero-phase bandpass frequency filtering using Cascaded Integrator-Comb (CIC) filter is implemented [21]. The filter performs forward-backward filtering successively. It has absolutely zero-phase for even filter order values and a phase-lag equal to a single stage CIC filter for odd filter order values. Therefore, for this application we recommend an even value to be chosen as the filter order to prevent any subsequent phase distortion.

An important issue here is to select the frequency band of interest for showing the time course of ERP [1]. This concern is very crucial due to the fact that different brain states can modulate different frequency bands of EEG rhythms. Additionally, even for the same action, one might observe intra-trial variability on dominant frequencies [1]. Accordingly, estimating the time-course of ERP signal related to the task for further analysis requires additional cautions to be taken into account.

Due to the subjective process of choosing a proper frequency band for ERP analysis, herein, the functions are designed to accept frequency ranges in a way that users can customize their frequency range of interest and experiment various bands of frequency easily. To read further in these regards, refer to [1, 2, 3].

Steps #2, #3 and #4:

Steps 2 through 4 are closely related (e.g. even could be considered as one step) where after squaring the samples of the filtered signal, two successive averaging steps are performed across trials and then time [1]. Ensemble averaging is the most important one among these steps since EEG signals are required to be averaged over synchronized data trials according to their respective trigger time. Based on this, first the onset of movement needs to be estimated (i.e. can be performed either through hardware settings (e.g. cue based recordings) or using our provided emg_onset.m function) as the trigger time. Then, the EEG trials have to be sorted and aligned with respect to the trigger time and eventually averaged across data trials. Afterwards, a simple temporal averaging is performed on obtained signal to extract the ERP sequence.

Interested individuals can refer to either the Reference Manual (i.e. Section3) or [1, 2, 3] for details of each averaging step and implemented signal trend estimation function.

Steps #5:

For this step, normally a 1-sec long reference period before the movement onset is chosen. Afterwards, the relative decrease intervals in power after movement onset are called ERD events and vice versa, the relative increase intervals in power are named ERS [22, 23].

According to discussions presented in [22, 23, 1], reference period is chosen some seconds before the task onset, completely based on application of interest and data and thus is empirically set. Reference value is calculated as the average of power samples within reference period. An axis according to the reference value is set to 100% and the corresponding crosses of ERP curve are named respectively, time step ”0”, start of ERD1/ERS1, end of ERD1/ERS1 and so on. Afterwards, standard deviation of the reference period is calculated and the corresponding confidence intervals are defined (mean plus/minus 2 or 3 times STD). Finally, the area enclosed by ERP curve and “mean - confidence interval” is calculated and normalized by the length of that segment as ERD area and vice versa; the area enclosed by ERP curve and “mean + confidence interval” is calculated and normalized by the length as ERS area.

In this toolbox, the function trigger_avg_erp.m is provided for this purpose which utilizes trigger_ synch.m function to synchronize cerebral signals according to their trigger time and BPFilter5.m function to perform the aforementioned CIC filtering. To operate this function you can use the demo m-file test_ERPanalysis_main.m or simply call function below with proper input parameters:

[erp, synch_eeg, synch_emg, trigger_time_sec, time_vec] = trigger_avg_erp(eeg, emg, fs, emg_onset_sampl, freq_band, duration)

where the inputs and outputs are described later in Section 3. Fig. 3 is an elaboration of how how this function works.

Moreover, function erp_quantification.m is provided to quantify the calculated ERP signal as the output of trigger_avg_erp.m according to the exact procedure described above. To operate this function you can simply call function below with proper input parameters:

[erd_area, ers_area, quant_erp] = erp_quantification(erp, fs, trigger_time_sec, ref_per, cof_intv)

where the inputs and outputs are described later in Section 3. Fig. 4 is an elaboration of how how this function works.

2.2.4 Time-frequency ERP Analysis

Time-frequency representations (TFRs) can monitor ERPs in a wide range of frequency components alongside with their temporal variations. Using TFRs, not only the risk of missing important information is lower, but also the interpretations provided can be more comprehensive and generic. Herein, as a part of ERP analysis package, a MATLAB function is dedicated to this purpose in which three different approaches for time/frequency analysis, namely Short-Time Fourier Transform (STFT), Continues Wavelet Transform (CWT) and Narrow-Band Channelization method (NBCH) as introduced in [1], are implemented.

Theoretical concepts of each TFR method is out of scope of the current document, however, interested users can refer to [24] for a detailed and comprehensive discussion on time/frequency analysis concepts and its techniques. Nevertheless, as a short description, in STFT, the TF map of each trial is calculated separately through Fourier Transform and averaged over all. Similarly, for CWT the scalogram displaying the squared and over-all-trials-averaged wavelet coefficients are used. Eventually, for NBCH, ERP maps are generated according to ERP behaviors within very narrow frequency bands. The ERP maps are matrices, the rows of which correspond to frequency-specific ERP estimations.

The implemented function enables variable frequency range. The procedures utilized for all three TFR methods here are in accordance with [1]. The function trigger_avg_TF_erp.m is provided for this purpose which utilizes trigger_ synch.m function to synchronize cerebral signals according to their trigger time and BPFilter5.m function to perform the aforementioned CIC filtering for NBCH method. Moreover, for STFT and CWT methods, MATLAB’s inner spectrogram and cwt functions are used. To operate this function you can use the demo m-file test_ERPanalysis_main.m or simply call function below with proper input parameters:

[erp_tf, synch_eeg, synch_emg, trigger_time_sec, freq_vec, time_vec] = trigger_avg_TF_erp(eeg, emg, fs, emg_onset_sampl, duration, method)

where the inputs and outputs are described later in Section 3. Fig. 5 is an elaboration of how how this function works. It is worth noting that the CWT produces similar results. Note that several benchmark studies such as [25], [26], [27] and [28] are research show-cases of the real-world applications how these TF analysis data can be put to use. For instance, as investigated and shown by [29] and also [27] and [30] inspecting the irregularities of the TF patterns through entropy quantities can be very informative regarding the underlaying motor-task [25], depth of sleep and Neurons activity rate [29, 27], levels of Alzheimer’s disease [31, 28] and even the cognitive skills and their respective performances in subjects [32, 26, 13].

2.2.5 Time-course Cortical Functional Connectivity

Time-course cortical functional connectivity is usually used to investigate the local-scale brain connectivity across different electrodes within same brain regions (i.e. primary motor cortex). To this end, Phase Locking Value (PLV) as the common quantity for this purpose [7, 9, 33] are implemented, measuring both temporal dynamics of local-scale and heat-maps of large-scale functional connectivity. The details of each method is described below.

Phase-locking Value (PLV):

PLV is a measure for quantifying how constant the phase difference between two signals is. In order to calculate the PLV in frequency $f$ for two signals (or channels) $x(t)$ and $y(t)$, the following steps are required [7, 9]:

• •

Using narrow-band filters centered at $f$, calculate the instantaneous frequency-specific phase values $\phi_{x}(t,f)$ and $\phi_{y}(t,f)$.

• •

Calculate the instantaneous phase-difference between $x(t)$ and $y(t)$ and quantify the local stability of this phase-difference over time as follows:

[TABLE]

where $T$ is the signal length and the summation is over all temporal samples of the instantaneous phases.

PLV varies between 0 and 1, corresponding to completely non-synchronized signals and complete synchronization, respectively [7, 9]. In simple wolds, considering the difference between firing rates of neurons as phase difference between their electrical potentials, PLV measures how ”connected” different neurons in far regions of brain are through quantifying the difference between their firing rates [7, 9, 33].

One very important consideration is the phase estimation approach that one needs to take into account. For phase estimation in first step, the Transfer Function Perturbation (TFP) method introduced in [34] and [35] has been adopted. TFP is a statistical Monte-Carlo based phase estimation approach which uses infinitesimal perturbations on filter and other signal estimation parameters in order to account for stochastic properties of EEG signals. According to the discussion presented in Section 2.2.2, using TFP is specifically required due to presence of background EEG noise (which in this case cannot be rejected as simple). It has been shown in [34] and [35] that without using TFP, estimated phase of EEG may contain variations and spikes which do not have any physiological origins within brain or are not due to brain activity and are merely side-effects of estimation method. For this purpose, cerebral signal phase estimation toolbox [13] has been used.

In this toolbox, function TCPLV.m is provided for this purpose which first utilizes phase_est.m function to estimate the instantaneous phase (IP) sequence of input signals through introduced TFP method and then uses PLV_PhaseSeq.m function to quantify the phase difference between calculated IPs. Function TCPLV.m measures time-course PLV for 1sec time-steps either between two specified electrode channels or between one important channel and all other channels (resulting in a heat-map across scalp). As recent examples of using such framework and implementation (time-course PLV) we can mention [35, 27, 25]. To operate this function you can use the demo m-file test_connectivity_main.m or simply call function below with proper input parameters:

tcplv = TCPLV(eeg, fs, freq_rng, pairofint, emg_onset_sampl, duration, pertnum)

where the inputs and outputs are described later in Section 3. Figures 6, 7 and 8 are elaborations of how this function works in either of aforementioned cases. Read the figure captions for more details.

2.2.6 Pair-wise Cortical Functional Connectivity

Pair-wise cortical functional connectivity is usually used to investigate the large-scale brain connectivity across scalp (i.e. between different brain regions). To this end, the Magnitude Squared Coherence (MSC) and Phase Locking Value (PLV) as the common quantities for this purpose [6, 36, 7, 9, 33] are implemented, measuring heat-maps of large-scale functional connectivity.

Magnitude Squared Coherence (MSC):

The conventional method presented in [6] for measuring coherence is based on weighted windowing of the Fourier transform of signals [33]. Considering $x(t)$ and $y(t)$ as two channels of recorded EEG signals, MSC can be measured in frequency $f$ as below:

[TABLE]

The MSC here, is calculated using a non-overlapping hamming window using FFT analysis for EEG, in a frequency range of interest for all trials, and then averaged across trials. The mean MSC then is analyzed within 1 second temporal windows covering -3 to -2, -2 to -1, -1 to 0, 0 to 1 and 1 to 2 seconds, with respect to 0 set as trigger time. MSC measures are calculated for all possible electrode combinations, resulting in N$\times$N MSC maps (for N-channel EEG recording), representing the functional connectivity across the whole scalp. Similar to PLV, MSC varies between 0 to 1 [6, 33, 37].

In this toolbox, function PWCoherence.m is provided for this purpose which first utilizes BPFilter5.m function to perform the aforementioned CIC filtering and then uses MATLAB’s internal mscohere.m function to quantify the frequency-specific MSC values. As mentioned, function PWCoherence.m measures MSC maps for 1sec time-steps between all possible pairs of electrode. As recent examples of using such framework and implementation (pair-wise MSC) we can mention [26, 27, 25, 35]. To operate this function you can use the demo m-file test_connectivity_main.m or simply call function below with proper input parameters:

pwcoher = PWCoherence(eeg, fs, freq_rng, emg_onset_sampl, duration, plot_flag)

where the inputs and outputs are described later in Section 3. Figures 9 are elaborations of how how this function works in either of aforementioned cases. Read the figure captions for more details.

Phase Locking Value (PLV) for Large-scale Functional Connectivity:

Similar to MSC, PLV can be measured between all possible electrode pairs, resulting in PLV maps across entire scalp. We provided this option in this toolbox by implementing function PWPLV.m which first utilizes phase_est.m function to estimate the instantaneous phase (IP) sequence of all channels of input data through introduced TFP method and then uses PLV_PhaseSeq.m function to quantify the phase difference between calculated IPs. Function TCPLV.m measures pair-wise PLVs for 1sec time-steps between all possible electrode pairs resulting in N$\times$N PLV maps (for N-channel EEG recording), representing the functional connectivity across entire scalp. As recent examples of using such framework and implementation (time-course PLV) we can mention [34, 35, 27, 25, 26]. To operate this function you can use the demo m-file test_connectivity_main.m or simply call function below with proper input parameters:

pwplv = PWPLV(eeg, fs, freq_rng, emg_onset_sampl, duration, pertnum, plot_flag)

where the inputs and outputs are described later in Section 3. This function also outputs figures similar to Figures 9.

2.2.7 EMG Quantification

The purpose of this section of toolbox is to estimate and quantify electromyogram (EMG) signals recorded from active muscles during performing motor tasks. The EMG quantification is a common practice within Neuroscience community in order to evaluate the work-load of brain (through EMG curve’s immediate slope after movement onset) and also for the rehabilitation purposes (stronger muscles after brain stimulation or training). The basic steps for EMG quantification are as below [10]:

• •

Removing Electrocardiogram (ECG) artifacts from EMG signals of each trial

• •

Form the full-rectified signal from the remaining EMG signal of each trial

• •

Trigger-average the rectified signals across trials

• •

Extract the trend of the averaged signal as the quantified EMG

The first step of this process, i.e. removing the ECG artifact, can be processed through several different simple or more complex approaches [38]. High-pass frequency filters (HPF) [39], ECG subtraction through QRS-complex detection (FAS) [40], adaptive filtering (AF) approaches [41, 42], ICA based approaches [43, 44] and also combined AF-ICA based methods [45]. Here in this toolbox however, we promise to stick with the simple and fast methods and thus, implement a modified high-pass frequency filtering approach not only to gain from the simplicity but to also improve the performance accuracy of removing ECG artifact from EMG signals. Accordingly, our modified ECG removal approach includes low-pass filter combined with median filtering (LPF+MF). HPF is normally used only to remove the major ECG frequency components which are restricted by the (1-30Hz) [46]. The main drawback of this approach is due to the fact that HPF will remove the frequency components of EMG too [46]. Nevertheless, to overcome this limitation, we implemented and used a LPF+MF to first detect the ECG signal (instead of blindly removing its components) and then perform a temporal subtraction. To this end, first we track the major ECG frequency components (1-30Hz) with a fourth order (can be modified within function parameters) elliptic IIR low-pass filter with pass-band ripple and stop-band attenuations of 0.1dB and -50dB, respectively. As before, the filtering is processed in a forward-backward zero-phase manner to avoid any phase distortion due to non-linear phase response of IIR filters. Afterwards, a median temporal filter with small window-length is used to estimate and extract the ECG patterns more accurately [10, 46, 47, 48].

In this toolbox, function emg_quantification.m is provided for this purpose which first utilizes MATLAB’s internal ellip function to generate an elliptic IIR filter and then uses BaseLine2.m function twice; first to track and remove the ECG components and then to quantify the EMG signal as described before. To operate this function you can use the demo m-file test_EMGanalysis_main.m or simply call function below with proper input parameters:

[emg_bl, synch_emg2, ecg_estimate2_bl, time_vec] = emg_quantification(emg_data, fs, emg_onset_sampl, duration)

where the inputs and outputs are described later in Section 3. Fig 10 is an elaboration of how this function works. Performance of simple HPF is also shown for comparison and elaborating the reliability of the propose approach (i.e. LPF+MF). Maximum magnitude (muscle activation strength), activation slope (muscle activation speed) and immediate activation slope after movement onset (muscle-brain work load) are some potentially important information captured from quantified EMG through our provided function. See Fig 11 for graphical elaboration of these concepts.

EMG Movement Onset Detection:

As promised before in Section 2.2.3, here, we provide users with a simple but fast and reliable EMG onset estimation function. The onset detection in EMG signals has a broad literature in which various methods aimed to satisfy different purposes and applications are proposed and discussed. The simplest proposed approach for this purpose is the single thresholding method [49] which later on has been modified into double-threshold methods [50, 51]. Although there are several complex highly sensitive methods for onset detection, their utility really depends on users’ application of interest and thus to reduce computational cost, here we designed and implemented a simple two-stage thresholding method which is highly reliable for ERP application. The proposed two-stage thresholding method is basically formed based on double-thresholding framework (which includes three steps) [50, 51]. The three steps are (1) signal conditioning (which is done in preprocessing phase), (2) detection of an event (performed by calculating a windowed standard deviation over time, named as STD vector) and (3) exact detection of the movement time (performed by using a threshold on STD vector’s temporal trend). As reported previously in [50] and [51], moving temporal filtering and thresholding methods are very common and could be satisfactory based on application requirements. Additionally, according to discussions in [50] and [51], parameters and thresholds in such approaches need to be set empirically based on the data and application.

In this toolbox, function emg_onset.m is provided for this purpose which follows the exact procedure as described above. To operate this function you can use the demo m-file test_EMGanalysis_main.m or simply call function below with proper input parameters:

[onset_sampl, onset_time] = emg_onset(emg, fs, W, th_coeff, Trl)

where the inputs and outputs are described later in Section 3. Fig 12 is an elaboration of how this function works. It worth noting that, although the implemented procedure is not as sensitive as complex new methods, it still satisfies our requirements for the application of ERP analysis [47, 48]. This is due to the fact that ERP is calculated by averaging over a large number of trials and thus here we are interested in intra-trial variations and information rather than inter-trial ones.

3 Reference Manual

3.1 bdf2mat_main.m

bdf2mat_main.m

Purpose:

Reading EEG/EMG signal from bdf files, preprocessing and conditioning.

Synopsis (global mode):

[eeg_data, emg_data, fs] = bdf2mat_main(trl_num, elec_num, emg_ch, eeg_ch, filename)

Synopsis (local mode):

[eeg_data, emg_data, fs, emg_onset_sampl, emg_onset_time] = bdf2mat_main(trl_num, elec_num, emg_ch, eeg_ch, filename, drift_flag, onset_flag)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

An empty bracket [.] Must be assigned to not-specified values.

• •

EMG and EEG channels of interest have to be a doubles (i.e. either an integer or a vector of indices).

• •

filename needs to be defined in a way that trl_num can be used as an index to track and load all of the stored data. For instance, data file names may be defined as filename_i where $i=1,2,...,N$ in which N represents the total number of data files. Refer to our sample data included in this toolbox to see another example of this, if needed.

• •

One should expect the emg_onset_sampl and emg_onset_time within output arguments in case they flagged the onset detection as 1 within input.

3.2 trigger_synch.m

trigger_synch.m

Purpose:

Synchronizing EEG/EMG signals according to movement onset time.

Synopsis (global mode):

[ensemble_eeg, synch_eeg, trigger_time_sec, time_vec] = trigger_synch(eeg, fs, onset_time)

Synopsis (local mode):

[ensemble_eeg, synch_eeg, trigger_time_sec, time_vec, synch_emg] = trigger_synch(eeg, fs, onset_time, duration, emg)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

An empty bracket [.] Must be assigned to not-specified values.

• •

The provided function ’bdf2mat_main.m’ outputs the compatible EEG/EMG data for this function; however, you can simply store your data in a cell array in which each element includes one data trial(s) of interest.

• •

The output argument ’ensemble_eeg’ is an essential variable to be used within the provided function ’trigger_avg_erp.m’ and thus you can simply replace it with ~ symbol if using this function independently.

• •

One should expect the synch_emg (e.g. synchronized EMG data) within output arguments in case they inputed EMG data.

3.3 trigger_avg_erp.m

trigger_avg_erp.m

Purpose:

Trigger-averaged ERP time-course estimation function to extract the temporal dynamics of event-related potentials.

Synopsis (global mode):

[erp, synch_eeg, trigger_time_sec, time_vec] = trigger_avg_erp(eeg, fs, freq_band, onset_time)

Synopsis (local mode):

[erp, synch_eeg, trigger_time_sec, time_vec, synch_emg] = trigger_avg_erp(eeg, fs, freq_band, onset_time, duration, emg)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

An empty bracket [.] Must be assigned to not-specified values.

• •

Available options for freq_band are: ’delta’, ’theta’, ’alpha’, ’beta’ or ’gamma’ where are defined in ranges 1-4(Hz), 4-8(Hz), 8-12(Hz), 12-32(Hz) and 32-80Hz, respectively. In case you need to change these ranges, open up the function script and change f0 and bw values for the band you wish to alter (i.e. lines 105 to 115). Note that f0 and bw are center and bandwidth of the frequency range, respectively.

• •

The provided function ’bdf2mat_main.m’ outputs the compatible EEG/EMG data for this function; however, you can simply store your data in a cell array in which each element includes one data trial(s) of interest.

• •

One should expect the synch_emg (e.g. synchronized EMG data) within output arguments in case they inputed EMG data.

3.4 trigger_avg_TF_erp.m

trigger_avg_TF_erp.m

Purpose:

Trigger-averaged ERP T/F representation to track the time-frequency dynamics of event-related potentials.

Synopsis (global mode):

[erp_tf, synch_eeg, trigger_time_sec, time_vec, freq_vec] = trigger_avg_TF_erp(eeg, fs, onset_time)

Synopsis (local mode):

[erp_tf, synch_eeg, trigger_time_sec, time_vec, freq_vec, synch_emg] = trigger_avg_TF_erp(eeg, fs, onset_time, duration, method, emg)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

An empty bracket [.] Must be assigned to not-specified values.

• •

The provided function ’bdf2mat_main.m’ outputs the compatible EEG/EMG data for this function; however, you can simply store your data in a cell array in which each element includes one data trial(s) of interest.

• •

One should expect the synch_emg (e.g. synchronized EMG data) within output arguments in case they inputed EMG data.

3.5 erp_quantification.m

erp_quantification.m

Purpose:

ERP area (ERD and ERS events’ area) quantification.

Synopsis (global mode):

[erd_area, ers_area, quant_erp] = erp_quantification(erp, fs, trigger_time)

Synopsis (local mode):

[erd_area, ers_area, quant_erp] = erp_quantification(erp, fs, trigger_time, ref_per, cof_intv)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

An empty bracket [.] Must be assigned to not-specified values.

• •

Reference period variable ’ref_per’ have to be a double vector in [-a, -b] form where -a and -b are the edges of reference segment in Seconds. Minus sign shows that this period belongs to before movement onset (i.e. trigger time).

• •

erd_area and ers_area are stored in a variable-size cell array the size of which depend on the duration of ERP signal after movement onset. For instance, in figure 4 there are 2 ERD events and only one ERS event detected. As a results, the corresponding area cell arrays will be of size 1$\times$2 and 1$\times$1, respectively.

3.6 TCPLV.m

TCPLV.m

Purpose:

PLV temporal dynamics estimated within 1 Second time-steps for any arbitrary time range before and after movement onset.

Synopsis (global mode):

tcplv = TCPLV(eeg, fs, onset_time)

Synopsis (local mode):

tcplv = TCPLV(eeg, fs, onset_time, freq_rng, duration, pairofint, pertnum)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

An empty bracket [.] Must be assigned to not-specified values.

• •

freq_rng has to be in [a, b] form double vector where ’a’ and ’b’ are edges of the frequency band of interest in Hz.

• •

duration has to be in [-a, b] form double vector where ’a’ is the required duration (Seconds) prior to movement onset and ’b’ is the required duration (Seconds) after the movement onset. The specified duration will be segmented into 1-sec-long windows to calculate PLV. Refer to Section 2.2.5 for more details.

• •

pairofint can either be a double vector or ’all’ string. As a double vector it has to be in [a, b] form resulting in PLV measures between electrode number a and electrode number b. Moreover, if you use ’all’, a PLV map will be calculated between C3 electrode and all other electrodes to cover large-scale functional connectivity dynamics between motor cortex and other brain regions. In this case, C3 is chosen as the default electrode; however, if you wish to change this for your data, open up the script, scroll to line 143 of code and change the value 6 to you electrode number of interest.

3.7 PWPLV.m

PWPLV.m

Purpose:

Pair-wise Phase-Locking Value (PLV) dynamics estimated within 1-Second time-steps for any arbitrary time range and all electrode pairs.

Synopsis (global mode):

pwplv = PWPLV(eeg, fs, onset_time)

Synopsis (local mode):

pwplv = PWPLV(eeg, fs, onset_time, freq_rng, duration, pertnum, plot_flag)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

An empty bracket [.] Must be assigned to not-specified values.

• •

freq_rng has to be in [a, b] form double vector where ’a’ and ’b’ are edges of the frequency band of interest in Hz.

• •

duration has to be in [-a, b] form double vector where ’a’ is the required duration (Seconds) prior to movement onset and ’b’ is the required duration (Seconds) after the movement onset. The specified duration will be segmented into 1-sec-long windows to calculate PLV. Refer to Section 2.2.5 for more details.

• •

This function utilizes all of the electrodes (i.e. not just motor cortex electrode C3) and all possible combinations to calculate PLV values and generate connectivity maps. Refer to Section 2.2.6 for more details.

3.8 PWCoherence.m

PWCoherence.m

Purpose:

Pair-wise Magnitude-Squared Coherence (MSC) dynamics estimated within 1-Second time-steps for any arbitrary time range and all electrode pairs.

Synopsis (global mode):

pwcoher = PWCoherence(eeg, fs, onset_time)

Synopsis (local mode):

pwcoher = PWCoherence(eeg, fs, onset_time, freq_rng, duration, plot_flag)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

An empty bracket [.] Must be assigned to not-specified values.

• •

freq_rng has to be in [a, b] form double vector where ’a’ and ’b’ are edges of the frequency band of interest in Hz.

• •

duration has to be in [-a, b] form double vector where ’a’ is the required duration (Seconds) prior to movement onset and ’b’ is the required duration (Seconds) after the movement onset. The specified duration will be segmented into 1-sec-long windows to calculate PLV. Refer to Section 2.2.5 for more details.

• •

This function utilizes all of the electrodes (i.e. not just motor cortex electrode C3) and all possible combinations to calculate MSC values and generate connectivity maps. Refer to Section 2.2.6 for more details.

3.9 emg_quantification.m

emg_quantification.m

Purpose:

EMG signal analysis and quantification.

Synopsis (global mode):

[emg_quant, synch_emg2, ecg_estimate2_bl, time_vec] = emg_quantification(emg_data, fs, emg_onset_sampl)

Synopsis (local mode):

[emg_quant, synch_emg2, ecg_estimate2_bl, time_vec] = emg_quantification(emg_data, fs, emg_onset_sampl, duration)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

3.10 drift_reject.m

drift_reject.m

Purpose:

Drift (baseline wander) cancellation from biosignal recordings.

Synopsis (global mode):

sig = drift_reject(raw_sig, L1)

Synopsis (local mode):

sig = drift_reject(raw_sig, L1, L2, approach)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

An empty bracket [.] Must be assigned to not-specified values.

3.11 emg_onset.m

emg_onset.m

Purpose:

EMG onset detection function based on introduced two-stage thresholding approach in Section 2.2.7.

Synopsis (global mode):

[onset_sampl, onset_time] = emg_onset(emg, fs, W)

Synopsis (local mode):

[onset_sampl, onset_time] = emg_onset(emg, fs, W, th_coeff, Trl)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

th_coeff is a coefficient that will be multiplied by the standard deviation of the EMG baseline and by default is set to 1 so that the estimation threshold will be equal to one-standard-deviation of baseline (i.e. 1$\times$STD{emg}).

• •

An empty bracket [.] Must be assigned to not-specified values.

3.12 phase_est.m

phase_est.m

Purpose:

Instantaneous phase estimation by the Transfer Function Perturbation (TFP) method [34, 35, 13].

Synopsis (global mode):

[phase_avg, freq_avg, amp_avg, analytic_sig_avg] = phase_est(sig, fs, f0, bw_base)

Synopsis (local mode):

[phase_avg, freq_avg, amp_avg, analytic_sig_avg] = phase_est(sig, fs, f0, bw_base, pertnum)

Inputs:

[TABLE]

Defaults:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

pertnum is by default set to 100, however the value mostly depends on the application and also the required level of reliability. Nevertheless, 100 has been tested in several applications before and is considered enough [35, 27, 25, 26].

• •

TFP perturbation parameters are set according to the findings of [34] and are chosen in a way that are physiologically irrelevant. In case one needs to change these values due to different application specs, refer to lines 83 to 88 of code.

• •

An empty bracket [.] Must be assigned to not-specified values.

3.13 PLV_PhaseSeq.m

PLV_PhaseSeq.m

Purpose:

Calculating Phase Locking Value (PLV) matrix (Pairwise PLV) using phase sequences [13].

Synopsis:

PLV = PLV_PhaseSeq(phase_sig)

PLV = PLV_PhaseSeq(phase_sig1, phase_sig2, phase_sig3, ...)

Inputs:

[TABLE]

Outputs:

[TABLE]

Notes:

• •

While using the first case, the phase_sig have to be a matrix with at least two rows where each row represents a phase signal. In second case, each of the phase_sig1...phase_sign are row vectors of phase sequences. This option is provided in case that someone needs to calculate the PLV matrix between separate phase signals.

3.14 sig_trend.m

sig_trend.m

Purpose:

Calculating the trend of a signal using local minima.

Synopsis:

[Tr_sig, loc] = sig_trend(sig)

Inputs:

[TABLE]

Outputs:

[TABLE]

3.15 task_separator.m

task_separator.m

Purpose:

Task based separation of data to avoid memory overuse.

Synopsis (global mode):

sep_file_names = task_separator(filename)

Inputs:

[TABLE]

Outputs:

[TABLE]

4 Acknowledgment

The authors would like to thank Dr. Maysam Ghovanloo222GT-Bionics Lab, School of Electrical and Computer Engineering, Georgia Tech, Atlanta, GA, USA, Dr. Minoru Shinohara333Human Neuromuscular Physiology Lab, School of Biological Sciences, Georgia Tech, Atlanta, GA, USA, Dr. Boris I. Prilutsky444Biomechanics and Motor Control Lab, School of Biological Sciences, Georgia Tech, Atlanta, GA, USA, Dr. Andrew J. Butler555B.F. Lewis School of Nursing & Health Professions, Georgia State University, Atlanta, GA, USA and Zhenxuan Zhang666GT-Bionics Lab, School of Electrical and Computer Engineering, Georgia Tech, Atlanta, GA, USA for their insightful discussions and comments.

5 References