Demystifying signal processing techniques to extract task-related EEG responses for psychologists

3.1.1 Measuring amplitudes

There are at least two ways of measuring the amplitude of an ERP component. The first one is to measure the voltage at the maximum point of a waveform within a specific time window, which can be called the “peak amplitude method” [5, 7]. The other is to calculate the mean voltage of all points in a defined time window, which can be termed the “mean amplitude method” [5].

The peak amplitude method is the oldest quantification of amplitudes and is used frequently [8 –10]. However, the amplitude measured by this method is sometimes counter‐intuitive, especially when the process is implemented automatically. Fig. 2(B) illustrates clearly why the most negative point within a time window, Ps, may not be properly used to define the negative peak. Normally, we would call Pl as the negative peak instead of Ps. As Luck pointed out, we can define Pl as the “local peak”, that is, the largest point in a measurement window with upward‐going and downward‐going deflections at its left and right side [5]. In contrast, absolute maximum or minimum points can be defined as the “simple peak”.

Generally, we want to measure the local peak amplitude, but not the simple peak amplitude. It is, therefore, important that we do not confuse one with the other. Mathematically, if the maximum (or minimum) point is not at the edges of a given time window, the simple peak amplitude is identical to the local peak amplitude. In practice, we can plot the ERP waveform to identify whether the maximum point occurs at the edges. If not, there is no need to distinguish them at all, and the traditional maximum algorithm produces the local peak amplitude.

The local peak method, however, is not without its own problems [5]. First and foremost, the local peak amplitude is positively biased by the noise level. That is, the measure is generally larger than the true peak amplitude when there are noises; the larger the noise, the larger the discrepancy. Unfortunately, EEG data are noisy, even all the preprocessing procedures have been performed. The noise level partly relies on the number of trials and participants, as well as the health conditions of participants. As a result, comparing the local peak amplitudes between markedly different conditions or groups can be problematic. Second, the local peak amplitude is attenuated by the latency jitter. If the latency varies across trials or participants, which is often the case, the average local peak amplitude would be smaller than the single‐trial amplitude [as illustrated in Fig. 2(A)]. This problem renders amplitude comparisons invalid when there are significant differences in latency variability and the number of trials between conditions or groups. Third, the local peak amplitude in the average ERP waveform is not a faithful representation of the average of local peak amplitude at the single‐participant level, unless different peak latencies are aligned to the same time point [11]. Without peak latency alignment, the average waveform may not reach its peak at the same time as the single-participant waveforms do; thus, the peak amplitude obtained from the average waveform is smaller than the average of single‐participant amplitudes in most cases.

By contrast, the mean amplitude is not (at least not easily) subject to any of the three problems. However, the mean amplitude method also has its own issues. The biggest drawback is that the time window, if too wide, may contain responses of other peaks or components, and thus the resulting mean amplitude is inaccurate. Fig. 2(B) illustrates a possible scenario where the mean P2 peak amplitude (area A) is partially canceled out by the N2 response (area B). One variant of the mean amplitude, termed the “mean signed area amplitude” [5], can circumvent the problem. To computing the signed area amplitude, we simply sum together either all the positive or negative points relative to the baseline within the time window. Since only positive or negative values are summed, the problem of cancelation disappears. In the figure, the mean signed area amplitude for the P2 component would be A/L instead of (A – B)/L (the mean amplitude), where L is the length of window. The price, however, is that the sign area amplitude tends to be positively biased, which is similar to the local peak amplitude [5]. In short, none of these methods is perfect, and we need to think carefully when choosing the appropriate method to extract ERP amplitudes.

3.1.2 Measuring latencies

Latency is the time interval between the starting time point of the targeted event (e.g., stimulus onset) and a specific ERP component. There is also more than one way of measuring latency. We can simply measure the time point where the waveform reaches its peak, or estimate the midpoint of a component by locating the time point which divides the area within the corresponding time window into two halves. The former latency is called the “peak latency”, while the latter is called the “50% area latency” [5].

These two measures of latency are based on the corresponding measures of amplitude, and thus inherit the strengths and weaknesses of the amplitude measures (see above). One thing to notice is that the peak latency is more sensitive to noises than the 50% area latency. The consequence is that the peak latency shows a larger variance, leading to lower statistical power. However, the 50% area latency is sensitive to the definition of the time window for extracting the ERP component, which would result in the inaccuracy of the measurement of ERP latency.

3.1.3 Practical issues

(1) How many trials do we need?

To derive a reliable ERP waveform, we need a sufficient number of trials [3, 5]. But how can we define “sufficient”? Unfortunately, we do not have any magic number for this question. The effect size, SNR, type of analysis, as well as many other factors determine the minimal number of trials needed to obtain a reliable and robust ERP waveform [5, 12 –14]. For a large ERP component, tens of trials would be enough; for a small ERP component, we may need hundreds of trials. One thing to remember is that the SNR improves in proportion to the square root of the trial number [5]. That is, to double the SNR, we should quadruple the number of trials. Therefore, if we already have hundreds of trials, it would be unnecessary to increase the number to a thousand. Instead, we should attempt to lower the noise level by other means.

(2) Do the baseline data reflect “pure” spontaneous activity?

No. To minimize voltage offsets and drifts due to skin hydration, skin potentials, and static charges in the electrodes, we correct the poststimulus data by subtracting the baseline mean from them, immediately after epoching the data [5]. However, the prestimulus baseline may be contaminated by the ERPs from the preceding trial or participant’s anticipations of the onset of the stimulus. Since the baseline‐corrected data are the difference between the uncorrected data and the baseline mean, the poststimulus responses are influenced by the baseline activity. We, therefore, should look at the baseline activity to ensure that the effect of interest is not an artifact caused by unwanted group or condition differences in the baseline. In practice, any differences that begin in the baseline or within 100 ms of stimulus onset may be spurious [5].

(3) What reference electrode should we choose?

The reference electrode has a huge impact on the ERP waveform. But there is no universally agreed‐upon reference electrode that fits well with all research questions. Cz, Fz, linkedears, linked‐mastoids, the ipsilateral‐ear, the contralateral‐ear, the tip of the nose, the average of all electrodes, and a point at infinite have been reported in the literature as the reference [15 –17]. In general, we should avoid the reference sites that are heavily biased toward one hemisphere, near the electrode of interest, or extremely noisy [5]. A practically‐accepted option is to choose the conventional one used in the relevant research field, thus allowing for comparisons of the results with those in previous studies.

(4) How can we visualize the results?

In the time domain, we have information about electric potentials across time and electrodes. We thus can depict graphs showing how electrical potentials vary across time points (i.e., ERP waveform) and how electrical potentials vary across electrodes (i.e., ERP topography). Providing these two kinds of graphs are mandatory for ERP studies [7].

ERP waveforms display the temporal characteristics of voltage without spatial characteristics [see the waveform in Fig. 2(B)]. Oftentimes we only choose one electrode of interest or compute the average waveform of several surrounding electrodes. We can, of course, draw the ERP waveform for every electrode and put them in the same graph (e.g., butterfly plots), but that might not be helpful as the reader may not be able to see the details like the difference of ERP waveforms between different conditions or groups. An oft‐used graph derived from ERP waveforms is the differential waveform, which could be created by subtracting one ERP waveform from another, showing how the difference changes over time.

ERP topographies, on the other hand, depict how the potentials are distributed across the scalp electrodes [see the topography in Fig. 2(B)]. We need to control another variable in this sort of graph, that is, time. The time point chosen usually corresponds to the latency of the component we are interested in. Alternatively, we can calculate the mean potentials within a time window and plot the corresponding scalp topography. One rule to remember is that the time variable (e.g., the interval for computing the average) must be identical for all the electrodes [7]; otherwise, the generated scalp topography would be confusing and difficult to interpret. However, there is no restriction on how many topographies we can draw so long as every one of them follows the rule. If necessary, we could provide a series of topographies for several time points (or several time intervals) for the ERP waveform.

3.2.1 Basic ideas of the STFT

The intuition of the STFT is straightforward. It cuts a single trial into numerous brief overlaping segments (i.e., time windows) where the assumption of stationarity supposedly holds. The STFT extracts the spectrum in every time window using the Fourier transform and assign the spectrum to the midpoint in the entire window, thereby providing temporal information. The STFT thus successfully addresses — or assumes away — the problems of unrealistic stationarity assumption and the lack of temporal information of the Fourier transform.

Mathematically, the STFT for a discrete signal x(n) (n = 1, 2,…, N) is defined as

X ( k , n ) = ∑ m = 0 N − 1 x ( m ) w ( m − n ) e − i 2 π k m N

where k and n are the frequency and time localizations respectively, w is the window function, and i is the imaginary unit [20]. In practice, we can follow the six steps below to run the STFT [21]:

“1. Select a window function of finite length.

Note that what we actually select is a window function rather than a simple window to diminish edge artifacts. Typical window functions include Hann (sometimes called Hanning), Hamming, and Gaussian windows. The Hann window is preferred since it, by definition, tapers the data to zero at its both sides [13]. The process of windowing and sliding the window along time indicates a vital issue in time‐frequency analysis: The spectra of contiguous time points would be strikingly similar, meaning the time precision of time‐frequency analysis is inherently by far poorer than the time domain analysis. Surely we can improve time resolution by narrowing the time window, but by doing so, we would have to compromise the frequency resolution, a phenomenon known as the uncertainty principle in signal processing [20].

3.2.2 Power‐based measures

The resulting X[k, n] from the STFT is a complex series. To quantify and test the time‐frequency results statistically, we need to derive more meaningful, power‐based measures, such as magnitude, power, and power spectral density (PSD). The only difference between these measures and those used in the spectral analysis is that we have the information for both time and frequency domains in task‐related EEG, whereas the information about time is irreverent in resting‐state EEG.

Magnitude is conceptualized as the absolute value of X(k, n). However, we need to consider the fact that a window function is applied, and results from the negative frequencies are redundant for real‐valued EEG signals. Consequently, magnitude is computed as

Magnitude = 2 | X ( k , n ) | S 1

where S ₁ represents the sum of window values w_j , that is, S 1 = ∑ j = 0 N − 1 w j . Power is the squared magnitude and defined as

Power = 2 | X ( k , n ) | 2 S 1 2

PSD is a derived measure from power and calculated as

PSD = 2 | X ( k , n ) | 2 F s S 2

where F _s denotes the sampling frequency, and S ₂ is the sum of squared window values ∑ j = 0 N − 1 w j 2 . Note that the window function plays an important part in all these three formulas.

Since all three measures are, by definition, non‐negative, we solve the non‐phase locked problem in time domain analysis, which results from the cancellation of positive and negative potentials. This problem is successfully solved only if we do the STFT at the single‐trial level.

With these measures at hand, we can define two important concepts describing task‐related power changes relative to the baseline activity: event‐related synchronization (ERS) and eventrelated desynchronization (ERD). The former means the poststimulus oscillatory power is greater than the baseline power; the latter, on the other hand, means the poststimulus oscillatory power is smaller than the baseline power. The exact formulas computing them rely on the baseline correction method, a topic we discuss later in 3.2.4 Practical issues.

3.2.3 Phase‐based measures

One particularly important measure we can extract from time‐frequency results is phase, which is the angle between the real part and the imaginary part of the complex series X(k, n) resulting from the STFT. Phase angle represents the timing of certain frequency‐band activity. To test the consistency of such timing across trials, we could calculate the phase‐locking value (PLV) [22], also known as phase resetting, phase coherence, or intertrial phase clustering [13]. Formally,

PLV ( t , f ) = | 1 M ∑ m = 1 M e i φ t , f , m |

where M is the trial number, i is the imaginary unit, φ_{t , f ,m} is the phase angle on trial m at certain time‐frequency point (t, f).

PLV can also be used to quantify connectivity between electrodes. The formula calculating PLV is

PLV ( t , f ) = | 1 M ∑ m = 1 M e i ( φ a − φ b ) |

where M is the number of time points (for PLV over time) or trials (for PLV over trial), φ_a and φ_b are the phase angles from electrode a and b, respectively. Note that M can be the number of time points OR trials here, whereas it can only be the number of trials in spectral analysis. This implies that we can compute PLV over time or trial. Intuitively, PLV is the absolute value of the averaged complex polar representation of phase angle differences. PLV ranges over [0, 1] and indicates the distribution range of the phase difference series over [0, 2π). Larger PLV implies higher connectivity.

PLV suffers from the problem of a common source, meaning that spurious connectivity occurs if activities from electrodes a and b are generated by a common source. To overcome this problem, we can use the phase lag index (PLI) [23], which is defined as

PLI ( t , f ) = 1 M ∑ m = 1 M sign ( φ a − φ b )

where M is also the number of time points (for PLI over time) or trials (for PLI over trial), sign() is the sign function, which equals −1 for negative inputs, 0 for zero, 1 for positive inputs. The range of PLI is also from 0 to 1, with higher values indicating higher connectivity. An improved measure, weighted PLI (wPLI) [24], is defined as

wPLI ( t , f ) = M − 1 ∑ m = 1 M | imag [ X a ( f ) X b * ( f ) ] | sign ( φ a , m − φ b , m ) M − 1 ∑ m = 1 M | imag [ X a ( f ) X b * ( f ) ] |

where X_a (f) is the Fourier coefficients for signals at electrode a, X_a* (f) is the complex conjugate of X_a (k), and imag() takes the imaginary part of a complex number. The complex conjugate of a complex number a + bi is a – bi. For example, 2 − 3i is the complex conjugate of 2 + 3i.

In practice, these phase‐based measures (e.g., PLV, PLI, and wPLI) were computed on the band‐passed filtered data. Note that these measures range from 0 to 1 and do not follow the normal distribution. When we compare phase‐based connectivity measures between conditions or groups, the results are bounded between −1 and 1 and are not normal-distributed, either. To facilitate meaningful statistical tests, we can Fisher‐z transform the condition or group connectivity differences to make them approximately normal‐distributed [13]. The Fisher‐z transformation is defined as

Fisher- z = 1 2 ln ( 1 + C 1 − C )

where ln() is the natural logarithm function, and C is some connectivity differences. Alternatively, we can use non‐parametric tests like the permutation test since they do not depend on the normal distribution assumption [13].

3.2.4 Practical issues

(1) How long should the time window for analysis be?

The length of the time window determines the trade‐off between time and frequency resolution. Shorter windows improve time resolution, but worsen frequency resolution; and vice versa. The window length also restricts the lowest frequency we can analyze. In theory, we can never extract the frequency information whose period is longer than the time window [13]. For example, the lowest frequency is 5 Hz if the window is 200 ms. To achieve a better SNR, we would want the window to contain two or three periods of certain frequency information, which further limits the lowest frequency we extract. In practice, we generally try out several different window lengths and choose the one that displays the best visually and statistically reasonable results [21].

(2) One STFT, one window function?

No. We can apply multiple window functions to the same data segment and average the spectrums, a method called “multitaper” [25]. In practice, to do multitaper analysis, we employ Slepian tapers, which are orthogonal to each other, rather than any arbitrarily chosen window functions [26]. The advantage of this method is that it improves the SNR, making it useful for analyzing high‐frequency activity and single‐trial data [13]. However, a higher SNR is achieved at the expense of worse time precision. The multitaper method thus is not particularly helpful for analysis that focuses on lower frequencies (e.g., below 30 Hz) or demands precise timing.

(3) How can we do baseline correction?

As in time domain analysis, baseline correction is necessary for time‐frequency analysis. One crucial reason is that raw time‐frequency power values follow the so‐called “1/f phenomenon”, meaning that the power decreases as the frequency increases [13]. This phenomenon renders it difficult to compare results from different frequency bands. Frequency‐wise baseline correction offers a partial solution.

There are four major approaches to perform the baseline correction: subtraction, percentage, decibel, and z‐score [13, 27], which are defined as follows:

Subtraction : P o w e r s u b = P o w e r − M e a n ( P o w e r b a s e l i n e )

Percentage : Power per = Power − Mean ( Power baseline ) Mean(Power baseline )

Decibel : Power dB = 10 log 10 Power Mean(Power baseline )

z - score : Power z = Power poststimulus − Mean ( Power baseline ) SD ( Power baseline )

where Mean(Power_baseline) is the average baseline power and SD(Power_baseline ) is the standard deviation of baseline power. It is noteworthy that the baseline need not be the entire prestimulus period. In fact, we had better choose a narrower baseline. The algorithm of the STFT implies that the power at the stimulus onset reflects the EEG activity of a time window whose midpoint is the stimulus onset. As a result, to avoid any contamination from the poststimulus activity, we should choose a baseline such that the distance between its ending and the stimulus onset is not larger than one‐half of the time window [27]. For example, if the window length is 200 ms, the baseline needs to end before −100 ms.

As for the approach per se, methods correcting the data by taking relative rather than absolute changes, that is, the percentage, decibel, and z‐score approaches, are superior to the subtraction approach in mitigating the 1/f phenomenon, making post‐corrected data more comparable between frequencies [27]. However, they all introduce a positive bias in estimating ERS and a negative bias in estimating ERD, and the resulting data are unproportionally influenced by the baseline power; the subtraction approach is free from these two problems [28]. A further advantage of the subtraction approach is that the order of baseline correction relative to averaging at the single‐participant level does not matter so long as the STFT is done at the single‐trial level. That is to say, performing the baseline correction either before or after averaging across participants provides exactly identical results. For the other three approaches, the order matters a lot. The bias in estimating ERS and ERD might be larger when relative‐changebased correction is performed at the single‐trial level than at the participant or group level [28].

(4) Should we compute connectivity measures over time or trial?

As we mentioned before, phase‐based measures, such as PLV, PLI, and wPLI, can be computed in two ways: over time or over trial. We are then forced to choose one of them. These two ways both have their own advantages and limitations [13]. Connectivity over time is less sensitive to trial‐to‐trial jitter of event timing and phase angle differences between trials. However, the epoch should be long enough for this kind of connectivity. Consequently, connectivity over time is more appropriate if we focus on high‐frequency connectivity. On the other hand, connectivity over trial puts no constraint on epoch length and achieves a higher level of evidence. Moreover, since connectivity is computed for every time point, the temporal precision is by far higher for connectivity over trial. Indeed, the STFT is not well suited for computing this connectivity; other methods like wavelet transform and Hilbert transform are preferred. The limitation is that connectivity over trial is highly sensitive to the event timing. Therefore, any uncertainty of timing may invalidate this measure.

(5) How can we visualize the results?

Time‐frequency analysis performed on a single trial produces three‐dimensional results: time, frequency, and certain measure (e.g., magnitude). A simple two‐dimensional plot thus is not adequate for visualization. To overcome this challenge, we could color the plot to add an extra dimension representing the measure we have computed [Fig. 3(B)]. Coloring the time‐frequency plot, however, is not easy. Please note that power follows the 1/f phenomenon. When using a single‐color scale, high‐frequency results would be distorted visually. That is, highly significant differences at high‐frequency bands would appear negligible compared with the results at lower frequency bands [21]. To tackle this nontrivial issue, we can use two color scales instead: one color for low frequency (delta, theta, alpha, and beta) results, the other for high frequency (gamma) results. However, we must keep in mind that the use of two color scales is for the purpose of visualization only, and that presented data must remain as they are.

In addition, we need to average the data across trials or participants. Although time-frequency results are available for every electrode, we may not want to show the results from all electrodes. Instead, we could choose one representative electrode or average several surrounding electrodes to plot time‐frequency results. To provide spatial information, we can draw a scalp topography by averaging power within the time‐frequency region of interest for each electrode [Fig. 3(C)].

3.3.1 Source analysis

One key feature of EEG data is that their spatial resolution is inherently poor. EEG data provide only the topographic distribution of electric potentials, not the spatial distribution of neural sources. Source analysis could be used to infer the latter from the former, which is a problem known as the inverse problem [29]. The inverse problem is notoriously difficult to solve. In actuality, this problem has no unique solution, meaning that more than one configuration of neural activity can produce the same topographic distribution and that we are theoretically unable to identify which configuration is absolutely correct for an EEG dataset. Localization of brain generators for specific electrical potentials is thus a challenging, if not impossible, task. Nevertheless, if conducted and interpreted carefully, source analysis provides valuable information about the brain regions generating the EEG activities recorded at the scalp, which is very helpful in a lot of basic and clinical applications. In practice, there are many methods to estimate the source, which are implemented in several freely- available softwares, for example, FieldTrip, sLORETA, and RECOR [29].

3.3.2 Single‐trial analysis

Traditionally, we average all single trials to boost the SNR and extract more accurate and consist brain features. However, the operation of averaging depends on an implicit and unrealistic assumption: The single‐trial response is largely invariant across trials [30]. We have mentioned in 3.1.1 Measuring amplitudes that latency jitter often occurs and, when they are extraordinarily larger, may invalidate the local peak amplitude. Nonetheless, trial variability includes at least some physiologically relevant information [31]. Therefore, across‐trial averaging is not desirable sometimes. Single‐trial measures do not confront these problems, yet its SNR is rarely satisfactory. Single‐trial analysis attempts to enhance the SNR and accurately estimates measures at the single‐trial level. There are several categories of single‐trial analysis methods, which include temporal filtering, spatial filtering, time‐frequency filtering, and model‐fitting approaches [30]. When properly conducted, single‐trial analysis enables us to extract plentiful dynamic information, such as single‐trial latency, amplitude, and morphology. This dynamic information allows us to reveal some additional information by performing a within‐subject correlation [32] and provide useful features for machine learning [33].