diff --git a/doc/conf.py b/doc/conf.py index 173c86876..0537f9793 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -126,6 +126,7 @@ 'examples_dirs': ['../examples', '../tutorials', '../motivations'], 'gallery_dirs': ['auto_examples', 'auto_tutorials', 'auto_motivations'], 'subsection_order' : ExplicitOrder(['../examples/manage', + '../examples/processing', '../examples/plots', '../examples/sims', '../examples/analyses', diff --git a/examples/processing/README.txt b/examples/processing/README.txt new file mode 100644 index 000000000..43ac437e1 --- /dev/null +++ b/examples/processing/README.txt @@ -0,0 +1,4 @@ +Processing +---------- + +Examples on how to process data related to spectral parameterization. \ No newline at end of file diff --git a/examples/processing/plot_line_noise.py b/examples/processing/plot_line_noise.py new file mode 100644 index 000000000..b1300f8d0 --- /dev/null +++ b/examples/processing/plot_line_noise.py @@ -0,0 +1,219 @@ +""" +Dealing with Line Noise +======================= + +This example covers strategies for dealing with line noise. +""" + +################################################################################################### + +# sphinx_gallery_thumbnail_number = 2 + +# Import the spectral parameterization object and utilities +from fooof import FOOOF +from fooof.plts import plot_spectrum, plot_spectra +from fooof.utils import trim_spectrum, interpolate_spectrum + +# Import simulation functions to create some example data +from fooof.sim.gen import gen_power_spectrum + +# Import NeuroDSP functions for simulating & processing time series +from neurodsp.sim import sim_combined +from neurodsp.filt import filter_signal +from neurodsp.spectral import compute_spectrum + +################################################################################################### +# Line Noise Peaks +# ---------------- +# +# Neural recordings typically have power line artifacts, at either 50 or 60 Hz, depending on +# where the data were collected, which can impact spectral parameterization. +# +# In this example, we explore some options for dealing with line noise artifacts. +# +# Interpolating Line Noise Peaks +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# One approach is to interpolate away line noise peaks, in the frequency domain. This +# approach simply gets rid of the peaks, interpolating the data to maintain the 1/f +# character of the data, allowing for subsequent fitting. +# +# The :func:`~fooof.utils.interpolate_spectrum` function allows for doing simple +# interpolation. Given a narrow frequency region, this function interpolates the spectrum, +# such that the 'peak' of the line noise is removed. +# + +################################################################################################### + +# Generate an example power spectrum, with line noise +freqs1, powers1 = gen_power_spectrum([3, 75], [1, 1], + [[10, 0.75, 2], [60, 1, 0.5]]) + +# Visualize the generated power spectrum +plot_spectrum(freqs1, powers1, log_powers=True) + +################################################################################################### +# +# In the plot above, we have an example spectrum, with some power line noise. +# +# To prepare this data for fitting, we can interpolate away the line noise region. +# + +################################################################################################### + +# Interpolate away the line noise region +interp_range = [58, 62] +freqs_int1, powers_int1 = interpolate_spectrum(freqs1, powers1, interp_range) + +################################################################################################### + +# Plot the spectra for the power spectra before and after interpolation +plot_spectra(freqs1, [powers1, powers_int1], log_powers=True, + labels=['Original Spectrum', 'Interpolated Spectrum']) + +################################################################################################### +# +# As we can see in the above, the interpolation removed the peak from the data. +# +# We can now go ahead and parameterize the spectrum. +# + +################################################################################################### + +# Initialize a power spectrum model +fm1 = FOOOF(verbose=False) +fm1.report(freqs_int1, powers_int1) + +################################################################################################### +# Multiple Interpolation Regions +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Line noise artifacts often also display harmonics, such that when analyzing broader +# frequency ranges, there may be multiple peaks that need to be interpolated. +# +# This can be done by passing in multiple interpolation regions to +# :func:`~fooof.utils.interpolate_spectrum`, which we will do in the next example. +# + +################################################################################################### + +# Generate an example power spectrum, with line noise & harmonics +freqs2, powers2 = gen_power_spectrum([1, 150], [1, 500, 1.5], + [[10, 0.5, 2], [60, 0.75, 0.5], [120, 0.5, 0.5]]) + +# Interpolate away the line noise region & harmonics +interp_ranges = [[58, 62], [118, 122]] +freqs_int2, powers_int2 = interpolate_spectrum(freqs2, powers2, interp_ranges) + +################################################################################################### + +# Plot the power spectrum before and after interpolation +plot_spectra(freqs2, [powers2, powers_int2], log_powers=True, + labels=['Original Spectrum', 'Interpolated Spectrum']) + +################################################################################################### + +# Parameterize the interpolated power spectrum +fm2 = FOOOF(aperiodic_mode='knee', verbose=False) +fm2.report(freqs2, powers_int2) + +################################################################################################### +# Fitting Line Noise as Peaks +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# In some cases, you may also be able to simply allow the parameterization to peaks to the +# line noise and harmonics. By simply fitting the line noise as peaks, the model can deal +# with the peaks in order to accurately fit the aperiodic component. +# +# These peaks are of course not to be analyzed, but once the model has been fit, you can +# simply ignore them. There should generally be no issue with fitting and having them in +# the model, and allowing the model to account for these peaks typically helps the model +# better fit the rest of the data. +# +# Below we can see that the model does indeed work when fitting data with line noise peaks. +# + +################################################################################################### + +# Fit power spectrum models to original spectra +fm1.report(freqs1, powers1) +fm2.report(freqs2, powers2) + +################################################################################################### +# The Problem with Bandstop Filtering +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# A common approach for getting rid of line noise activity is to use bandstop filtering to +# remove activity at the line noise frequencies. Such a filter effectively set the power +# of these frequencies to be approximately zero. +# +# Unfortunately, this doesn't work very well with spectral parameterization, since the +# parameterization algorithm tries to fit each power value as either part of the aperiodic +# component, or as an overlying peak. Frequencies that have filtered out are neither, and +# the model has trouble, as it and has no concept of power values below the aperiodic component. +# +# In practice, this means that the "missing" power will impact the fit, and pull down the +# aperiodic component. One way to think of this is that the power spectrum model can deal with, +# and even expects, 'positive outliers' above the aperiodic (these are considered 'peaks'), but +# not 'negative outliers', or values below the aperiodic, as there is no expectation of this +# happening in the model. +# +# In the following example, we can see how bandstop filtering negatively impacts fitting. +# Because of this, for the purposes of spectral parameterization, bandstop filters are not +# recommended as a way to remove line noise. +# +# Note that if one has already applied a bandstop filter, then you can still +# apply the interpolation from above. +# + +################################################################################################### + +# General settings for the simulation +n_seconds = 30 +fs = 1000 + +# Define the settings for the simulated signal +components = {'sim_powerlaw' : {'exponent' : -1.5}, + 'sim_oscillation' : [{'freq' : 10}, {'freq' : 60}]} +comp_vars = [0.5, 1, 1] + +# Simulate a time series +sig = sim_combined(n_seconds, fs, components, comp_vars) + +################################################################################################### + +# Bandstop filter the signal to remove line noise frequencies +sig_filt = filter_signal(sig, fs, 'bandstop', (57, 63), + n_seconds=2, remove_edges=False) + +################################################################################################### + +# Compute a power spectrum of the simulated signal +freqs, powers_pre = trim_spectrum(*compute_spectrum(sig, fs), [3, 75]) +freqs, powers_post = trim_spectrum(*compute_spectrum(sig_filt, fs), [3, 75]) + +################################################################################################### + +# Plot the spectrum of the data, pre and post bandstop filtering +plot_spectra(freqs, [powers_pre, powers_post], log_powers=True, + labels=['Pre-Filter', 'Post-Filter']) + +################################################################################################### +# +# In the above, we can see that the the bandstop filter removes power in the filtered range, +# leaving a "dip" in the power spectrum. This dip causes issues with subsequent fitting. +# + +################################################################################################### + +# Initialize and fit a power spectrum model +fm = FOOOF() +fm.report(freqs, powers_post) + +################################################################################################### +# +# In order to try and capture the data points in the "dip", the power spectrum model +# gets 'pulled' down, leading to an inaccurate fit of the aperiodic component. This is +# why fitting frequency regions that included frequency regions that have been filtered +# out is not recommended. +# \ No newline at end of file diff --git a/fooof/tests/utils/test_data.py b/fooof/tests/utils/test_data.py index fb4d7ad67..d16b4c27c 100644 --- a/fooof/tests/utils/test_data.py +++ b/fooof/tests/utils/test_data.py @@ -21,11 +21,31 @@ def test_trim_spectrum(): def test_interpolate_spectrum(): + # Test with single buffer exclusion zone freqs, powers = gen_power_spectrum(\ [1, 75], [1, 1], [[10, 0.5, 1.0], [60, 2, 0.1]]) - freqs_out, powers_out = interpolate_spectrum(freqs, powers, [58, 62]) + exclude = [58, 62] + + freqs_out, powers_out = interpolate_spectrum(freqs, powers, exclude) assert np.array_equal(freqs, freqs_out) assert np.all(powers) assert powers.shape == powers_out.shape + mask = np.logical_and(freqs >= exclude[0], freqs <= exclude[1]) + assert powers[mask].sum() > powers_out[mask].sum() + + # Test with multiple buffer exclusion zones + freqs, powers = gen_power_spectrum(\ + [1, 150], [1, 100, 1], [[10, 0.5, 1.0], [60, 1, 0.1], [120, 0.5, 0.1]]) + + exclude = [[58, 62], [118, 122]] + + freqs_out, powers_out = interpolate_spectrum(freqs, powers, exclude) + assert np.array_equal(freqs, freqs_out) + assert np.all(powers) + assert powers.shape == powers_out.shape + + for f_range in exclude: + mask = np.logical_and(freqs >= f_range[0], freqs <= f_range[1]) + assert powers[mask].sum() > powers_out[mask].sum() diff --git a/fooof/utils/data.py b/fooof/utils/data.py index 8410b1985..97b6c738f 100644 --- a/fooof/utils/data.py +++ b/fooof/utils/data.py @@ -1,5 +1,7 @@ """Utilities for working with data and models.""" +from itertools import repeat + import numpy as np ################################################################################################### @@ -60,9 +62,10 @@ def interpolate_spectrum(freqs, powers, interp_range, buffer=3): Frequency values for the power spectrum. powers : 1d array Power values for the power spectrum. - interp_range : list of float + interp_range : list of float or list of list of float Frequency range to interpolate, as [lowest_freq, highest_freq]. - buffer : int + If a list of lists, applies each as it's own interpolation range. + buffer : int or list of int The number of samples to use on either side of the interpolation range, that are then averaged and used to calculate the interpolation. @@ -101,23 +104,35 @@ def interpolate_spectrum(freqs, powers, interp_range, buffer=3): >>> freqs, powers = interpolate_spectrum(freqs, powers, [58, 62]) """ - # Get the set of frequency values that need to be interpolated - interp_mask = np.logical_and(freqs >= interp_range[0], freqs <= interp_range[1]) - interp_freqs = freqs[interp_mask] + # If given a list of interpolation zones, recurse to apply each one + if isinstance(interp_range[0], list): + buffer = repeat(buffer) if isinstance(buffer, int) else buffer + for interp_zone, cur_buffer in zip(interp_range, buffer): + freqs, powers = interpolate_spectrum(freqs, powers, interp_zone, cur_buffer) + + # Assuming list of two floats, interpolate a single frequency range + else: + + # Take a copy of the array, to not change original array + powers = np.copy(powers) + + # Get the set of frequency values that need to be interpolated + interp_mask = np.logical_and(freqs >= interp_range[0], freqs <= interp_range[1]) + interp_freqs = freqs[interp_mask] - # Get the indices of the interpolation range - ii1, ii2 = np.flatnonzero(interp_mask)[[0, -1]] + # Get the indices of the interpolation range + ii1, ii2 = np.flatnonzero(interp_mask)[[0, -1]] - # Extract & log the requested range of data to use around interpolated range - xs1 = np.log10(freqs[ii1-buffer:ii1]) - xs2 = np.log10(freqs[ii2:ii2+buffer]) - ys1 = np.log10(powers[ii1-buffer:ii1]) - ys2 = np.log10(powers[ii2:ii2+buffer]) + # Extract & log the requested range of data to use around interpolated range + xs1 = np.log10(freqs[ii1-buffer:ii1]) + xs2 = np.log10(freqs[ii2:ii2+buffer]) + ys1 = np.log10(powers[ii1-buffer:ii1]) + ys2 = np.log10(powers[ii2:ii2+buffer]) - # Linearly interpolate, in log-log space, between averages of the extracted points - vals = np.interp(np.log10(interp_freqs), - [np.median(xs1), np.median(xs2)], - [np.median(ys1), np.median(ys2)]) - powers[interp_mask] = np.power(10, vals) + # Linearly interpolate, in log-log space, between averages of the extracted points + vals = np.interp(np.log10(interp_freqs), + [np.median(xs1), np.median(xs2)], + [np.median(ys1), np.median(ys2)]) + powers[interp_mask] = np.power(10, vals) return freqs, powers