SPHARA – The theoretical background in a nutshell

Motivation

The discrete Fourier analysis of two-dimensional data, defined on a flat surface and represented by a Cartesian or regular grid, is a fundamental tool in digital image processing. For such data, the basis functions (BF) of the Fourier transformation are implicitly specified by the transformation rule, compare [Rao et al., 2010].

In many applications, however, sensor positions are not arranged on a flat surface and cannot be described by Cartesian or regular grids. A typical example in biomedical engineering is electroencephalography (EEG), in which sensors are placed at predetermined positions on the head surface in \(\mathbb{R}^3\). Such sensor setups can be represented by triangular meshes. Because of the irregular spatial arrangement, standard 2D Fourier analysis is not applicable. Nevertheless, a spatial harmonic analysis remains useful even for irregularly sampled data.

This Python package implements a method for SPatial HARmonic Analysis (SPHARA) of multisensor data using the eigenbasis of the Laplace–Beltrami operator of the mesh describing the sensor positions. With this approach, spatial harmonic basis functions can be constructed for arbitrary sensor arrangements. Recorded multisensor data can be decomposed by projection onto these basis functions. For a more comprehensive introduction see [Graichen et al., 2015, Graichen et al., 2019].

Discrete representation of surfaces

For the discrete setting we assume that the sensors lie on an arbitrary surface represented by a triangular mesh in \(\mathbb{R}^3\). The mesh \(M = \{V,E,T\}\) consists of vertices \(v \in V\), edges \(e \in E\), and triangles \(t \in T\). Each vertex \(v_i \in \mathbb{R}^3\) corresponds to a sensor position. The neighbourhood \(i^{\star}\) of a vertex \(v_i\) is defined by

\[i^{\star} = \{ v_x \in V : e_{ix} \in E \}.\]

Angles \(\alpha_{ij}\) and \(\beta_{ij}\) are opposite the edge \(e_{ij}\), and the set of triangles sharing a vertex \(v_i\) is

\[i^{\triangledown} = \{ t_x \in T : v_i \in t_x \}.\]

Figure Fig. 1 illustrates these components.

Fig. 1 Approximation of the Laplace–Beltrami operator. (a) continuous representation; (b) discrete representation; (c) area of the barycell \(A_i^\mathrm{B}\) for vertex \(v_i\).

Discrete Laplace–Beltrami operators

A function \(\boldsymbol{f}\) defined on the vertices \(v_i\) associates a real value with each vertex. A discrete approximation \(\Delta_D\) of the Laplace–Beltrami operator is given by

(1)\[\Delta_D \boldsymbol{f}_i = b_i^{-1} \sum_{x \in i^{\star}} w(i,x)\,\bigl( \boldsymbol{f}_i - \boldsymbol{f}_x \bigr),\]

with edge weights \(w(i,x)\) and normalization coefficients \(b_i\).

Using the coefficients of equation (1), matrices \(\boldsymbol{B}^{-1}\) and \(\boldsymbol{S}\) are constructed as

(2)\[\begin{split}\boldsymbol{B}^{-1}_{ij} = \begin{cases} b_i^{-1} & i = j, \\ 0 & \text{otherwise}, \end{cases}\end{split}\]

and

(3)\[\begin{split}\boldsymbol{S}_{ij} = \begin{cases} \sum_{x \in i^{\star}} w(i,x) & i = j, \\ -w(i,x) & e_{ix} \in E, \\ 0 & \text{otherwise}. \end{cases}\end{split}\]

The Laplacian matrix \(\boldsymbol{L}\) is then

(4)\[\boldsymbol{L} = \boldsymbol{B}^{-1}\boldsymbol{S}.\]

The operator \(\Delta_D\) applied to \(\boldsymbol{f}\) can be written as

(5)\[\Delta_D \boldsymbol{f} = - \boldsymbol{L} \boldsymbol{f}.\]

Graph-theoretical and geometric discretization approaches

Four approaches are commonly used to discretize the Laplace–Beltrami operator:

1. Topological Laplacian (TL) No geometric information is used. \(w(i,x) = 1\), \(b_i^{-1}=1\), see [Taubin, 1995, Chung, 1997].

2. Inverse Euclidean (IE) Geometric distances are considered via \(w(i,x) = \lVert e_{ix} \rVert^\alpha\) with \(\alpha=-1\).

3. Cotangent-weighted Laplacian (COT) Based on minimizing the Dirichlet energy [Pinkall and Polthier, 1993, Polthier, 2002], using weighted cotangents of opposite angles.

4. Finite Element Method (FEM) The most faithful geometric discretization, based on hat functions \(\psi_i\) and integrals of these basis functions over the mesh, resulting in a mass matrix \(\boldsymbol{B}\) and stiffness matrix \(\boldsymbol{S}\). See [Dyer et al., 2007, Zhang et al., 2007, Vallet and Levy, 2007].

Eigensystems of discrete Laplace–Beltrami operators

Desirable properties for discrete Laplacians include symmetry, locality, positive weights, linear precision, and convergence. No discretization fulfils all properties simultaneously; see [Wardetzky et al., 2007].

For symmetric Laplacians (TL, IE, COT), the eigenvalue problem reads

(6)\[\boldsymbol{L}\,\boldsymbol{x}_i = \lambda_i\,\boldsymbol{x}_i,\]

with eigenvalues \(\lambda_i \ge 0\) and orthonormal eigenvectors \(\boldsymbol{x}_i\). Collecting all eigenvectors yields the matrix

(7)\[\boldsymbol{X} = \bigl[ \boldsymbol{x}_1\;\boldsymbol{x}_2\;\cdots\;\boldsymbol{x}_n \bigr], \qquad \boldsymbol{X}^\top \boldsymbol{X} = \boldsymbol{I}.\]

These eigenvectors form a complete orthonormal basis and can be used for a spectral decomposition of signals on the mesh.

Spatial frequencies and eigenvalues

For a geometric FEM discretization, a SPHARA basis is obtained by solving

(8)\[\boldsymbol{S} \,\boldsymbol{\phi}_i = \tau_i \,\boldsymbol{B}\,\boldsymbol{\phi}_i,\]

with \(\tau_i \ge 0\) and \(\boldsymbol{\phi}_i\) orthonormal in the \(\boldsymbol{B}\)-weighted inner product. Collecting these eigenvectors yields

\[\boldsymbol{\Phi} = \bigl[ \boldsymbol{\phi}_1 \;\boldsymbol{\phi}_2 \;\dots\;\boldsymbol{\phi}_n \bigr].\]

Here, \(\tau_i\) approximate the eigenvalues of the continuous Laplace–Beltrami operator. They act as squared spatial angular frequencies, i.e.,

\[\omega_i = \sqrt{\tau_i}, \qquad f_i = \frac{\sqrt{\tau_i}}{2\pi}, \qquad \lambda_i = \frac{2\pi}{\sqrt{\tau_i}}.\]

Small \(\tau_i\) correspond to smooth spatial patterns (large wavelengths), large \(\tau_i\) to oscillatory patterns (small wavelengths). The columns of \(\boldsymbol{\Phi}\) are thus ordered from low to high spatial frequency.

Note

This quantitative interpretation of \(\tau_i\) in terms of spatial angular frequencies, frequencies and wavelengths only holds for the FEM discretization. For TL, IE and COT the eigenvalues still provide an ordering from smooth to oscillatory eigenmodes, but they do not correspond to geometric spatial frequencies. Accordingly, absolute cut-off frequencies for SPHARA filters are meaningful only for FEM-based bases (\(\boldsymbol{\Phi}\)).

SPHARA as a signal processing framework

Basis functions

For symmetric Laplacians the eigenvectors \(\boldsymbol{x}_i\) ((6)) are orthonormal. For FEM discretization, the \(\boldsymbol{B}\)-orthonormality of \(\boldsymbol{\phi}_i\) is ensured (see (8)). In both cases, the number of basis functions equals the number of vertices, ensuring completeness.

Analysis and synthesis

For symmetric Laplacians, coefficients for a signal \(\boldsymbol{f}\) are obtained via the standard inner product

(9)\[c_i = \langle \boldsymbol{f}, \boldsymbol{x}_i \rangle,\]

leading to

(10)\[\boldsymbol{c}^\top = \boldsymbol{f}^\top \boldsymbol{X}.\]

For FEM discretization, the \(\boldsymbol{B}\)-inner product is used:

(11)\[c_i = \boldsymbol{f}^\top \boldsymbol{B}\,\boldsymbol{\phi}_i,\]

and

(12)\[\boldsymbol{c}^\top = \boldsymbol{f}^\top \boldsymbol{B}\,\boldsymbol{\Phi}.\]

Synthesis uses

(13)\[\boldsymbol{f} = \sum_{i=1}^n c_i \,\boldsymbol{x}_i,\]

or

(14)\[\boldsymbol{f}^\top = \boldsymbol{c}^\top \boldsymbol{\Phi}^\top.\]

Spatial filtering using SPHARA

For filtering, a subset of basis functions is selected. Using the matrix \(\boldsymbol{X}\) of eigenvectors (symmetric Laplacians), and a selection matrix \(\boldsymbol{R}\), the filter matrix becomes

\[\boldsymbol{F} = \boldsymbol{X}\,\boldsymbol{R}\,(\boldsymbol{X}\,\boldsymbol{R})^\top.\]

For FEM discretization, the mass matrix \(\boldsymbol{B}\) enters:

\[\boldsymbol{F}_{\mathrm{FEM}} = \boldsymbol{B}\,\boldsymbol{X}\,\boldsymbol{R}\,(\boldsymbol{X}\,\boldsymbol{R})^\top.\]

Filtering a data matrix \(\boldsymbol{D}\) (time × sensors) yields

\[\tilde{\boldsymbol{D}} = \boldsymbol{D}\,\boldsymbol{F}.\]

Frequency-domain filter design

A more flexible approach replaces \(\boldsymbol{R}\) with a diagonal matrix \(\boldsymbol{H} = \mathrm{diag}(h_1,\dots,h_n)\), where \(h_i = H(\lambda_i)\) for a transfer function \(H\). Then

\[\boldsymbol{F}_H = \boldsymbol{X}\,\boldsymbol{H}\,\boldsymbol{X}^\top,\]

or for FEM,

\[\boldsymbol{F}_{H,\mathrm{FEM}} = \boldsymbol{B}\,\boldsymbol{X}\,\boldsymbol{H}\,\boldsymbol{X}^\top.\]

Typical choices include:

• Ideal filters (low-/high-/band-pass),

• Gaussian filters (smooth roll-off, amplitude specified at cut-off),

• Butterworth filters (maximally flat magnitude response).

SpharaPy implements these filters in the module spharapy.spectral_filters, where \(H(\lambda)\) is evaluated at the SPHARA eigenvalues and passed to spharapy.spharafilter.SpharaFilter.