FOOOF

using NeuroAnalyzer
using CairoMakie
eeg = load("files/eeg.hdf")
e10 = epoch(eeg, ep_len = 10)

FOOOF is the algorithm for parameterizing neural power spectra into aperiodic and periodic components.

NeuroAnalyzer re-implements the fitting algorithm of fooof 1.1.0 / specparam 2.0.0rc (the core algorithm and its default settings are identical in both), with these differences:

• curve fitting uses an internal box-constrained Levenberg-Marquardt solver instead of scipy.optimize.curve_fit; results agree with the Python reference to within fitting tolerance, but are not bit-identical
• for aperiodic_mode = :fixed the knee is reported as 0 (the fixed model is the knee model with knee = 0), so that the output always has 3 aperiodic parameters

For details of the FOOOF method, please visit the project documentation.

FOOOF analysis in NeuroAnalyzer, analyze channel F3, 10th epoch, in the 1 to 40 Hz range; detect up to 2 peaks:

ep = 10
psd_data = psd(e10, ch = "F3", db = false)
frq = psd_data.f
f_idx1 = vsearch(1, frq)
f_idx2 = vsearch(40, frq)
frq = frq[f_idx1:f_idx2]
pow = log10.(psd_data.p[1, f_idx1:f_idx2, ep])
f = fooof(e10, ch = "F3", flim=(1, 40), max_n_peaks = 2)

(!) For both operations the default method is :welch.

Aperiodic parameters:

println("Offset: ", round(f.ap[1, :, ep][1], digits = 3))
println("Knee: ", round(f.ap[1, :, ep][2], digits = 3))
println("Exponent: ", round(f.ap[1, :, ep][3], digits = 3))

Offset: 0.093
Knee: 0.0
Exponent: 0.99

Peaks:

peaks = f.pe[ep]
println("Number of peaks: ", size(peaks, 1))
for idx in axes(peaks, 1)
    println("Peak $idx: CF: ", round(peaks[idx, 1], digits = 3),
     " PW: ", round(peaks[idx, 2], digits = 3), 
     " BW: ", round(peaks[idx, 3], digits = 3))
end

Number of peaks: 2
Peak 1: CF: 4.559 PW: 0.526 BW: 3.463
Peak 2: CF: 22.102 PW: 0.307 BW: 3.764

Goodness of fit metrics:

println("R^2 of the model fit: ", round(f.r2[1, ep], digits = 3))
println("Error of the fit: ", round(f.err[1, ep], digits = 3))

R^2 of the model fit: 0.959
Error of the fit: 0.069

Visual presentation:

fig = Figure()
ax = Axis(fig[1, 1],
          xlabel = "Frequency",
          ylabel = "log10(Power)")
lines!(frq, pow, color = :black, label = "original spectrum")
lines!(frq, f.ap_fit[1, :, ep], color = :blue, linestyle = :dash, label = "aperiodic fit")
lines!(frq, f.model[1, :, ep], color = :red, label = "full model fit")
axislegend(ax)
fig