diff --git a/doc/conf.py b/doc/conf.py index 0537f979..253382ed 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -100,6 +100,7 @@ ("Motivations", "auto_motivations/index"), ("Tutorials", "auto_tutorials/index"), ("Examples", "auto_examples/index"), + ("Visualizers", "visualizers"), ("Reference", "reference"), ("GitHub", "https://github.com/fooof-tools/fooof", True), ], @@ -125,8 +126,9 @@ sphinx_gallery_conf = { 'examples_dirs': ['../examples', '../tutorials', '../motivations'], 'gallery_dirs': ['auto_examples', 'auto_tutorials', 'auto_motivations'], - 'subsection_order' : ExplicitOrder(['../examples/manage', - '../examples/processing', + 'subsection_order' : ExplicitOrder(['../examples/processing', + '../examples/manage', + '../examples/models', '../examples/plots', '../examples/sims', '../examples/analyses', diff --git a/doc/contents.rst b/doc/contents.rst index aece81e7..343ee95b 100644 --- a/doc/contents.rst +++ b/doc/contents.rst @@ -9,6 +9,7 @@ Table of Contents glossary.rst reference.rst changelog.rst + visualizers.rst auto_tutorials/index.rst auto_examples/index.rst auto_motivations/index.rst diff --git a/doc/reference.rst b/doc/reference.rst index aae1fb0c..24725d7e 100644 --- a/doc/reference.rst +++ b/doc/reference.rst @@ -30,7 +30,7 @@ The code for the simulations and applications is indexed and available in the `Paper repository `_. You can also find the full list of articles that cite the `Parameterizing Neural Power Spectra` paper at this -`Google scholar link `_. +`Google scholar link `_. Methods Reporting ~~~~~~~~~~~~~~~~~ @@ -54,7 +54,7 @@ In addition, we recommend that reports should include information on: Reporting Template & Example ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -To assist in reporting on using FOOOF, we have created some templates for reporting on spectral parameterization methods. There are also some utilities included in the code to collected the required information. +To assist in reporting on using FOOOF, we have created some templates for reporting on spectral parameterization methods. There are also some utilities included in the code to collect the required information. The following box is an example of what a methods report might look like (where all of the *X*'s should be filled in with the relevant information). @@ -65,6 +65,18 @@ The following box is an example of what a methods report might look like (where peak threshold : *XX*; and aperiodic mode : *XX*. Power spectra were parameterized across the frequency range *XX* to *XX* Hz. +Checking module version +~~~~~~~~~~~~~~~~~~~~~~~ + +If you are not sure which version of the module you have installed, you can +check the `__version__` from Python, using the following code: + +.. code-block:: python + + # Check the version of the tool + from fooof import __version__ as fooof_version + print('Current fooof version:', fooof_version) + Generating Methods Reports ~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/doc/visualizers.rst b/doc/visualizers.rst new file mode 100644 index 00000000..d4d9f3cf --- /dev/null +++ b/doc/visualizers.rst @@ -0,0 +1,14 @@ +Visualizers +=========== + +This page includes animated visualizers to explore topics related to spectral parameterization. + +The source code to create these visualizations is in the +`visualizers repository `_. + +Spectral Rotation +----------------- + +Animated visualizer showing spectral rotation, emphasizing why offset and exponent are often correlated: + +.. image:: https://raw.githubusercontent.com/fooof-tools/Visualizers/main/gifs/specrot.gif diff --git a/examples/manage/plot_data_exporting.py b/examples/manage/plot_data_exporting.py new file mode 100644 index 00000000..9a3125a2 --- /dev/null +++ b/examples/manage/plot_data_exporting.py @@ -0,0 +1,174 @@ +""" +Exporting Model Results +======================= + +This example covers utilities for extracting and exporting model fit results. + +Note that the functionality to export to pandas DataFrames was added in version 1.1, +and requires the optional dependency `pandas` to be installed. +""" + +################################################################################################### + +# Import model objects, and Bands object to define bands of interest +from fooof import FOOOF, FOOOFGroup, Bands + +# Import simulation functions to create some example data +from fooof.sim import gen_power_spectrum, gen_group_power_spectra + +################################################################################################### +# Exporting Results +# ----------------- +# +# In this example we will explore available functionality for extracting model fit results, +# specifically the options available for exporting results to pandas objects. +# +# Note that the main use case of exporting models to pandas DataFrames is for +# analysis across models. If you are just trying to access the model fit results from +# a fit model, you may want the :meth:`~fooof.FOOOF.get_results` and/or +# :meth:`~fooof.FOOOF.get_params` methods. +# +# Defining Oscillation Bands +# ~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# A difficulty with exporting and collecting model results across model fits is that the +# models may be different sizes - most notably, they may contain different numbers of peaks. +# +# This means that we need to define some kind of strategy to organize the peak +# parameters across different models. Across these examples, this will include using the +# :class:`~fooof.Bands` object to define oscillations bands of interest. +# + +################################################################################################### + +# Initialize bands object, defining alpha band +bands1 = Bands({'alpha' : [7, 14]}) + +# Initialize model object +fm = FOOOF() + +################################################################################################### + +# Simulate example power spectrum +freqs, powers = gen_power_spectrum([1, 50], [0, 10, 1], [10, 0.25, 2], freq_res=0.25) + +# Fit model to power spectrum +fm.fit(freqs, powers) + +################################################################################################### +# +# The :meth:`~fooof.FOOOF.to_df` method supports exporting model fit results to pandas objects. +# + +################################################################################################### + +# Export results +fm.to_df(None) + +################################################################################################### +# +# In the above, we export the model fit results to a pandas Series. +# +# Note that we explicitly passed in `None` to the `peak_org` argument, meaning we did not define +# a strategy for managing peaks. Because of this, the export did not include any peak parameters. +# +# Next, let's can try exporting again, this time passing in an integer to define the number +# of peaks to extract. +# + +################################################################################################### + +# Export results - extracting 1 peak +fm.to_df(1) + +################################################################################################### +# Using Band Definitions +# ~~~~~~~~~~~~~~~~~~~~~~ +# +# In the above, we extract the results specifying to extract 1 peak. By default, this approach +# will extract the dominant (highest power) peak. +# +# Notably, specifying a set number of peaks to extract does allow for combining across peaks +# (in terms of enforcing the same model size), but may not be the ideal way to examine across +# peaks (since the dominant extract peak may vary across model fits). +# +# Therefore, we may often want to instead define a set of band definitions to organize peaks, +# as can be done by passing a `Bands` object in to the `to_df` method. +# + +################################################################################################### + +# Export results - extracting peaks based on bands object +fm.to_df(bands1) + +################################################################################################### +# +# Note that there are limitations to using the bands definitions - notably that doing so +# necessarily requires extracting at most 1 peak per band. In doing so, the extraction will +# select the dominant peak per band. However, this may not fully reflect the full model fit +# if there are, for example, multiple peaks fit within a band. +# + +################################################################################################### +# Example on Group Object +# ~~~~~~~~~~~~~~~~~~~~~~~ +# +# In the above, we used the model object to show the basic exporting functionalities. +# +# This functionality is more useful when considering multiple model fits together, such +# as can be done using the :meth:`~fooof.FOOOFGroup.to_df` method from the Group object, +# which allows for exporting DataFrames of organized model fit parameters across power spectra. +# +# As with the above, keep in mind that for some cases you may want the +# :meth:`~fooof.FOOOFGroup.get_results` and/or :meth:`~fooof.FOOOFGroup.get_params` methods +# instead of doing a DataFrame export. +# + +################################################################################################### + +# Simulate an example group of power spectra +freqs, powers = gen_group_power_spectra(5, [1, 50], [0, 1], [10, 0.25, 2]) + +# Initialize a group model object and fit power spectra +fg = FOOOFGroup(verbose=False) +fg.fit(freqs, powers) + +################################################################################################### + +# Export results to dataframe, with no peak definition +fg.to_df(None) + +################################################################################################### + +# Export results to dataframe, specifying to export a single peak +fg.to_df(1) + +################################################################################################### + +# Export results to dataframe, using bands definition with defines the alpha band +fg.to_df(bands1) + +################################################################################################### +# +# In the above examples, we showed the same exports on the Group object as we previously +# explored on the single spectrum in the model object. +# +# Note that we can also define new bands objects to change the peak output organization, +# as demonstrated in the following example. +# + +################################################################################################### + +# Define a new bands object, specifying both the alpha and beta band +bands2 = Bands({'alpha' : [7, 14], + 'beta' : [15, 30]}) + +################################################################################################### + +# Export results to dataframe, using bands object with alpha & beta +fg.to_df(bands2) + +################################################################################################### +# +# That covers the pandas export functionality available from the model objects. +# diff --git a/examples/models/README.txt b/examples/models/README.txt new file mode 100644 index 00000000..168deb76 --- /dev/null +++ b/examples/models/README.txt @@ -0,0 +1,4 @@ +Models & Parameters +------------------- + +Examples of evaluating and exploring model fits and parameters. \ No newline at end of file diff --git a/examples/models/plot_aperiodic_params.py b/examples/models/plot_aperiodic_params.py new file mode 100644 index 00000000..96f018e4 --- /dev/null +++ b/examples/models/plot_aperiodic_params.py @@ -0,0 +1,177 @@ +""" +Aperiodic Parameters +==================== + +Exploring properties and topics related to aperiodic parameters. +""" + +################################################################################################### + +from scipy.stats import spearmanr + +from fooof import FOOOF, FOOOFGroup +from fooof.plts.spectra import plot_spectra +from fooof.plts.annotate import plot_annotated_model +from fooof.plts.aperiodic import plot_aperiodic_params +from fooof.sim.params import Stepper, param_iter +from fooof.sim import gen_power_spectrum, gen_group_power_spectra +from fooof.utils.params import compute_time_constant, compute_knee_frequency + +################################################################################################### +# 'Fixed' Model +# ------------- +# +# First, we will explore the 'fixed' model, which fits an offset and exponent value +# to characterize the 1/f-like aperiodic component of the data. +# + +################################################################################################### + +# Simulate an example power spectrum +freqs, powers = gen_power_spectrum([1, 50], [0, 1], [10, 0.25, 2], freq_res=0.25) + +################################################################################################### + +# Initialize model object and fit power spectrum +fm = FOOOF(min_peak_height=0.1) +fm.fit(freqs, powers) + +################################################################################################### + +# Check the aperiodic parameters +fm.aperiodic_params_ + +################################################################################################### + +# Plot annotated model of aperiodic parameters +plot_annotated_model(fm, annotate_peaks=False, annotate_aperiodic=True, plt_log=True) + +################################################################################################### +# Comparing Offset & Exponent +# --------------------------- +# +# A common analysis of model fit parameters includes examining and comparing changes +# in the offset and/or exponent values of a set of models, which we will now explore. +# +# To do so, we will start by simulating a set of power spectra with different exponent values. +# + +################################################################################################### + +# Define simulation parameters, stepping across different exponent values +exp_steps = Stepper(0, 2, 0.25) +ap_params = param_iter([1, exp_steps]) + +################################################################################################### + +# Simulate a group of power spectra +freqs, powers = gen_group_power_spectra(\ + len(exp_steps), [3, 40], ap_params, [10, 0.25, 1], freq_res=0.25, f_rotation=10) + +################################################################################################### + +# Plot the set of example power spectra +plot_spectra(freqs, powers, log_powers=True) + +################################################################################################### + +# Initialize a group mode object and parameterize the power spectra +fg = FOOOFGroup() +fg.fit(freqs, powers) + +################################################################################################### + +# Extract the aperiodic values of the model fits +ap_values = fg.get_params('aperiodic') + +################################################################################################### + +# Plot the aperiodic parameters +plot_aperiodic_params(fg.get_params('aperiodic')) + +################################################################################################### + +# Compute the correlation between the offset and exponent +spearmanr(ap_values[0, :], ap_values[1, :]) + +################################################################################################### +# +# What we see above matches the common finding that that the offset and exponent are +# often highly correlated! This is because if you imagine a change in exponent as +# the spectrum 'rotating' around some frequency value, then (almost) any change in exponent +# has a corresponding change in offset value! If you note in the above, we actually specified +# a rotation point around which the exponent changes. +# +# This can also be seen in this animation showing this effect across different rotation points: +# +# .. image:: https://raw.githubusercontent.com/fooof-tools/Visualizers/main/gifs/specrot.gif +# +# Notably this means that while the offset and exponent can change independently (there can be +# offset changes over and above exponent changes), the baseline expectation is that these +# two parameters are highly correlated and likely reflect the same change in the data! +# + +################################################################################################### +# Knee Model +# ---------- +# +# Next, let's explore the knee model, which parameterizes the aperiodic component with +# an offset, knee, and exponent value. +# + +################################################################################################### + +# Generate a power spectrum with a knee +freqs2, powers2 = gen_power_spectrum([1, 50], [0, 15, 1], [8, 0.125, 0.75], freq_res=0.25) + +################################################################################################### + +# Initialize model object and fit power spectrum +fm = FOOOF(min_peak_height=0.05, aperiodic_mode='knee') +fm.fit(freqs2, powers2) + +################################################################################################### + +# Plot annotated knee model +plot_annotated_model(fm, annotate_peaks=False, annotate_aperiodic=True, plt_log=True) + +################################################################################################### + +# Check the measured aperiodic parameters +fm.aperiodic_params_ + +################################################################################################### +# Knee Frequency +# ~~~~~~~~~~~~~~ +# +# You might notice that the knee *parameter* is not an obvious value. Notably, this parameter +# value as extracted from the model is something of an abstract quantify based on the +# formalization of the underlying fit function. A more intuitive measure that we may +# be interested in is the 'knee frequency', which is an estimate of the frequency value +# at which the knee occurs. +# +# The :func:`~.compute_knee_frequency` function can be used to compute the knee frequency. +# + +################################################################################################### + +# Compute the knee frequency from aperiodic parameters +knee_frequency = compute_knee_frequency(*fm.aperiodic_params_[1:]) +print('Knee frequency: ', knee_frequency) + +################################################################################################### +# Time Constant +# ~~~~~~~~~~~~~ +# +# Another interesting property of the knee parameter is that it has a direct relationship +# to the auto-correlation function, and from there to the empirical time constant of the data. +# +# The :func:`~.compute_time_constant` function can be used to compute the knee-derived +# time constant. +# + +################################################################################################### + +# Compute the characteristic time constant of a knee value +time_constant = compute_time_constant(fm.get_params('aperiodic', 'knee')) +print('Knee derived time constant: ', time_constant) diff --git a/examples/models/plot_data_components.py b/examples/models/plot_data_components.py new file mode 100644 index 00000000..768434e4 --- /dev/null +++ b/examples/models/plot_data_components.py @@ -0,0 +1,123 @@ +""" +Exploring Data Components +========================= + +This example explores the different data components, exploring the isolated aperiodic +and periodic components as they are extracted from the data. +""" + +################################################################################################### + +# sphinx_gallery_thumbnail_number = 3 + +# Import FOOOF model objects +from fooof import FOOOF, FOOOFGroup + +# Import function to plot power spectra +from fooof.plts.spectra import plot_spectra + +# Import simulation functions to create some example data +from fooof.sim import gen_power_spectrum, gen_group_power_spectra + +################################################################################################### + +# Simulate example power spectrum +freqs, powers = gen_power_spectrum([1, 50], [0, 10, 1], [10, 0.25, 2], freq_res=0.25) + +# Initialize model object and fit power spectrum +fm = FOOOF() +fm.fit(freqs, powers) + +################################################################################################### +# Data Components +# ~~~~~~~~~~~~~~~ +# +# The model fit process includes procedures for isolating aperiodic and periodic components in +# the data, fitting each of these components separately, and then combining the model components +# as the final fit (see the Tutorials for further details on these procedures). +# +# In doing this process, the model fit procedure computes and stores isolated data components, +# which are available in the model. +# +# Before diving into the isolated data components, let's check the data (`power_spectrum`) +# and full model fit of a model object (`fooofed_spectrum`). +# + +################################################################################################### + +# Plot the original power spectrum data from the object +plot_spectra(fm.freqs, fm.power_spectrum, color='black') + +################################################################################################### + +# Plot the power spectrum model from the object +plot_spectra(fm.freqs, fm.fooofed_spectrum_, color='red') + +################################################################################################### +# Aperiodic Component +# ~~~~~~~~~~~~~~~~~~~ +# +# To fit the aperiodic component, the model fit procedure includes a peak removal process. +# +# The resulting 'peak-removed' data component is stored in the model object, in the +# `_spectrum_peak_rm` attribute. +# + +################################################################################################### + +# Plot the peak removed spectrum data component +plot_spectra(fm.freqs, fm._spectrum_peak_rm, color='black') + +################################################################################################### + +# Plot the peak removed spectrum, with the model aperiodic fit +plot_spectra(fm.freqs, [fm._spectrum_peak_rm, fm._ap_fit], + colors=['black', 'blue'], linestyle=['-', '--']) + +################################################################################################### +# Periodic Component +# ~~~~~~~~~~~~~~~~~~ +# +# To fit the periodic component, the model fit procedure removes the fit peaks from the power +# spectrum. +# +# The resulting 'flattened' data component is stored in the model object, in the +# `_spectrum_flat` attribute. +# + +################################################################################################### + +# Plot the flattened spectrum data component +plot_spectra(fm.freqs, fm._spectrum_flat, color='black') + +################################################################################################### + +# Plot the flattened spectrum data with the model peak fit +plot_spectra(fm.freqs, [fm._spectrum_flat, fm._peak_fit], colors=['black', 'green']) + +################################################################################################### +# Full Model Fit +# ~~~~~~~~~~~~~~ +# +# The full model fit, which we explored earlier, is calculated as the combination of the +# aperiodic and peak fit, which we can check by plotting these combined components. +# + +################################################################################################### + +# Plot the full model fit, as the combination of the aperiodic and peak model components +plot_spectra(fm.freqs, [fm._ap_fit + fm._peak_fit], color='red') + +################################################################################################### +# Notes on Analyzing Data Components +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The above shows data components as they are available on the model object, and used in +# the fitting process. Some analyses may aim to use these isolated components to compute +# certain measures of interest on the data. Note that these data components are stored in +# 'private' attributes (indicated by a leading underscore), meaning in normal function they +# are not expected to be accessed by the user, but as we've seen above they can still be accessed. +# However, analyses derived from these isolated data components is not currently officially +# supported by the module, and so users who wish to do so should consider the benefits and +# limitations of any such analyses. +# diff --git a/examples/manage/plot_freq_by_freq_error.py b/examples/models/plot_freq_by_freq_error.py similarity index 100% rename from examples/manage/plot_freq_by_freq_error.py rename to examples/models/plot_freq_by_freq_error.py diff --git a/examples/models/plot_peak_params.py b/examples/models/plot_peak_params.py new file mode 100644 index 00000000..5498e0bb --- /dev/null +++ b/examples/models/plot_peak_params.py @@ -0,0 +1,102 @@ +""" +Periodic Parameters +=================== + +Exploring properties and topics related to peak parameters. +""" + +################################################################################################### + +from fooof import FOOOF, FOOOFGroup +from fooof.plts.spectra import plot_spectra +from fooof.plts.periodic import plot_peak_params +from fooof.sim.utils import set_random_seed +from fooof.sim.params import Stepper, param_iter +from fooof.sim import gen_power_spectrum, gen_group_power_spectra +from fooof.plts.annotate import plot_annotated_model +from fooof.utils.params import compute_time_constant, compute_knee_frequency + +################################################################################################### +# Gaussian Peak Model +# ------------------- +# +# By default, the spectral parameterization model fits Gaussians to any detected peaks. +# +# These Gaussians are defined by a mean, height, and standard deviation, which we in turn +# interpret as the center frequency (CF), power (PW), and bandwidth (BW) of the putative +# oscillation. +# +# In this example, we will further explore these peak parameters and some topics and +# limitations related to their use and interpretations. +# + +################################################################################################### + +# Simulate an example power spectrum +freqs, powers = gen_power_spectrum([3, 40], [0, 1], [10, 0.3, 1.], freq_res=0.25) + +################################################################################################### + +# Initialize model object and fit power spectrum +fm = FOOOF(min_peak_height=0.1) +fm.fit(freqs, powers) + +################################################################################################### + +# Plot annotated model labelling the peaks +plot_annotated_model(fm, annotate_peaks=True, annotate_aperiodic=False, plt_log=True) + +################################################################################################### +# +# In the above we can see an example of the fit peak parameters +# + +################################################################################################### +# Overlapping Peaks +# ----------------- +# +# Let's now consider some features of fitting Gaussian peaks, and how this relates to +# the neural data under study. In particular, Gaussian's are symmetric and while they do +# seem to approximate the peaks we observe in emprical data quite well, not all peaks +# in empirical power spectra are quite symmetric. +# +# To deal with this, the model sometimes fits overlapping peaks, whereby two or more peaks +# are used by the model to capture the shape of what otherwise looks like a single peak. +# +# We can explore this in a simplified simulation. +# + +################################################################################################### + +# Set the random seed +set_random_seed(10) + +################################################################################################### + +# Simulate an example power spectrum created with an asymmetric peak +freqs, powers = gen_power_spectrum([3, 40], [0, 1], [[10, 0.3, 1.], [11.25, 0.175, 0.5]], freq_res=0.25) + +################################################################################################### + +# Initialize model object and fit power spectrum +fm = FOOOF(min_peak_height=0.1) +fm.report(freqs, powers) + +################################################################################################### +# +# As we can see in the above model solution, in the data, it looks like there is a single +# oscillatory peak, and yet the model has captured this peak with two Gaussians. +# +# This example serves to demonstrate two key points. First, not all Gaussians fit in the model +# necessary reflect separate peaks, as some may overlap. Second, when peaks overlap, +# the parameters of each individually may accurately capture a peak in the data, as the +# overall shape of the peak may be captured as the interaction across multiple Gaussians +# (this is most common / notable for the bandwidth measure, whereby the width of the peak is +# best described as the combined width of the two adjacent peaks). +# +# Note that, by construction, this simulated example was created by simulating two overlapping +# peaks, and so in that sense the model is actually correct in it's solution. In empirical +# data, we do not know if a power spectrum that looks like this does reflect two underlying +# oscillatory processes, or perhaps a single oscillatory process that happens to be asymmetric +# in the frequency domain. +# diff --git a/fooof/core/strings.py b/fooof/core/strings.py index 15ef6137..a118cca8 100644 --- a/fooof/core/strings.py +++ b/fooof/core/strings.py @@ -209,7 +209,7 @@ def gen_methods_report_str(concise=False): # Methods report information 'To report on using FOOOF, you should report (at minimum):', '', - '- the code version that was used used', + '- the code version that was used', '- the algorithm settings that were used', '- the frequency range that was fit', diff --git a/motivations/concepts/plot_WigglyPeaks.py b/motivations/concepts/plot_WigglyPeaks.py new file mode 100644 index 00000000..2a605c33 --- /dev/null +++ b/motivations/concepts/plot_WigglyPeaks.py @@ -0,0 +1,187 @@ +""" +'Oscillations' as Peaks +======================= + +Exploring the idea of oscillations as peaks of power. +""" + +################################################################################################### + +# Imports from NeuroDSP to simulate & plot time series +from neurodsp.sim import sim_powerlaw, sim_oscillation, sim_combined, set_random_seed +from neurodsp.spectral import compute_spectrum +from neurodsp.plts import plot_time_series, plot_power_spectra +from neurodsp.utils import create_times + +################################################################################################### + +# Define simulation settings +n_seconds = 30 +fs = 1000 + +# Create a times vector +times = create_times(n_seconds, fs) + +################################################################################################### +# Frequency Specific Power +# ------------------------ +# +# Part of the motivation behind spectral parameterization is dissociating aperiodic +# activity, with no characteristic frequency, to periodic power, which is defined as +# having frequency specific power. This leads to the idea of oscillations as 'peaks' +# of power in the power spectrum, which can be detected and measured. +# +# In this exploration, we will use simulated time series to examine how rhythmic signals +# do display as 'peaks' of power in frequency domain representations. We will also +# explore some limitations of this idea. +# + +################################################################################################### + +# Simulate an ongoing sine wave +sine_wave = sim_oscillation(n_seconds, fs, 10) + +# Visualize the sine wave +plot_time_series(times, sine_wave, xlim=[0, 5]) + +################################################################################################### +# +# In the above, we simulated a pure sinusoid, at 10 Hz. +# + +################################################################################################### + +# Compute the power spectrum of the sine wave +freqs, powers = compute_spectrum(sine_wave, fs) + +# Visualize the power spectrum +plot_power_spectra(freqs, powers) + +################################################################################################### +# +# The power spectrum of the sine wave shows a clear peak of power at 10 Hz, reflecting +# the simulated rhythm, with close to zero power at all other frequencies. +# +# This is characteristic of a sine wave, and demonstrates the basic idea of considering +# oscillations as peaks of power in the power spectrum. +# +# Next lets examine a more complicated signal, with multiple components. +# + +################################################################################################### + +# Define components for a combined signal +components = { + 'sim_oscillation' : {'freq' : 10}, + 'sim_powerlaw' : {'exponent' : -1}, +} + +# Simulate a combined signal +sig = sim_combined(n_seconds, fs, components) + +# Visualize the time series +plot_time_series(times, sig, xlim=[0, 5]) + +################################################################################################### +# +# Here we see a simply simulation meant to be closer to neural data, reflecting an oscillatory +# component, as well as an aperiodic component, which contributes power across all frequencies. +# + +################################################################################################### + +# Compute the power spectrum of the sine wave +freqs, powers = compute_spectrum(sig, fs) + +# Visualize the power spectrum +plot_power_spectra(freqs, powers) + +################################################################################################### +# Interim Conclusion +# ~~~~~~~~~~~~~~~~~~ +# +# In the power spectrum of the combined signal, we can still the peak of power at 10 Hz, as +# well as the pattern of power across all frequencies contributed by the aperiodic component. +# +# This basic example serves as the basic motivation for spectral parameterization. In this +# simulated example, we know there are two components, and so a procedure for detecting the +# peaks and measuring the pattern of aperiodic power (as is done in spectral parameterization) +# provides a method to measuring these components in the data. +# + +################################################################################################### +# Harmonic Peaks +# -------------- +# +# The above seeks to demonstrate the basic idea whereby a peak of power _may_ reflect +# an oscillation at that frequency, where as patterns of power across all frequencies +# likely reflect aperiodic activity. +# +# In the this section, we will further explore peaks of power in the frequency domain, +# showing that not every peak necessarily reflect an independent oscillation. +# +# To do so, we will start by simulating a non-sinusoidal oscillation. +# + +################################################################################################### + +# Simulate a sawtooth wave +sawtooth = sim_oscillation(n_seconds, fs, 10, 'sawtooth', width=0.5) + +# Visualize the sine wave +plot_time_series(times, sawtooth, xlim=[0, 5]) + +################################################################################################### +# +# In the above, we can see that there is again an oscillation, although it is not sinusoidal. +# + +################################################################################################### + +# Compute the power spectrum of the sine wave +freqs, powers = compute_spectrum(sawtooth, fs) + +# Visualize the power spectrum +plot_power_spectra(freqs, powers) + +################################################################################################### +# +# Note the 10 Hz peak, as well as the additional peaks in the frequency domain. +# +# Before further discussing this, let's create an example with an aperiodic component. +# + +################################################################################################### + +# Define components for a combined signal +components = { + 'sim_oscillation' : {'freq' : 10, 'cycle' : 'sawtooth', 'width' : 0.25}, + 'sim_powerlaw' : {'exponent' : -1.}, +} + +# Simulate a combined signal +sig = sim_combined(n_seconds, fs, components) + +# Visualize the time series +plot_time_series(times, sig, xlim=[0, 5]) + +################################################################################################### + +# Compute the power spectrum of the sine wave +freqs, powers = compute_spectrum(sig, fs) + +# Visualize the power spectrum +plot_power_spectra(freqs, powers) + +################################################################################################### +# +# In the power spectrum above, we see that there is a peak of power at the fundamental frequency +# of the rhythm (10 Hz), but there are also numerous other peaks. These additional peaks are +# 'harmonics', and that are components of the frequency domain representation that reflect +# the non-sinusoidality of the time domain signal. +# +# This serves as the basic motivation for the claim that although a peak _may_ reflect an +# independent oscillation, this need not be the case, as a given peak could simply be the harmonic +# of an asymmetric oscillation at a different frequency. For this reason, the number of peaks in +# a model can not be interpreted as the number of oscillations in a signal. +# diff --git a/tutorials/plot_09-Reporting.py b/tutorials/plot_09-Reporting.py new file mode 100644 index 00000000..8d530da9 --- /dev/null +++ b/tutorials/plot_09-Reporting.py @@ -0,0 +1,170 @@ +""" +09: Reporting & Referencing +=========================== + +This section covers how to access reporting info and reference use of the module. + +This page is a hands-on example of the reporting and referencing information on the +`Reference page `_. +""" + +################################################################################################### + +# Import FOOOF model objects +from fooof import FOOOF, FOOOFGroup + +# Import simulation functions to create some example data +from fooof.sim import gen_power_spectrum, gen_group_power_spectra + +# Import utilities to print out information for reporting +from fooof.utils.reports import methods_report_info, methods_report_text + +# sphinx_gallery_start_ignore +# Note: this code gets hidden, but serves to create the text plot for the icon +from fooof.core.strings import gen_methods_report_str +from fooof.core.reports import REPORT_FONT +import matplotlib.pyplot as plt +text = gen_methods_report_str(concise=True) +text = text[0:142] + '\n' + text[142:] +_, ax = plt.subplots(figsize=(8, 3)) +ax.text(0.5, 0.5, text, REPORT_FONT, ha='center', va='center') +ax.set_frame_on(False) +_ = ax.set(xticks=[], yticks=[]) +# sphinx_gallery_end_ignore + +################################################################################################### +# Checking Module Version +# ----------------------- +# +# It's useful and important to keep track of which version of the module you are using, +# and to report this information when referencing use of the tool. +# +# There are several ways to check the version of the the module that you are using, +# including checking what is installed in the Python environment you are using. +# +# From within Python, you can also check the version of the module by checking +# `__version__`, as shown below: +# + +################################################################################################### + +# Check the version of the module +from fooof import __version__ as fooof_version +print('Current fooof version:', fooof_version) + +################################################################################################### +# Getting Model Reporting Information +# ----------------------------------- +# +# To assist with reporting information, the module has the following utilities: +# +# - the :func:`~.methods_report_info`, which prints out methods reporting information +# - the :func:`~.methods_report_text`, which prints out a template of methods reporting text +# +# Both of the these functions take as input a model object, and use the object to +# collect and return information related to the model fits. +# +# Note that not all information may be returned by the model - these methods only have access +# to the current object. This means it is also important that if you use these functions, +# you use them with objects that match the settings and data used in the analysis to be reported. +# +# In the following, we will explore an example of using these functions for an example model fit. +# + +################################################################################################### + +# Initialize model object +fooof_obj = FOOOF() + +################################################################################################### + +# Print out all the methods information for reporting +methods_report_info(fooof_obj) + +################################################################################################### +# +# You might notice that the above function prints out several different sections, +# some of which might look familiar. +# +# The settings information, for example, is the same as printed using the +# # - :meth:`~fooof.FOOOF.print_settings` method. +# +# Next, let's check out the text version of the methods report. +# + +################################################################################################### + +# Generate methods text, with methods information inserted +methods_report_text(fooof_obj) + +################################################################################################### +# Additional Examples +# ~~~~~~~~~~~~~~~~~~~ +# +# In the above examples, not all the information was printed, as not all information was +# available in the example object (for example, it had no data). +# +# In the next example, let's see how the outputs look for an example object with model fit results. +# + +################################################################################################### + +# Simulate an example power spectrum +freqs, powers = gen_power_spectrum([1, 50], [0, 10, 1], [10, 0.25, 2], freq_res=0.25) + +# Initialize model object +fm = FOOOF(min_peak_height=0.1, peak_width_limits=[1, 6], aperiodic_mode='knee') +fm.fit(freqs, powers) + +################################################################################################### + +# Generate methods text, with methods information inserted +methods_report_info(fm) + +################################################################################################### + +# Generate methods text, with methods information inserted +methods_report_text(fm) + +################################################################################################### +# +# The report text is meant as a template / example that could be added to a methods section. +# +# Note that there may be missing information that needs to be filled in (indicated by 'XX's), +# and you can and should consider this a template to be edited as needed. +# + +################################################################################################### +# Other Model Objects +# ~~~~~~~~~~~~~~~~~~~ +# +# Note that the reporting functions work with any model object. +# +# For example, next we will use them on a :class:`~fooof.FOOOFGroup` object. +# + +################################################################################################### + +# Simulate an example group of power spectra +freqs, powers = gen_group_power_spectra(10, [1, 75], [0, 1], [10, 0.25, 2]) + +################################################################################################### + +# Initialize and fit group model object +fg = FOOOFGroup(max_n_peaks=4, peak_threshold=1.75) +fg.fit(freqs, powers) + +################################################################################################### + +# Generate the methods report information +methods_report_info(fg) + +################################################################################################### + +# Generate the methods report text +methods_report_text(fg) + +################################################################################################### +# +# That concludes the example of using the helper utilities for generating methods reports! +#