fooof-tools / fooof

Parameterizing neural power spectra into periodic & aperiodic components.
https://fooof-tools.github.io/
Apache License 2.0
359 stars 97 forks source link

difficulty understanding aperiodic fits #231

Closed shenweiss closed 2 years ago

shenweiss commented 2 years ago

Hi,

I am a new fooof user and I am trying to understand the meaning of the aperiodic offset and exponent. I have two sets of data and the average power spectrum and fit demonstrates a large difference in the offset. However, from the averages the lower offset visually has aperiodic parameters of [2.76 1.66] and the higher offset has aperiodic parameters of [2.77 1.38]. I know the exponent is lower in the average psd that appears to have a higher offset. My Python code is appended below. Please let me know if I did something wrong, or maybe have a fundamental misunderstanding of the methodology.

Thanks, Shennan

fooofdebug1 fooofdebug2

import numpy as np import matplotlib.pyplot as plt from matplotlib import cm, colors, colorbar

import mne from mne import io from mne.datasets import sample from mne.viz import plot_topomap from mne.time_frequency import psd_welch

from fooof import FOOOFGroup from fooof.bands import Bands from fooof.analysis import get_band_peak_fg

from fooof.plts.spectra import plot_spectra, plot_spectra_shading

Import simulation utilities for creating test data

from fooof.sim.gen import gen_power_spectrum, gen_group_power_spectra from fooof.sim.params import param_iter, Stepper import scipy.io from fooof.objs.utils import average_fg, combine_fooofs, compare_info

Non-propagating

aperiodic=[] peak=[] gaussian=[] rsquared=[]

file0='/data/downstate/wip/completed_analyses/part 2/4145fooofdemo.mat' matfile_vars = scipy.io.loadmat(file0) eegdata=matfile_vars['events_out_np'] n_channels = 1 sampling_freq = 2000 ch_names =['outnp'] ch_types='eeg' info = mne.create_info(ch_names=ch_names, ch_types=ch_types, sfreq=sampling_freq) eegdata = np.expand_dims(eegdata, axis=1) epochs=mne.EpochsArray(eegdata,info)

Calculate power spectra across the the continuous data

spectra, freqs = psd_welch(epochs, fmin=1, fmax=60, tmin=0, tmax=1000, n_overlap=150, n_fft=2000)

Initialize a FOOOFGroup object, with desired settings

fooofobjs=[]; for i in range(spectra.shape[0]): fg = FOOOFGroup(peak_width_limits=[1, 6], min_peak_height=0.15, peak_threshold=2., max_n_peaks=6, verbose=False)

Define the frequency range to fit

freq_range = [1, 60]
fg.fit(freqs, spectra[i], freq_range)
fooofobjs.append(fg)

fm=combine_fooofs(fooofobjs) bands=Bands({'delta': [2, 4]}) afm = average_fg(fm, bands, avg_method='median') afm.plot() np_aperiodic=afm.get_params('aperiodic_params') np_peak=afm.get_params('peak_params') np_gauss=afm.get_params('gaussian_params') np_rsqr=afm.get_params('r_squared')

out p

eegdata=matfile_vars['events_out_p'] n_channels = 1 sampling_freq = 2000 ch_names =['outp'] ch_types='eeg' info = mne.create_info(ch_names=ch_names, ch_types=ch_types, sfreq=sampling_freq) eegdata = np.expand_dims(eegdata, axis=1) epochs=mne.EpochsArray(eegdata,info)

Calculate power spectra across the the continuous data

spectra, freqs = psd_welch(epochs, fmin=1, fmax=60, tmin=0, tmax=1000, n_overlap=150, n_fft=2000)

Initialize a FOOOFGroup object, with desired settings

fooofobjs=[]; for i in range(spectra.shape[0]): fg = FOOOFGroup(peak_width_limits=[1, 6], min_peak_height=0.15, peak_threshold=2., max_n_peaks=6, verbose=False)

Define the frequency range to fit

freq_range = [1, 60]
fg.fit(freqs, spectra[i], freq_range)
fooofobjs.append(fg)

fm=combine_fooofs(fooofobjs) plt.figure() bands=Bands({'delta': [2, 4]}) afm = average_fg(fm, bands, avg_method='median') afm.plot() p_aperiodic=afm.get_params('aperiodic_params') p_peak=afm.get_params('peak_params') p_gauss=afm.get_params('gaussian_params') p_rsqr=afm.get_params('r_squared') print("done")

shenweiss commented 2 years ago

Sorry just realized it is a Lorentzian function so offset works in reverse, never mind. Cool!

TomDonoghue commented 2 years ago

Hey @shenweiss - I'm not sure I quite understand the original question here. Note that when you extract aperiodic parameters (in fixed model), you get [offset, exponent], so the parameters you printed imply that there is a difference in exponent, not offset, in the data you fit.

The fit function is a Lorentzian, but I'm not sure what you mean by the offset "working in reverse". Anyways, it seems you might have figured out what is going on here - but let me know if there are follow up questions.