# Decomposing an EEG Signal from Scratch

**URL:** https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904
**Category:** General Usage
**Tags:** visualization, juliahealth, dsp
**Created:** [June 20, 2024, 6:05am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904 "2024-06-20T06:05:53Z")
**Posts on this page:** 20
**Page:** 1

<div class="post-metadata">

### Author: ![TheCedarPrince](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thecedarprince/32/17323_2.png) [@TheCedarPrince](https://discourse.julialang.org/u/TheCedarPrince)
#### Post date: [June 20, 2024, 6:05am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/1 "2024-06-20T06:05:54Z")

</div>

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](https://openneuro.org/datasets/ds004902/versions/1.0.5)). 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**
>
> ```julia
> 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**
>
> ```julia
> 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):

 ![display](https://global.discourse-cdn.com/julialang/original/3X/8/0/80b3169b93d26d4b923f0ae81bcd08520008c23e.png)

> **Topoplot Code**
>
> ```julia
> 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):

 ![display](https://global.discourse-cdn.com/julialang/original/3X/6/6/661bc30df90a044a98984c339fa715db3d7fd35d.png)

> **Raw Voltage Plot**
>
> ```julia
> 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:

```julia
s = welch_pgram(channel_data[1, :], fs = 500)

f = s.freq
p = s.power

```

Then, taking liberal usage of [Prof. Makoto Miyakoshi’s](https://sccn.ucsd.edu/wiki/Makoto%27s_useful_EEGLAB_code) 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:

```julia
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:

 ![display](https://global.discourse-cdn.com/julialang/original/3X/8/e/8e67d928d948f52419d7a1577a9b11cdc9945593.png)

> **Segmented Wave Bands**
>
> ```julia
> 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):

 ![display](https://global.discourse-cdn.com/julialang/original/3X/2/3/23b5224fd0a454b5f434275c6cd2b0aac7d77c68.png)

> **Abbreviated, Banded PSD Plot**
>
> ```julia
> 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](https://eeglab.org) or Python’s [MNE-Python](https://mne.tools/stable/index.html)). 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:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/1/d1b93b5a43423d83414d29c270f63c2cbac99b8d.png)

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 🌳

---

<div class="post-metadata">

### Author: ![TheCedarPrince](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thecedarprince/32/17323_2.png) [@TheCedarPrince](https://discourse.julialang.org/u/TheCedarPrince)
#### Post date: [June 20, 2024, 6:10am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/2 "2024-06-20T06:10:44Z")

</div>

P.S. I do not mean to be a bother but am going to CC some folks who I’d be curious to hear from and to hear how you tackle problems like these: @Zach_Christensen @Jakub_Mitura @tim.holy @jcorream @AdamWysokinski

---

<div class="post-metadata">

### Author: ![Jakub\_Mitura](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_mitura/32/19496_2.png) [@Jakub\_Mitura](https://discourse.julialang.org/u/Jakub_Mitura)
#### Post date: [June 20, 2024, 6:23am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/3 "2024-06-20T06:23:19Z")

</div>

I had not tried but knowing the location of the electrodes one should be able to triangulate the position to some area of the brain (not precisely without MRI of the patient) still relative amplitude of the extracted event from diffrent spots should give the ability to show on a template brain model where is the source of given electrical activity

---

<div class="post-metadata">

### Author: ![mkoculak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkoculak/32/28310_2.png) [@mkoculak](https://discourse.julialang.org/u/mkoculak)
#### Post date: [June 20, 2024, 11:30am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/4 "2024-06-20T11:30:31Z")

</div>

Hey,  
what you did looks fine, but I think it is not exactly what you wanted to get. 😉  
You have extracted the spectral representation, however from the graphic in the first point it seems you want timeseries.  
So just as you have mentioned - you need to bandpass filter the data with bounds matched to different types of waves. Can you share the code for filtering that you used and the error you got?

As for the third question, I am not sure what are you referring to. Do you mean taking an average of the raw signal? Or average of some metric? Typically people calculate metrics over smaller segments of data (e.g. couple of seconds) and then average the values from segments. It will however vary depending on the chosen measure. Averaging raw signal is not a typical practice.

---

<div class="post-metadata">

### Author: ![jcorream](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jcorream/32/46276_2.png) [@jcorream](https://discourse.julialang.org/u/jcorream)
#### Post date: [June 20, 2024, 1:05pm UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/5 "2024-06-20T13:05:58Z")

</div>

Hi! Happy to chime in.  
Yes, the frequency ranges that you’ve used are consistent with the literature. See Herrmann, C.S., Strüber, D., Helfrich, R.F. and Engel, A.K., 2016. EEG oscillations: from correlation to causality. _International Journal of Psychophysiology_, _103_, pp.12-21.

The extent to which the oscillations are extracted accurately in the time-domain depends on the properties of the filter that you design. One way to assess this is to look at the frequency response of the filter and see whether there’s good attenuation outside of the passband of interest. See the DSP documentation here and [tutorial](https://docs.juliadsp.org/v0.4/filters.html) at the end. There are of course many other criteria. Oppenheim’s Digital Signal Processing book has more information on this. The [MIT OCW](https://ocw.mit.edu/courses/res-6-008-digital-signal-processing-spring-2011/pages/study-materials/) on his DSP class has excellent resources as well.

To compute the spectrum, I recommend using Multitaper methods. It’ll reduce the variance of your spectral estimate and uses tapers that maximize the average power within a passband of interest. There’s a Julia [package](https://docs.juliahub.com/Multitaper/OT9LO/0.2.0/) for doing this.

---

<div class="post-metadata">

### Author: ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [June 20, 2024, 1:36pm UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/6 "2024-06-20T13:36:51Z")

</div>

Just to add [DSP.jl has multitaper methods](https://docs.juliadsp.org/stable/periodograms/#Multitaper-periodogram-estimation) as well, and they are reference-tested against MATLAB and PyMNE. I actually hadn’t heard of Multitaper.jl so I’m not sure how DSP’s methods compare to those.

---

<div class="post-metadata">

### Author: ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [June 20, 2024, 1:36pm UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/7 "2024-06-20T13:36:58Z")

</div>

shameless plug here, sorry for being off-topic! If some new methods for decomposing signals are discussed here, perhaps they can be added via PR to [GitHub - JuliaDynamics/SignalDecomposition.jl: Decompose a signal/timeseries into structure and noise or seasonal and residual components](https://github.com/JuliaDynamics/SignalDecomposition.jl) or open an issue mentioning the process so that at least it is there in the library sense!

edit: right, after reading the original post in more detail, I can see that I am _way too off topic_. This post is about decomposition in frequency domain while what I mentioned above is decomposition in time (albeit they are conjugate of each other but oh well).

---

<div class="post-metadata">

### Author: ![jcorream](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jcorream/32/46276_2.png) [@jcorream](https://discourse.julialang.org/u/jcorream)
#### Post date: [June 20, 2024, 1:49pm UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/8 "2024-06-20T13:49:30Z")

</div>

Right! Thanks. Forgot to mention that DPS.jl also has multitaper methods. 🙂 I am also unsure as to how they compare to one another.

---

<div class="post-metadata">

### Author: ![AdamWysokinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adamwysokinski/32/206422_2.png) [@AdamWysokinski](https://discourse.julialang.org/u/AdamWysokinski)
#### Post date: [June 20, 2024, 6:36pm UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/9 "2024-06-20T18:36:22Z")

</div>

Hi all,  
I’m very busy for the next few days, so please excuse me for my very brief reply. All these operations (and much more 😀) can be done using [NeuroAnalyzer.jl](https://codeberg.org/AdamWysokinski/NeuroAnalyzer.jl). I’ll prepare a brief tutorial on this problem, just give me day or two 🙂  
PS 1. Your plot is the power spectrum, the reference plot is a 5-second signal segment split into 5 frequency bands, e.g. using a band-pass filter for each band.  
PS 2. Delta band is usually defined from 0.1 Hz or 0.5 Hz (e.g. in sleep studies).  
ATB,  
Adam

---

<div class="post-metadata">

### Author: ![AdamWysokinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adamwysokinski/32/206422_2.png) [@AdamWysokinski](https://discourse.julialang.org/u/AdamWysokinski)
#### Post date: [June 20, 2024, 6:46pm UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/10 "2024-06-20T18:46:06Z")

</div>

Unfortunately it’s much more difficult since source localization is an undetermined problem. For each set of channel recordings there is an infinite number of solutions for their sources.

---

<div class="post-metadata">

### Author: ![TheCedarPrince](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thecedarprince/32/17323_2.png) [@TheCedarPrince](https://discourse.julialang.org/u/TheCedarPrince)
#### Post date: [June 21, 2024, 4:07am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/11 "2024-06-21T04:07:28Z")

</div>

Hey @mkoculak,

> [@mkoculak](#):
>
> So just as you have mentioned - you need to bandpass filter the data with bounds matched to different types of waves. Can you share the code for filtering that you used and the error you got?

Thanks for the comments! Sure thing, here is the code I was trying to use:

```julia
using DSP

#=

...
Code for loading data here
...

=#

fs = 500 # Sampling frequency

Wp = (4, 8) ./ (fs / 2);
Ws = (3.5, 8.5) ./ (fs / 2);
Rp = 3
Rs = 40

n, Wn = buttord(Wp, Ws, Rp, Rs)

theta_band = Bandpass(4, 8; fs = fs);
theta_filter = digitalfilter(theta_band, Butterworth(n));

z = theta_filter.z
p = theta_filter.p
k = theta_filter.k

```

After this point, I was a bit confused on how to properly apply the filter to my data as well as what form I should have my data in (i.e. \mu V or something else). Do you know what I might want to do after this point if what I did seemed reasonable? I was using this discussion to try to track down how to do this: [matlab - EEG bandpass filter in mat lab - Stack Overflow](https://stackoverflow.com/questions/23664631/eeg-bandpass-filter-in-mat-lab)

> [@mkoculak](#):
>
> Averaging raw signal is not a typical practice.

Yea, I was more referring to a processed signal (i.e. after some denoising and associated filtering had been done).

P.S. You may want to see the following comment now as I retried using DSP.jl and I think I may have gotten something more like what I wanted for question 1.

---

<div class="post-metadata">

### Author: ![TheCedarPrince](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thecedarprince/32/17323_2.png) [@TheCedarPrince](https://discourse.julialang.org/u/TheCedarPrince)
#### Post date: [June 21, 2024, 4:52am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/12 "2024-06-21T04:52:26Z")

</div>

> [@jcorream](#):
>
> See the DSP documentation here and [tutorial](https://docs.juliadsp.org/v0.4/filters.html) at the end. There are of course many other criteria.

After refreshing myself with a bit more EEG DSP material, I took another look at this tutorial you mentioned and things started clicking! After playing around with the tutorial some more, I _think_ I may be getting closer to what I am after. Here’s the latest diagram:

 ![display](https://global.discourse-cdn.com/julialang/original/3X/e/e/ee8624775cd5fe34609e59ad6ee79c033c332767.jpeg)

Here’s my associated code where I am using DSP.jl and processing my raw channel data (in \mu V):

> **Brain Wave Extraction Attempt Code**
>
> ```julia
> delta_range = [1, 4]
> theta_range = [4, 8]
> alpha_range = [8, 12]
> beta_range = [12, 30]
> gamma_range = [30, 200]
> 
> fig = Figure(
> size = (1200, 1400), 
> fontsize = 20
> );
> 
> for (idx, r) in enumerate([delta_range, theta_range, alpha_range, beta_range, gamma_range])
> if r[1] > 0
> Ws = (r[1] - 0.5, r[2] + .5)./(fs/2)
> else
> Ws = (0.5, r[2] + .5)./(fs/2)
> end
> Wp = (r[1], r[2])./(fs/2)
> Rp = 3 # Passband ripple
> Rs = 40 # Stopband attenuation
> 
> n, ωn = buttord(Wp, Ws, Rp, Rs)
> 
> ax = CairoMakie.Axis(
> fig[idx, 1], 
> )
> 
> btr_band_filt = digitalfilter(Bandpass(r[1], r[2]; fs = fs), Butterworth(n))
> signal = filt(btr_band_filt, channel_data[1, :]);
> 
> lines!(ax, 1:5000, signal[1001:6000]);
> ax.title = L"\%$(labels[idx])\text{-Band}"
> ax.ylabel = "Frequency"
> if idx == 5
> ax.xlabel = "Time"
> end
> end
> 
> ```

Additionally, here is information about my estimated order (n) for my Butterworth Filter and other information about the signal I am showing here:

> **Variable Information**
>
> ```plaintext
> delta Band Information:
> Wp = (0.004, 0.016) 
> Ws = (0.002, 0.018) 
> n = 25 
> 
>     
> theta Band Information:
> Wp = (0.016, 0.032) 
> Ws = (0.014, 0.034) 
> n = 28 
> 
>     
> alpha Band Information:
> Wp = (0.032, 0.048) 
> Ws = (0.03, 0.05) 
> n = 25 
> 
>     
> beta Band Information:
> Wp = (0.048, 0.12) 
> Ws = (0.046, 0.122) 
> n = 120 
> 
>     
> gamma Band Information:
> Wp = (0.12, 0.8) 
> Ws = (0.118, 0.802) 
> n = 380 
> 
> ```

The signals for the brain waves look almost like what I’d expect but then gamma is just obliterated by the high order estimate of n = 380. I am just not sure if this is exactly right but many thanks @jcorream – I might be getting… Closer?

---

<div class="post-metadata">

### Author: ![TheCedarPrince](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thecedarprince/32/17323_2.png) [@TheCedarPrince](https://discourse.julialang.org/u/TheCedarPrince)
#### Post date: [June 21, 2024, 4:55am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/13 "2024-06-21T04:55:53Z")

</div>

I had no idea this was the case and makes using DSP.jl even more salient in my eyes! Could I add a quick docs PR to highlight this fact? I didn’t come across this fact anywhere outside of your reference.

---

<div class="post-metadata">

### Author: ![TheCedarPrince](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thecedarprince/32/17323_2.png) [@TheCedarPrince](https://discourse.julialang.org/u/TheCedarPrince)
#### Post date: [June 21, 2024, 4:57am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/14 "2024-06-21T04:57:53Z")

</div>

Oh thanks for the comments – looking forward to seeing more about NeuroAnalyzer.jl 's usage! I think it is a lovely package and am eager to hear more soon! Hope you are doing great Adam (and hope to see you at a JuliaHealth meeting again sometime! 😃 )! Take care!

---

<div class="post-metadata">

### Author: ![Jakub\_Mitura](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_mitura/32/19496_2.png) [@Jakub\_Mitura](https://discourse.julialang.org/u/Jakub_Mitura)
#### Post date: [June 21, 2024, 6:10am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/15 "2024-06-21T06:10:57Z")

</div>

thanks for correcting me here ! I was extrapolating from ECG, where you can read on the basis of electrode tell quite a lot about spatial localisation of anomaly.

---

<div class="post-metadata">

### Author: ![AdamWysokinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adamwysokinski/32/206422_2.png) [@AdamWysokinski](https://discourse.julialang.org/u/AdamWysokinski)
#### Post date: [June 21, 2024, 7:46am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/16 "2024-06-21T07:46:33Z")

</div>

I’ve discovered a nasty bug in the .FDT import routine, it is fixed, but please use the current NeuroAnalyzer version (0.24.7-dev):

```sh
git clone https://codeberg.org/AdamWysokinski/NeuroAnalyzer.jl
cd NeuroAnalyzer.jl
julia

```

Ok, briefly, here’s the code:

```julia
using Pkg
using Plots

Pkg.activate(@ __DIR__ )
Pkg.instantiate()
using NeuroAnalyzer

# load the data
eeg = import_set("sub-01_ses-1_task-eyesopen_eeg.set")

# there is a huge low-frequency drift, we should remove it
NeuroAnalyzer.filter!(eeg, fprototype=:fir, ftype=:hp, cutoff=0.5)

# split into bands
# s is the split signal (band × channel × samples × epochs), bn is the bands names and bf is their boundary frequencies
s, bn, bf = bpsplit(eeg)

# let's plot the second channel
# we will plot the first 20 seconds (`1:20*sr(eeg)`)
ch = 2
# time points of the segment
t = eeg.time_pts[1:20*sr(eeg)]
# plot individual bands
p_d = Plots.plot(t, s[1, ch, 1:20*sr(eeg), 1], ylims=(-15, 15), title="Band: $(string(bn[1]))", legend=false, xlabel="[s]", ylabel="[μV]")
p_t = Plots.plot(t, s[2, ch, 1:20*sr(eeg), 1], ylims=(-15, 15), title="Band: $(string(bn[2]))", legend=false, xlabel="[s]", ylabel="[μV]")
p_a = Plots.plot(t, s[3, ch, 1:20*sr(eeg), 1], ylims=(-15, 15), title="Band: $(string(bn[3]))", legend=false, xlabel="[s]", ylabel="[μV]")
p_b = Plots.plot(t, s[6, ch, 1:20*sr(eeg), 1], ylims=(-15, 15), title="Band: $(string(bn[6]))", legend=false, xlabel="[s]", ylabel="[μV]")
p_g = Plots.plot(t, s[9, ch, 1:20*sr(eeg), 1], ylims=(-15, 15), title="Band: $(string(bn[9]))", legend=false, xlabel="[s]", ylabel="[μV]")
# combine all plots
Plots.plot(p_d, p_t, p_a, p_b, p_g, layout=(5, 1), labelfontsize=5, titlefontsize=6, ytickfontsize=4, xtickfontsize=5, lw=0.25)

```

To plot the topographical map:

```julia
# first, load channel locations
load_locs!(eeg, file_name="sub-01_ses-1_electrodes.tsv")
# plot the topomap: averge amplitude of the first 20 seconds
plot_topo(eeg, seg=(0, 20))

```

 ![topo](https://global.discourse-cdn.com/julialang/original/3X/1/8/18359fd08a5cc5026a04a38a36ab7a89f7234d6e.png)

ATB, Adam

---

<div class="post-metadata">

### Author: ![AdamWysokinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adamwysokinski/32/206422_2.png) [@AdamWysokinski](https://discourse.julialang.org/u/AdamWysokinski)
#### Post date: [June 21, 2024, 7:53am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/17 "2024-06-21T07:53:25Z")

</div>

The major difference in source localization between ECG and EEG is that in the heart you have a periodical, sequential activity (P peak P-Q segment, QRS peak, S-T segment, T peak and back to P peak). Each part represents electrophysiological activity of a distinctive part of the heart (atria, ventricles, etc.) Therefore, you may map these elements onto the heart structures when you have them recorded at various positions. The brain activity is at the macroscopic level aperiodical and non-sequential, hence the problem.

---

<div class="post-metadata">

### Author: ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [June 21, 2024, 1:14pm UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/18 "2024-06-21T13:14:11Z")

</div>

Sure, doc improvements are of course always welcome! I only know about it because I helped upstream it from some private code. Multitaper methods seem somewhat niche outside of EEG analysis though so it probably shouldn’t really be much more prominent than the standard stuff most folks are probably looking for.

---

<div class="post-metadata">

### Author: ![TheCedarPrince](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thecedarprince/32/17323_2.png) [@TheCedarPrince](https://discourse.julialang.org/u/TheCedarPrince)
#### Post date: [June 21, 2024, 3:17pm UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/19 "2024-06-21T15:17:10Z")

</div>

Added a wee little doc PR here: [[DOCS] Add Small Note about Method Validations by TheCedarPrince · Pull Request #561 · JuliaDSP/DSP.jl · GitHub](https://github.com/JuliaDSP/DSP.jl/pull/561)

Might not make sense in that location so I could move this notice around depending.

---

<div class="post-metadata">

### Author: ![Emmanuel-R8](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/emmanuel-r8/32/11839_2.png) [@Emmanuel-R8](https://discourse.julialang.org/u/Emmanuel-R8)
#### Post date: [June 22, 2024, 10:31am UTC](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904/20 "2024-06-22T10:31:36Z")

</div>

Out of curiosity, I guess Independent Component Analysis has been tried in the literature. Is that useful at all? What would be its interpretation?

[Next page](https://discourse.julialang.org/t/decomposing-an-eeg-signal-from-scratch/115904.md?page=2)
