Schrödinger filtering: a precise EEG despiking technique for EEG-fMRI gradient artifact

In EEG data acquired in the presence of fMRI, gradient-related spike artifacts contaminate the signal following the common preprocessing step of average artifact subtraction. Spike artifacts compromise EEG data quality since they overlap with the EEG signal in frequency, thereby confounding frequency-based inferences on activity. As well, spike artifacts can inflate or deflate correlations among time series, thereby confounding inferences on functional connectivity. We present Schrödinger filtering, which uses the Schrödinger equation to decompose the spike-containing input. The basis functions of the decomposition are localized and pulse-shaped, and selectively capture the various input peaks, with the spike components clustered at the beginning of the spectrum. Schrödinger filtering automatically subtracts the spike components from the data. On real and simulated data, we show that Schrödinger filtering (1) simultaneously accomplishes high spike removal and high signal preservation without affecting evoked activity, and (2) reduces spurious pairwise correlations in spontaneous activity. In these regards, Schrödinger filtering was significantly better than three other despiking techniques: median filtering, amplitude thresholding, and wavelet denoising. These results encourage the use of Schrödinger filtering in future EEG-fMRI pipelines, as well as in other spike-related applications (e.g., fMRI motion artifact removal or action potential extraction).


Introduction
Simultaneous acquisition of electroencephalography and functional magnetic resonance imaging (EEG-fMRI) has been popular for over two decades ( Abreu et al., 2018 ;Gotman et al., 2006 ;Huang-Hellinger et al., 1995 ;Ives et al., 2001 ;Lemieux et al., 2001 ) due to each modality's complementary strengths. EEG, using scalp electrodes, measures changes in electric potential associated with ion fluxes from action potentials and postsynaptic potentials at the cortical surface ( Beres, 2017 ). As the electric potential passes through the head's tissues before reaching the scalp electrodes, it gets blurred in both space and time to the extent that spatial resolution is poor while temporal resolution is high enough to study the classical frequency bands ranging up to the order of 100 Hz ( Burle et al., 2015 ;Niedermeyer, 2005 ). In fMRI, a time series of images is acquired, with the voxel size ranging from microns to centimeters ( Glover, 2011 ). The voxel intensity is spatially specific to the change in concentration of deoxygenated hemoglobin, providing an indirect measure of neural activity ( Ogawa et al., 1990 ) via hemodynamic processes ( Logothetis et al., 2001 ). The hemodynamic response is on the order of 1 s, giving fMRI poor temporal resolution regardless of how fast each image is acquired ( Glover, 2011 ). Thus, the relatively high temporal resolution and low spatial resolution of EEG complements the relatively Following the popular artifact removal technique of average artifact subtraction (AAS) (red), there is clearly residual artifact, with sharp spikes at the ends of the artifact-containing regions. Mild low-pass (400 Hz) filtering was applied here following AAS to accentuate the spikes. (C) Power spectra, computed using Welch's method, for the time series depicted in panel B, but including all epochs. The gradient artifact peaks roughly every 7 Hz, which corresponds to the fMRI fundamental slice acquisition frequency (21 slices per volume / 3 s per volume) and harmonics. The persistence of these peaks after AAS is largely attributed to the persistent time-domain spikes.

Spike artifacts as confounds in EEG data
A common feature of the residual post-AAS gradient artifact is timedomain spikes. A spike is a high-amplitude transient artifact in a time series. Spike artifacts are a major nuisance in EEG data for two reasons.
Firstly, since spikes overlap with the EEG signal in frequency, they potentially confound interpretations of frequency-based EEG activity. Spikes are highly localized in time, and, as a consequence, their Fourier transform is highly delocalized. Treating a spike as a narrow Gaussian, it is clear that the Fourier transform amplitude spectrum is a wide Gaussian centered at 0 Hz and spanning the EEG frequency bands. This is true whether the spike is positive or negative, and regardless of the spike's position in the time series ( Fig. 2 A). The more spikes there are, the higher-amplitude their collective Fourier transform. The more heterogeneous their spacing in time, the less smooth their collective Fourier transform, making it difficult to distinguish from the coloured noise is the same wide Gaussian centered at 0 Hz and spanning the EEG frequency bands, up to roughly 150 Hz in this example. (E) Simulated EEG, modeled with evoked sinusoidal activity, as well as pink noise, alpha-band noise, and white noise, with narrow Gaussian spikes of random position and polarity superimposed. (F) The Fourier transforms of the separate simulated EEG (blue) and spikes (red) from panel E. The EEG spectrum is dominated by pink noise, while the spectrum from the spikes is a noisy-looking Gaussian comparable in magnitude and bandwidth. structure of EEG signal ( Fig. 2 B). Additional Fourier components come from any periodicity of the spikes' occurrence in the time series, which also overlaps with the EEG signal. Thus, frequency-domain bandpass filtering is inadequate for spike removal ( Boroujeni et al., 2020 ).
Secondly, spikes are detrimental to EEG correlation analyses commonly done for explorations of functional connectivity. For spikes found at the same locations in different time series due to a common confounding source, which, in this case, is the post-AAS gradient artifact, correlations are inflated. Conversely, correlations are deflated if spike artifacts are present in one time series but not others.
For these reasons, a robust time-domain despiking technique is needed.

Schrödinger filtering
We present a new signal processing technique, called Schrödinger filtering, which accurately removes gradient artifact spikes and preserves EEG signal. Schrödinger filtering uses semi-classical signal analysis (SCSA) ( Laleg-Kirati et al., 2013 ), in which an input signal is treated as an attractive potential within the Schrödinger operator. From the discrete spectrum of the operator, the input is decomposed into a series of smooth pulse-shaped functions called Schrödinger components. Schrödinger components are squared eigenfunctions scaled by their respective eigenvalues, making SCSA analogous to the Fourier transform, which uses complex sinusoids to decompose an input. Low-order Schrödinger components capture energy-dense peaks of the input, making this technique well-suited to applications of peak retention or removal. Spike removal by Schrödinger filtering is the subtraction of these low-order components from the spike-contaminated input.
The implemented Schrödinger filtering algorithm is automatic in determining the boundary between these low-order spike components and the rest of the spectrum. The algorithm does not use any templates or references of the artifact or signal. It is therefore not a closed-loop adaptive filter, unlike many existing gradient artifact removal algorithms ( Allen et al., 2000 ;Niazy et al., 2005 ). Although closed-loop adaptive filters are very powerful if the reference signal is known with a high degree of accuracy, a means of acquiring an accurate reference of the artifact ( Chowdhury et al., 2014 ;Steyrl et al., 2017 ) is not yet widespread. Instead, the proposed algorithm uses kernel density estima-tion on the SCSA eigenspectrum, which robustly and automatically clusters the components into spike-containing and non-spike-containing.

Theory
The one-dimensional semi-classical Schrödinger operator is given by where ℎ ∈ ℝ > 0 is a scalar parameter and ( ) ∈ ℝ ≥ 0 is a one-dimensional input signal of time . The Schrödinger equation is an eigenvalue problem: where is the eigenvalue and ( ) is the corresponding eigenfunction on which ℎ ( ) operates (Laleg-Kirati et al., 2013) . In general, there is a continuous spectrum for ≥ 0 and a discrete spectrum for negative such that where ,ℎ > 0 and ℎ is the number of negative eigenvalues. After disregarding the nonnegative , ( ) is represented by the series (Helffer and Laleg-Kirati, 2011 ;Laleg-Kirati et al., 2013 ), where each term 4 ℎ ℎ, 2 ℎ, ( ) is a Schrödinger component. Schrödinger components are pulse-shaped, giving them affinity in capturing energy-dense input peaks. A given component contains one or more pulses symmetrically distributed about a particular input peak. Furthermore, the more pulses belonging to a component, the less energy-dense it is, and the higher its order in the series of Eq. (1) . This negative gradient of energy density is represented in the spectrum of eigenvalues { − 2 ℎ, } ℎ =1 , which is a monotonically decreasing function ( Fig. 3 A).
The value of ℎ has a critical impact on the manner in which SCSA reconstructs ( ) . For a high ℎ -value, few Schrödinger components are produced in the series of Eq. (1) . These components capture the largest The series of Schrödinger components begins with a single-pulsed, energy-dense component centered at the input peak. This component has the greatest eigenvalue in the series. The next component has one more pulse, with both of its pulses distributed symmetrically about the first component, and is slightly less energy-dense, with a slightly lesser eigenvalue. Progressively higher-order components are progressively less energy-dense with progressively more pulses and progressively lesser eigenvalues (C) , while remaining symmetric about the first component. (D) The first four components, shown for detail. (E) A noisy signal contains three signal peaks, the identities of which are not readily apparent. (F) From SCSA decomposition, the first three components (red) correctly locate the signal peaks (green). The eigenspectrum (G) shows a distinct separation between the first three components and the rest of the spectrum, highlighting SCSA's ability to detect peaks. (H) A noisy input with a single signal peak (black) is reconstructed using SCSA with an h -value of 36 (red), corresponding to 15 components (I) . Though the fine details of the noise are not captured, the prominent signal peak is. (J) The same noisy input signal (black) from panel H, consisting of one signal peak, is reconstructed with the lower h -value of 1 (red). As a result, many more components are generated (roughly 460 in total) (K) , and the input is nearly entirely reconstructed, including the fine details. peaks -i.e. , those with the densest energy -while finer details of the input are not captured ( Fig. 3 H). For increasingly lower values of ℎ , more components are generated in the series, with the higher-order components capturing the fine details of ( ) , enabling high-fidelity reconstruction ( Fig. 3 J). For an input of infinite resolution, SCSA reconstruction is perfect in the limit of ℎ tending to zero (Laleg-Kirati et al., 2013) . For an input of finite resolution, such as real-world data, this limit does not hold, and the fidelity of the reconstruction is limited by the data resolution.
SCSA can be used as a smoothing technique, with smoothing parameter ℎ in Eq. (1) , that well-preserves input peak structure ( Li and Laleg-Kirati, 2019 ) . If the input peaks correspond to signal, the SCSAsmoothed output, using the appropriate ℎ -value, is a denoised version ( Fig. 3 H). SCSA has been used to denoise MRI and magnetic resonance spectroscopy data ( Chahid et al., 2019 ;Laleg-Kirati et al., 2016 ). Alternatively, if the input peaks correspond to artifact, the SCSA-smoothed output is a precise representation of the artifact that may be subtracted from the input. SCSA has been used in this manner for magnetic resonance spectroscopy water peak suppression (Chahid et al., 2019) . We coin the term Schrödinger filtering for the general class of data filtering methods that use SCSA. In this work, Schrödinger filtering constitutes estimating any existing spike content and subtracting it from the data.

Despiking via Schrödinger filtering
With the appropriate choice of ℎ in Eq. (1) , the spike-containing input ( Fig. 4 A) is processed so the spike content is nearly unaffected while the rest of the input is almost entirely smoothed away ( Fig. 4 E,F). Let this ℎ -value be denoted by ℎ * . ℎ * is found by maximizing the kurtosis of the output. Kurtosis is a measure of the tailedness of a function ( Westfall, 2014 ), and is a suitable proxy for this purpose since spikecontaining inputs appear heavily tailed.  4. The Schrödinger filtering algorithm. Following preprocessing, which includes detrending and epoching, the spike-containing input (A) is processed in two separate streams. In one stream, the input is half-wave-rectified (B) , and in the other stream, it is inverted (C) and then half-wave-rectified (D) . These intermediate outputs are separately smoothed by semiclassical signal analysis (SCSA) to form a template of the spikes that contains very little signal (E, F) . The spikes (K, L) are contained in the first few components of each template's eigenspectrum (G, H) , the identities of which are clustered via kernel density estimation, a technique that estimates the probability density functions of the spectra (I, J) . The spike components of the negative stream are inverted (M) , and the spike components from both branches are subtracted from the original input, giving the spike-free output (N) .
Following SCSA smoothing with ℎ = ℎ * , the spikes are clearly distinguished from the rest of the output. Thanks to a few more unique characteristics of SCSA, this output may be processed to further extract the spikes. For a Schrödinger component of arbitrary index , there is uniquely one eigenvalue , which is proportional to the component's energy. Since Schrödinger components are localized, pulse-shaped func-tions, a high -value corresponds to an energy-dense peak such as an artifact spike. Moreover, since the SCSA eigenspectrum is a monotonically decreasing function, the spike content is found in the first components ( Fig. 4 G,H). Despiking entails subtraction of these components from the unsmoothed spike-containing input.
Determination of the value of is a clustering problem between the first eigenvalues and the rest of the spectrum ( Fig. 4 G,H). For the SCSA spectra in the present work, which are one-dimensional, the probability density function is expected to be multimodal, with one mode representing the first components, which capture the spikes, and one or more separate modes representing the remaining components, which do not capture the spikes. Thus, kernel density estimation, whereby the probability density function is estimated, enables automatic identification of ( Fig. 4 I,J). The algorithm is outlined in Section 3.1 , and a schematic is provided in Fig. 4 .

Schrödinger filtering algorithm
1. Per channel: a. Detrend : Low-frequency ( < 30 Hz) fluctuations in the data can be of considerable magnitude, and can result in spike underestimation if the spike is on a trough of a fluctuation, or overestimation if the spike is on a crest. A high-pass filter robustly flattens such data while preserving spike structure. Additionally, high-frequency ( > 700 Hz) fluctuations, in extreme cases, are of comparable magnitude to the spikes, potentially causing spike underestimation. A low-pass filter cleans away these fluctuations while also preserving spike structure. Lastly, high mains noise can compromise despiking if it is not notch-filtered. Note that the detrended data are only used for spike definition, and that spike components are subtracted from the nondetrended data. b. Chunk : The runtime of the present SCSA code grows significantly with input length. Chunking the data into epochs reduces computational load. Note that the resolution of the SCSA decomposition is limited by the data resolution but not the chunk size. 2. Per epoch: a. Rectify : Since Schrödinger components are exclusively nonnegative, the positive and negative peaks are analyzed in two parallel branches. In the positive branch, the data are half-wave rectified -i.e. , nulled for all values less than zero ( Fig. 4 B). In the negative branch, the data are inverted -i.e. , reflected about the x-axis ( Fig. 4 C) -and then half-wave rectified ( Fig. 4 D). 3. Per branch: a. SCSA-smooth : SCSA reconstruction using ℎ = ℎ * in Eq. (1) preserves spike structure while smoothing away most of the rest of the input ( Section 2.1 ) ( Fig. 4 E,F), and ℎ * is determined by maximizing the kurtosis of the SCSA reconstruction. b. Kernel density estimation : The estimated probability density function ( Fig. 4 I,J) of the SCSA eigenspectrum ( Fig. 4 G,H) depends on the bandwidth of the smoothing kernel, which determines the smoothness of the estimated function. The bandwidth is constrained according to Silverman's rule of thumb for Gaussian distributions ( Silverman, 1986 ). c. Select components : In a spike-containing probability density function, the boundary between the spike and non-spike eigenvalues is constrained as the local minimum between the major mode on the left and the first mode immediately right ( Fig. 4 I,J). All components whose eigenvalue is greater than this boundary are the spikecontaining components ( Fig. 4 K,L). Note that no such boundary exists for data with no spikes. d. If negative branch, invert the selected components ( Fig. 4 M). e. Subtract the selected components from the non-detrended, unsmoothed, un-rectified input.

Fig. 5. (A)
Example of simulated EEG data components. A weighted sum of white noise (orange), pink noise (blue), and alpha-band noise (yellow) is referred to as the noise floor and represents the EEG data in the absence of evoked activity or artifact, where the white noise represents hardware noise. Evoked activity is simulated as ten-second spurts of alpha-band activity purple). (B) The root mean square of this evoked activity relative to the root mean square of the noise floor is topographically distributed so that activity is greatest near visual cortex. Fig. 4 is a schematic of the algorithm.

Data
Real data 1. A freely available dataset from the FASTR (Niazy et al., 2005) toolbox website (fsl.fmrib.ox.ac.uk/eeglab/fmribplugin/#tutorial) was used. During scanning, the subject opened and closed their eyes in consecutive alternating 10-s intervals. The sampling rate was 2048 Hz and EEG-MRI clock synchronization was not implemented. This dataset will be referred to as the FMRIB dataset. 2. The dataset described by Georgie et al. (2018) , available through the Open Science Framework (https://osf.io), was used. Data were sampled at 5 kHz and include runs from 17 subjects acquired outside the MRI scanner as well as in the scanner during fMRI. The data acquired outside the scanner comprised one run per subject, and the data acquired in the scanner during fMRI comprised up to five runs per subject. EEG-MRI clock synchronization was implemented. The subjects engaged in tasks of decision-making in response to presented images ( Georgie et al., 2018Ostwald et al., 2012.

Simulated data
The code for all simulated data is available at github.com/gbenigno/schrodinger_filtering. 1. For a distilled demonstration of despiking, white noise was simulated with narrow Gaussian spikes superimposed. Here, the white noise represents the signal. To give the spikes heterogeneity, their amplitudes and locations were randomly generated. One-thousand trials were performed. 2. In a model of spontaneous activity and noise closely following that of the EEGSourceSim toolbox ( Barzegaran et al., 2019 ), broadband and alpha-band activity, which, to a first order, appear stochastic in EEG data, were modeled as pink noise and alpha noise, respectively. Additionally, white noise represented hardware noise. These three components, collectively referred to as the noise floor, were weighted per-channel according to the resting-state EEG data in the FMRIB dataset recorded before scanning began but while the subject was in the scanner. The EEG was simulated in the same 30 channels as the FMRIB dataset. Evoked activity was simulated as spurts of increased alpha-band (8-12 Hz) activity, with this activity strongest in the electrodes closest to visual cortex ( Fig. 5 ). Post-AAS gradient artifact spikes were modeled as narrow Gaussians. A pair of one positive and one negative spike was placed at the beginning and end of each slice epoch in the simulated data since spikes were common at these locations in the post-AAS FMRIB data. The spike amplitude and width were given by the group statistics of the outlying peaks of the post-AAS FMRIB data. Spikes were superimposed in three separate realizations -one with high-amplitude spikes, one with mediumamplitude spikes, and one with low-amplitude spikes -to represent differing degrees of AAS effectiveness due to a factor such as the rate of sampling of the raw gradient artifact. Here, high-, medium-, and low-amplitude spikes correspond to 100%, 75%, and 50% of the spike amplitudes observed in the FMRIB data, respectively. 3. Two EEG time series were simulated, also following EEGSourceSim ( Barzegaran et al., 2019 ), and were made correlated with one another. An appropriate correlation matrix was synthesized and the data were transformed by the Cholesky decomposition of the correlation matrix. This was done for correlation values of 0.25, 0.5, and 0.75. Spikes of differing amplitude were then superimposed as in simulated dataset 2.

Processing steps
The code for these steps is available at github.com/gbenigno/schrodinger_filtering. AAS was applied on the raw gradient-artifact containing data using the FACET ( Glaser et al., 2013 ) toolbox, except for the decision-making dataset, for which efficient and equivalent custom code was used. The following four steps were applied in parallel: (1) Schrödinger filtering as described in Section 3.1 ; (2) amplitude thresholding using the MATLAB (MathWorks, Natick, MA, USA) filloutliers function on the detrended input, which was then added to the trend; (3) median filtering using the MATLAB medfilt function on the detrended input, which was then added to the trend; and (4) wavelet denoising using the MATLAB function wdenoise on the detrended input, which was then subtracted from the non-detrended input. For the decision-making dataset, preprocessing following AAS and preceding despiking comprised re-referencing to the average channel time series, baseline correction, and peri-stimulus window extraction. Low-pass filtering was not performed.

Performance metric
High despiking performance corresponds to simultaneous high spike removal and high signal preservation, and is directly measured using the ground truths. Let the EEG and spike artifact content in the processed output be denoted by ̃ ( ) and ̃ ( ) , respectively. The signal-toartifact ratio, analogous to the canonical formula for signal-to-noise ratio ( Johnson, 2006 ) where is the number of samples.

Schrödinger filtering robustly despikes
In the first simulation, a simple model of spikes and signal as narrow Gaussians and white noise, respectively, facilitated a distilled demonstration of Schrödinger filtering of an arbitrary spike-containing signal that is without trends or evoked activity. Of the three compared techniques (median filtering, amplitude thresholding, and wavelet denoising), median filtering performed poorest, removing lots of signal and little artifact ( Fig. 6 D). Amplitude thresholding was good at signal preservation in general, yet was unable to remove the portions of the spikes below the threshold, and was unable to detect spikes of amplitude below threshold. As well, amplitude thresholding is prone to false positives where the signal surpasses the threshold ( Fig. 6 E). Wavelet denoising ( Fig. 6 C), among the three compared techniques, performed best, indicated by its signal-to-artifact ratio scores ( Fig. 6 F). Schrödinger filtering outperformed all compared techniques, most closely resembling the ground-truth signal ( Fig. 6 A,B), and having a signal-to-artifact ratio significantly higher than the other techniques in general ( Fig. 6 F). In Fig. 6 F, the pairwise scatter plot is shown of the signal-to-artifact ratio of Schrödinger filtering as a function of the signal-to-artifact ratio of a compared technique. One-thousand unique trials were simulated. Thus, there is one unique point for each group, and one-thousand points per group. Points along the line of zero intercept and unit slope correspond to trials for which the performance of Schrödinger filtering equals that of the compared technique. For wavelet denoising, which performed best among the three compared techniques, the points were clustered about a zero-intercept line of roughly slope three, indicating that Schrödinger filtering performed roughly three times better than wavelet denoising, and that Schrödinger filtering exceeded the other two techniques in performance by an even greater amount.
The FMRIB data, following AAS, have many spikes, and the spike magnitude is very high on average. This could be partly due to high subject motion. Additionally, this is likely due in some part to the fact that these data were acquired without synchronization of the EEG and MRI clocks, causing phase jitter in the gradient artifact waveform, and thereby hampering the performance of AAS. Thus, this dataset serves as an extreme test case. The goal of processing this dataset was to qualitatively show Schrödinger filtering's ability to selectively remove spikes and retain signal. In the time domain , it is clear that Schrödinger filtering despikes without affecting any other portions of the input. This corresponds to a cleaned power spectrum ( Fig. 7 D). The ten-second bursts of visual evoked alpha-band (8-12 Hz) activity are visible in the seven channels closest to the visual cortex for the post-AAS data, and even more so following Schrödinger filtering ( Fig. 7 A). These bursts are clear in the spectrograms both before and after Schrödinger filtering ( Fig. 7 C). Lastly, retention of the topographic distribution of this activity is also evident ( Fig. 7 B).
To complement the FMRIB dataset, a similar dataset was simulated. EEG signal comprised various noise components and alpha-band bursts, the latter of which were strongest in visual cortex. Simulated gradient artifact spikes were superimposed, with their locations, widths, amplitudes, and polarities derived from the FMRIB data. Schrödinger filtering performed best in general, most closely resembling the ground truth ( Fig. 8 A,D,G) and according to the signal-to-artifact scores ( Fig. 8 B,E,H). In particular, the logarithmically-spaced pairwise scatter plots clearly exhibit two clusters: one cluster about the zero-intercept line of unit slope, corresponding to comparable performance, and one cluster above the zero-intercept line of slope three, corresponding to out-performance by Schrödinger filtering by a factor greater than three. The majority of trials correspond to the latter cluster. For the third simulation, which featured low-amplitude spikes, Schrödinger filtering performed comparably to the other techniques on the majority of trials ( Fig. 8 H). This is likely because spikes were so low-amplitude in comparison to the EEG fluctuations that they were undetectable. Furthermore, the spikes were so low-amplitude that they do not appear as spikes and do not bear spike artifacts' spectral properties of meaningfully contaminating the EEG data. In these cases, the signal-to-artifact ratio underestimates performance. Lastly, the evoked alpha-band activity was not altered following Schrödinger filtering or the other techniques in any of the simulation variants ( Fig. 8 C,F,I; Fig. S1) The decision-making dataset contained few post-AAS spike artifacts, likely because the data were acquired with EEG-MRI clock synchronization. Schrödinger filtering exhibited high-quality despiking in these data ( Fig. 9 A) and did not alter the event-related potentials in terms of morphology ( Fig. 9 B) or topographic distribution ( Fig. 9 C).
Preservation of low-frequency evoked activity, such as the alphaband activity in this report, is not of concern with the presently implemented algorithm due to the bandpass (30-700 Hz) filtering performed in the detrending preprocessing step ( Section 3.1 ). Thus, the spike-containing Schrödinger components were generated from a filtered version of the input that was missing most of this evoked activity. This filtering was also applied for the other three despiking techniques as a matter of methodological congruence. Thus, all four techniques successfully did not alter the evoked activity in any of the aforementioned data.
These results, which span data from a range of acquisition-related parameters, support that Schrödinger filtering (1) is robust at simultaneous artifact removal and signal preservation and (2) does not affect evoked activity.

Schrödinger filtering reduces spurious correlations in broadband spontaneous activity
In the common case of the spike artifacts being coincident across electrode channel time series, correlations are inflated, and spike reduction should lower correlation measurements closer to the ground truth. As well, given the superficially stochastic nature of spontaneous EEG activity, correlation measurements of spontaneous EEG time series may be sensitive to small changes, such as signal loss following despiking. These two points motivate the hypothesis that correlations of spontaneous EEG activity contained in spike-contaminated data are more accurately measured following despiking that has simultaneously high spike removal and high signal preservation. Such despiking performance was quantified earlier using the signal-to-artifact ratio ( Figs. 6 and 8 ), for which Schrödinger filtering scored highest.
To explore this hypothesis directly, pairs of correlated spontaneous EEG time series were simulated. Realizations of different ground-truth correlations (0.25, 0.5, and 0.75) and spike amplitudes (low, medium, and high) were performed. The error in Pearson's correlation, relative to the ground-truth spontaneous signal, was calculated as the absolute difference following despiking. The results ( Fig. 10 ) show overall superior reduction of spurious correlation by Schrödinger filtering, represented by the lowest error score on each occasion. These results suggest that Schrödinger filtering's propensity to reduce spurious spike-related correlation is insensitive to the spike amplitude and the ground-truth correlation, which may be because of the technique's robustness in preserving signal outside spike regions.
These results provide evidence that Schrödinger filtering enables more accurate correlation-based functional connectivity analyses of broadband spontaneous activity that is initially contaminated with spikes. Fig. 8. Top row: high-amplitude spikes; m iddle row: medium-amplitude spikes; bottom row: low-amplitude spikes. Left column: time series of example epochs of the spike-containing signal (black), as well as the ground-truth signal and the output followed by processing by SF, WD, AT, or MF (red). Middle column: pairwise scatter plot of the trial-wise signal-to-artifact ratio of the spike-contaminated epochs from all channels following Schrödinger filtering versus that of the other techniques. With 30 channels and 300 epochs per channel, there are 9,000 trials per plot. Dashed lines of slopes unity and three are shown for reference. Right column: Topographic plots of the alpha-band activity. SF, Schrödinger filtering; WD, wavelet denoising; AT, amplitude thresholding; MF, median filtering.

Conclusions and perspectives
The success of Schrödinger filtering with high-spike artifact datasets like the FMRIB dataset and the corresponding simulated datasets encourages the use of EEG-fMRI in scenarios where EEG-MRI clock synchronization or low subject motion is not possible, and enables a greater volume of retained EEG-fMRI data in general thanks to Schrödinger filtering's spike-cleaning capabilities.
Since despiking by Schrödinger filtering robustly removes spike artifacts while preserving signal, and since spike artifacts overlap in frequency with the EEG bands, Schrödinger filtering relaxes the need to aggressively low-pass filter the data, which is presently done to remove artifacts present in the beta and gamma bands. Thus, Schrödinger filter-ing can potentially broaden the horizon for EEG-fMRI research, enabling studies of higher-frequency activity.
Since Schrödinger filtering is so precise at removing spurious spikes while leaving underlying signal unaffected, it is especially suited for analysis of spontaneous activity, where functional connectivity measures like correlation are sensitive to the functional form of the time series.
The dataset for which the EEG and MRI clocks were synchronized had the fewest post-AAS spike artifacts, further supporting that one of the best acquisition practices is to perform such synchronization. Given that post-AAS spike artifacts are inevitable even with synchronized acquisition, and that they are expected to increase in number and severity since gradient coils are being made faster and stronger, a robust despik- Fig. 9. (A) One electrode channel time series following average artifact subtraction (AAS) (blue) from the decision-making dataset, which was recorded with synchronization of the EEG and MRI clocks. A closeup of a post-AAS spike-containing region is shown, for which Schrödinger filtering is able to despike (orange). (B) The average event-related potential across trials, channels, and subjects for one of the conditions following AAS (blue) and further following Schrödinger filtering, applied to each event-related potential separately (orange). The curves are virtually identical, indicating that Schrödinger filtering did not compromise this activity. (C) Preservation of topographic distribution is also evident. Shown are topographic plots for the event-related potential averaged across trials and five subjects for one condition, per channel, and shown at 0.36 s post-stimulus. Fig. 10. Schrödinger filtering reduces spurious correlation in spontaneous spike-containing EEG data. Pairs of correlated spontaneous EEG timeseries contaminated with spikes were despiked by Schrödinger filtering (SF) or one of the compared techniques (MF, median filtering; AT; amplitude thresholding; WD; wavelet denoising). The error in Pearson's correlation relative to the ground truth was computed for the pre-despiked data as well as the despiked data. Panels are ordered according to the values of the ground-truth correlation and spike amplitude.
ing technique like Schrödinger filtering may be an ideal complement that could potentially be implemented in real time and with hardware feedback.
Lastly, given that the parent technique -semi-classical signal analysis (SCSA) -is young and actively studied, Schrödinger filtering is expected to improve at a rapid pace. Along this vein, the authors see promise in notch filtering applications of Schrödinger filtering, in which the SCSA basis functions, expressed in terms of frequency instead of time, are used to design high-fidelity notch filters. Furthermore, since SCSA is robust at detecting, preserving, or removing spikes, Schrödinger filtering should find various uses in the realm of processing of other spike-containing data, such as motion spike removal in fMRI data or action potential segmentation in electrophysiological data.

Data and code availability statement
All code is available at https://github.com/gbenigno/schrodinger_ filtering.