Merge pull request #203 from fooof-tools/linen

[DOC] - Add a line noise example

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',

examples/processing/README.txt

@@ -0,0 +1,4 @@
+Processing
+----------
+
+Examples on how to process data related to spectral parameterization.

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.
+#

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()

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