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_spectrum , 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_spectrum (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