Source reconstruction of broadband EEG/MEG data using the frequency-adaptive broadband (FAB) beamformer

A beamformer enhances the signal from a voxel of interest by minimising interference from all other locations represented in the sensor covariance matrix. However, the presence of narrowband oscillations in EEG/MEG implies that the spatial structure of the covariance matrix, and hence also the optimal beamformer, depends on the frequency. The frequency-adaptive broadband (FAB) beamformer introduced here exploits this fact in the Fourier domain by partitioning the covariance matrix into cross-spectra corresponding to different frequencies. For each frequency bin, an individual spatial filter is constructed. This assures optimal noise suppression across the frequency spectrum. After applying the spatial filters in the frequency domain, the broadband source signal is recovered using the inverse Fourier transform. MEG simulations using artificial data and real resting-state measurements were used to compare the FAB beamformer to the LCMV beamformer and MNE. The FAB beamformer significantly outperforms both methods in terms of the quality of the reconstructed time series. To our knowledge, the FAB beamformer is the first beamforming approach tailored for the analysis of broadband neuroimaging data. Due to its frequency-adaptive noise suppression, the reconstructed source time series is suited for further time-frequency or connectivity analysis in source space.

Electroencephalography (EEG) and magnetoencephalography (MEG), jointly 2 abbreviated as MEEG here, measure the electrical currents or magnetic fields associated 3 with synchronised activity of neuronal populations [1][2][3]. Compared to hemodynamic 4 methods such as functional magnetic resonance imaging (fMRI), the Achilles' heel of 5 MEEG is the low accuracy in localising the sources. Since the sensors are located at a 6 distance to the brain, the activity of brain sources needs to be statistically 7 reconstructed. The inverse problem of estimating the sources from the sensor 8 measurements is underconstrained because the number of sources exceeds the number of 9 sensors. Furthermore, neurons in some subcortical areas are aligned such that they 10 constitute a closed field, producing little measurable electrical activity [1]. Consequently, 11 in order to obtain a unique solution, additional assumptions about the sources (such as 12 their uncorrelatedness) are required to constrain the inverse model [4,5]. The linear 13 forward model of MEEG, specifying how activity in the brain projects into the sensors, 14 can be formulated as x(t) = L s(t) + (1) where x(t) ∈ R N is activity measured in N sensors at time t, s(t) ∈ R M is the 16 source activity in M brain sources at time t, L ∈ R M ×N is the leadfield matrix relating 17 source activity to the EEG/MEG sensors, and ∈ R N is a noise term representing 18 sensor noise and non-neural artefacts such as muscle activity. The i-th row of the 19 leadfield matrix is called the spatial pattern of the i-th source. The leadfield matrix L is 20 typically derived from forward calculations using a subsampled MRI or a cortical 21 surface model [1]. Here, it is assumed that the dipole orientations are fixed, but the 22 model is very similar for free dipole orientations (see below). 23 Source localisation versus source reconstruction 24 Often, the primary purpose of source imaging is pinpointing the neural location that 25 gives rise to observed EEG/MEG activity, i.e. source localisation. Dipole fitting models 26 such as Multiple Signal Classification (MUSIC) and their extensions [6][7][8] search for a 27 small number of sources that explain the data with sufficient accuracy. 28 Some imaging approaches such as minimum-norm estimates [9][10][11] and beamformers 29 also provide linear functionals called spatial filters which allow for the estimation of the 30 time series at a given voxel, i.e. source reconstruction. It has been shown that these 31 methods can be cast into the generalised model 32 s(t) = C s L C −1 x(t) =: w x(t) (2) where C ∈ R N ×N is the sensor covariance matrix, C s ∈ R M ×M is the source 33 covariance matrix,ŝ is the estimated source amplitude (a.k.a. virtual channel), and 34 w ∈ R N is the spatial filter [12]. 35 Time-domain analysis: the LCMV beamformer 36 The purpose of beamforming is to optimally enhance the neuronal signal from a specific 37 location in the brain with a known spatial pattern [3][4][5]13]. Several beamforming 38 approaches have been shown to be equivalent to the Linearly Constrained Minimum 39 Variance (LCMV) beamformer [14,15]. They differ mainly in the way that neural 40 activity is quantified [13]. The spatial filter is given as where l ∈ R N is the row of the leadfield matrix corresponding to the target source.

42
The principal assumption for beamforming is that all sources are uncorrelated and 43 hence the source covariance matrix C s is diagonal. Indeed, inserting a diagonal source 44 covariance matrix into Eq 2 leads to a scaled version of the beamformer solution in Eq 3. 45 Recently, it has been shown that LCMV beamforming is formally equivalent to Linear 46 Discriminant Analysis (LDA) [16].

47
For source localisation, source amplitude or variance can be calculated for all source 48 locations. Peaks in the source power map indicate an active source. However, the 49 relatively small norm of leadfields for deep sources causes a depth bias by amplifying 50 noise. To account for this, the source power estimate is usually divided by noise power. 51 The resulting quantity was originally introduced as neural activity index by Van Veen et 52 al. [15] and is given here for fixed dipoles [13,17] as For free dipoles, the spatial pattern l = [l x l y l z ] consists of 3 components along the 54 x, y, and z axes. In the original approach by Veen et al., unreliable directions with a 55 small norm can dominate the resulting quantity. The modified index by Huang et 56 al. [13] remedies the problem by normalising each direction separately: The noise covariance matrix Q can be estimated e.g. by using pre-stimulus activity 58 in an event-related design or empty room recordings in the case of MEG. Coherent Sources (DICS) approach [18,19]. The cross-spectral density matrix C(ω) is 62 determined for a particular frequency ω. The real part of C(ω) then takes the role of 63 the covariance matrix, leading to a frequency-dependent solution analogous to LCMV 64 beamforming: where (·) designates the real part. Beamforming in the time-frequency domain was 66 introduced by Dalal et al. [20]. Sensor time series are passed through filter banks 67 enabling frequency-specific beamformers. Additionally, cross-spectra are calculated for 68 specific time segments, accounting for nonstationarity. As a result, beamformer weights 69 are obtained for each time-frequency bin where the spatial filter and the cross-spectral density are now functions of both time 71 t and frequency ω.

73
What is the optimal spatial filter for reconstructing a broadband signal? As illustrated 74 in Fig 1, the spatial structure of the background activity varies as a function of 75 frequency. Clearly, a spatial filter that is optimal for a given frequency band is most 76 likely not optimal for a different frequency band.

77
The LCMV beamformer, when applied to broadband data, is biased towards 78 cancelling noise in low frequencies because they contribute most of the variance [20], 79 producing a flattened frequency spectrum of the estimated source [21]. DICS, on the 80 other hand, has an excellent signal-to-noise ratio because beamforming is restricted to a 81 narrow frequency band. The disadvantage is that DICS cannot be directly applied to 82 broadband signals. The time-frequency beamformer solves this problem by splitting the 83 signal into time-frequency bins [20]. However, it is computationally expensive and the 84 different frequencies are analysed separately. Power spectra. Power spectra for different electrodes (coloured lines) for an example EEG dataset, depicted on a logarithmic frequency axis. For the frequencies 4, 8, and 80 Hz, the real parts of the cross-spectra are shown as small rectangles. To better illustrate the spatial structure of the cross-spectra, each topography depicts the eigenvector of the covariance matrix corresponding to the largest eigenvalue with magnitude being colour-coded. Clearly, the spatial structure of the cross-spectrum changes as a function of frequency.
Here, a different approach is taken by considering a broadband signal in the 86 frequency domain. Due to the linearity of the Fourier transform, applying a spatial 87 filter in the time domain is equivalent to applying it in the Fourier domain. However, 88 instead of using the same spatial filter for each frequency, one can divide the spectrum 89 up into frequency bins and calculate the cross-spectrum corresponding to each bin. This 90 allows for the construction of a different spatial filter for each frequency bin. After 91 applying these spatial filters, the Fourier transform can then be inverted to recover the 92 broadband signal in source space. We call this approach frequency-adaptive broadband oscillations in MEEG assures that this partitioning is non-trivial, that is, the 104 cross-spectra have different spatial structures (see Fig 1). Analogous to minimising the 105 covariance matrix by a single spatial filter (as in LCMV beamforming), one can instead 106 perform beamforming by minimising the cross-spectra.

107
Let X ∈ R N ×T be the sensors × time data matrix, and let x i ∈ R T be the i-th row 108 of X corresponding to the i-th MEEG channel. Then its Fourier transform is defined as 109 f i (ω) = T −1 t=0 x i (t) exp(−i ω t) and x i (t) can be expressed using the inverse Fourier where the frequency index k is related 111 to the angular frequency as ω k = 2πk/T for k = 0, 1, . . . , T − 1.

112
Assuming that all channels have zero mean, the covariance between the i-th and the 113 j-th channel, corresponding to the (i,j)-th entry of the covariance matrix C, and their 114 cross-spectral density, corresponding to the (i,j)-th entry of the cross-spectral density 115 matrix C(ω) for frequency ω can be defined as where · denotes the complex dot product. Using the orthogonality of the Discrete not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available The copyright holder for this preprint (which was this version posted December 20, 2018. ; https://doi.org/10.1101/502690 doi: bioRxiv preprint Together, Eq 8 and Eq 9 show that the cross-terms vanish and hence the broadband 121 covariance matrix can be partitioned into cross-spectral densities: Now, the conjugate symmetry of the Fourier spectrum of real signals can be 123 exploited, implying that C(ω k ) = C(ω T −k ) * . Hence, the covariance matrix can be 124 expressed by the real parts of the cross-spectra as

126
The periodogram is an unbiased but inconsistent estimator of the power spectral density. 127 In Welch's method, the signal is divided up into overlapping segments. Each segment is 128 multiplied by a window function and then DFTs are calculated. The cross-spectra 129 derived from the individual segments are then averaged. This significantly improves the 130 estimate at the expense of frequency resolution. However, a resolution of e.g. less than 131 0.5 Hz is rarely needed for neural oscillations in MEEG, justifying this trade-off. An 132 alternative way of spectrum estimation is the use of multitapers. Here, a Slepian 133 sequence of mutually orthogonal tapers is used to obtain independent power 134 estimates [22].

135
Regularisation 136 EEG and MEG signals are low-dimensional, often leading to ill-conditioned covariance 137 and cross-spectral density matrices. Moreover, if a subspace has been removed in the 138 course of artefact correction (e.g. using maxfiltering or Independent Components 139 Analysis), the cross-spectral matrices are singular. To guarantee invertibility, they can 140 be slightly manipulated using a regularisation approach. In analogy to shrinkage 141 estimation in covariance matrices [23][24][25], the shrinkage regularisation approach is 142 adapted here to the cross-spectra: Here, I ∈ R N ×N is the identity matrix and ν ω = trace(C(ω))/N scales the identity 144 matrix such that its trace matches the trace of C(ω). λ is the regularisation parameter 145 that blends between the empirical cross-spectrum (λ = 0) and a diagonal matrix with 146 equal total variance (λ = 1). Typically, a small value of λ in the range of 10 −3 is 147 sufficient to ensure invertibility without distorting the empirical cross-spectrum too 148 much. Note that regularisation does not need to come at the expense of a bad estimate. 149 In some cases, the bias introduced by regularisation can even increase the accuracy of 150 the inverse model [24,26].  Simulation data. Graphical illustration of how the data for the simulations was created. In simulation 1, real source spaces were used together with artificially generated signals consisting of a narrowband or 1/f signal (yellow power spectra) and narrowband noise sources (red). The signals were mapped to sensor space using forward models. In simulation 2, real resting-state measurements were used as background noise. In order to have a ground truth, additionally a source time series from another subject was used as the target signal. It was mapped to sensor space using the first subject's source space and then added to the resting-state data.
Frequency-adaptive broadband (FAB) beamforming 152 Using the regularised cross-spectra, the frequency-adaptive spatial filter for frequency ω 153 is given by Let F ∈ R N ×T be the matrix containing the [1 × T ] Fast Fourier Transforms (FFTs) 155 of the sensor time series as rows, and let F(ω) ∈ R N be the column of F corresponding 156 to the frequency ω. The source-space Fourier coefficient for frequency ω is then given by 157 Note that in practice the cross-spectra are calculated across multiple segments, 158 leading to a lower frequency resolution but better estimate. For the purpose of 159 obtaining the Fourier coefficients of the source, the spatial filter w(ω) can therefore be 160 applied in the frequency neighbourhood [ω − /2, ω + /2], where is the width of the 161 frequency bins. Once Fourier coefficients for all frequencies have been calculated, the 162 inverse Fourier transform is calculated to recover the source time series.

163
In order to use the FAB beamformer for source localisation, information needs to be 164 accumulated across the different frequency bins. In Eq 4, the neural activity index was 165 introduced as a way to represent signal power normalised by noise power. As a 166 straightforward extension of this idea, the summed neural activity index can be 167 introduced, which essentially sums the neural activity scores across the frequency bins 168 to measure the collective evidence for the presence of a signal. Here, the formula is 169 given for fixed dipoles, but the extension to free dipoles is straightforward: It is assumed that the signal is analysed in K different frequency bins, and Q(ω k ) is 171 the noise cross-spectrum at frequency ω k . Due to the 1/f shape of the MEEG power 172 spectrum, it can be expedient to use a logarithmic spacing of the frequency bins. This is 173 because with a linear spacing, higher frequency bands such as the gamma band would 174 dominate the outcome by contributing more frequency bins than the lower frequency To evaluate the efficacy of the FAB beamformer in reconstructing sources, simulations 178 were performed wherein it was compared to LCMV beamforming and minimum-norm estimates (MNE) [9,10,27,28]. Simulations were based on artificial signals in source 180 space (simulation 1) or real MEG data from a resting-state measurement (simulation 2). 181 To build source models, individual MRI scans of 582 subjects that took part in the 182 Cambridge Center for Ageing and Neuroscience (Cam-CAN) study on healthy 183 ageing [29][30][31] were used. See [29] for a detailed protocol of the study and [30] for the 184 processing pipeline. Helsinki, Finland). It uses spherical basis functions to reject environmental noise based 207 on magnetic fields outside a sphere enclosing the sensor array, by removing temporally 208 correlated sources inside this sphere but outside a sphere enclosing the brain using the 209 temporal extension of signal space separation (tSSS) [32], with a correlation threshold of 210 0.98 and a 10 s sliding window. The tSSS was further used to perform movement 211 correction every 200 ms using the HPI coils data, and to align all participants' data in a 212 common space (as if the centre of their heads were in the same position relative to the 213 sensors). Finally, MaxFilter was also used to detect and reconstruct bad channels, and 214 to notch-filter the line noise at 50 Hz and its harmonics. ICA was then used to identify 215 physiological artefacts from blinks, eye-movements, and cardiac pulse. This was done by 216 identifying those components whose time courses and spatial topographies correlated 217 highly with reference time courses (correlation greater than three standard deviations 218 from mean) and spatial topographies (correlation greater than two standard deviations 219 from mean), respectively, for each of the above artefact types. (162 per hemisphere) [35]. Surfaces were coregistered with the MEG sensors using 226 FieldTrip [36]. An initial alignment was performed using the fiducial markers (nasion, 227 left and right preauricular points). The coregistration was then refined using the not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available The copyright holder for this preprint (which was this version posted December 20, 2018. ; https://doi.org/10.1101/502690 doi: bioRxiv preprint headshape points. Forward modelling for 204 planar gradiometers was performed using 229 a realistic single-shell volume conductor [37]. This yielded a set of spatial patterns in x, 230 y, and z directions for each vertex collected in a N × 3 matrix. To obtain sources with 231 fixed dipoles, a singular value decomposition (SVD) was performed on the matrix for 232 each vertex. The principal component corresponding to the maximum eigenvalue was 233 then chosen as fixed dipole. There was a total of 582 subjects, aged 18-88, for which 234 source spaces had been created in this way. 235 Inverse modelling 236 For the LCMV beamformer, Eq 3 was used along with a regularised covariance matrix 237 using shrinkage. For MNE, the regularised inverse operator can be defined as where C n ∈ R N,N is the regularised noise covariance matrix and λ 1 is a 239 regularisation parameter. For the simulations, the noise covariance matrix C n was 240 extracted from empty room recordings and then regularised as  Separate simulations were performed for every of the 582 subjects, using the source 260 spaces based on their individual MRIs. Furthermore, simulations were repeated 20 times 261 for every subject, with a randomly chosen signal vertex and randomly generated signals 262 in each iteration, yielding a total of 582 · 20 = 14, 550 MEG simulations. Signal-to-noise 263 ratio (SNR), defined in sensor space as the trace of the signal covariance divided by the 264 trace of the noise covariance, was increased from 10 −6 to 10 0 in 25 logarithmic steps. 265 Finally, to ensure that any difference between the inverse models is not due to a 266 specific choice of the regularisation parameters (the λ's), grid searches were conducted 267 for optimising the parameters. Ten different values for lambda were tested in 268 logarithmic steps, i.e. λ = 10 −9 , 10 −8 , ..., 10 0 . For the FAB beamformer, λ determined 269 the regularisation of the cross-spectral matrices. For the LCMV beamformer, λ 270 determined the regularisation of the covariance matrix. For the MNE, a 271 two-dimensional search grid for λ 1 , λ 2 as introduced above was used.

272
In each iteration, and for each SNR and λ value, the time series in the signal vertex 273 was reconstructed using the three inverse models. For the FAB beamformer, 274 cross-spectra were calculated by splitting the MEG into 1s segments with 50% overlap. 275 Each segment was multiplied by a Hanning window and an FFT was performed, yielding 276 a spectral resolution of 1 Hz. The cross-spectra were then averaged across segments.

277
This way, cross-spectra were obtained for all frequencies between 1 and 100 Hz in 1 Hz 278 steps. For the LCMV beamformer, the covariance matrix was calculated from the  To investigate whether the FAB beamformer outperforms the other methods in real 287 MEG data, resting-state data downsampled to 200 Hz was used as described above. In 288 order to be able to quantify reconstruction quality, a 'ground truth', i.e. a known source 289 time series was required. To this end, a source time series obtained with MNE from a 290 second, randomly selected subject and a randomly selected vertex was used. This time 291 series was assigned to a randomly chosen vertex and then forward projected into the  To reduce computational load, simulations were performed for a medium SNR value 295 of 10 −3 for which FAB, LCMV and MNE differed well in simulation 1. SNR was defined 296 as the trace of the multivariate signal (source signal of second subject mapped to sensor 297 space) divided by the trace of the resting-state data. Since SNR was defined on basis of 298 the broadband data and the power varies substantially as a function of frequency, the 299 effective SNR was somewhat different for each frequency band. However, this was not a 300 problem since we contrasted the different methods within the frequency bands.

301
For FAB and LCMV, the regularisation parameter was fixed to λ = 10 −3 because it 302 was the largest value that had been selected in simulation 1. For MNE, which 303 performed worse than the other two methods in simulation 1, a grid search was again 304 performed in order to maximise MNE performance. Grid search is less costly for MNE 305 since the spatial filters are independent of the MEG data and thus have to be computed 306 only once for every subject. For every subject, the simulation was repeated 25 times. In 307 each iteration, a different source time series from a randomly selected second subject 308 was used as signal. 309 Additionally, to investigate in how far the performance of the FAB beamformer 310 depends on the cross-spectral estimation method, we compared FFT-based estimation 311 using Hanning windows (used in simulation 1) to spectral estimation based on 312 multitapers and Morlet wavelets. To this end, the MEG was split into 2s segments with 313 50% overlap. Multitaper-based estimation was performed using discrete prolate 314 spheroidal sequences (Slepian sequences) with 5 tapers. Cross-spectra were obtained for 315 all frequencies between 1 and 100 Hz in 1 Hz steps and then averaged across tapers and 316 across segments. Wavelet-based estimation was performed using Morlet wavelets in the 317 Fourier domain. The FFT of each segment was tapered with a Gaussian window, the 318 Fourier domain representation of a Morlet wavelet. The time-frequency trade-off was 319 determined by setting wavelet width to 5. For instance, at a frequency of 10 Hz, this 320 trade-off leads to a bandwidth of 10/5*2 = 4 Hz. The cross-spectrum was then 321 calculated by summing across all non-zero frequency bins. Cross-spectra were obtained 322 for all frequencies between 1 and 100 Hz in logarithmic steps and then averaged across 323 Fig 3. Results of simulation 1. (a) Reconstruction quality, measured as correlation between true source signal and the source reconstructed signal a for narrowband signal (top) and a 1/f signal (bottom). The standard error is too small to be displayed. Note that the results for LCMV beamformer (green) and MNE (blue) are almost identical, which is why the MNE data is not easily seen. (b) Bivariate distribution of best regularisation parameters λ 1 and λ 2 across subjects for MNE. (c) Distributions of best regularisation parameter λ for FAB and LCMV beamformer. segments.

324
The resulting broadband source signals were bandpass filtered in the delta (1-4 Hz), 325 theta (4-8 Hz), alpha (8-12 Hz), beta (12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30), lower gamma (30-60 Hz), and upper 326 gamma (60-99 Hz) band. Reconstruction quality was again quantified using Pearson The regularisation parameters have different meanings for different methods, so their 354 distributions were not subjected to hypothesis testing. However, one can qualitatively 355 appreciate that for MNE, the bivariate distribution of λ 1 and λ 2 shows multiple modes. 356 Mostly, large regularisation values were selected for both parameters. In general, the λ's 357 tend to be either large or small, with a low occurrence of intermediate values. The 358 parameters for LCMV and FAB beamforming are comparable to some extent since they 359 both refer to the data covariance. For both beamformers, regularisation tends to be 360 stronger for narrowband signals (blue line) than for 1/f signals (orange). At least for the 361 FAB beamformer, this may be partially explained by the fact that 1/f signals contribute 362 variance at all frequencies, whereas for narrowband signals the information is In medium to high frequencies, the FAB beamformer with wavelet-based estimation outperforms both the LCMV beamformer and MNE. (b) Comparing different ways of calculating cross-spectra in the FAB beamformer. The wavelet-based method outperforms multitapers, although it underperforms for low frequencies. Multitapers in turn outperform estimation based on FFT. (c) Comparing the FAB beamformer using wavelet-based estimation to separate LCMV beamformers that have been applied to bandpass-filtered data. Except for the theta band, where FAB underperforms, a single FAB beamformer is as good as separate LCMV beamformers.
concentrated at a specific frequency, requiring more regularisation for the other 364 frequency bins. Furthermore, the strength of regularisation in LCMV beamforming 365 tends to be an order of magnitude smaller than in FAB beamformer. This could be due 366 to the fact that the LCMV covariance matrix might be slightly better conditioned 367 because it combines information across all frequency bands.  It was then tested whether the FAB beamformer performs better than the other two 372 methods in terms of reconstruction quality across different frequency bands (Fig 4a). To 373 test statistical significance, two different 2-way analyses of variance (ANOVA) were

390
Finally, we compared the wavelet-based FAB beamformer to LCMV beamformers 391 applied to the bandpass-filtered data instead of the broadband data using a 2-way 392 ANOVA. The FAB beamformer tended to be worse, but this was not significant (p = 393 .06). The effect was mostly caused by differences in the theta band (Fig 4c), whereas 394 FAB beamformer and bandpass LCMV beamformer were nearly on par for the other 395 frequency bands.

397
A frequency-adaptive broadband (FAB) beamforming approach was introduced for the 398 source reconstruction of broadband MEEG data. The inverse model was benchmarked 399 against two state of the art competitors, namely the LCMV beamformer and MNE. The 400 Pearson correlation between the known true time series and the reconstructed source 401 was used as reconstruction quality measure.

402
In simulation 1, using narrowband noise sources and a single narrowband or 1/f 403 target signal, the FAB beamformer consistently outperformed both the LCMV 404 beamformer and MNE for a large range of SNR values (Fig 3). The methods converged 405 only for very high or very low SNR values. At low SNR values, reconstruction quality in 406 terms of Fisher-z values increased by up to 90% from FAB as compared to MNE/LCMV. 407 In simulation 2, each subject's resting-state MEG served as real background activity. 408 A source time series from a different subject was used as the target signal. Applied to 409 broadband data, the LCMV beamformer reconstructed lower frequency bands very well 410 but failed to efficiently reconstruct gamma activity (Fig 4). This is not surprising, since 411 the broadband covariance matrix is dominated by lower frequencies. Vice versa, MNE 412 performed relatively poor for low frequency bands but outperformed the LCMV 413 beamformer in the gamma band. Crucially, the FAB beamformer showed a better 414 overall performance than the other two methods, although it performed on par with the 415 LCMV beamformer in the theta band and tended to be worse in the delta band. 416 Furthermore, the effect of spectral estimation procedure on FAB beamformer 417 performance was investigated (Fig 4b). In terms of performance, it was found that 418 wavelet > multitaper > FFT, although the wavelet-based estimation was 419 underperforming for the delta band. This is possibly due to the logarithmic scaling of 420 the frequency bandwidth of Morlet wavelets that does not match the true bandwidth of 421 neural oscillations. Smoothing across frequency bins for low frequencies could increase 422 frequency bandwidth and thereby alleviate this problem. Alternatively, the 423 time-frequency trade-off could be changed towards a more coarse frequency resolution. 424 Apart from this, the overall superior performance of the wavelet-based approach 425 suggests that the distribution of information is broader at higher frequencies. 426 Crucially, except for the delta band, the FAB beamformer performed equally well as 427 separate LCMV beamformers tested on the standard frequency bands (Fig 4c). This 428 implies that the broadband time series obtained with the FAB beamformer has a high 429 signal-to-noise ratio at every frequency.

431
The FAB beamformer reconstructs the broadband time series of a vertex that has a 432 high SNR across the frequency spectrum, unlike DICS and LCMV which are optimal 433 only for a frequency band. Therefore, FAB beamforming can serve as a first step for 434 further analysis in source space, for instance decoding of broadband ERPs or 435 connectivity analysis. It might also be useful for representing neural events with 436 complex neural signatures. For instance, during sleep K-complexes are often followed by 437 sleep spindles. A single frequency band might not be able to capture the full dynamics 438 of such a complex signal. Finally, in this paper we focused on source reconstruction.