Skip to content

Commit fd78ad7

Browse files
committed
add draft of aperiodic params example
1 parent ee0ea5a commit fd78ad7

File tree

1 file changed

+177
-0
lines changed

1 file changed

+177
-0
lines changed
Lines changed: 177 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,177 @@
1+
"""
2+
Aperiodic Parameters
3+
====================
4+
5+
Exploring properties and topics related to aperiodic parameters.
6+
"""
7+
8+
###################################################################################################
9+
10+
from scipy.stats import spearmanr
11+
12+
from fooof import FOOOF, FOOOFGroup
13+
from fooof.plts.spectra import plot_spectra
14+
from fooof.plts.annotate import plot_annotated_model
15+
from fooof.plts.aperiodic import plot_aperiodic_params
16+
from fooof.sim.params import Stepper, param_iter
17+
from fooof.sim import gen_power_spectrum, gen_group_power_spectra
18+
from fooof.utils.params import compute_time_constant, compute_knee_frequency
19+
20+
###################################################################################################
21+
# 'Fixed' Model
22+
# -------------
23+
#
24+
# First, we will explore the 'fixed' model, which fits an offset and exponent value
25+
# to characterize the 1/f-like aperiodic component of the data.
26+
#
27+
28+
###################################################################################################
29+
30+
# Simulate an example power spectrum
31+
freqs, powers = gen_power_spectrum([1, 50], [0, 1], [10, 0.25, 2], freq_res=0.25)
32+
33+
###################################################################################################
34+
35+
# Initialize model object and fit power spectrum
36+
fm = FOOOF(min_peak_height=0.1)
37+
fm.fit(freqs, powers)
38+
39+
###################################################################################################
40+
41+
# Check the aperiodic parameters
42+
fm.aperiodic_params_
43+
44+
###################################################################################################
45+
46+
# Plot annotated model of aperiodic parameters
47+
plot_annotated_model(fm, annotate_peaks=False, annotate_aperiodic=True, plt_log=True)
48+
49+
###################################################################################################
50+
# Comparing Offset & Exponent
51+
# ---------------------------
52+
#
53+
# A common analysis of model fit parameters includes examining and comparing changes
54+
# in the offset and/or exponent values of a set of models, which we will now explore.
55+
#
56+
# To do so, we will start by simulating a set of power spectra with different exponent values.
57+
#
58+
59+
###################################################################################################
60+
61+
# Define simulation parameters, stepping across different exponent values
62+
exp_steps = Stepper(0, 2, 0.25)
63+
ap_params = param_iter([1, exp_steps])
64+
65+
###################################################################################################
66+
67+
# Simulate a group of power spectra
68+
freqs, powers = gen_group_power_spectra(\
69+
len(exp_steps), [3, 40], ap_params, [10, 0.25, 1], freq_res=0.25, f_rotation=10)
70+
71+
###################################################################################################
72+
73+
# Plot the set of example power spectra
74+
plot_spectra(freqs, powers, log_powers=True)
75+
76+
###################################################################################################
77+
78+
# Initialize a group mode object and parameterize the power spectra
79+
fg = FOOOFGroup()
80+
fg.fit(freqs, powers)
81+
82+
###################################################################################################
83+
84+
# Extract the aperiodic values of the model fits
85+
ap_values = fg.get_params('aperiodic')
86+
87+
###################################################################################################
88+
89+
# Plot the aperiodic parameters
90+
plot_aperiodic_params(fg.get_params('aperiodic'))
91+
92+
###################################################################################################
93+
94+
# Compute the correlation between the offset and exponent
95+
spearmanr(ap_values[0, :], ap_values[1, :])
96+
97+
###################################################################################################
98+
#
99+
# What we see above matches the common finding that that the offset and exponent are
100+
# often highly correlated! This is because if you imagine a change in exponent as
101+
# the spectrum 'rotating' around some frequency value, then (almost) any change in exponent
102+
# has a corresponding change in offset value! If you note in the above, we actually specified
103+
# a rotation point around which the exponent changes.
104+
#
105+
# This can also be seen in this animation showing this effect across different rotation points:
106+
#
107+
# ![spectralrotation](https://raw.githubusercontent.com/fooof-tools/Visualizers/main/gifs/specrot.gif "specrot")
108+
#
109+
# Notably this means that while the offset and exponent can change independently (there can be
110+
# offset changes over and above exponent changes), the baseline expectation is that these
111+
# two parameters are highly correlated and likely reflect the same change in the data!
112+
#
113+
114+
###################################################################################################
115+
# Knee Model
116+
# ----------
117+
#
118+
# Next, let's explore the knee model, which parameterizes the aperiodic component with
119+
# an offset, knee, and exponent value.
120+
#
121+
122+
###################################################################################################
123+
124+
# Generate a power spectrum with a knee
125+
freqs2, powers2 = gen_power_spectrum([1, 50], [0, 15, 1], [8, 0.125, 0.75], freq_res=0.25)
126+
127+
###################################################################################################
128+
129+
# Initialize model object and fit power spectrum
130+
fm = FOOOF(min_peak_height=0.05, aperiodic_mode='knee')
131+
fm.fit(freqs2, powers2)
132+
133+
###################################################################################################
134+
135+
# Plot annotated knee model
136+
plot_annotated_model(fm, annotate_peaks=False, annotate_aperiodic=True, plt_log=True)
137+
138+
###################################################################################################
139+
140+
# Check the measured aperiodic parameters
141+
fm.aperiodic_params_
142+
143+
###################################################################################################
144+
# Knee Frequency
145+
# ~~~~~~~~~~~~~~
146+
#
147+
# You might notice that the knee _parameter_ is not an obvious value. Notably, this parameter
148+
# value as extracted from the model is something of an abstract quantify based on the
149+
# formalization of the underlying fit function. A more intuitive measure that we may
150+
# be interested in is the 'knee frequency', which is an estimate of the frequency value
151+
# at which the knee occurs.
152+
#
153+
# The `:func:`~.compute_knee_frequency` function can be used to compute the knee frequency.
154+
#
155+
156+
###################################################################################################
157+
158+
# Compute the knee frequency from aperiodic parameters
159+
knee_frequency = compute_knee_frequency(*fm.aperiodic_params_[1:])
160+
print('Knee frequency: ', knee_frequency)
161+
162+
###################################################################################################
163+
# Time Constant
164+
# ~~~~~~~~~~~~~
165+
#
166+
# Another interesting property of the knee parameter is that it has a direct relationship
167+
# to the auto-correlation function, and from there to the empirical time constant of the data.
168+
#
169+
# The `:func:`~.compute_time_constant` function can be used to compute the knee-derived
170+
# time constant.
171+
#
172+
173+
###################################################################################################
174+
175+
# Compute the characteristic time constant of a knee value
176+
time_constant = compute_time_constant(fm.get_params('aperiodic', 'knee'))
177+
print('Knee derived time constant: ', time_constant)

0 commit comments

Comments
 (0)