Hello everyone,
I am working on estimating the 1/f (aperiodic) component of my power spectrum for source space power mapping with Dynamic Imaging of Coherent Sources (DICS).
DICS implementation as in:
MNE DICS Tutorial
FLUX MNE DICS TUTORIAL
I am facing the issue that the fooof algorithm does not currently support complex values.
# Check if power values are complex
if np.iscomplexobj(power_spectrum):
raise DataError("Input power spectra are complex values. "
"FOOOF does not currently support complex inputs.")
Has anyone tried this before?
When and how would you implement the fooofing in the pipeline?
My first intuition would be to estimate the aperiodic and periodic components before creating and applying the DICS filters. However, this would require the complex power spectrum to be converted to real values, which complicates/hampers the later DICS power mapping…
Alternatively, one could estimate the power per frequency via DICS and do the fooofing at the very end. This could be very time- and computation-consuming… Any suggestions?
If you need more information about the code or some sample data, just let me know 
Thanks for your support!
Cheers,
Katharina
#########################
- MNE version: e.g. 1.3.0
- fooof version: 1.0.0
- operating system: Ubuntu 20.04
#########################
Example Code:
import numpy as np
from fooof import FOOOF
from joblib import Parallel, delayed
from tqdm import tqdm
# Fooofing Algorithm
def fooofing_power_single_vertex(vertex, freqs, power):
fm = FOOOF(peak_width_limits=[2, 6], min_peak_height=0.1, aperiodic_mode='fixed', verbose=True)
fm.fit(freqs, power[vertex, :])
r2, error = fm.r_squared_, fm.error_
offset, slope = fm.aperiodic_params_
L = 10 ** (offset - slope * np.log10(freqs))
fooofed_power = power[vertex, :] - L
fooof_meta = np.array([offset, slope, r2, error])
return fooofed_power, fooof_meta
def fooofing_power(power, freqs, n_jobs):
results = Parallel(n_jobs=n_jobs)(
delayed(fooofing_power_single_vertex)(vertex, freqs, power) for vertex in tqdm(range(power.shape[0])))
fooofed_power = np.zeros(power.shape)
fooof_meta = np.zeros((power.shape[0], 4))
for vertex, (fp, fm) in enumerate(results):
fooofed_power[vertex, :] = fp
fooof_meta[vertex, :] = fm
return fooofed_power, fooof_meta
#####################################################
from mne.time_frequency import csd_morlet
from mne.beamformer import make_dics, apply_dics_csd
cond = 'HW' # High working memory load
csd = csd_morlet(epochs[cond], np.arange(fmin, fmax, 1),
tmin=epochs[cond].tmin,
tmax=epochs[cond].tmax, n_cycles=3)
csd_fooofed = csd.copy()
# Traceback due to complex data
csd_fooof_cond._data, fooof_meta = fooofing_power(csd._data, csd.frequencies, n_jobs)
# Converting the csd to a real power spectrum (np.abs()) leads unsurprisingly to a RuntimeWarning during make_dics():
# data covariance does not appear to be positive semidefinite, results will likely be incorrect