SPHARA - A Generalized Spatial Fourier Analysis for Multi-Sensor Systems with Non-Uniformly Arranged Sensors: Application to EEG

Important requirements for the analysis of multichannel EEG data are efficient techniques for signal enhancement, signal decomposition, feature extraction, and dimensionality reduction. We propose a new approach for spatial harmonic analysis (SPHARA) that extends the classical spatial Fourier analysis to EEG sensors positioned non-uniformly on the surface of the head. The proposed method is based on the eigenanalysis of the discrete Laplace-Beltrami operator defined on a triangular mesh. We present several ways to discretize the continuous Laplace-Beltrami operator and compare the properties of the resulting basis functions computed using these discretization methods. We apply SPHARA to somatosensory evoked potential data from eleven volunteers and demonstrate the ability of the method for spatial data decomposition, dimensionality reduction and noise suppression. When employing SPHARA for dimensionality reduction, a significantly more compact representation can be achieved using the FEM approach, compared to the other discretization methods. Using FEM, to recover 95% and 99% of the total energy of the EEG data, on average only 35% and 58% of the coefficients are necessary. The capability of SPHARA for noise suppression is shown using artificial data. We conclude that SPHARA can be used for spatial harmonic analysis of multi-sensor data at arbitrary positions and can be utilized in a variety of other applications.


Introduction
The discrete Fourier analysis of 2D data defined on a flat surface and represented by a Cartesian or a regular grid is very common in digital image processing and a fundamental tool in many applications. For such data, the basis functions (BF) for the Fourier transformation are usually implicitly specified in the transformation rule, compare [1].
In many applications, however, the sensors for data acquisition are not located on a flat surface and can not be represented by Cartesian or regular grids. One important example for nonregular sensor positions is the application of electroencephalography (EEG). In EEG, the sensors are placed at predetermined positions on a surface in space R 3 , the head shape. The positions of the sensors in these systems can be described by means of triangular meshes. Because of the particular sensor arrangement, the spatial analysis of multi-sensor data can not be performed using a standard 2D Fourier analysis.
Several methodologies have been proposed previously that are applicable for signal enhancement, feature selection and dimensionality reduction of multichannel EEG data. Examples are the principal component analysis (PCA), the independent component analysis (ICA), the parallel factor analysis (PARAFAC) and matching pursuit (MP) using Bessel or multichannel Gabor atoms.
PCA, also known as Karhunen-Loève transform, is a statistical procedure commonly used to reduce a large number of variables to a smaller number of principal components (PC). These components are pairwise orthogonal and describe the maximum amount of variance in the original multivariate data. PC are determined by computing the eigensystem of the covariance matrix of the measured data. PCA is used to decompose multichannel EEG data into orthogonal components with respect to their spatial distribution. PCA is applied for signal enhancement, feature selection and dimensionality reduction in BCI [2][3][4], for source localization [5,6] and for artifact recognition and removal [7]. However, the basis functions of the PCA are determined by using the measured EEG data. Thus, PCA can be used only after a certain number of time samples are recorded.
Another linear decomposition technique is ICA. In contrast to PCA, the ICA approach maximizes the statistical independence between the signal components (independent components). Furthermore, basis vectors of the ICA are in general not orthogonal and there is no closed form expression to calculate the components. Therefore, iterative approaches are used to determine the individual components, which results in higher computational costs. One application of ICA is the separation of brain activity from artifacts. It is assumed that brain activity and artifacts are statistically independent processes, which result in statistically independent components of the measured signal [8][9][10][11]. Another application of ICA is the decomposition of evoked fields or potentials [12,13]. Furthermore, ICA is employed to separate spatially overlapping EEG activities [10]. In BCI applications, ICA is used for signal enhancement and feature selection [3,14,15]. Like PCA, ICA can be used only after a certain number of time samples is recorded.
Both PCA and ICA are used to decompose two-way data. Usually they are employed to perform a space-time decomposition of multichannel EEG signals. Using parallel factor analysis (PARAFAC), also known as canonical (polyadic) decomposition, a multi-way decomposition of multichannel EEG data is possible. By employing PARAFAC, wavelet transformed multichannel EEG signals can be decomposed into channel, frequency and time components [16]. In another application of PARAFAC, a five-way analysis of EEG data is performed and the recorded data are decomposed into channel, frequency, time, subject and condition components [17]. PARAFAC requires a high computational effort. Further, preprocessing steps, such as a time-frequency decomposition, are often necessary. In BCI, PARAFAC is used for feature selection [18].
The Matching Pursuit method is used to decompose single channel EEG signals into timefrequency components using an over-complete dictionary [19]. There also exists an extended version of the matching pursuit method for the application to multichannel EEG data using multichannel Gabor atoms [20,21]. In addition, Bessel atoms are used for the spatial decomposition of MEG data [22]. When using Bessel atoms, the basis functions are determined prior to data recording. However, discretized Bessel functions are only applicable for the spatial decomposition of the data of certain sensor configurations; otherwise, the required orthogonality of the Bessel functions is lost. As part of the data decomposition, a very large overdetermined dictionary with basis functions has to be searched repeatedly.
Another approach for spatial analysis of EEG data is the use of spherical harmonics, where the head shape is approximated by a spherical surface [23]. The effect of inevitable mapping errors that arise in this approach on the spatial spectrum is difficult to estimate. Moreover, due to the anatomy of the head, only parts of the spherical surface can be sampled in EEG, therefore, they must be combined with window functions.
Laplace spectra, which we use for the spatial decomposition, are widely applied in the field of graph theory [24] and computer graphics [25][26][27]. A similar approach, called Spherical Cap Harmonic Analysis is used in geosciences to analyze the earth magnetic field [28,29].
Surface Laplacian approaches are applied in EEG to improve the spatial resolution [30][31][32]. Laplacian eigenmaps are also deployed for EEG studies [33,34] to investigate EEG time series and topographies. However, the Laplacian eigenbasis of these two approaches is not taking into account the spatial arrangement of the EEG electrodes in R 3 .
In this publication we present a new method for spatial harmonic analysis (SPHARA) of EEG data using the eigenbasis of the Laplace-Beltrami operator of the meshed surface of electrode positions. Using this approach, BF of spatial harmonics for arbitrary arrangements of EEG electrodes can be generated. The recorded EEG data are decomposed by projection into the space of the BF. In addition, a number of possible applications are presented, such as analysis of the spatial properties of SEP data, dimensionality reduction and noise reduction.
We presented preliminary research results in the field of spatial harmonic analysis of multisensor in a talk at the 8th International Symposium on Noninvasive Functional Source Imaging of the Brain and Heart & International Conference on Bioelectromagnetism -NFSI & ICBEM, 2011, Banff [35].

Overview
In this section, the new method for the spatial harmonic analysis of EEG data is derived. Starting from the well-investigated continuous case, we introduce several discretization schemes, which can be applied to EEG sensor arrangements. The BF, which are used for the spatial signal decomposition, are determined by solving the eigenvalue problem for the Laplace-Beltrami operator, which is also referred to as the Helmholtz equation. We show relevant properties of the BF that are necessary for the application in a signal processing frameworks. Finally, computational complexity and algorithmic aspects are considered.

Continuous case
Continuous Laplace-Beltrami Operator. We assume a compact Riemannian manifold (M, g) of dimension m, where M is a connected manifold and real-differentiable of class C 1 . The function g defines for each point p 2 M the inner product of the tangent space T p M. The union of all tangent spaces of M is TM. An m-dimensional manifold is a topological space, that locally resembles the Euclidean space R m . Within this paper, we focus on two-dimensional manifolds representing surfaces in R 3 .
If the manifold M has a boundary B = @M, it is assumed that M is oriented and that C 1 also applies to B. The outward unit normal vector field on B is denoted by n.
We consider a real-valued function f with f 2 L 2 (M, g) and f is C k with k ! 2. The directional derivative of f at p 2 M for each ξ 2 T p M is denoted by ξf. The gradient of f is the vector Let X, Y be vector fields, which are C k with k ! 1, and let r ξ X be the covariant derivative of X with respect to ξ for all p 2 M and ξ 2 T p M. And r ξ additionally fulfills r ξ (X + Y) = r ξ X + r ξ Y and r ξ (fX) = (ξf)X(p) + f(p)r ξ X. The divergence of X is defined by where ξ ranges over T p M. The continuous Laplace-Beltrami operator Δ is defined by For a more extensive derivation of the continuous Laplace-Beltrami operator see for example [36,37]. If the manifold M possesses a boundary B, several boundary conditions (BC) can be applied with the nabla operator r and η 6 ¼ 0, z 6 ¼ 0. The utilized boundary condition affects the solution of Eq (7) and thus the properties of the determined basic functions. When using Dirichlet BC, the values of functions f on the boundary B of the manifold M are set to 0. For Neumann BC, the normal derivative of the function f on the boundary B of the manifold M is fixed to 0. The Robin BC is a linear combination of Dirichlet and Neumann BC using the coefficients η and z. In our application, the manifold M describes the area of the head, which is covered by the EEG sensor system. The potential field of the head surface extends beyond the range that is covered by M. The potentials on the boundary B of the manifold is generally not 0 or constant. Therefore, the Neumann BC is more appropriate to compute the basis functions using the presented approach.
Properties of Continuous Laplace-Beltrami Operators and their Eigensystems. The operator Δ possesses a number of properties that enable the computation of harmonic BF by solving the Laplacian eigenvalue problem These properties are: Δf = 0 for constant f, symmetry, local support, linear precision, maximum principle, and positive semi-definiteness [38]. A basis for a harmonic analysis of functions, which are defined on the Riemannian manifold, can be determined by finding all eigenvalues λ i 2 R and the associated non-trivial eigenfunctions ϕ i , with ϕ i 2 C 2 (M). The set of all eigenvalues fl i g 1 i¼1 defines the spectrum of M specðMÞ with lim i ! 1 λ i ! 1. The set of the eigenfunctions f i g 1 i¼1 forms an orthonormal basis and can be used for the spectral analysis of functions defined on M. The eigensystem of the Laplace-Beltrami operator can be considered as basis of a generalized Fourier analysis on M, see also [36,37,39]. The continuous approach to determine the BF for M will be referred to as CONT later in the text.
For any point p i on M, an approximation Δ A for Δ is given by the curvilinear integral where γ is a closed simple curve on M surrounding the point p i , and jγj is the length of γ [40,41], see Fig 1(a). In practical applications, only a few discrete points of γ are usually known. In the recording of EEG data, the continuous function f describing the potential distribution on the surface of the head is sampled only at the electrode positions. For this reason, discretization approaches for Δ are necessary. If a geometric discretization approach is used, the discrete Laplace-Beltrami operator converges to the continuous form Δ D ! Δ when the grid is refined [38].

Discrete case
Discrete representation of surfaces. For the discrete case we assume that the sensors are located on an arbitrary surface, which is represented by a triangular mesh in share the edge e ij . The set of triangles sharing the vertex v i is given by The area of a triangle t is given by jtj. An example for these mesh components is illustrated in Fig 1  (b).
A discretization Δ D of the continuous Laplace-Beltrami operator in Eq (9) forf is given by with the weighting function w(i, x) for edges e ix 2 E and the normalization coefficient b i for the vertex v i . For practical applications it is convenient to transform Eq (10) into matrix notation. The elements of the matrices B −1 and S are determined using the coefficients b i and w(i, x) of Eq (10). B −1 is a diagonal matrix, the elements are and the entries of S are A Laplacian matrix L can be expressed as product of a diagonal matrix B −1 and a matrix S compare also [42]. The size of the Laplacian matrix L for a mesh M is n × n, with n = jVj. Using the Laplacian matrix L, Δ D applied tof can be written as In the following we present four different approaches to discretize the Laplace-Beltrami operator, a graph-theoretical approach and three geometric approaches. First, we look at the graph-theoretical approach, where the coordinates of the positions of the vertices are not considered. The topological Laplacian results from Eq (10) by using [24,26,40]. The graph-theoretical approach will be referred to as TL later in the text.
Second, for inhomogeneous triangular meshes, where the distances between vertices and the sizes of angles and triangles are different, the weighting function w has to be adapted according to the mesh geometry. In these approaches, the positions of the vertices are also considered. They are referred to as geometric approaches. There are different approaches to treat inhomogeneous meshes.
The first possibility is to use the Euclidean distance of adjacent vertices raised to the power of a value α. For Eq (10) the coefficients b À1 i ¼ 1 and w(i, x) = ke ix k α are chosen. A common choice is to use the inverse of the Euclidean distance with α = −1 [40,43]. This approach will be referred to later as IE.
The second approach for a geometric discretization of the Laplace-Beltrami operator is derived by minimizing the Dirichlet energy for a triangulated mesh [44,45]. It uses cotangent weights with with the two angles α ix and β ix opposed to the edge e ix , see Fig 1(b). For edges on the boundary of the mesh, the term cot(β ix ) is omitted, which leads to Neumann BC. A drawback of using the cotangent weights is that the value representing the integral of the Laplacian over a 1-ring neighborhood (area of the i ? -neighborhood) is assigned to a point sample [26]. To resolve this issue and to guarantee the correspondence between the continuous and the discrete approaches, the weights in Eq (15) are divided by the area A B i of the barycell for the vertex v i [46], resulting in The barycell for a vertex v i is framed by a polygonal line that connects the geometric centroids of triangles in i ! and the midpoints of the adjoined edges e ix , see Fig 1(c). The area of the i ?neighborhood for a vertex v i , which is the area of the triangles that are enclosed by the vertices v x 2 i ? , is referred to as For the discretizations using the cotangent weighted formulation, the parameter b À1 i in Eq (10) is set to b À1 i ¼ 1. This approach, using cotangent weights will be referred to as COT later in the manuscript.
The third geometric approach to discretize the Laplace-Beltrami operator is the finite element method (FEM), which is related to the approach using cotangent weights. Assuming that the function f is piecewise linear and defined by its values f i on the vertices v i of a triangular mesh, f can be interpolated using nodal basis functions We use the hat function for ψ i , with For two functions f and g defined on M, a scalar product is given by with the area element da on M and the mass matrix B. The sparse mass matrix B is given by For the FEM approach using hat functions, the elements of B can be calculated by where t a and t b are the two triangles adjacent to the edge e ij , see Fig 1(b). For the FEM discretization of the Laplace-Beltrami operator also a stiffness matrix S has to be calculated. The elements of S ij can be estimated using the Eqs (12) and (15), compare [26,47,48].
To simplify the FEM approach, the mass matrix B can be replaced by a lumped mass matrix B, which is defined asB The simplified FEM approach, which uses a lumped mass matrixB, is equivalent to the method using cotangent weights presented in Eq (16), see [48].
Properties of Discrete Laplace-Beltrami Operators and their Eigensystems. Desirable properties of the discrete Laplacian L are symmetry, positive weights, positive semi-definiteness, locality, linear precision and convergence [38]. The symmetry L ij = L ji leads to real eigenvalues and orthogonal eigenvectors. Positive weights w(i, j) ! 0 assure, together with the symmetry, the positive semi-definiteness of L. The locality of the discrete Laplace-Beltrami operator enables the determination of weights w(i, j) using the i ? -neighborhood of a vertex v i , with w(i, j) = 0, if e ij = 2 E. The linear precision implies for a linear function f defined on vertices v i that D Df i ¼ 0 applies, which ensures the exact recovery of f from the samples. The convergence property provides the convergence from the discrete to the continuous Laplace-Beltrami operator Δ D ! Δ for a sufficient refinement of the mesh. The Laplacian matrix L for the TL and the IE approach are positive semi-definite, symmetric and use positive weights. The COT and the FEM approach do not fulfill the positive weight property, if the mesh contains triangles with interior angles in the interval (π/2, π), for which the cotangent is negative. The TL approach is no geometric discretization, because it violates the linear precision and the convergence property. In contrast, the COT and the FEM approach are geometric discretizations as they fulfill the linear precision and the convergence property, but they violate the symmetry property. None of the presented discretization methods fulfill all desirable properties, see also [38].
The discrete Laplacian eigenvalue problem for a real and symmetric Laplacian matrix is given by with eigenvectorsx i and eigenvalues λ i of L. The Laplacian matrix L is real and symmetric for the TL, the IE and the COT approach. Because L is real and symmetric, eigenvalues λ i 2 R with λ i ! 0 are obtained. The eigenvectorsx i are real-valued and form a harmonic orthonormal basis. The corresponding eigenvalues λ i can be considered as spatial frequencies. The eigenvectorsx i can be used for a spectral analysis of functions defined on the mesh M. The projection of a discrete functionf defined on M onto the basis of spatial harmonic functions is performed by the inner product for Euclidean n-spaces hf ;x i i. For the matrix X, where the eigenvectorsx i represent the columns, it applies with the identity matrix I. For the FEM formulation, the BFỹ i are computed by solving the generalized symmetric definite eigenproblem Thus, the inversion of the mass matrix B is avoided. Because B −1 S is not symmetric, the eigenvectorsỹ i are real-valued, but not orthonormal with respect to the inner product for Euclidean n-spaces h.i. To use these eigenvectors as BF, the inner product, defined in Eq (19), has to be used which assures the B-orthogonality, compare also [48]. The eigenvectors computed by the FEM approach can be normalized by using the B-relative norm For a matrixỸ , where the normalized eigenvectorsỹ i represent the columns it appliesỸ Discrete spatial harmonic signal processing framework Requirements. To use the eigenvectors of the discrete Laplacian Beltrami operator in the context of a signal processing framework, it is necessary that they exhibit certain properties. The eigenvectors have to form a set of BF. An inner product for the decomposition and for the reconstruction of the data has to be defined, which is used for the transformation into the domain of the spatial frequencies and for the back-transformation into the spatial domain. To be utilized for practical applications, the transformation from the spatial domain into the spatial frequency domain has to have linear properties and should fulfill Parseval's theorem.
Basis functions. A complete set of linearly independent vectors can be used as basis. For real and symmetric Laplacian matrices L the orthonormality of the eigenvectors is given inherently, see Eqs (23) and (25). For the FEM approach, the orthonormality of the eigenvectors is assured explicitly, see Eqs (26)- (30). The property of orthogonality includes the linear independence. To use the eigenvectors as BF, they must further fulfill the property of completeness. The completeness can be shown by the dimension theorem for vector spaces. The dimensionality is equal for both the spatial representation and the representation in the spatial frequency domain. For a mesh with n vertices, n unit impulse functions are used as BF for the spatial representation. For the same mesh, we obtain n discrete spatial harmonic functions (eigenvectors) for the representation using spatial frequencies. The calculated eigenvectors are orthonormal and complete; therefore, they can be used as orthonormal BF.
Decomposition and reconstruction. For the decomposition of discrete data defined on the vertices of the triangular mesh, the inner product is used (transformation from spatial domain to spatial frequency domain). For a decomposition using the eigenvectors of a symmetric Laplacian matrix L (TL, IE and COT), the vector space inner product is applied. The coefficient c i for a single spatial harmonic BFx i can be determined by The transformation from the spatial into the spatial frequency domain is computed bỹ For a decomposition using eigenvectors computed by the FEM approach, the inner product that assures the B-orthogonality needs to be applied The transformation from the spatial into the spatial frequency domain is then be computed bỹ Discrete data are reconstructed using the linear combination of the coefficients c i and the Linearity. The property of linearity is required for many applications in which a spatial harmonic decomposition is performed. Eq (32) describes a linear transformation from the spatial into the spatial frequency domain by using the transformation matrix X. For two transformationsf with the functions in the spatial domain,f i , and functions in spatial frequency domain,c i , and scalar values a 1 and a 2 , it follows The property of linearity can be shown analogously for transformations corresponding to Eq (34) using the operator BỸ .
Parseval's theorem. The energy representations of a signal in the spatial domain and in the spatial frequency domain are equivalent; they must possess the same energy (Parseval's theorem). For a discrete real-valued functionf defined in the spatial domain and for its representationc in the spatial frequency domain, which was calculated by the eigenvectors of a real and symmetric Laplacian matrix L (TL, IE and COT), it applies hf ;f i ¼ hc;ci : When using the FEM approach, where the B-relative inner product has to be employed to compute the signal energy in the spatial domain, it applies

Numerical considerations
The Laplacian operator L is a sparse square matrix n × n, where n is the number of vertices. In the case of a symmetric Laplacian matrix (TL, IE and COT), the eigenvalues and eigenvectors of L are computed using a parallel implementation of the LAPACK routine DSYEVR [49]. DSYEVR consists of the two subroutines DSYTRD and DSTEMR. The DSYTRD-routine performs a Householder tridiagonalization of L and has a complexity of O(n 3 ) [50]. The second subroutine DSTEMR computes the eigenvalues and eigenvectors of the resulting real symmetric tridiagonal matrix and has a complexity of O(n 2 ) [51]. Thus, the computational complexity of DSYEVR is O(n 3 ). For asymmetric Laplacian matrices L (FEM approach) the eigenvalues and eigenvectors are computed using the LAPACK routine DGEEV, which has a complexity of O(n 3 ). To compute the eigenvalues and eigenvectors, this method performs the reduction to an upper Hessenberg form (DGEHRD) and a Schur factorization (DHSEQR) [49].
For the FEM approach, a symmetric definite eigenproblem has to be solved. Thus, the LAPACK routine DSYGV can be applied. By using a Cholesky factorization of the Laplacian matrix, the problem is reduced to a standard symmetric eigenproblem [49]. The complexity of DSYGV is also O(n 3 ). When considering the partial eigenproblem only and by using multigrid approaches and preconditioning, the computational complexity of the FEM approach can be reduced to O(n log n) [52]. It should be noted that two matrices have to be computed in the FEM approach for the data decomposition, since the mass matrix B is also required for the inner product h.i B .
The inner product is used for the decomposition of the data defined on n spatial sample points. For a single BF, the complexity of computing the inner product is O(n); n multiplications and n − 1 additions of floating point numbers are carried out. The multiplications are independent; the additions can be cascaded. The complete decomposition of the data using all n normalized BF is performed by the multiplication of a vector containing the data and a matrix, where the columns are the normalized BF. The complexity is O(n 2 ).

Spatial decomposition methods used for comparison with SPHARA
As mentioned in the introduction, several methods for spatial decomposition of multichannel data are already available. In order to assess the performance, we have compared SPHARA with PCA and ICA. When using SPHARA, EEG data are decomposed by projection into the space of basis functions. This is similar to PCA or ICA, where projections into the space of principal or independent components are used. For the comparison, we utilized an own implementation of PCA and the Infomax-based ICA [8].

Data
We illustrate the application of our proposed method to data sets that originate from a previously performed EEG experiment addressing the cortical activation related to somatosensoryevoked potentials (SEP). Eleven healthy male volunteers, aged between 26 and 40 years (mean 30.7 ± 4.3 years), ten right-handed and one left-handed, participated in the experiment. No volunteer suffered from neurological disease, and all participants were free of medication. The study was approved by the Ethics Committee of the medical faculty of the Friedrich Schiller University Jena and written informed consent was obtained. The median nerve of the right forearm was stimulated by bipolar electrodes (stimulation rate: 3.7 Hz, interstimulus interval: 270 ms, stimulation strength: motor plus sensor threshold [53,54], constant current rectangular pulse wave impulses with a length of 50 μs, number of stimulations: 6000).
The positions of the EEG electrodes were digitized using an optical ANT Neuro BV (antneuro.com) xensor 3D electrode digitizer system, see Fig 2. EEG signals were recorded using an ANT Neuro BV waveguard 256-channel EEG cap with an equidistant electrode layout and two cascaded RefaExt 128 channel amplifiers and transformed to common average reference. Data were sampled at 2048 Hz and software high-pass (24 dB/oct, cutoff-frequency 2 Hz) and notch (50 Hz and two harmonics) filtered. All trials were manually checked for artifacts, the remaining trials were averaged, see S1 Dataset. The mean global field power (MGFP) and topographic maps of the averaged and filtered SEP data of a single volunteer at selected points in time are shown in Fig 3. In summary, the grand average (GA) and the range of ± the standard deviation of the MGFP of all subjects is shown in Fig 4. The exact points in time of the maximum potentials vary between volunteers. For each volunteer, the specific time of the P14, N20 and N30 potentials are manually estimated.
To demonstrate the ability of the noise suppression of SPHARA, Gaussian white noise was added to averaged SEP data. The following signal-noise-ratios (SNR) were used (dB): 10, 6, 3, 0, -3, -6, -10. In the simulation, the average SEP data were used as reference (signal). For each SNR 100 000 noise realizations were analyzed.

Results
The SEP data of the eleven volunteers were decomposed by using the base functions for each of the three discretization approaches. The procedure of the spatial harmonic decomposition of the SEP data for one volunteer is illustrated in Fig 5. The SEP data consist of 256 EEG-channels, which are shown in Fig 5(a). The decomposition was performed by calculating the inner product of the multichannel data and the BF, see Fig 5(d). The power contribution of the Generalized Spatial Fourier Analysis for Multi-Sensor Systems spatial harmonic BF to the SEP data at each time step is shown in Fig 5(b), which shows the time versus spatial harmonic frequency representation of the decomposed data.
We investigated how many of the spatial-harmonic coefficients are required to achieve 90%, 95% and 99% of the signal energy of the samples in the time range from 50 ms before to 130 ms after stimulation. We compared BF, computed by the TL, IE and FEM approaches. In the first case, the individually tracked sensor positions of the volunteers were used to determine  Generalized Spatial Fourier Analysis for Multi-Sensor Systems the BF using the IE and FEM approach. The BF of the TL approach can be determined without information of the sensor positions. The BF with the largest energy contribution are used for signal reconstruction. The results are shown in Fig 6. The available compression ratios in the best, in the median and in the worst case are listed in Table 1. The compression ratio is the quotient of the maximum number of available BF and the number of BF that are used to restore the signal at a given quality. A maximum of 256 basis functions are available to reconstruct the entire signal energy. The significance of the influence of the three discretization methods on the achieved compression ratio was analyzed using the Mann-Whitney U test. This test was chosen because not all examined data are normally distributed, tested by the Pearson χ 2 test. None of the methods is significantly superior in reconstructing 90% of the signal energy. There is no significant difference between the TL and IE approach in the ability to present EEG data compactly. FEM provides significantly better results than TL and IE in the reconstruction of 95% and 99% of the signal energy, see Fig 6 and Table 2. The main contribution to the signal power of the SEP data is provided by the low frequency spatial-harmonic BF as exemplarily shown in Fig 5(b). For the FEM approach, we have calculated two sets of BF, first using individually tracked sensor positions (ITSP) and second using standard sensor positions (SSP) provided by the manufacturer of the EEG cap. There is no difference between the two methods in the compact representation of the analyzed EEG data. In the median case, the same number of BF are necessary to achieve 90%, 95% and 99% of the total signal energy, see Fig 7. We compared the FEM approach with PCA and ICA regarding the ability of a compact representation of EEG data in a spatial aspect. In PCA, the PC are sorted in decreasing order by there variances, which represents their energy contribution to the investigated data. In the Infomax-based ICA, the independent components (IC) are sorted in descending order of mean  projected variance. The maximum number of IC that are to be calculated is determined by the rank of the covariance matrix of the analyzed data. The number of BF, PC and IC to achieve 90%, 95% and 99% of the signal energy of the investigated EEG data are evaluated. The results are shown in Fig 8. PCA comprises clearly the best performance on the compact representation of the EEG data. The second best compression ratio achieved the SPHARA approach. ICA revealed the worst performance for the investigated EEG data. The capability of the new method for noise suppression is demonstrated by implementing a spatial low pass filter. For this purpose, we designed a filter with 25 low-frequency basis functions. Depending on the data of the volunteers, between 96.4 and 98.6% of the signal energy of the EEG data can be reconstructed with 25 low-frequency basis functions. In the simulation,  Generalized Spatial Fourier Analysis for Multi-Sensor Systems performed to investigate the noise reduction, Gaussian white noise with different noise ratios was added to the averaged SEP data. Subsequently, the noisy data were low-pass filtered with SPHARA. An example of the averaged SEP data, with additional Gaussian white noise and the results of the spatially harmonic low-pass filtering are shown in Fig 9. In a simulation with 100 000 repetitions, the ability of the spatial harmonic analysis for noise reduction could be shown. The improvement of the SNR by spatially harmonic low-pass filtering is in the median case be-

Discussion
The eigensystem of the discrete Laplace-Beltrami operator defined on a triangular mesh forms a set of harmonic orthonormal BF, which can be adaptively computed using the topology of the sensor setup and the sensor positions in R 3 . The set of BF enables the spatial Fourier analysis of arbitrarily arranged multi-sensor data. Example of noise suppression using spatial harmonic low pass filter. Left column unfiltered data, right column spatial harmonic low-pass filtered data; in the first row averaged SEP data without additional noise; in the following rows SEP data with additional Gaussian white noise with 3, 0, -3 dB signal to noise ratio. Activations in the human brain cause electrical potentials around the entire head. For anatomical reasons, it is difficult or impossible to observe the electrical potential around the whole head. In EEG measurements, the head surface is only partially covered by sensors, therefore only a part of the electrical potentials that emerge at the surface of the human head are captured. The potential on the boundary of the examined area of the surface are not zero or constant. For this reason, Neumann boundary conditions are more appropriate to compute the BF.
For the data analysis described in this article, three discretization schemes are used to compute the BF. The three discretization methods require different information about the triangular mesh for the calculation of the BF. Also the effort of calculating the BF and of the decomposition of the EEG data is different. For the determination of BF by means of TL, only the layout of the triangle mesh must be known. The particular information about coordinates of the sensor positions are not necessary. For different EEG cap sizes with the same triangulation, the same BF can be used. These BF can be already determined during the design process of the EEG cap. For the calculation of the BF using the IE and the FEM approach, additional information about the sensor positions are necessary. The eigensystem of the discrete Laplace-Beltrami operator needs to be computed only once for a EEG electrode setup, if the electrode positions are not individually determined. If the electrode positions are individually tracked for each volunteer, the BF have to be determined for each measurement session. In contrast to TL and IE, the Laplace-Beltrami operator using the FEM approach represents a geometric discretization and fulfills the convergence property. It converges to the continuous solution for a sufficient refinement of the mesh. Subsequently, the decomposition of the recorded multi-sensor Fig 10. Noise suppression using spatial harmonic low pass filter. Improvement of the SNR (shown on the Y-axis) using a spatial harmonic low-pass filter applied to data with different noise ratios (shown on the Xaxis). data for a single point in time is carried out by a multiplication of the matrix with the BF and a vector with the measurement data. For the data decomposition using the FEM approach, in addition to the matrix of the BF the mass matrix is required and the computational effort is slightly increased compared to TL and IE.
We compared TL, IE and FEM with regard to their applicability for a compact representation of EEG data. When reconstructing 90% of the signal energy of the EEG data, none of the discretization methods is superior. Using the FEM approach, a significantly more compact representation of data can be achieved, in the case of 95% and 99% reconstructed signal energy. However, compared to TL and IE it requires a slightly higher effort for the determination of the BF and the data decomposition.
We determined the BF for IE and FEM using for each volunteer individually tracked and standard sensor positions, and we compared the influence on the compact representation of the EEG data. For the considered data, no differences were observed in the median case. However, there are applications conceivable in which individually tracked sensor positions are advantageous, e. g. approaches that incorporate individual volume conductor models of the heads of volunteers. This can be the field of future investigations.
In terms of compact representation of the EEG data, PCA has performed better and ICA worse than SPHARA. The better performance of PCA can be explained by the fact, that the maximization of the energy, represented by the principal components, is a design criterion of PCA. To achieve this, the measured time series need to be analyzed. ICA is a decomposition method, which splits the multivariate data into statistically independent components. These statistically independent components do not necessarily have a high energy content. For most multivariate decomposition methods, such as PCA and ICA, the recorded multivariate time series are employed to generate the components for a spatial decomposition. In contrast to these methods, in SPHARA the BF for the spatial decomposition are generated using topological information of the electrode arrangement and the sensor positions. Therefore, the BF can be determined prior to the recording of the time series. This is particularly beneficial for online data processing and for other time critical applications.
The ability of the new approach to denoise EEG data is demonstrated by averaged SEP data with additional artificial noise. Utilizing the SPHARA BF a spatial low-pass filter was implemented. Depending on the amount of added artificial noise, an median improvement of SNR between 4.31 dB and 9.74 dB could be achieved by the spatial low-pass filtering. The filtering is done only in the spatial domain. Different time samples of the EEG data are considered separately. This facilitates the online denoising of EEG data.
The SPHARA approach can also be combined with frequency and time-frequency analysis methods. Furthermore, the spatial-harmonic approach can be applied to various fields where multi-sensor data are measured on surfaces and manifolds in R 3 .
In most multivariate decomposition methods, such as PCA, ICA and PARAFAC, the recorded multivariate time series are employed to generated components for a spatial decomposition. In contrast to these methods, in SPHARA the BF for the spatial decomposition are generated using topological information of the electrode arrangement and the sensor positions solely. Therefore, the BF can be determined prior to the recording of the time series.
Supporting Information S1 Dataset. SEP data that has been analyzed. The data set contains averaged and filtered SEP data of the volunteers as well as the standard electrode positions and triangle list of the ANT Neuro BV waveguard 256-channel EEG cap. (ZIP)