Characterizing peaks in the EEG power spectrum

Despite the mixing of periodic and nonperiodic components, EEG as a field has come to confound these aspects, referring to broad bands of the aggregate power spectrum or power spectral density (PSD) as oscillations. Such spectral decomposition historically gained popularity in the field of radio transmission and was used in teasing apart periodic components. However, when applied to complex nonperiodic signals like the EEG, interpretation must be more careful.

The convention that has emerged in the EEG field is to divide the power spectrum into broad bands of different frequency ranges and refer to each of these bands as 'waves' or 'oscillations' within that range. This definition has its origins in its early history in the 1930s and 40s when the absence of computers necessitated simplifications based on visually recognizable aspects of the signal [1, 9–11]. The tedium of performing manual computation on a signal obtained on a film or ink writing oscillograph gave rise to photomechanical frequency analyzers (e.g. the Walter and Kaiser analyzers) that produced outputs of particular signal bands directly on the oscillograph [10, 12]. Due to cost and technical constraints these analyzers had to limit to broad band analysis with pre-specified ranges. While the early researchers who built these analyzers warned of its limitations, particularly in application to more complex nonperiodic signals like the EEG [9, 10, 12], it has nonetheless become convention over the subsequent decades leading to use of language with an implicit supposition that the signal entails propagating waves of a particular frequency range where one frequency band is in essence independent of another.

The PSD of the EEG signal is typically a decreasing function with lower power at higher frequencies [13]. This decay in the PSD approximates, but does not strictly conform to, a 1/f^γ pattern (example shown in figure 1(A) red). Such 1/f structure implies underlying temporal correlations between the frequencies and reinforces a complex relationship between them [14], although the origin of this structure is still very much a subject of debate [8, 15, 16]. In human EEG the exponent γ that defines the steepness of the 1/f decay in the 2–25 Hz range can have a broad range of values from 0.5 to 2.25 across individuals and conditions [17]. However, under certain circumstances, periodicity or oscillations become evident in the signal. This periodicity, when it emerges and is strong enough, reflects as a peak in the power spectrum of the EEG signal over and above the 1/f-like spectral envelope (for examples see figures 4(D)–(F)). Periodicity is most evident in the signal when visual activity is blocked (eyes closed) and is typically around 10 Hz (within the alpha band, defined as the range from 7 to 15 Hz or 7.5 or 8 to 12 or 13 Hz depending on research group since the band delineation is not always agreed upon) [18–20]. Less frequently, it can arise in other ranges as well or at multiple frequencies as we demonstrate below.

Zoom In Zoom Out Reset image size

However, spectral analysis methods typically used in EEG analysis do not separate such embedded oscillations from the background spectral activity. Researchers interchangeably use the term 'alpha waves' and 'alpha oscillations' to refer to the sum of the PSD across the alpha range of frequencies. This is so regardless of whether there is indeed an evident oscillatory component or not and irrespective of where the peak occurs within this range. However, as we will demonstrate here, the power of the alpha band is largely a reflection of the overall decay of the nonperiodic 1/f spectral envelope than the emergent periodic component. Alpha power as a measure is therefore not particularly sensitive to the oscillatory component and therefore misses important distinguishing information.

Here we define a method that can better isolate the peak or oscillatory component from the overall spectral envelope, defining a new measure which we term the oscillation energy or E_o. Other methods [21] have approached this challenge of estimation of the peak in isolation by subtraction of the background activity using sophisticated curve fitting. Our method represents an alternative approach that is more sensitive to detection and discrimination of small peaks that may be overwhelmed by the background and is relatively insensitive to the shape of the background activity.

We demonstrate our method using simulated colored noise with a 1/f^γ, structure in the power spectrum (i.e. gaussian distributed noise with a power law spectrum), embedded with a 10 Hz sine wave. We generated this colored noise with γ varying from 1 to 2, which is the typical range for the slope of the PSD between 2 and 24 Hz for resting state EEGs [17]. To do so we used the python library colored noise 1.0.0 (https://pypi.org/project/colorednoise/) that is based on the algorithm defined by Timmer and Koenig [22]. Briefly, this method generates pseudo-random numbers, multiplies these numbers with a curve typical for that decay or γ value and then normalizes them. We then applied a band pass filter to the generated signals between 1 and 50 Hz as is typical in EEG studies, utilizing a total signal length of 3 minutes.

The power spectral density or PSD of these simulated signals is shown in figure 1(A) (grey to black) alongside the PSD of an actual EEG signal (also three minutes in duration, 14 channel recording) that shows no evidence of an oscillation (red). In this case the PSD was computed using the standard Welch method (8 segments with hamming window and 50% overlap). The example EEG signal approximates a decay exponent γ of 1.75. (Note the log scales for both axes in contrast to the semi-log scale in other panels in this figure).

Each of these signals was then added together with a 10 Hz sine wave with a range of amplitudes. Specifically, the amplitude of the sine wave was varied from 0.05 to 0.2 times the standard deviation (SD) of the amplitude distribution of the original signal. Figure 1(B) shows an example of a power noise signal with γ = 1.75 overlaid with the sine wave (above) along with the aggregate signal overlaid on the original signal (below). The two examples shown are for sine waves of amplitude 0.1 (top) and 0.2 (bottom). As can be seen, at these amplitudes, the added oscillatory activity is barely visible to the eye in the aggregate signal. However, it is evident in the PSD as a sharp spike at 10 Hz in both the simulated noise signals and the base EEG where the sine wave has been added (figures 2(D)–(F)).

Zoom In Zoom Out Reset image size

Periodicity embedded within a signal can best be visualized by looking at its autocorrelation function. Figures 2(A)–(C) show the autocorrelation for the simulated signals with varying amplitudes of the embedded sine wave (2(A), (B) power noise with γ = 1.75 and 10 Hz sine wave amplitude of 0.1SD and 0.2SD the signal amplitude; C: sine waves at 10 and 15 Hz at 0.2SD amplitude). Our method therefore involves first computing the standard autocorrelation (R) of the signal.

We next compute the fast fourier transform or FFT of the autocorrelation which provides an estimate of the power spectrum or PSD. We note that one could use any spectral transformation such including the Welch method (pwelch), the FFT² or the autoregressive spectral analysis of the signal in its raw or pre-processed forms Fundamentally these are all intended as estimates of the power spectrum of the signal, but practically, since they are only estimates, they differ considerably due to various parameters within the algorithm such as the type of window and window overlap used. The Welch method uses greater window overlap and averaging which smooths the spectrum, broadens any peaks and also results in shallower slopes overall, while the FFT² has more noise. Autoregressive spectral analysis is a different algorithm that is more suited for short signal lengths. These aspects are extensively discussed in signal processing textbooks so we do not address this in further detail here. We have chosen to use the FFT of the autocorrelation (R) since the autocorrelation provides an intermediate step that provides a useful visualization of the presence of periodicity. We compute this using a standard FFT function that utilizes a hanning window and represent it with the notation S_R so as to distinguish from the PSD computed using other methods.

$$\begin{eqnarray}&&{S}_{R}=\left|\mathrm{FFT}(R)\right|\end{eqnarray} \tag{ 1 }$$

Like the PSD shown in figure 1, this produces peaks corresponding to the frequency of the embedded oscillation (figures 2(D)–(F)). However, when the amplitude of the embedded oscillation is small (for example as in 1(C)), the peak is often barely visible. Thus, to detect small oscillations the peak must be amplified.

To amplify detection of the oscillation and suppress the background decay we use a second order differential of the function. This is a method commonly used in spectroscopy and other fields to amplify the detection of small signals embedded within decaying backgrounds. The higher order the differential, the greater the amplification of small signals. Here we compute $S{{\prime\prime} }_{R}$ as the second order difference function of the FFT of the autocorrelation (figures 2(G)–(I), grey line).

$$\begin{eqnarray}&&S{{\prime\prime} }_{R}=\displaystyle \frac{{d}^{2}{S}_{R}}{{{df}}^{2}}\end{eqnarray} \tag{ 2a }$$

for which the numerical derivative is equivalent to

$$\begin{eqnarray}&&{S}_{R}^{{\prime\prime} }=\displaystyle \frac{{X}_{j+1}-2{X}_{j}+{X}_{j-1}}{{\rm{\Delta }}{f}^{2}}\end{eqnarray} \tag{ 2b }$$

We then compute the Hilbert amplitude envelope (H) of this second order difference function (i.e. $H(S{{\prime\prime} }_{R})$ (figures 2(G)–(I), black line). Essentially this represents the absolute value of the analytic function obtained by the Hilbert transform. We do this using the env function in R within the seewave package (https://cran.r-project.org/web/packages/seewave/seewave.pdf) using a mean sliding window of length 4 points and an overlap of 50% between successive windows. The length of the window is the primary determinant of the degree of smoothing of the envelope.

The next step is to identify the peaks in the envelope. To do so we use a method that deliberately avoids the use of a fixed window or range (e.g. 7.5 to 12.5 Hz for the alpha band) within which to identify peaks, since fixed windows can bias the detection and miss peaks on or outside the boundaries. Instead, we use a threshold of one standard deviation (SD) of the envelope function to identify those segments with values continuously greater than 1SD. To accomplish this all values of the envelope below the 1SD threshold are set to 0. The thresholded Hilbert envelope H_T is thus

$$\begin{eqnarray}{H}_{T}:= \left\{\begin{array}{ll}0 & H\lt {SD}(H)\\ H & H\gt {SD}(H)\end{array}\right.;\end{eqnarray} \tag{ 3 }$$

We then find the index values for the nonzero values of H_T and compute the difference between successive index values. Differences greater than 1 indicate the end of one segment and start of the other. This can then be used to create vectors of the index values of H_T that represent the start and end values of each segment. For example if H_T = [0, 5, 7, 9, 0, 0, 6, 8, 6, 0] then the index vector I_S representing the start values would be [2, 7] and the Index vector I_E representing the end values would be [4, 9]. Each peak segment P₁, P₂, ..., P_n would thus be represented as:

$$\begin{eqnarray}\begin{array}{rcl}{P}_{1} & = & {H}_{T}\left({I}_{S}(1)\right):{H}_{T}\left({I}_{E}(1)\right)\\ {P}_{2} & = & {H}_{T}\left({I}_{S}(2)\right):{H}_{T}\left({I}_{E}(2)\right){\rm{etc.}}\end{array}\end{eqnarray} \tag{ 4 }$$

We then exclude all segments that have a length or width of less than 0.1 Hz (∼12 data points) to eliminate short noise periods that are not real peaks. Each segment greater than this length cut-off is then considered a peak segment. In our simulated example the width is about 1 Hz while peaks in real EEG data typically span 2–5 Hz.

To characterize the peaks we first take each segment of the envelope that meets the above criteria and compute the maximal value to identify the peak position. The peak frequency (P_f) is then computed as the position or index on the frequency axis corresponding to the peak value shifted by two positions to the left since the envelope is computed on the second order derivative or difference function and is therefore shifted by two data points.

$$\begin{eqnarray}&&{P}_{f}=\displaystyle \frac{\left({\rm{index}}(\max P)-2\right)}{{F}_{S}}\end{eqnarray} \tag{ 5 }$$

Where F_S is the sampling rate.

We then compute the energy of the oscillation (E₀) for each peak as the area under the Hilbert envelope for each segment identified above. as:

$$\begin{eqnarray}&&{E}_{0}=\displaystyle \sum _{i=1}^{\max }P(i)\end{eqnarray} \tag{ 6 }$$

Since the embedded oscillation in our first two examples have a peak at 10 Hz which is within the alpha band, we call the computed metric E_α. In the third case the second peak is at 15 Hz (within the beta band) so is defined as E_β.

We note that the bandpass or low pass filters applied at 1 Hz create a large peak in the power spectrum around the filter cut-off. Given this, any oscillations at frequencies less than 3 to 5 Hz are unlikely to be accurately detected as they will be indistinguishable from peaks arising from the filter as well as measurement artifacts such as baseline drifts. However, if low pass filters applied are closer to 0.1, it's possible that peaks at 1 Hz could be detected.

We next demonstrate the sensitivity of this measure in detecting and discerning differences in the amplitude of the embedded oscillation. In figure 3(A) we show E_α for different amplitudes of embedded 10 Hz sine waves. As can be seen, the curves are fairly similar for signals of different background decays (grey to black). Conversely, when visualized as E_α vs background decay for different amplitudes of the embedded sine wave (figure 3(B)), we see that the lines (corresponding to a particular embedded sine wave amplitude) are fairly flat allowing estimation of the embedded amplitude independent of the characteristics of the background decay

Zoom In Zoom Out Reset image size

For comparison we show some commonly used alpha power measures (both relative and absolute from 7.5 to 12.5 Hz) computed as a function of the amplitude of the embedded oscillation for different decays, and correspondingly as a function of decay for different amplitudes (figures 3(C)–(E)). Relative power was computed as the power within the band of interest (here 7.5 to 12.5 Hz) as a fraction of the summed power across the entire computed spectrum. Figures 3(C) and (D) show the relative alpha power computed using the square of the FFT (FFT²) as a function of sine wave amplitude and 1/f decay (i.e. γ value) respectively. As can be seen, in contrast to E_α, the alpha power was more separated for each value of γ for the same sine wave amplitude (figure 3(C)) and conversely declined more substantially as a function of γ for each embedded sine wave amplitude value (figure 3(D)).

We next show in figure 3(E) how a shift in decay (i.e. γ) of the 1/f spectrum from 1.5 to 1.75 distorts the estimation of the same embedded alpha oscillation for various methods of characterization of this feature. We use percentage difference in the metric value between γ = 1.5 and 1.75 since the scales of each measure are different and show the values for two different amplitudes (0.15 and 0.2). While E_α resulted in a difference of 8% between the two decays for an amplitude of 0.15 and 6% of an amplitude of 0.2, alpha power measures resulted in differences of anywhere from 13% to 78% depending on the method used. In general, using the pwelch method resulted in less decay dependent distortion relative to using the FFT². This is most likely due to the additional averaging of windows in the welch method which can result in broader peaks in the estimated power spectrum. Further both relative and absolute alpha power gave similar results for pwelch but relative performed better for FFT².

Thus, studies that seek to study alpha oscillations using alpha power (relative or absolute) are prone to misinterpret changes in decay of the 1/f spectrum as changes on the alpha oscillation. Indeed alpha changes reported using power spectrum could easily arise from a small shift in the decay rather than any change in the periodicity or oscillation, particularly when the oscillations in the signal are weak. Further, these comparisons demonstrate the substantial variability in results depending on the method used for computing the alpha band power. This suggests caution in the interpretation of results using these traditional metrics.

Our measure of oscillation energy (E₀) will allow researchers to detect smaller oscillations with greater sensitivity and specifically study changes in the peaks, which thus describes the oscillatory component more independent of the decay of the background envelope.

Evidence suggests that the background 1/f signal is largely composed of complex waveforms rather than multiple independent periodic components. However, the conventional method of assessing oscillations are based on summing the power spectral density or PSD within a particular frequency band that (i) include both the oscillation and 1/f background component in estimation and (ii) using fixed frequency boundaries. Both of these aspects introduce significant error into estimation of the periodic or oscillatory component. While oscillations most frequently arise in the alpha range (7.5 to 12.5 Hz), the peak values and the width of the power spectrum that they span can vary based on both intrinsic aspects of the signal and parameter choices used in computing the PSD (e.g. window overlap). In some cases (as in the example shown in figure 4(A)), the peak may spread across two bands resulting in an underestimation when considering one band alone. Further, oscillations can arise simultaneously at multiple frequencies (i.e. two peaks) and can be missed.

Using colored noise simulations as well as an example of a real EEG signal embedded with either a single 10 Hz sine wave or multiple sine waves at 10 and 15 Hz, we have defined a method to characterize multiple peaks emerging from embedded oscillations in a manner that is relatively independent of the nonperiodic background spectral envelope. We have termed this measure Oscillation Energy, to distinguish it from power, and abbreviated to E₀, it provides an aggregate measure of the strength and fidelity of the oscillation. Our measure is less sensitive to changes in the 1/f decay of the overall spectrum and therefore performs better than the traditionally used estimations of power within any a spectral band such that differences between individuals or conditions can be more directly attributed to the periodic component. Further, since it amplifies the peaks it has a higher sensitivity to smaller peaks relative to other methods to separate out the peak [21]. We use E_α to indicate the energy of the oscillation whose peak is detected within the alpha band, E_β to indicate the energy of the oscillation in the beta band and so on.

The fundamental premise of this method rests on the property of derivatives of peak-like signals whereby the amplitude of the nth derivative is inversely proportional to the nth power of the width of the peak (for a peak of the same shape and amplitude). Consequently, when a peak arises out of a slower decaying background, representing a relatively narrow period of faster change, each higher order derivative will amplify the peak while suppressing the background which is effectively a stronger but broad peak.

The key advantage of this measure is that it specifically amplifies the periodicity and minimizes the impact of the decay of the PSD, allowing both identification of the presence of small oscillations as well as greater discrimination of specific differences in the oscillatory component. It is also aided by defining the period of interest based on a threshold rather than by an arbitrary frequency range, and therefore provides a 'band' independent metric. Since the second derivative with respect to a position is acceleration or momentum, the nomenclature of energy might be thought of in terms of its relation to the momentum with which the oscillation exceeds the background.

We do note however that the method does have some limitations in that it amplifies differences in a nonlinear manner and also cannot fully eliminate the impact of background decay. It therefore does not allow for a specific estimate of the oscillation amplitude. Second, while the evidence suggests that the background spectral decay is largely composed of nonperiodic waveforms, this may not always be the case. Third, very broad peaks and peaks close to the resolution of the power spectrum cannot be easily detected, making it a challenge to detect low frequency peaks arising at less than 1 Hz. Finally, the signal is non stationary and peaks that arise may not always reflect periodic activity. Therefore, a more generalizable interpretation of this measure would be the degree to which any peak exceeds the overall background activity, analogous to a signal to noise measure.

Finally, we note that the metric is impacted by several parameter choices such as the window type and overlap used in the FFT computation as well as the smoothing parameter of the Hilbert transform and the threshold of the Hilbert envelope used to define the peak segments and therefore oscillation energy. Consequently, unless the parameters used are identical (as is the case for other measures such as the power as well), the absolute values will not be informative. Rather this serves more as a relative measure to compare across two records.

In conclusion, E₀ provides a more precise metric for periodicity that will enable researchers in the field to develop a clearer understanding of the function of this feature of the signal. In this context, it could ultimately prove to be useful in understanding and predicting behaviors that rely on this feature. For example, attention has been shown to have a relationship already to the overall power in the alpha band [24, 25], although at low correlation levels, which possibly could be improved by use of this metric alone or in combination with overall alpha band activity. Further, since it is not yet known which features of the EEG signal, alone or in combination, are fundamentally important to each different behavior or outcome, there is a need in the field to explore many different features. For example, features of the EEG such as different frequency bands generally have on average 0.3 to 0.4 correlations with different mental health or behavioral outcomes [23]. Features extracted by this algorithm represents a new dimension that may be better able to predict different outcomes than other methods, or even perform better for some outcomes over others, eventually finding application in various areas.