|
| 1 | +""" |
| 2 | +Dealing with Line Noise |
| 3 | +======================= |
| 4 | +
|
| 5 | +This example covers strategies for dealing with line noise. |
| 6 | +""" |
| 7 | + |
| 8 | +################################################################################################### |
| 9 | + |
| 10 | +# sphinx_gallery_thumbnail_number = 2 |
| 11 | + |
| 12 | +# Import the spectral parameterization object and utilities |
| 13 | +from fooof import FOOOF |
| 14 | +from fooof.plts import plot_spectra |
| 15 | +from fooof.utils import trim_spectrum, interpolate_spectrum |
| 16 | + |
| 17 | +# Import simulation functions to create some example data |
| 18 | +from fooof.sim.gen import gen_power_spectrum |
| 19 | + |
| 20 | +# Import NeuroDSP functions for simulating & processing time series |
| 21 | +from neurodsp.sim import sim_combined |
| 22 | +from neurodsp.filt import filter_signal |
| 23 | +from neurodsp.spectral import compute_spectrum |
| 24 | + |
| 25 | +################################################################################################### |
| 26 | +# Line Noise Peaks |
| 27 | +# ---------------- |
| 28 | +# |
| 29 | +# Neural recordings typically have power line artifacts, at either 50 or 60 Hz, depending on |
| 30 | +# where the data were collected, which can impact spectral parameterization. |
| 31 | +# |
| 32 | +# In this example, we explore some options for dealing with line noise artifacts. |
| 33 | +# |
| 34 | +# Interpolating Line Noise Peaks |
| 35 | +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 36 | +# |
| 37 | +# One approach is to interpolate away line noise peaks, in the frequency domain. This |
| 38 | +# approach simply gets rid of the peaks, interpolating the data to maintain the 1/f |
| 39 | +# character of the data, allowing for subsequent fitting. |
| 40 | +# |
| 41 | +# The :func:`~fooof.utils.interpolate_spectrum` function allows for doing simple |
| 42 | +# interpolation. Given a narrow frequency region, this function interpolates the spectrum, |
| 43 | +# such that the 'peak' of the line noise is removed. |
| 44 | +# |
| 45 | + |
| 46 | +################################################################################################### |
| 47 | + |
| 48 | +# Generate an example power spectrum, with line noise |
| 49 | +freqs1, powers1 = gen_power_spectrum([3, 75], [1, 1], |
| 50 | + [[10, 0.75, 2], [60, 1, 0.5]]) |
| 51 | + |
| 52 | +# Visualize the generated power spectrum |
| 53 | +plot_spectra(freqs1, powers1, log_powers=True) |
| 54 | + |
| 55 | +################################################################################################### |
| 56 | +# |
| 57 | +# In the plot above, we have an example spectrum, with some power line noise. |
| 58 | +# |
| 59 | +# To prepare this data for fitting, we can interpolate away the line noise region. |
| 60 | +# |
| 61 | + |
| 62 | +################################################################################################### |
| 63 | + |
| 64 | +# Interpolate away the line noise region |
| 65 | +interp_range = [58, 62] |
| 66 | +freqs_int1, powers_int1 = interpolate_spectrum(freqs1, powers1, interp_range) |
| 67 | + |
| 68 | +################################################################################################### |
| 69 | + |
| 70 | +# Plot the spectra for the power spectra before and after interpolation |
| 71 | +plot_spectra(freqs1, [powers1, powers_int1], log_powers=True, |
| 72 | + labels=['Original Spectrum', 'Interpolated Spectrum']) |
| 73 | + |
| 74 | +################################################################################################### |
| 75 | +# |
| 76 | +# As we can see in the above, the interpolation removed the peak from the data. |
| 77 | +# |
| 78 | +# We can now go ahead and parameterize the spectrum. |
| 79 | +# |
| 80 | + |
| 81 | +################################################################################################### |
| 82 | + |
| 83 | +# Initialize a power spectrum model |
| 84 | +fm1 = FOOOF(verbose=False) |
| 85 | +fm1.report(freqs_int1, powers_int1) |
| 86 | + |
| 87 | +################################################################################################### |
| 88 | +# Multiple Interpolation Regions |
| 89 | +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 90 | +# |
| 91 | +# Line noise artifacts often also display harmonics, such that when analyzing broader |
| 92 | +# frequency ranges, there may be multiple peaks that need to be interpolated. |
| 93 | +# |
| 94 | +# This can be done by passing in multiple interpolation regions to |
| 95 | +# :func:`~fooof.utils.interpolate_spectrum`, which we will do in the next example. |
| 96 | +# |
| 97 | + |
| 98 | +################################################################################################### |
| 99 | + |
| 100 | +# Generate an example power spectrum, with line noise & harmonics |
| 101 | +freqs2, powers2 = gen_power_spectrum([1, 150], [1, 500, 1.5], |
| 102 | + [[10, 0.5, 2], [60, 0.75, 0.5], [120, 0.5, 0.5]]) |
| 103 | + |
| 104 | +# Interpolate away the line noise region & harmonics |
| 105 | +interp_ranges = [[58, 62], [118, 122]] |
| 106 | +freqs_int2, powers_int2 = interpolate_spectrum(freqs2, powers2, interp_ranges) |
| 107 | + |
| 108 | +################################################################################################### |
| 109 | + |
| 110 | +# Plot the power spectrum before and after interpolation |
| 111 | +plot_spectra(freqs2, [powers2, powers_int2], log_powers=True, |
| 112 | + labels=['Original Spectrum', 'Interpolated Spectrum']) |
| 113 | + |
| 114 | +################################################################################################### |
| 115 | + |
| 116 | +# Parameterize the interpolated power spectrum |
| 117 | +fm2 = FOOOF(aperiodic_mode='knee', verbose=False) |
| 118 | +fm2.report(freqs2, powers_int2) |
| 119 | + |
| 120 | +################################################################################################### |
| 121 | +# Fitting Line Noise as Peaks |
| 122 | +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 123 | +# |
| 124 | +# In some cases, you may also be able to simply allow the parameterization to peaks to the |
| 125 | +# line noise and harmonics. By simply fitting the line noise as peaks, the model can deal |
| 126 | +# with the peaks in order to accurately fit the aperiodic component. |
| 127 | +# |
| 128 | +# These peaks are of course not to be analyzed, but once the model has been fit, you can |
| 129 | +# simply ignore them. There should generally be no issue with fitting and having them in |
| 130 | +# the model, and allowing the model to account for these peaks typically helps the model |
| 131 | +# better fit the rest of the data. |
| 132 | +# |
| 133 | +# Below we can see that the model does indeed work when fitting data with line noise peaks. |
| 134 | +# |
| 135 | + |
| 136 | +################################################################################################### |
| 137 | + |
| 138 | +# Fit power spectrum models to original spectra |
| 139 | +fm1.report(freqs1, powers1) |
| 140 | +fm2.report(freqs2, powers2) |
| 141 | + |
| 142 | +################################################################################################### |
| 143 | +# The Problem with Bandstop Filtering |
| 144 | +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 145 | +# |
| 146 | +# A common approach for getting rid of line noise activity is to use bandstop filtering to |
| 147 | +# remove activity at the line noise frequencies. Such a filter effectively set the power |
| 148 | +# of these frequencies to be approximately zero. |
| 149 | +# |
| 150 | +# Unfortunately, this doesn't work very well with spectral parameterization, since the |
| 151 | +# parameterization algorithm tries to fit each power value as either part of the aperiodic |
| 152 | +# component, or as an overlying peak. Frequencies that have filtered out are neither, and |
| 153 | +# the model has trouble, as it and has no concept of power values below the aperiodic component. |
| 154 | +# |
| 155 | +# In practice, this means that the "missing" power will impact the fit, and pull down the |
| 156 | +# aperiodic component. One way to think of this is that the power spectrum model can deal with, |
| 157 | +# and even expects, 'positive outliers' above the aperiodic (these are considered 'peaks'), but |
| 158 | +# not 'negative outliers', or values below the aperiodic, as there is no expectation of this |
| 159 | +# happening in the model. |
| 160 | +# |
| 161 | +# In the following example, we can see how bandstop filtering negatively impacts fitting. |
| 162 | +# Because of this, for the purposes of spectral parameterization, bandstop filters are not |
| 163 | +# recommended as a way to remove line noise. |
| 164 | +# |
| 165 | +# Note that if one has already applied a bandstop filter, then you can still |
| 166 | +# apply the interpolation from above. |
| 167 | +# |
| 168 | + |
| 169 | +################################################################################################### |
| 170 | + |
| 171 | +# General settings for the simulation |
| 172 | +n_seconds = 30 |
| 173 | +fs = 1000 |
| 174 | + |
| 175 | +# Define the settings for the simulated signal |
| 176 | +components = {'sim_powerlaw' : {'exponent' : -1.5}, |
| 177 | + 'sim_oscillation' : [{'freq' : 10}, {'freq' : 60}]} |
| 178 | +comp_vars = [0.5, 1, 1] |
| 179 | + |
| 180 | +# Simulate a time series |
| 181 | +sig = sim_combined(n_seconds, fs, components, comp_vars) |
| 182 | + |
| 183 | +################################################################################################### |
| 184 | + |
| 185 | +# Bandstop filter the signal to remove line noise frequencies |
| 186 | +sig_filt = filter_signal(sig, fs, 'bandstop', (57, 63), |
| 187 | + n_seconds=2, remove_edges=False) |
| 188 | + |
| 189 | +################################################################################################### |
| 190 | + |
| 191 | +# Compute a power spectrum of the simulated signal |
| 192 | +freqs, powers_pre = trim_spectrum(*compute_spectrum(sig, fs), [3, 75]) |
| 193 | +freqs, powers_post = trim_spectrum(*compute_spectrum(sig_filt, fs), [3, 75]) |
| 194 | + |
| 195 | +################################################################################################### |
| 196 | + |
| 197 | +# Plot the spectrum of the data, pre and post bandstop filtering |
| 198 | +plot_spectra(freqs, [powers_pre, powers_post], log_powers=True, |
| 199 | + labels=['Pre-Filter', 'Post-Filter']) |
| 200 | + |
| 201 | +################################################################################################### |
| 202 | +# |
| 203 | +# In the above, we can see that the the bandstop filter removes power in the filtered range, |
| 204 | +# leaving a "dip" in the power spectrum. This dip causes issues with subsequent fitting. |
| 205 | +# |
| 206 | + |
| 207 | +################################################################################################### |
| 208 | + |
| 209 | +# Initialize and fit a power spectrum model |
| 210 | +fm = FOOOF() |
| 211 | +fm.report(freqs, powers_post) |
| 212 | + |
| 213 | +################################################################################################### |
| 214 | +# |
| 215 | +# In order to try and capture the data points in the "dip", the power spectrum model |
| 216 | +# gets 'pulled' down, leading to an inaccurate fit of the aperiodic component. This is |
| 217 | +# why fitting frequency regions that included frequency regions that have been filtered |
| 218 | +# out is not recommended. |
| 219 | +# |
0 commit comments