Hey folks! 
I’ve been recently tinkering with decomposing some EEG data into various visualizations and want to make sure I am going about things the right way (the data using is openly available from here). Here is how I approached processing some of this EEG data as well as questions I had at the end.
Minor Data Prep
First, as the data was in fdt files, I wrote my own parser to read it (this is where raw channel level data comes from):
fdt Parser
fdt_parser(fdt_file, mat_file) =
open(fdt_file) do io
read!(
io,
Array{Float32}(undef, (Int(mat_file["nbchan"]), Int(mat_file["pnts"]))),
)
end
Then, I loaded up the .set and .fdt files to have as a reference for the session recording that I was analyzing. This was pretty straightforward:
Loading Data
using MAT
mat_file = matread(
joinpath(
session_1_path,
"sub-01_ses-1_task-eyesclosed_eeg.set"
)
)
channel_data = fdt_parser(
joinpath(
session_1_path,
"sub-01_ses-1_task-eyesclosed_eeg.fdt"
),
mat_file)
Visualizing Electrode Positions
This session recorded EEG signal across 61 different channels (i.e. electrodes). For a quick visualization of that, I used TopoPlot.jl by @behinger which gave the following photo (code is below the photo):
Topoplot Code
fig = eeg_topoplot(
data[1][1:61, 1, 1],
electrodes.name,
positions = pts,
interpolation = NullInterpolator(),
label_scatter =
(
markersize=8,
color = :black
),
axis =
(
type=Axis,
title = "61-Lead Electrode Array Layout",
aspect= DataAspect(),
leftspinevisible = false,
rightspinevisible = false,
topspinevisible = false,
bottomspinevisible = false,
xgridvisible = false,
ygridvisible = false,
xminorgridvisible = false,
yminorgridvisible = false,
xminorticksvisible = false,
yminorticksvisible = false,
xticksvisible = false,
yticksvisible = false,
xticklabelsvisible = false,
yticklabelsvisible = false,
),
labels = electrodes.name,
label_text =
(;
fontsize = 6,
offset = (-3, 0)
),
)
Quick Data Preview
For the next parts of this work, I decided to only use 1 channel for my analysis as I thought I could scale up complexity later (i.e. analyzing and averages all channels together later). To get a sense of the signal on this channel, I did a quick plot of the voltage amplitude recorded (measured in \mu V); here’s that plot (with code below):
Raw Voltage Plot
fig = Figure(
size = (1200, 500),
fontsize = 20
);
ax = CairoMakie.Axis(
fig[1, 1],
)
lines!(ax, 1:length(channel_data[1, :]), channel_data[1, :])
ax.xlabel = "Samples"
ax.ylabel = "Voltage"
ax.title = "Raw Signal (Freq: 500 Hz)"
Signal Decomposition
After this, I get to the part that is the most challenging for me – decomposing the signal into different bands (i.e. delta, theta, alpha, beta, and gamma). This part is where I could use the most checking. How I first went about this was to use DSP.jl to compute a Welsh Periodogram of the time series:
s = welch_pgram(channel_data[1, :], fs = 500)
f = s.freq
p = s.power
Then, taking liberal usage of Prof. Makoto Miyakoshi’s notes, I think I somewhat figured out how to extract the bands after translating some MATLAB ideas into Julia. Here is how I approached it:
delta = findall(x -> x>=0 && x<4, f);
theta = findall(x -> x>=4 && x<8, f);
alpha = findall(x -> x>=8 && x<12, f);
beta = findall(x -> x>=12 && x<30, f);
gamma = findall(x -> x>=30 && x<200, f);
And with these parts segmented out, I created the following plot here:
Segmented Wave Bands
fig = Figure(
size = (1000, 1200),
fontsize = 20
);
colors = ["red", "green", "blue", "orange", "black"]
labels = ["delta", "theta", "alpha", "beta", "gamma"]
bands = [delta[50:end], theta, alpha, beta, gamma]
for (idx, band) in enumerate(bands)
ax = CairoMakie.Axis(
fig[idx, 1],
)
lines!(ax, f[band], p[band], color = colors[idx])
ax.title = L"\%$(labels[idx])\text{-Band}"
ax.ylabel = "Power"
if idx == 5
ax.xlabel = "Frequency (Hz)"
end
end
Additionally, I created a Power Spectral Density plot as follows for additional visualization (with some aesthetic banding for easier viewing):
Abbreviated, Banded PSD Plot
fig = Figure(
size = (1600, 400),
fontsize = 20
);
ax = CairoMakie.Axis(
fig[1, 1],
)
colors = ["red", "green", "blue", "orange", "black"]
labels = ["delta", "theta", "alpha", "beta", "gamma"]
bands = [delta[50:end], theta, alpha, beta, gamma[1:800]]
band_ymax = p[vcat(bands...)] |> maximum
for (idx, band) in enumerate(bands)
band!(ax, f[band], fill(0, length(band)), fill(band_ymax, length(band)), color = (colors[idx], 0.1))
lines!(ax, f[band], p[band], color = colors[idx], label = L"\%$(labels[idx])")
ax.title = "EEG Recording Session"
ax.ylabel = "Power"
if idx == 5
ax.xlabel = "Frequency (Hz)"
axislegend()
end
end
Note: I cut off the gamma band in this particular diagram prematurely as it resulted in a rather “scrunched up” graphic. This was purely for presentation purposes.
Outstanding Questions
I hope that more or less made sense so far! It was somewhat tough tracking down these methods outside of established software (like MATLAB’s EEGLAB or Python’s MNE-Python). I definitely could have made errors so far and that leads me to outstanding questions I had from this work so far:
First, in the plot where I segmented out the different regions, I am unsure how to make a visualization like this:
As one can see, the different bands of one signal are all somehow harmonized against the same x-axis. My segmented one is not. How can I make a plot like this?
Second, are the methods I have used so far to extract the bands correct? I tried creating a nice bandpass filter with DSP.jl but couldn’t get it to work so just opted to directly use the Welch Periodogram function that comes with DSP.jl instead. Would there be a better way to do things here?
Third, is it often standard practice to take the mean of signals across all channels over time when computing brainwaves? In my readings, it seems like that is one way to do it to aggregate occurrences across the entire recording during a session but was curious what else folks do.
Thanks all!
~ tcp 
I am also unsure as to how they compare to one another.
) can be done using 
)! Take care!