Effectively combining temporal projection noise suppression methods in magnetoencephalography

BACKGROUND
Magnetoencephalography (MEG) is an excellent non-invasive tool to study the brain. However, measurements often suffer from the contribution of external interference, including noise from the sensors. Suppression of noise from the data is critical for an accurate representation of brain signals. Due to MEG's limited spatial resolution and superior temporal resolution, noise suppression methods that operate in the temporal domain can be favorable.


NEW METHOD
We examined the independent and joint effects of two temporal projection noise suppression algorithms for MEG measurements: One commonly used algorithm which suppresses correlated noise; temporal signal space separation (tSSS) and one new method which suppresses uncorrelated sensor noise; oversampled temporal projection (OTP).


RESULTS
We found that both OTP and tSSS effectively suppress noise in raw MEG data and have the greatest effect of joint operation in cases where SNR is low, or when detecting higher SNR single-trial responses from raw data. We additionally demonstrate how the combination of OTP and tSSS is useful for the detectability of high-frequency brain oscillations (HFO).


COMPARISON WITH EXISTING METHODS
Although the mathematical description of OTP has been described before (Larson & Taulu, 2017), OTP's effect on HFOs in MEG data is novel. Additionally, the combination of OTP and commonly used temporal noise suppression algorithms (i.e., tSSS) has not been shown.


CONCLUSIONS
This finding is applicable to clinical populations such as epilepsy, where HFO signals are thought to be important markers for areas of seizure onset and are typically difficult to detect with non-invasive neuroimaging methods.


Introduction
Magnetoencephalography (MEG) recordings are performed with multi-channel sensor arrays that are optimized to detect magnetic fields produced by neuronal currents in the brain. The measured signal is a combination of brain activity, physiological interference (e.g., heart, muscle activity, eye blinks), environmental interference (e.g., power lines, electronics), and sensor noise (e.g., transducer or electronic noise) (Taulu, 2014). The brain-related magnetic fields are extremely weak (10-100 fT)  at typical measurement distances, and the amplitude of the interfering signals is often significantly larger than that of the brain activity of interest (Taulu, 2014). Therefore, noise suppression is critical for accurate data analysis and interpretation of results.
In the spatial domain, the number of degrees of freedom or effective rank of the MEG data has been described in the past (Ahonen et al., 1993;. According to the spatial sampling theory, neuromagnetic signals detectable by MEG are limited to the low end of the spatial frequency spectrum because the amplitude in higher spatial frequencies decays faster as a function of distance from the source and there is always a considerable gap between the brain and the MEG sensors (Taulu, 2014). Although the number of sensors and spatial resolution is limited, the temporal sampling rate for MEG is high. For this reason, the application of noise suppression methods that operate in the temporal domain have the potential to significantly enhance the MEG signal, especially for high frequency signals. Several methods exist for suppressing physiological artifacts as well as external environmental interference. Spatial and temporal filtering are robust suppression methods under a large variety of interference conditions. These methods rely on a priori assumptions about the temporal or spatial structure of the noise. Several methods exist to suppress interference that is spatially correlated between sensors, including signal space separation (SSS), signal space projection (SSP), independent component analysis (ICA), and principal component analysis (PCA). In particular, the SSS method spatially distinguishes signals produced by magnetic sources inside the MEG sensor array from contributions of sources located outside by applying quasistatic Maxwell's equations to the conditions of the MEG recording, making SSS applicable to any multichannel MEG system. In the spatial domain, brain and interference signal vectors (N-dimensional vectors, composed of signals from all N channels for a given time point) are linearly independent, which enables their unique separation with an oblique projection , even though the brain and interference signals are not orthogonal in the spatial domain. The tSSS method (Taulu and Simola, 2006), which is the temporal extension of SSS, uses temporal information to suppress signals from sources that are too close to the sensors for the associated signals to be clearly classified as internal or external. Because the signals produced by nearby artifact sources are statistically independent of brain signals in the time domain, temporal orthogonal projectors can be used to suppress their contribution in the data.
In contrast to external interference, sensor noise fluctuations are random processes without specific temporal or spatial patterns . Without an accurate representation of the spatiotemporal structure of the sensor noise, spatial filtering can spread the noise between channels, and temporal filtering can result in residual noise in the data. Although environmental interference signals are typically strong and tend to dominate MEG data, random sensor noise is typically weak and therefore not always taken into account in the data processing pipeline (Volegov et al., 2004). Sensor noise becomes problematic when the signal of interest is weak and occupies a frequency band dominated by uncorrelated sensor noise (Larson and Taulu, 2017). A prime example of brain responses that suffer from these conditions is high-frequency oscillations (HFO). These high-frequency brain signals are potential biomarkers for several brain functions and disorders (Ozaki and Hashimoto, 2011) and typically have a much smaller peak-to-peak amplitude than responses at lower frequencies (Zijlmans et al., 2017). For example, HFOs are thought to be a reliable indicator for the epileptic seizure onset zone in patients with epilepsy (Graef et al., 2013). Localization of these signals could represent epileptogenic regions better than localization of interictal spikes and optimize the diagnosis and treatment of epilepsy. However, these events can be challenging to detect due to low SNR compared to other epileptiform discharges (Zijlmans et al., 2017). Environmental and physiological noise signals are typically concentrated to lower frequencies whereas random sensor noise typically covers the entire frequency spectrum. Therefore, suppression of sensor noise would significantly increase SNR for high-frequency oscillations, generated for example, by inter-ictal spikes.
Oversampled temporal projection (OTP) (Larson and Taulu, 2017) is a cross-validation (CV) method designed to suppress sensor noise in multi-channel systems. OTP does not require visual inspection of channels by the user, or a priori assumptions about the temporal or spatial structures of the noise. It requires only that the multichannel measurement provides oversampling of the signal of interest, as is the case in multichannel MEG measurements. Unlike other CV methods for sensor noise suppression, such as the "sensor noise suppression" (SNS) algorithm (de Cheveigné and Simon, 2008), which uses the spatial correlation between channels to build a spatial denoising operator, the CV-based operator in OTP is derived and applied in a time-varying manner in the temporal domain. The current implementation of OTP uses singular value decomposition (SVD) to extract temporal patterns in the data to create an equivalent representation as a spatial domain operator. By calculating and applying operators in a time-varying manner across overlapping temporal windows, each operator is optimized for the given window. This results in a significant reduction of induced signal localization bias as compared to using a single spatial operator method, which might reduce the amplitude of the signals of interest (see Larson and Taulu, 2017). The effectiveness of using temporal projection operators was investigated either in isolation or in combination to suppress spatially correlated noise sources (tSSS) and spatially uncorrelated noise sources (OTP). Although the performance of tSSS alone has been previously described (Taulu and Simola, 2006;Taulu, 2008;Medvedovsky et al., 2007), it is unclear how tSSS is affected by the level of random noise in the data. The goal of this paper is to evaluate the effectiveness of tSSS and OTP for noise suppression in MEG by quantifying the effects on data quality, SNR, and source localization bias. The effects were compared using recordings of empty room data, a calibrated phantom, and a dataset containing high-frequency oscillations from the early somatosensory evoked response. The combination of OTP and tSSS is shown to be most effective in reducing noise in raw and processed MEG data and that combining these methods has the greatest effect of joint operation on data when the signal to noise ratio is low (e.g., few trials or weak signals), or on the detection of higher SNR signals in raw data (e.g. epileptic spikes). The results highlight the importance of readjusting tSSS parameters to operate properly on data that have already been denoised using OTP.

Materials and methods
To quantify the effects of OTP and tSSS on multichannel MEG data, the following datasets were used: (a) empty room data to compare spectra, (b) recordings from a calibrated phantom head to compare source localization and SNR, and (c) human median nerve stimulation data to compare quality of high frequency evoked responses. Each dataset was compared alone and in combination with OTP: (a) unprocessed (raw) data, (b) data processed with SSS, and (c) data processed with tSSS.

Oversampled temporal projection (OTP)
Oversampled temporal projection (Larson and Taulu, 2017) is a sensor noise suppression algorithm that does not require spatial or temporal assumptions about the structure of the noise and requires minimal intervention from the user. It only requires that the noise is spatially uncorrelated. OTP works on any multi-channel measurement system with N independent measurement channels that provides oversampling of the signal of interest. In a multi-channel system with independent channels, each channel is expected to have some independent sensor noise as well as a combination of the signal of interest and other sources of external noise. OTP works as a leave-one-out cross validation method, when one channel at a time is left out and its signal is reconstructed from all remaining channels. Overlapping windowed segments of data are used; within each window, singular value decomposition (SVD) is applied to create a basis of temporal components formed by all other channels. The window length is used to determine a trade-off between noise suppression and adaptation to sensor characteristics which change over time. This window length can be changed based on the data of interest. For example, longer windows have a larger noise suppression factor, but may not be suitable for data where the noise characteristics can change over time (i.e., not perfectly stationary). The original channel signal is then projected onto the temporal basis spanned by the right singular vectors of the SVD. In this projection, the brain signal of interest will be mostly explained by the aforementioned temporal basis, whereas the random noise of the channel under investigation should correspond to a vector that is orthogonal to the temporal basis. Thus, its projection approaches zero with long data segments. For further details and mathematical description, see Larson & Taulu, (2018).

Signal space separation (SSS)
Signal space separation (SSS) ) is a method that compensates for external interference artifacts and is commonly used in MEG. SSS is derived from quasistatic Maxwell's equations, which enables one to apply Laplace equation to uniquely separate magnetic fields arising from inside and outside of the sensor array. In practice, SSS separates multi-channel magnetic signals into two subspaces: signals generated in the brain (internal subspace) and signals generated by sources outside the sensor array (external subspace). Thus, SSS is able to suppress external interference signals while preserving the signals of interest.

Temporal signal space separation (tSSS)
Temporal signal space separation (tSSS) (Taulu and Simola, 2006;Taulu and Hari, 2009), the temporal extension of SSS, additionally suppresses the contribution of nearby artifact sources (e.g., metal implants, dental work) by utilizing temporal information in addition to Maxwell's equations. Nearby signals that exist in both the internal and external subspaces are recognized by comparing the two in the temporal domain and suppressing the part that is common for both. This method assumes that the brain signals and any artifact signals are uncorrelated in time.

tSSS correlation limit (CL)
The tSSS subspace correlation limit (CL) is a value between 0 and 1 that is used to quantify the degree of synchrony in both the internal and external subspaces above which magnetic signals will be projected out. In the case when CL = 0, no synchrony is required at all (i.e., all components will be projected out), and in the case when CL = 1, absolute synchrony is required (i.e., in practice, no components will be projected out). Reducing the CL will promote more liberal acceptance of the signal components as behaving synchronously and thus more components will be projected out. In other words, the rejection of artifacts can become more efficient. In principle, this does not pose a significant risk of accidental removal of brain signals as their contribution should be minimal in the residual signal that results when the internal SSS-reconstructed signal has been subtracted from the data. However, the residual might contain a small fraction of the brain signal due to truncation or calibration errors, which might be detected by tSSS as an artifact, in terms of a given CL, when the amplitude of random noise in the data decreases. Therefore, CL should be reduced when there is decreased random noise in the data, and increased when the data have less noise. The effect of CL has been studied in detail previously (Medvedovsky et al., 2009) on data with typical sensor noise levels, and they generally recommend CL values in the range of 0.8 to 0.98.

Fig. 2.
OTP reduces noise floor in data from an empty room recording. The power spectral density across each channel type calculated from empty room data before and after OTP processing. The mean +/-1 STD (across channels) is shown. Dashed grey lines indicate a notch filter (60 Hz) and low pass filter (330 Hz). The top row shows raw data, the middle row shows data processed with SSS and bottom row shows data processed with tSSS. A tSSS CL of 0.999 was used for data processed with OTP and a CL of 0.98 was used for unprocessed data. All data were recorded on using a whole-head Elekta Neuromag Vectorview system (Elekta, Oy, Helsinki, Finland) with 204 axial gradiometers and 102 magnetometers located inside a magnetically shielded room (IMEDCO, Noblesville, IN) at the University of Washington's Institute for Learning & Brain Sciences. Recruitment and experimental procedures for resting state data were approved by the university Institutional Review Board (IRB) and written informed consent was acquired from the participant. The participant received financial reimbursement for their time. Empty room data were collected to compare frequency spectra. To compare the SNR and source localization accuracy, a calibrated phantom (Elekta Oy, Helsinki, Finland), designed to mimic the neural activity of the human brain, was used. The phantom contained 32 current dipoles at known locations, oriented tangentially relative to the phantom's origin. Each of the 32 dipoles was fed with a current magnitude of 500 nA m peak-to-peak and was activated for two cycles of a 20 Hz sinusoid for 100 repetitions. A 500 nA m phantom was chosen in order to mimic larger, high-SNR single trial responses such as epileptic spikes. Five fixed head-position-indicator (HPI) coils were attached to the surface of the phantom and digitized to register the location of the phantom within the MEG helmet. Signals were recorded with analog band-pass filtering (0.03-650 Hz) and sampled at 2 kHz.

Human median nerve data
Data were collected on one adult participant with the Elekta Neuromag Triux MEG system at the BioMag Laboratory, Helsinki University Central Hospital, Finland, with 204 axial gradiometers and 102 magnetometers. Written informed consent was acquired from the participant and they received financial reimbursement for their time. The study is in line with the Declaration of Helsinki and was approved by the ethics committee of the Helsinki University Central Hospital. Data were filtered in hardware with a combination of high pass (0.1 Hz), notch (50 Hz) and lowpass (1 kHz) filters and acquisition at a rate of 3 kHz. The participant's right median nerve was stimulated with a Lucius & Baer Konstant-Strom Stimulator. A total of 2852 squarewave electric pulses (200 μs duration) were delivered at 5.5 mA with an ISI of 300-350 ms. A simultaneous trigger from the stimulator was sent to the MEG acquisition system to mark the onset of each pulse for signal averaging.

Phantom data
The phantom recording was processed with OTP, SSS, tSSS or no environmental noise suppression. OTP window lengths of 1, 10 and 60 s were used. For both SSS and tSSS, an inside expansion order of 8 was used, and an outside expansion order of 3. For tSSS, a window length of 20 s was used. Because the phantom has a higher SNR than realistic brain data, and because OTP reduces the amount of independent sensor noise in the raw data thereby increasing the SNR, higher tSSS CLs were required. A CL of 0.9999 was used for OTP processed data and a CL of 0.995 was used on data that was not processed with OTP. The choice of these two correlation limits led to roughly the same number of noise components being projected out from the data when processing both raw and OTP data with tSSS. All datasets were epoched around the onset of the phantom signal (−100 to 250 ms). All 100 trials for each current source were averaged together to generate an inter-trial average. MEG sources were estimated using a single equivalent current dipole (ECD) model. Dipoles were localized to the first largest peak in the sinusoid (latency 40 ms after stimulus onset). A sphere was used as a conductor model where the origin was (0, 0, 0) in the head coordinate system, and the noise covariance was calculated from the epoch baseline (−100 to −25 ms). The localization of each dipole was compared with the nominal data provided by the MEG manufacturer. The  displacement (mm) of each dipole from its nominal location was calculated and the mean was taken as a measure of the localization error. This value was compared across conditions. In order to observe the effects of tSSS and OTP in realistic noise conditions, the phantom data were combined with a brain noise resting state measurement, when the participant was asked to rest with his eyes closed.

Human median nerve & resting state data
For both resting state and median nerve data, OTP was applied using the default window length of 10 s. SSS and tSSS were applied separately to compare their individual performances. For both SSS and tSSS, an inside expansion order of 8 was used, and an outside expansion order of 3. For tSSS, a window length of 20 s was used. A tSSS CL of 0.999 was used for OTP processed data. In the case of raw data that were not processed with OTP, the default CL of 0.98 was used. The choice of these two correlation limits led to roughly the same number of noise components being projected out from the data when processing both raw and OTP data with tSSS.
In order to see the high-frequency somatosensory response, a high pass filter of 500 Hz was applied to the median nerve data. All event related responses were extracted as epochs with a time window of −50 to 200 ms averaged with respect to the median nerve stimulus trigger. The trials were baseline-corrected by applying a baseline of −50 to 0 ms before the trigger onset.

Median nerve data: single-trial responses
In order to see the effects of OTP and tSSS on single-trial responses, brain responses were compared from 10 s of raw median nerve data before and after processing with OTP and tSSS. Visual inspection showed that single-trial responses become more easily identified after processing with OTP (Fig. 1). Fig. 1. Also shows the joint combination of OTP and tSSS and the effect of altering the tSSS correlation limit. When tSSS is used on OTP processed data, a higher CL is required. With a lower CL (e.g., 0.98), responses can become flattened due to the removal of signal components from the data. When using a higher CL (e.g., .999), tSSS removes less components and the responses become more visible.

Empty room data: spectra
To compare overall noise reduction, empty room data spectra across conditions for each sensor type were examined. Empty room spectra showed a ∼10 dB reduction in the noise floor in all conditions after processing with OTP (see Fig. 2). In addition, applying tSSS reduced the noise floor in gradiometers compared to raw data.

Phantom data: localization and SNR
Simulated phantom data were compared to observe whether OTP and method type (MT) have an effect on source localization and SNR.

Fig. 5.
Performance of tSSS is affected by the noise level. Averaged phantom signal for a superficial dipole (n = 100), processed with tSSS is shown before and after combining 'brain noise'. Top row is phantom data alone, and bottom row is phantom data combined with brain noise. The right column is data processed with OTP. It is observed that the application of tSSS on phantom data with no 'brain noise' leads to worse results than one would expect with realistic brain signals.
M. Clarke, et al. Journal of Neuroscience Methods 341 (2020) 108700 Data from a single trial (N = 1) and data from multiple trials (N = 100) from each of the 32 dipoles were analyzed using two-way repeated measures ANOVAs manipulating OTP (with and without OTP processing) and MT (raw, SSS, and tSSS). For single trial data (Fig. 3a), there was a significant effect of OTP, F(1, 31) = 32.54, p = < 0.001, with significantly lower localization error when OTP was applied. There was also a significant effect of MT, F(2, 62) = 12.11, p = < .001, however, there was no significant interaction between MT and OTP, F(2, 62) = 0.87, p = 0.42. For N = 100 trials (Fig. 3b), there was no significant effect of OTP, F(1, 31) = 1.88, p = 0.17, MT, F(2, 62), p = 0.12 0 , or interaction between them, F(2, 62) = 0.38, p = 0.68. To examine the significant effect of MT on localization error, post- Fig. 6. Detectability of HFOs is improved after tSSS and OTP. High frequency somatosensory evoked field shown with raw data processed with tSSS (top) and raw data processed with OTP + tSSS (bottom). Data were high-pass filtered at 500 Hz, and the time scale is 5-50 ms. SEF topography becomes clear after the elimination of artifacts after OTP processing.
M. Clarke, et al. Journal of Neuroscience Methods 341 (2020) 108700 hoc paired t tests on MT were performed: 1. For N = 1 trial: raw versus SSS processed data, t(31) = −3.40, p = .002, and raw versus tSSS processed data, t(31) = −3.32 p = .002. A slight increase in localization error after SSS/tSSS has been previously discussed by  however, this effect is seen for N = 1 trial and is not significant for N = 100 trials.

Phantom data: Performance of tSSS is affected by noise level
To analyze whether the performance of tSSS can be affected by the noise level, we compared phantom data alone and in combination with human resting state data, acting as surrogate brain noise. The phantom data were processed with both OTP and tSSS with a number of CLs, and then responses (N = 100) from a superficial dipole were averaged together. The performance of tSSS is affected by the level of noise in the data (see Fig. 5), whereby applying tSSS on data with high SNR (e.g., phantom) requires a higher CL than data with additional noise (here, more realistic brain signals). Fig. 7. HFO response is more easily detectable with a fewer number of trials after OTP in combination with tSSS. The top row shows the averaged waveform of a gradiometer with the highest peak amplitude in somatosensory cortex MEG0442 without (left) and with (right) OTP processing. The middle row shows amplitude over time as a function of the increasing number of trials in the averaged response. The bottom row shows the amplitude values as z-scores (relative to baseline).

High-frequency somatosensory data: improved quality
The high frequency somatosensory evoked field (SEF) topology becomes more easily visible after processing with OTP and tSSS from sensors over the somatosensory cortex (see Fig. 6). The gradiometer with the highest peak amplitude in the somatosensory cortex (MEG0442) was used to compare amplitude as a function of increasing number of trials in the averaged response (Fig. 7) and time-frequency representation (Fig. 8). The HFO response starting at ∼ 20 ms is more clearly visible after noise suppression with OTP. Noise has been reduced across the time course after OTP processing, especially in averaged responses containing fewer trials (Fig. 7). The analysis of inter-trial coherence (ITC) values on sensor MEG0442 revealed increases across the 100−1000 Hz spectrum during the time of the HFO response (∼ 20 ms) after OTP processing (Fig. 8), appearing roughly twice as many standard deviations above the mean noise level (based on the z-score relative to baseline levels).

Discussion
MEG data can be separated into the brain signal of interest, external interference sources, and intrinsic noise and artifacts arising from the MEG sensors. Suppression of noise is critical for successful analysis of the measured brain signals. Due to MEG's limited spatial resolution, and high temporal resolution, the application of noise suppression methods that operate in the temporal domain can be favorable.
We investigated the validity of two temporal projection noise suppression algorithms for MEG measurements: one which suppresses interference that is spatially correlated between sensors (tSSS), and one which suppresses interference that is independent across sensors (OTP). Previous work (Larson and Taulu, 2017), describes the potential application of OTP for noise suppression in MEG, but does not show the effect on brain data and when combined with existing noise suppression methods such as SSS and tSSS. Taulu and Hari, 2009 have also demonstrated tSSS as an effective external noise suppression method for MEG data, but do not describe the effect of noise level on performance.
We have shown that both OTP and tSSS effectively suppress noise in raw MEG data and have the greatest effect of joint operation in cases when SNR is low, such as when there is a low trial count or a low-amplitude signal of interest. The application of OTP to raw data decreases localization error in single trials and increases SNR in both single and multiple trials. Previous research shows the ability of tSSS to increase the detectability of single-trial responses (Taulu and Hari, 2009). However, we have shown that the combination of tSSS and OTP further increases the SNR of single-trial and multiple-trial responses. We have also demonstrated how tSSS performance is affected by the level of random noise in the data: tSSS becomes less sensitive with the addition of noise, and the correlation limit should be adjusted accordingly in order for optimal performance. In particular, when using a higher SNR signal (such as the phantom), a higher correlation limit should be used, especially when combined with OTP.
In addition, we demonstrate how the combination of OTP and tSSS is useful for the detectability of HFOs and higher SNR signals from raw data. This finding is important, and notably applicable to clinical populations such as epilepsy, when HFO signals are thought to be important markers for areas of seizure onset and are difficult to detect with non-invasive techniques, such as MEG, due to low signal-to-noise ratio (SNR) levels. Suppression of sensor noise would be especially beneficial to increase SNR for these types of high-frequency oscillations and could help improve the diagnosis and treatment of epilepsy. The OTP method is simple to use and can provide vast improvements to a range of MEG datasets alone and in combination with tSSS.
We note that there are three important caveats with using OTP and tSSS. First, OTP can introduce a small amount of amplitude bias, therefore we would not recommend using OTP on datasets when the SNR is already very high (e.g., evoked datasets with hundreds/thousands of trials); OTP should be used primarily on datasets with low SNR, or on single trial high SNR responses (as we have demonstrated with the phantom). The reason for the potential amplitude bias is as follows. For each channel, OTP computes the temporal subspace for the signals of interest via the singular value decomposition of all other channels. The idea is that, due to the oversampling provided by MEG arrays, the signal of any channel can be reconstructed based on the information of other channels. In OTP, the reconstruction is done by first determining the (temporal domain) right singular vectors of the aforementioned SVD operation and then projecting the original time domain signal of the channel under investigation onto this subspace. For relatively long data segments (at least several seconds), the random M. Clarke, et al. Journal of Neuroscience Methods 341 (2020) 108700 noise vector of the channel should have a large angle compared to the temporal subspace and therefore its projection will be short, corresponding to a reduction of random sensor noise. However, if the channel under investigation has a very high signal-to-noise ratio with respect to the signal of interest, it could be that not all of the details of the signal will be fully captured by the right singular vectors of the SVD. In such a case, significant amplitude distortion could occur. Second, the algorithm tends to be computationally expensive. Multiprocessing computers are expected to speed up the calculations significantly. Further optimization of the OTP algorithm is thereby an important future goal. Third, when using OTP in combination with tSSS, a higher tSSS correlation limit should be used. Considering OTP reduces the noise in raw data, tSSS becomes more sensitive and requires a higher correlation limit when used on cleaner data. This approach is also true for any instances of applying tSSS on data with very high SNR. As we have shown in tests with the phantom signal, tSSS requires a higher correlation limit to yield optimal results.

Declaration of Competing Interest
The authors report that Dr. Samu Taulu and Dr. Eric Larson are consultants for S3 Signal Processing, a company that has a commercial interest in the methodology presented in this manuscript.