Using MNE to Characterize ADHD.

27 min readDec 28, 2023

Attention Deficit Hyperactivity Disorder (ADHD) is a neurodevelopmental disorder, affecting millions of children and adults worldwide.

It’s typically characterized by the lack of ability to focus, hyperactivity and impulsivity, or both combined.

To be more specific, ADHD can be represented by an inability to carry out desired tasks, forgetfulness, lethargy, fluctuating levels of motivation, and much more.

In a life where a person needs to use their cognition to perform at a high level, in order to accomplish what they want and become who they’d love to be, ADHD turns out to be detrimental.

Too many people, on the order of millions, grapple with the challenges of ADHD, enduring its impact on their ability to achieve their aspirations and fulfill their potential.

According to a meta-analysis conducted by Song et al., out of the entire adult population, approximately 366.33 million (6.76%) adults have exhibited ADHD-like symptoms¹.

Another meta-analysis of 175 research studies worldwide, showed that in children aged 18 and under, approximately 7.2% have ADHD².

This means that about 185,760,000 children and adolescents worldwide are grappling with the challenges of ADHD³ ⁴.

Therefore, a total of 552,090,000 worldwide are under the influence of ADHD.

While many people learn to live and thrive with ADHD, there’s a significant portion of the ADHD population that continuously struggle in their day-to-day lives.

9.1% (~3,421,600) of students that have ADHD drop out after the second year of university⁵ ⁶
1 out of 3 people diagnosed with ADHD are jobless at any given time⁷
50 to 60 percent of children with ADHD have difficulty with peer relationships which can carry over into adulthood⁸

Ultimately, ADHD poses a significant barrier to achieving optimal performance in various areas of life.

A Supportive Approach: Neurofeedback.

I’ve been curious about the efficacy of using neurofeedback (NFB), to both characterize ADHD and enhance the lives of those with ADHD.

For context, NFB is a type of biofeedback that leverages neural signals recorded through electroencephalography (EEG) to allow for real-time regulation of a persons brain state.

A person that has ADHD could use NFB to mitigate their hyperactive or inattentive symptoms in order to promote more focused behaviors.

While doing some research, I landed on a meta-analysis conducted by Arns et al., which showed NFB to be “Efficacious and Specific” for the treatment of ADHD⁹.

To be more specific, the use of the theta/beta frequency band ratio through EEG NFB, has been found to effective for gauging the levels of cognitive processing in ADHD¹⁰.

This study presents an important implication; that the theta/beta ratio may indeed be an index of cognitive processing
— Picken et al (2019)

The Theta / Beta ratio can also serve as a reliable measure for ADHD prognosis and progression¹¹.

TBR has prognostic value in this sub-group, warranting its use as a prognostic measure rather than a diagnostic measure.
— Arns et al (2012)

This particular EEG ratio has been used in clinical NFB systems. But unfortunately, clinical NFB systems are riddled with inefficiencies proving to be impractical for ease in daily use.

If we somehow manage to pack an EEG NFB system into a small yet highly effective wearable device, then perhaps ADHD treatment through NFB can become increasingly accessible, effective, and cost-efficient.

Given this, I thought to myself, “What components are needed to make such a neurofeedback system, reality?”

The Importance of Signal Processing

Signal processing is the analysis, filtering, and cleaning of a recorded signal.

A Very Basic Overview of Signal Processing Analysis

A key component in building such a neurofeedback system is the use of effective signal processing methodologies in order to create a robust signal processing pipeline.

Neural signals recorded through EEG are typically VERY noisy and exhibit a low signal to noise ratio.

This can be due to a variety of factors such as volume conduction, ocular/muscular/cardiac artifacts, high impedance, electrical interference, and much more which lead to a low signal to noise ratio (not good).

In order to ensure that a recorded signal is reliable for use, effective signal processing methodologies need to be implemented to properly filter and clean a recorded signal.

Given the important role that signal processing plays, I decided to process, clean, and filter a signal myself in order to gain an understanding in how it works.

I downloaded an ADHD EEG dataset (ADHD_part2) from IEEEDataPort, cleaned the signal, and then calculated the Theta/Beta ratio to characterize ADHD.

If you want access to the dataset, check it out here

Let me show you how I did this.

Implementing MNE

MNE is an, “Open-source python package for exploring, visualizing, and analyzing human neurophysiological data”.

It provides a comprehensive set of tools and functions needed for the analysis of MEG, ECoG, NIRS, EEG, and other neurophysiological recording formats.

The MNE library allows us to conduct important steps and apply certain tools needed when processing a signal such as bandpass/stop filters, notch filters, epoching, principal/independent component analysis (PCA / ICA), time-domain analysis, spectral analysis, event related analysis, and much more.

MNE also allows for the seamless integration of machine learning through scikit-learn and keras in order to better classify EEG data based on specific features.

All in all, it’s pretty robust.

For my project, MNE served as the primary library alongside matplotlib for data visualization, NumPy for data manipulation, and SciPy for loading data.

If you decide to replicate this project, I’d check out the MNE documentation for guidance.

Feel free to check out my video explanation.

@vxnuaj on youtube

Also, if you want to check out the project repository, feel free to click the link below to do so.

MNEPROCESS

Check out how I processed an ADHD dataset and calculated the theta / beta ratio!

github.com

The Processing Pipeline

For better visualization, I used JupyterLab.
Each code snippet represents a single cell.
For clarity, feel free to check out the full notebook on my Github.

First things first, we need to import all our libraries and dependencies.

import glob

import numpy as np

import matplotlib
import matplotlib.pyplot as plt

from scipy.io import loadmat

import mne
from mne.preprocessing import ICA

We import the module glob which allows us to interact with the current file directory in order to identify and get files that match a specific requirement.

We import the library, NumPy under the alias np in order for easy manipulation of our data.

We import matplotlib and the function, pyplot under the alias plt to allow for data visualization.

From the module under the SciPy library, scipy.io , we import loadmat which essentially allows us to load our EEG dataset in the form of the MATLAB file format, .mat.

Finally, we import mne and the ICA function from the mne.preprocessing module to allow us to adequately filter and clean our signal.

Now, we locate our dataset in our directory and grab it using the glob.glob() function.

matfiles = glob.glob('ADHD_part2/*.mat') 

matfiles

ADHD_part2 is a folder in the current directory and within that folder, the glob.glob() function grabs any .mat file and assigns it to the the list, matfiles .

When matfiles is called in the Jupyter cell, the output should be a list of files ending in .mat .

#This is the output of the previous code snippet

['ADHD_part2/v177.mat', 'ADHD_part2/v215.mat', 'ADHD_part2/v200.mat',
'ADHD_part2/v213.mat', 'ADHD_part2/v206.mat', 'ADHD_part2/v204.mat',
'ADHD_part2/v238.mat', 'ADHD_part2/v198.mat', 'ADHD_part2/v263.mat',
'ADHD_part2/v288.mat', 'ADHD_part2/v274.mat', 'ADHD_part2/v270.mat',
'ADHD_part2/v265.mat', 'ADHD_part2/v254.mat', 'ADHD_part2/v286.mat',
'ADHD_part2/v279.mat', 'ADHD_part2/v250.mat', 'ADHD_part2/v244.mat',
'ADHD_part2/v246.mat', 'ADHD_part2/v284.mat', 'ADHD_part2/v181.mat',
'ADHD_part2/v234.mat', 'ADHD_part2/v209.mat', 'ADHD_part2/v196.mat',
'ADHD_part2/v236.mat', 'ADHD_part2/v183.mat', 'ADHD_part2/v227.mat',
'ADHD_part2/v179.mat', 'ADHD_part2/v190.mat', 'ADHD_part2/v231.mat', 
'ADHD_part2/v219.mat']

After doing so, the dataset we intend to process is grabbed from the matfiles list.

mat_data = loadmat(matfiles[10])
eeg_data = mat_data['v274']

I used 11th .mat file in our directory which happens to be a file named v274.mat. We then assign it to mat_data .

Then the data is extracted from the v274 file and assigns it to eeg_data .

ch_names = ['Fz', 'Cz', 'Pz', 'C3', 'T3', 'C4', 'T4', 'Fp1', 'Fp2', 'F3', 'F4', 'F7', 'F8', 'P3', 'P4', 'T5', 'T6', 'O1', 'O2']
sfreq = 128

info = mne.create_info(ch_names = ch_names, sfreq=sfreq, ch_types = 'eeg')

eeg_data = eeg_data.T

Here, we define the meta-data for our EEG recording. We define our channel names in list format and assign it to the variable, ch_names .
The sampling frequency of our dataset, 128 Hz, is assigned to sfreq .

Using the MNE function, create_info() , we create an info instance based on the aforementioned metadata and the channel type, eeg , and then assign it to the info variable.

Using the NumPy function, ndarray.T, we transpose eeg_data into a NumPy array making it readable for MNE.

To make sure that the code worked as intended, we can run info in the next cell. The output should be as follows.

Now, we prepare our EEG data for visualization in the time-domain.

raw = mne.io.RawArray(eeg_data, info)
raw.apply_function(lambda x: x*1e-6, channel_wise=False)

The function mne.io.RawArray() allows us to take our eeg_data and the metadata represented by info and turn it into what’s called a RawArray object by reading our data directly from the NumPy array eeg_data . This RawArray object is assigned to the variable raw.

Afterward, we use the function RawArray.apply_function to perform the mathematical operation, x*1e-6 onto all channels at once.

This is done because our original data is expressed in microvolts (µV), since neural signals are typically have a minuscule level of power.

But MNE expects our data to be in expressed in volts so we use the function x*1e-6 to convert our data from µV to V.

Now, we can finally visualize our data in plot format.

mne.viz.set_browser_backend("matplotlib")
raw.plot(scalings="auto", duration = 10.0)

We use mne.viz.set_browser_backend("matplotlib") which sets our 2D backend to matplotlib. Now, all of our plots and visualizations will be drawn using matplotlib.

So, we then use the RawArray.plot() function to plot our RawArray object, raw. We set the duration to 10 as we’re only going to examine the first 10 seconds of our dataset. ‘Scalings’ is set to auto , meaning our plot will be automatically scaled to make it readable.

After running that cell, we should be able to visualize our data in a plot.

If you take a look at the entirety of channels F4, F8, and then Cz within the 1st half-second you can tell that our dataset is pretty noisy. There also seems to be a some noise at the 2nd and 7th second of our data.

Hence, we need to clean it.

In order to do so, we first need to figure out what frequency bands are responsible for the excess noise in order to filter it out.

spectrum = raw.compute_psd()
spectrum.plot()

We use the function RawArray.compute_psd() to compute the power spectral density (PSD). It allows us to perform spectral analysis on our EEG data. The function returns our data object in the frequency domain. We then plot our data using spectrum.plot() .

We get the following output:

We can clearly see that the noise in our data comes from the frequencies of 50 and ~56.

This is because the EEG recording originates in Tehran, Iran.
There, the power line frequency is at 50 Hz.
The origin of a spike at 56 Hz is for unknown reasons.

So we apply bandpass and notch filters to get rid of this noise.

bp_low_cut_freq = 2
bp_high_cut_freq = 40

raw_filtered = raw.filter(bp_low_cut_freq, bp_high_cut_freq)
raw_filtered.compute_psd().plot()

Our bandpass filter is set to have a low cut at 2 Hz and a high cut at 40 Hz. Using the function RawArray.filter() , we apply the low cut and high cut frequencies to create a bandpass filter and then apply it onto our data.

Our filtered data is then assigned to the object raw_filtered.
To see the effects in action, we can once more use the compute_psd() function to get the PSD and then perform spectral analysis. We use plot() to plot our dataset.

As you can see, frequencies > 40 Hz begin to drop off until there’s essentially no signal at frequencies > 50 Hz.

Yet we still see a small spike at 50 Hz. To mitigate such spike, we now apply what’s called a notch filter.

raw_filtered = raw_filtered.notch_filter(50)
raw_filtered.compute_psd().plot()

We use the function RawArray.notch_filter() to apply the notch filter directly at a frequency of 50 Hz.

After we, compute the PSD and plot our data once more, we should get the following output.

Our filtered EEG data in the spectral domain| Bandpass @ 2–40 Hz, Notch @ 50 Hz

As you can see, there’s no more excess noise at unnecessary frequencies.
We can now plot our filtered data in the time domain.

raw_filtered.plot(n_channels = 19, scalings = "auto", duration = 10.0)

Now, we get the following output:

Our filtered EEG data in the time domain

If we take a closer look at this plot, we can see that channel Cz seems to be excessively noisy. If we leave it in, it’d damage the fidelity of our signal and negatively affect any further analysis.

So we select it as a bad channel and remove it.

bad_channels = ['Cz']
raw_filtered.info['bads']=bad_channels
raw_filtered.plot(n_channels = 19, scalings = "auto", duration = 10.0)

We add the channel Cz to a list and assign it to bad_channels .
With the function, RawArray.info we modify the metadata of our dataset, specifically, 'bads', and set it equal to bad_channels .

Now, we can once again plot our dataset to verify that Cz has been marked as bad.

As you can see, channel Cz is greyed out indicating that it has been marked as a bad channel.

At this point, we can either interpolate our data using the function RawArray.interpolate_bads() or remove the channel altogether.

In this case, I chose to remove the entire channel from the dataset rather than interpolate it. I still encountered a noisy Cz channel post-interpolation so it was best to remove it.

raw_filtered.drop_channels(ch_names = bad_channels)
raw_filtered.plot(scalings = "auto", n_channels = 18, duration = 10.0)

The function RawArray.drop_channels drops the data our channels assigned to bad_channels which is channel Cz.

To ensure channel Cz was completely taken out of the equation, we plot our data once more.

We should get the following output.

As you can see, channel Cz, which is typically between channel Pz and Fz, is no longer in the equation.

Now, taking a closer look at our plot, we can see some noisy segments at about the 2 and 6.65 second mark.

To visualize this better, here’s the same plot with a modified scaling factor.

Those segments of data may represent redundant artifacts which can contaminate our signal.

So, using the following code mark them as ‘bad’.

onset = [1.95, 6.65]
duration = [.2, .30]
description = ['bad'] * len(onset)

annotations = mne.Annotations(onset, duration, description)
raw_filtered.set_annotations(annotations)

We set the onset or beginning of our noisy segments to the 1.95 and 6.65 second mark. Then the duration if each noisy segment is assigned as a list to duration .

For each value in the list onset , we generate a string `bad` .

`bad` is the description for our noisy segment will let MNE know which segment we intend to remove.

Per MNE documentation. Will be useful later in our program.

Through the mne.Annotations() function, we input our onset , duration , and description. This function creates an annotation object that will allow us to then annotate our dataset. We then assign it to annotations.

To actually annotate our data we use the RawArray.set_annotations() function. The function takes in our annotations variable to mark our specific segments as `bad`.

To verify that the function worked properly we can now plot our data once more.

raw_filtered.plot(scalings = 'auto', duration = 10)

We then get the following output:

To remove those bad segments, we first epoch our data into different fixed length time windows.

epochs = mne.make_fixed_length_epochs(raw_filtered, duration=.3, preload=False, reject_by_annotation=True)

Through the mne.make_fixed_length_epochs() we input our raw_filtered dataset and create epochs with a duration of 3 seconds.

By setting reject_by_annotation to True, we essentially get rid of our channels that were marked as `bad`. We assign our newly created Epochs object to epochs.

By plotting our epoched data, we can verify that our noisy segments have been removed.

epochs.plot(scalings = 'auto', n_epochs=32)

We should get the following output.

If you take a closer look, you’ll see that the noisy segments aren’t where they once were. They’re gone!

Now, given that we’re missing some time points (epochs 6,7 and 22,23) in our epoched dataset, it’s best that we join them together to recreate our continuous data.

We do so with the following code.

raw_epoched = mne.concatenate_epochs([epochs])
raw_epoched = raw_epoched.get_data()
raw_epoched = np.concatenate(raw_epoched, axis=1)

We use the function mne.concatenate_epochs() to join our epoched data (epochs) into a singular Epochs object. This is assigned to raw_epoched.

Then, we get the EEG data from raw_epoched in the form of a 3D NumPy array using get_data() and reassign it to raw_epoched.

Finally, we concatenate this 3D NumPy array along axis 1, representing the channel axis, which effectively merges all our epochs into a single continuous time series. It’s then reassigned once more to raw_epoched.

But now, our data is in the format of an Epochs object. We want it back in RawArray format which’ll allow us for easier analysis. So we convert it back to RawArray using the following code.

ch_names = ['Fz', 'Pz', 'C3', 'T3', 'C4', 'T4', 'Fp1', 'Fp2', 'F3', 'F4', 'F7', 'F8', 'P3', 'P4', 'T5', 'T6', 'O1', 'O2']
info = mne.create_info(ch_names = ch_names, sfreq=sfreq, ch_types = 'eeg')
raw = mne.io.RawArray(raw_epoched, info)

We have to update our ch_names list as we’ve removed the Cz channel from our data. Our sampling frequency (sfreq) of 128 Hz can remain as is.

Using the function mne.create_info() we input ch_names, sfreq, and our ch_types to create the metadata for our dataset. This is assigned to info.

Then we create a RawArray object by inputting our raw_epoched data and info metadata into the function mne.io.RawArray(). This newly created RawArray object is assigned to raw.

To verity that it worked, we can plot our data once again.

raw.plot(n_channels=19, scalings="auto", duration = 10.0)

You should get the following output.

If you compare this to our original plot, we’ve removed a lot of unnecessary noise.

Now, it’s time we conduct a very useful analysis method called Independent Component Analysis or ICA for short.

ICA essentially takes our recorded signal and transforms it into separate independent components that represent singular sources of data.

Using ICA will allow for a more specific form of signal analysis allowing us to remove specific sources of noise such as eye blinks, cardiac artifacts, muscular movement, and more from our data.

If you’re curious and would like to learn more about ICA, check our this article. Or this lesson by the University of Washington.
I found them very useful.

We prepare our dataset and program for ICA through the following code.

raw.set_montage('standard_1020')

raw_ica = raw.crop(tmin= 0, tmax=10)
raw_ica.plot(scalings = dict(eeg=1e-4))

Using the function RawArray.set_montage('standard_1020') we set our EEG electrode positions and digitization points to the international 10–20 system.

In essence, the international 10–20 system is an universally recognized method to apply scalp electrodes when using EEG.

Here’s what a 21 channel EEG system would look like when using the 10–20 system:

A 21 Channel EEG system using the 10–20 system.

Anyways, back to our code.

We apply the function RawArray.crop() and input our start time of 0 seconds as well as our end time of 10 seconds. This crops our dataset to only the first 10 seconds as we’re only going to analyze the first 10 seconds of our data.

Our cropped data is assigned to raw_ica and we then plot raw_ica to ensure that our data was actually cropped.

We should get the following output.

As you can see, the scrollbar at the bottom ends at the 10 second mark indicating our cropped dataset is only 10 seconds in length.

Now, we define our components and pre-whiten our data.

Essentially, pre-whitening our data equalizes the variances in signal strength to make them more uniform. It allows for ICA to identify and separate components more accurately.

n_components = 18
ica = ICA(n_components=n_components, random_state=10, max_iter = 1000)
ica.fit(raw_ica, verbose = 'ERROR')

We set the number of ICA components we want to output to 18. It can’t be greater than 18 as we’re limited to only 19 indices given that we only have 19 channels.

We use the mne.preprocessing.ICA() function to create an object which’ll take in the information needed for our dataset to identify individual components.

Random_state can be any value but the max_iter should be intentionally chosen based on your computational resources. A higher max_iter value will take up more computational resources yet will provide higher accuracy. 1000 is typically a reasonable iteration value.

Afterward, we use the ICA.fit() function to run the ICA decomposition on our dataset, raw_ica. Feel free to set the verbosity as ERROR to limit the verbosity of comments or warnings as an output.

Now, we should get the following output:

As you can see, “1000 iterations on raw data” indicates that our dataset was “Fit” through 1000 iterations using the FastICA method outputting a total of 18 independent components (IC).

Now, let’s visualize our ICs.

ica.plot_components()

The function ICA.plot_components() writes out all our ICs in plot format.

We should get the following output:

The next step in our EEG analysis is to take a look at the properties of each IC in order to decide whether to exclude them or not.

To visualize the properties of each individual IC, we can use the function ICA.plot_properties()

ica.plot_properties(raw_ica, picks = 0)

We input our dataset, raw_ica, and the individual component we want to look at through picks.

In this example, I chose the to take a look at first component at index 0.

After running the cell, I got the following output:

When taking a look at the IC, it’s fairly straightforward to tell that any noise in this IC originates from lateral eye movement, particularly on the left side.

If you look at the ERP image on the top right, you can see the noise occurs at latter stages of the 5th segment which is at about 9.5 seconds in our dataset (Each segment is 2 seconds each, hence the 5th segment represents seconds 8 to 10. Given that there’s a higher level of noise at about 75% of the 5th segment, it amounts to around the 9.5 second mark).

Given the lateral eye movement and the fact that this is the first IC, which typically represents the more variable and noisy ICs, I identified this IC for exclusion (we’ll get into how to exclude ICs next).

To take a look at the properties of each individual component and accurately identify them for exclusion, you’d have to change the picks value between 0 and 17 in the ica.plot_properties(raw_ica, picks = 0) function.

If you’re going about this project, I’d recommend you take a look at each individual IC and their properties yourself.
You might either see something I didn’t see or get a good learning experience (at least if it’s valuable to you)!
If you do so, feel free to use the following resources.
They greatly assisted me in identifying components.
They regard the use of EEGLAB rather than MNE,
but they carry the same principles which is really all that matters.
Resource 1, Resource 2, Resource 3

After looking through each individual component and their properties, I identified the ICs 0, 10, 12, and 17 for exclusion.

To exclude these components we run the following code:

ica.exclude = [0, 10, 12, 17]
raw_cleaned = ica.apply(raw_ica)

We create a list for our bad ICs and assign it to ica.exclude. Then, through the ICA.apply() function,

The ICA unmixing matrix is applied
The specified ica.exclude components are removed
Our data is remixed using a mixing matrix

Under that cell, you should get the following message:

It confirms the aforementioned steps. The ICA unmixing matrix was applied, the 4 ICA components were excluded, and the mixing matrix was reapplied to return our data to the original format which in our case, was a RawArray object.

Our newly cleaned RawArray is assigned to raw_cleaned.

Afterward, we can plot our clean data to visualize it once more.

raw_cleaned.plot(scalings=dict(eeg=1e-4))

You should the following plot:

Though, if you excluded different ICs then your data will look different.

Now, finally approaching the pinnacle of our MNE pipeline, it’s time to begin preparing our data to visualize the spectral density and then calculate the Theta / Beta ratio.

We now must epoch our data once more. Our Theta / Beta ratio will be averaged over each 2-second interval.

Though, if you’d like to have a longer or shorter interval, you can most definitely change the duration of each epoch.

epochs = mne.make_fixed_length_epochs(raw_cleaned, duration=2, preload=False, reject_by_annotation=True)
epochs.plot(scalings = "auto")

Once again, we use the function mne.make_fixed_length_epochs() and then specify our duration of 2 seconds within the function in order to epochs our data into 2 second intervals.

We once again plot our data but this time using the function Epochs.plot().

You should get the following output.

You can see dotted the black lines that denote each 2 second epoch indicating that our function worked as intended.

Now, we compute the power spectral density (PSD) in order to take a look at our data in the frequency domain which will allow us to identify which frequencies are strongest in our signal.

epochs_spectrum = epochs.compute_psd()
epochs_spectrum.plot()

We use the function Epochs.compute_psd() to get our data in the frequency domain and assign it to epochs_spectrum.

We then plot our data once more and we should get the following output:

As you can see, the frequencies <10 Hz seem to be the strongest, indicating a dominance of the delta and theta bands.

Frequencies > 10 Hz but < 35 Hz are seemingly moderate.

Frequencies > 40 Hz drop off due to the previous application of our bandpass and notch filters.

Each individual color in our spectrum is color coded to each individual channel. You can see each channel in the top right corner of our plot.

In order to view the power of each frequency relative to each channel, we can plot our PSD in the format of a scalp topography.

epochs_spectrum.plot_topomap(ch_type = 'eeg', show_names = True)

The function Spectrum.plot_topomap() plots the PSD in a scalp topography. To visualize each individual channel name, we can set show_names to True in the function.

After running the code, you should get the following plot:

You can now see the power of each individual frequency band per brain region. You’ll notice, the delta band (0–4 Hz), seems to have the strongest power as it ranges from 61248.683 µV/Hz to 230953.283 µV/Hz.

The theta band is of similar power, ranging from 41662.315 µV/Hz to 177575.821 µV/Hz.

In contrast, the beta band is nowhere near the power of delta nor theta. It merely ranges from 6104.920 µV/Hz to 15134.243 µV/Hz.

The gamma band is even lower, ranging from 768.656 µV/Hz to 3744.635 µV/Hz.

Again, if you removed different ICs, segments, or filtered the data differently, for better purposes you’d see different values.

All in all, we can see that the delta and theta bands have a higher prominence & power than the higher frequency bands of beta and gamma.

But what’s the exact ratio?

Well that’s what we’re about to find out.

raw_epoched_2 = epochs.get_data()

We use the function Epochs.get_data() to get the data from our Epochs object and transform it into a 3D NumPy array. This makes it readable for our next function. The 3D NumPy array is assigned to raw_epoched_2.

Afterward, we use Welch’s method to estimate the PSD values of our signal.

Feel free to learn about Welch’s method, from Stanford, here.
I’m still yet to fully grasp the math behind it, lol.

psds, freqs = mne.time_frequency.psd_array_welch(raw_epoched_2, sfreq = 128, fmin = 3, fmax = 40, n_fft = 256)

The function mne.time_frequency.psd_array_welch() allows us to employ Welch’s method onto our dataset, raw_epoched_2 by specifying it’s sampling frequency, the minimum and maximum Hz we want, and the window size of each Fast Fourier Transform (FFT) computation or bin.

For context, FFT is the mathematical method used to transform a signal from the time domain into the frequency domain

In our case, our sampling frequency is 128, the minimum Hz we want is 4 Hz as that marks the lower end of the theta band, and the maximum Hz we want is 30 as that marks the upper end of the beta band.

We set the n_fft value to 256 as our window size can’t be more than twice our sampling frequency.

The PSD calculation then returns psds and freqs which is then assigned to their corresponding variables. psds will hold the power spectral densities and freqs will hold the frequencies.

So now, we can define the range of the theta and beta bands.

fmin_theta, fmax_theta = 4, 8
fmin_beta, fmax_beta = 13, 30 

theta_indices = (freqs >= fmin_theta) & (freqs <= fmax_theta)
beta_indices = (freqs >= fmin_beta) & (freqs <= fmax_beta)

The theta band ranges from 4–8 Hz and the beta band ranges from 13–30 Hz. Hence, in our code we assign the fmin_theta to 4, the fmax_theta to 8, the fmin_beta to 13, and the fmax_beta to 30.

So then the limits of the theta band is assigned to a boolean array, theta_indices, and the limits of the beta band is assigned to a boolean array, beta_indices.

Afterward, we get the average theta and beta values from our psds array per 2 second epoch.

psd_theta = psds[:, :, theta_indices].mean(axis=-1)
psd_beta = psds[:, :, beta_indices].mean(axis=-1)

We look at the 3rd dimension of our 3D psds array as that holds the frequencies. We use the boolean array theta_indices / beta_indices to then extract the theta and beta bands from the psds array.

Once done, .mean(axis=-1) calculates the average of of the extracted theta and beta bands to along the last axis or the 3rd dimension of the 3D array which holds the frequencies.

Both the average of the theta and beta bands are assigned to the corresponding psd_theta and psd_beta.

Now we calculate the EEG ratio.

theta_beta_ratio = psd_theta / psd_beta

It’s a simple division.

Once calculated, we print the data out.

print(theta_beta_ratio)

You should get the following output.

[[14.24585974 17.00854954  8.73644889 18.66439298 10.35378318 24.80291788
  15.53794867 23.18015227 10.19410801 22.81926067 18.47482058 15.85629424
   8.02090181 14.65767321  6.88291899 10.25030782 14.01618669 34.09462749]
 
[14.14268983 23.74512228 13.7848695  18.73838238 11.57862147 16.02754893
  10.77896311 15.59140046  6.46788619 12.07241129  3.1350443  16.57906325
   1.21217917 10.00192009  1.13746622 18.61028476 17.29716471 13.4303745 ]
 
[ 9.72703461 11.6778984  11.86698283  9.20870894  9.83779854 11.52738504
   9.04488704 17.76950501  5.22511413  8.59557069  8.53960467 11.44716499
   1.57528786  9.28832179  1.45063329 16.27778208 16.73691082 15.04533255]
 
[22.83171869 27.73154983 10.6697696  24.90448829 11.33849498 12.91807595
   8.56041358 10.12631801  5.96810591 18.87360354  6.27188118 14.98131042
   8.51126961  7.17636289  5.83429323 21.62826478  8.65682859 10.66859703]
 
[ 5.88553342 10.21934253  8.70449541  9.97584756  5.40288631 16.34475158
   8.43228802 18.31463182  5.86610279  4.97029056  7.12865796  9.07346259
   1.24249726  8.00679333  1.26999772  7.76765771  7.34158132 17.44056921]]

If you examine the output, the theta / beta ratio is split into 5 lists which indicate the five 2-second epochs. Within each list you can see the average theta / beta ratio for each channel in order per 2 second segments.

For example, 14.24585974 represents the average theta / beta for channel Fz.

To make our data more readable we can separate the data per epoch and into singular channels.

For my analysis, I wanted to extract the data from the mid-frontal electrodes, Fz, Fp1, Fp2, F4, F7, and F8 and average all channels out per epoch.

midfrontal_ch_indices = [0, 6, 7, 8, 10, 11, 12] 
midfrontal_ch_names = ['Fz', 'Fp1', 'Fp2', 'F3,' 'F4', 'F7', 'F8']

midfrontal_ch_data = theta_beta_ratio[:, selected_channels_indices]

for i, channel_data in enumerate(midfrontal_ch_data):
    print(f"EPOCH {i + 1} | {', '.join(f'{name}: {value:.5f}' for name, value in zip(midfrontal_ch_names, channel_data))}")

We set the indices of our channel to a new list, midfrontal_ch_indices in reference to our original ch_names list. Our channel names are put in a list and assigned to midfrontal_ch_names.

We then take our theta_beta_ratio 2D array and extract data from the selected channels.

Then, for each Epoch of data, we print the average theta / beta ratio of each channel.

We should get the following output:

EPOCH 1 | Fz: 14.24586, Fp1: 15.53795, Fp2: 23.18015, F3,F4: 18.47482, F7: 15.85629, F8: 8.02090
EPOCH 2 | Fz: 14.14269, Fp1: 10.77896, Fp2: 15.59140, F3,F4: 3.13504, F7: 16.57906, F8: 1.21218
EPOCH 3 | Fz: 9.72703, Fp1: 9.04489, Fp2: 17.76951, F3,F4: 8.53960, F7: 11.44716, F8: 1.57529
EPOCH 4 | Fz: 22.83172, Fp1: 8.56041, Fp2: 10.12632, F3,F4: 6.27188, F7: 14.98131, F8: 8.51127
EPOCH 5 | Fz: 5.88553, Fp1: 8.43229, Fp2: 18.31463, F3,F4: 7.12866, F7: 9.07346, F8: 1.24250

Now, we know the average of the midfrontal channels which will prove to be useful for interpreting our data.

I also wanted to know the average of all the mid-frontal channels combined per epoch.

So I used the following code to calculate it.

avg_midfrontal = midfrontal_ch_data.mean(axis=1)

for i, average_value in enumerate(avg_midfrontal):
    print(f"EPOCH {i + 1} | Average Midfrontal: {average_value:.5f}")

The mean of our midfrontal_ch_data is taken and assigned to avg_midfrontal. Then, for each Epoch of data, we print out the average theta / beta ratio of all midfrontal channels.

We should get the following output:

EPOCH 1 | Average Midfrontal: 15.88600
EPOCH 2 | Average Midfrontal: 10.23989
EPOCH 3 | Average Midfrontal: 9.68391
EPOCH 4 | Average Midfrontal: 11.88049
EPOCH 5 | Average Midfrontal: 8.34618

Now, this provides a more general overview of our data. While it may not be as specific, averaging the mid-frontal channels allows for a simplified way of conducting our analysis.

Both, the average and the individual theta / beta ratios are extremely useful.

So, regarding the actual code of our pipeline, we’re finished!

Now, it’s time we take a look at our theta / beta values and figure out what they really mean.

Characterizing Selective Attention and Cognitive Processing

There’s a very good reason I chose to extract data from the mid-frontal EEG electrodes.

A study conducted by Guo et al., found that an elevated theta band in the mid-frontal region of children, characterizes ADHD and diminished levels of selective attention¹².

The present study provides a novel evidence that children with ADHD show overelevated midfrontal theta ERS and abnormal posterior theta modulation. The elevated theta ERS may reflect a lower level of attention arousal and compensatory cognitive processing related to the poorer behavioral performance in children with ADHD.
The “adult-like” posterior theta modulation observed in children with ADHD may reflect compensatory maturation for attention deficits that promotes the stability of their behavioral performance.
— Guo et al.

What’s interesting, is that as a result of higher levels of theta in the midfrontal region indicating a lowered level selective attention, there appeared to be compensatory cognitive processing in the posterior region of the brain that reflected a level of cognitive processing of a grown adult¹². This allowed for their behavioral performance to remain “stable” ¹².

Interestingly, it appears that the cognitive functions of their brain have a certain degree of antifragility.

CC: Nassim Nicholas Taleb

Anyways…

If we take a look at our data,

EPOCH 1 | Fz: 14.24586, Fp1: 15.53795, Fp2: 23.18015, F3,F4: 18.47482, F7: 15.85629, F8: 8.02090
EPOCH 2 | Fz: 14.14269, Fp1: 10.77896, Fp2: 15.59140, F3,F4: 3.13504, F7: 16.57906, F8: 1.21218
EPOCH 3 | Fz: 9.72703, Fp1: 9.04489, Fp2: 17.76951, F3,F4: 8.53960, F7: 11.44716, F8: 1.57529
EPOCH 4 | Fz: 22.83172, Fp1: 8.56041, Fp2: 10.12632, F3,F4: 6.27188, F7: 14.98131, F8: 8.51127
EPOCH 5 | Fz: 5.88553, Fp1: 8.43229, Fp2: 18.31463, F3,F4: 7.12866, F7: 9.07346, F8: 1.24250

It appears that the first 2 Epochs of data have a slightly higher theta / beta ratio than the latter 3.

If we take a look at our averaged data,

EPOCH 1 | Average Midfrontal: 15.88600
EPOCH 2 | Average Midfrontal: 10.23989
EPOCH 3 | Average Midfrontal: 9.68391
EPOCH 4 | Average Midfrontal: 11.88049
EPOCH 5 | Average Midfrontal: 8.34618

we can confirm that Epoch 1 and 2 both have higher theta / beta ratios respectively at 15.886 and 10.239

It seems that the latter 3 Epochs begin to drop off to 8.346.

What does this mean?

Well, based on the 10 seconds of analyzed data and the current literature, we can infer that at the beginning of our data, in the first 4 seconds, the subject exhibits a lower level of selective attention.

Through the final 6 seconds, theta levels begin to drop off in favor of beta which indicate a higher level of selective attention.

So what now?

Now that I have a deeper understanding in signal processing and MNE, the next step would be to develop a similar pipeline that’s able to take in live EEG data, clean that EEG data, and output a metric to characterize ADHD-like electrophysiology.

All in real-time.

Perhaps maybe then, we’ll be a step closer to building a small yet highly effective wearable device that’ll let ADHD become more manageable on a day-to-day basis.

I’m also really curious to figure out and see for myself, how machine learning, deep learning, and / or neural networks can play a role in signal processing and neurofeedback.

It’s what I intend to explore very soon.