Abstract

Multichannel time series often arise from a finite number of latent oscillatory sources observed through multiple sensors. A fundamental inverse problem is to determine how many such sources are active and, when a propagation or attenuation model is available, where the sources are located. Previous work by Bernadotte formulated a Hankel-embedding approach to the inverse problem of estimating the number of oscillatory sources in EEG and brain-computer-interface recordings. In that framework, the effective rank of a trajectory embedding provides an estimate of the number of active oscillators. Here we extend this framework from source-number estimation to source-coordinate recovery for attenuating oscillatory systems.

We formulate Multichannel Singular Spectrum Analysis (MSSA) as a block-Hankel method for extracting latent oscillatory modes from multichannel signals. Each channel is embedded into a Hankel trajectory matrix, and all channel-wise embeddings are concatenated into a joint block-Hankel matrix. Singular value decomposition of this matrix separates reproducible lagged multichannel dynamics from residual fluctuations; its signal subspace yields the oscillatory modes, with their frequencies and damping rates, and a spatial loading vector across sensors for every mode. We show that when a latent oscillatory source undergoes distance-dependent amplitude attenuation, the MSSA-derived spatial loading vector is an amplitude fingerprint of the source coordinate, and that without noise it is recovered exactly whenever the sources differ in frequency or damping. If sensor coordinates are known and the attenuation law is known and monotone, source localization reduces to matching this empirical amplitude fingerprint to theoretical attenuation profiles generated by candidate source locations.

The conditions of identifiability are made explicit. The rank of the block-Hankel matrix counts distinct oscillatory modes, so that sources of equal frequency and damping are counted once; they are detected through the phase structure of their common loading and separated by a scan of a two-dimensional subspace, which does not need the statistical independence that blind separation requires. For power-law attenuation, two positions have the same normalized fingerprint exactly when all sensors lie on a sphere with respect to which the positions are mutually inverse, or on a hyperplane in which they are mirror images; $d+2$ sensors in general position therefore determine a source in $\mathbb{R}^d$, and an array on a sphere is unambiguous inside it. For exponential attenuation the corresponding surfaces are hyperboloids. A first-order formula gives the accuracy of the recovered position, and simulations attain it with the number of sources estimated from the data. With a gain matrix in place of the scalar attenuation law the same procedure localizes current dipoles, for which a distance-only model fails. The resulting framework applies to neurophysiology, wearable sensing, industrial monitoring, geophysics, wireless systems, and financial markets, where latent oscillatory regimes may be observed through multiple attenuated channels or indicators.

Keywords: multichannel singular spectrum analysis; MSSA; Hankel embedding; inverse problem; oscillatory sources; attenuation; source localization; intrinsic oscillatory dimensionality; multichannel time series.

1. Introduction

Many biological, physical, technological, and economic systems are observed through multichannel time series. In such systems, the measured channels are often mixtures of a smaller number of latent dynamical sources. These sources may be oscillatory, quasi-periodic, transient, damped, amplified, or locally coherent over finite windows. A central problem is therefore inverse: from observed multichannel signals, infer the latent structure that generated them.

Two inverse problems are particularly important. The first is structural: how many latent oscillatory sources are active [1, 2]? The second is geometric: where are these sources located, either in physical space or in an abstract sensor, feature, or factor space?

Previous work by Bernadotte addressed the first inverse problem: estimating the number of oscillatory sources using Hankel embeddings and singular-spectrum analysis [1, 2]. In the paper Estimating the Number of Sources in EEG with Hankel Embedding for Brain-computer Interface [2], the Polyharmonic Brain Signal Model was introduced to justify trajectory embedding and low-rank approximation. Each EEG or MEG trace was mapped to a Hankel matrix, and the effective rank of the unfolding was interpreted as the number of oscillatory components in the ideal noise-free case. The paper proposed an algorithm for evaluating the number of sources, or oscillators, from the singular or eigenvalue spectrum of unfolded data. The approach rests on the theory of the unfolding of time series and its geometry on Grassmann manifolds [3, 4]. It has since been developed into a classifier of time series [5] and into a marker of brain state, the number of oscillatory modes in a window of EEG [6].

A related work, A Lightweight Hankel-Embedded Pipeline for Real-Time EEG Filtering and Classification [7], developed a practical Hankel-based pipeline for oscillator extraction, adaptive filtering, jPCA-based dynamics, and lightweight EEG classification. This showed that Hankel-derived oscillatory structures are not only theoretically meaningful but also useful for real-time multichannel signal processing.

The present article extends this line of work. We retain the first result — that Hankel/MSSA structure can estimate the number of active oscillatory sources — and add a second result: if oscillatory sources undergo distance-dependent amplitude attenuation, then the spatial loading vector recovered by MSSA becomes an amplitude fingerprint of the source coordinate. Under appropriate identifiability conditions, this fingerprint can be used to recover the coordinates of the sources.

Thus, MSSA is used not only as a decomposition or denoising method, but as a regularizing operator for inverse problems. It first estimates the number of active oscillatory sources and their spatial loadings. Then, using a known attenuation model and known sensor coordinates, these loadings are matched to theoretical source-to-sensor profiles to recover source coordinates.

The paper states both results precisely and tests them.

1. The rank of the block-Hankel matrix equals the number of distinct exponential modes of the recording, whatever the number of channels (Theorem 1). It counts sources when their frequencies or damping rates differ, and counts once those that coincide in both (Corollary 1).
2. The spatial loading of a mode is defined through the poles and amplitudes of the signal subspace, and not through individual singular vectors; without noise it is recovered exactly (Theorem 2). Its phase structure distinguishes a single source from two sources that share a mode (Proposition 1).
3. The injectivity of the normalized attenuation map, on which localization depends, is characterized geometrically for power-law and exponential attenuation (Theorems 3 and 4), and the accuracy of the recovered position is given to first order (Proposition 3).
4. The scalar attenuation law is the simplest case of a forward model with a known gain matrix, which is the form that EEG and MEG require (Section 5.4).
5. Every statement is checked numerically, the complete procedure is run with the number of sources estimated from the data, and the spatial step is compared with blind separation of independent components (Section 11).

1.1 Relation to earlier work

The framework continues our earlier work on the unfolding of time series: the estimation of the number of sources from the rank of a Hankel embedding [1, 2], the real-time extraction of oscillators [7], classification by the subspace of the embedding [5], and the oscillatory dimension of brain states [6]. Singular spectrum analysis, on which it draws, decomposes a series through the singular value decomposition of its trajectory matrix [9, 10, 11, 12]; the multichannel form stacks the trajectory matrices of several channels [13, 14]. Three further bodies of work meet in the second inverse problem.

Array signal processing. In [1, 2] the number of sources is read from the rank of the unfolding of the recording. In array processing it is estimated from the eigenvalues of a covariance matrix [15, 16], and the parameters of the sources from its signal subspace: MUSIC scans a manifold of candidate array responses for vectors close to the subspace [17], and ESPRIT reads the parameters from a shift invariance [18]. In that literature the subspace is spatial, formed from the covariance between sensors, and the manifold is parametrized by the direction of arrival. Here the subspace is temporal — the trajectory space of the recording — and space enters afterwards, through the loading of each temporal mode.

Localization from amplitudes. The oscillatory components extracted in [7] were used for filtering and classification, and their distribution over the sensors was not interpreted. Locating a source from the energy received at distributed sensors is known as energy-based localization: ratios of energies at pairs of sensors confine the source to hyperspheres, which are intersected by least squares or fitted by maximum likelihood [19, 20, 21]. Gain ratios of arrival have been combined with time differences of arrival [22], and positioning from received signal strength is standard in wireless sensor networks [23]. The Apollonius and hyperboloid systems of Section 9 are the noise-free form of these estimators. Such methods measure the energy of the received signal in a band, so several sources have to be fitted jointly [20].

Electromagnetic source imaging. The source counts of [1, 2] and the oscillatory dimension of [6] say how many generators are active in an EEG, not where they are. In EEG and MEG the forward model is a lead field and localization from a spatial pattern is dipole fitting: spatio-temporal dipole models and MUSIC scans of the lead field [24], and the fitting of an equivalent dipole to the scalp map of an independent component [25].

Blind separation in space. The modes counted in [1, 2] are temporal, and sources that share a mode have to be separated in space. Independent component analysis [26, 27] and second-order blind identification [28] do this from the statistics of the mixture, without a forward model; in EEG the maps of the components are then fitted with dipoles [25]. These methods require statistically independent sources, at least as many sensors as sources and, for independent component analysis, non-Gaussian signals. The scan of Section 6.2 uses the forward model in place of independence, and Section 11.9 compares the two.

What the present framework adds is the temporal separation that precedes the spatial fit. Each oscillatory mode, defined by its frequency and damping, receives its own spatial fingerprint; the number of modes is decided by a calibrated test; and sources that overlap in space and in frequency band are fitted one at a time, provided their frequencies or damping rates differ by a fraction of the Fourier resolution. The identifiability analysis of Section 7 applies to any method that localizes from relative amplitudes, whatever its front end.

2. Multichannel Singular Spectrum Analysis

The embedding is the unfolding of a time series into a Hankel matrix, on which the source-number method of [1, 2] and the pipeline of [7] are built. This section writes it for several channels at once.

2.1 Multichannel time series

Let

$$X = \{x_j(t)\}_{j=1,\,t=1}^{m,\,N} \in \mathbb{R}^{m \times N}$$

be a multichannel time series with $m$ observed channels or sensors and $N$ time samples. The goal of MSSA is to extract reproducible lagged structures from $X$, identify the effective number of structured oscillatory modes, and estimate their spatial loading patterns across the sensors.

2.2 Channel-wise Hankel embedding

For each channel $x_j(t)$, choose an embedding length $L$ and define $K = N - L + 1$. The Hankel trajectory matrix of channel $j$ is

$$H_L(x_j) = \begin{bmatrix} x_j(1) & x_j(2) & \cdots & x_j(K) \\ x_j(2) & x_j(3) & \cdots & x_j(K+1) \\ \vdots & \vdots & \ddots & \vdots \\ x_j(L) & x_j(L+1) & \cdots & x_j(N) \end{bmatrix} \in \mathbb{R}^{L \times K}.$$

Each column of $H_L(x_j)$ is a delayed window of the original signal. The channel is therefore transformed from a one-dimensional time series into a trajectory in an $L$-dimensional delay-coordinate space.

2.3 Block-Hankel MSSA matrix

The multichannel trajectory matrix is obtained by stacking all channel-wise Hankel matrices:

$$\mathcal{H}_L(X) = \begin{bmatrix} H_L(x_1) \\ H_L(x_2) \\ \vdots \\ H_L(x_m) \end{bmatrix} \in \mathbb{R}^{mL \times K}.$$

This construction preserves time alignment across sensors: the $k$-th column of $\mathcal{H}_L(X)$ contains delayed segments of all channels beginning at the same time.

MSSA differs from independent SSA because the decomposition is performed jointly across sensors. A weak oscillatory source may not be visible in one channel alone, but if it is coherent across several channels, it can appear as a significant component in the block-Hankel spectrum. Appendix B gives the size of the effect: an oscillation shared by eight channels at $-17.5$ dB is detected in 96 % of joint analyses and in 11 % of single-channel ones. The gain belongs to shared structure only; an oscillation confined to one channel is detected better by that channel alone.

2.4 Singular value decomposition

The MSSA matrix is decomposed as

$$\mathcal{H}_L(X) = U \Sigma V^{\top}, \qquad \Sigma = \operatorname{diag}(\sigma_1, \dots, \sigma_p), \quad \sigma_1 \ge \sigma_2 \ge \dots \ge 0 .$$

The leading singular components represent reproducible lagged multichannel patterns. The singular values measure the strength of these patterns. The lower spectral tail is typically associated with noise, unstructured fluctuations, and residual activity.

The decomposition provides three objects: singular values $\sigma_i$, used to estimate the number of structured modes; temporal components, encoded in $V$; and spatial-lagged loading patterns, encoded in $U$. For inverse source reconstruction, the important object is the spatial loading profile across sensors associated with a recovered latent mode.

Two remarks prepare what follows. The same decomposition is obtained from the eigenvectors of the lagged scatter $\mathcal{H}_L(X)\mathcal{H}_L(X)^{\top}$ of size $mL \times mL$, which is the faster route when $mL \ll K$. And a single singular vector is in general not the pattern of one source: when two sources contribute components of similar strength, the singular vectors are mixtures of both (Sections 5.2 and 11.4). What the data determine is the subspace spanned by the leading singular vectors.

3. Oscillatory source model

The model is the polyharmonic signal model of [1, 2]: every source is a finite combination of exponentially modulated harmonics. The classifier of [5] rests on the same finite-rank property.

3.1 Latent oscillatory sources

Assume that the observed multichannel signal is generated by $r$ latent oscillatory sources $s_1(t), \dots, s_r(t)$. Each source may be represented locally as an exponentially modulated oscillation,

$$s_k(t) = e^{\lambda_k t} \sin(\omega_k t + \phi_k),$$

or, more generally, as a finite combination of exponential atoms. The observed signal at sensor $j$ is

$$x_j(t) = \sum_{k=1}^{r} a_{jk}\, s_k(t) + \varepsilon_j(t), \tag{1}$$

where $a_{jk}$ is the contribution of source $k$ to sensor $j$, $s_k(t)$ is the temporal signal of source $k$, and $\varepsilon_j(t)$ is noise. In matrix form,

$$X = AS + E,$$

where $A \in \mathbb{R}^{m \times r}$ is the spatial mixing matrix and $S \in \mathbb{R}^{r \times N}$ contains the temporal source signals. The complex number $z_k = e^{\lambda_k + i\omega_k}$ is the pole of source $k$; it carries its frequency and its damping rate.

3.2 Finite-rank Hankel structure

The key mathematical fact is that finite sums of exponential atoms generate finite-rank Hankel matrices. Let

$$x(t) = \sum_{i=1}^{n} c_i z_i^{\,t},$$

where $z_i \in \mathbb{C} \setminus \{0\}$ are distinct and $c_i \neq 0$. Then the Hankel trajectory matrix $H_L(x)$ admits a Vandermonde factorization

$$H_L(x) = V_L(z)\, \operatorname{diag}(c_1 z_1, \dots, c_n z_n)\, V_K(z)^{\top},$$

in which the columns of $V_L(z)$ are the vectors $v_L(z_i) = (1, z_i, \dots, z_i^{L-1})^{\top}$. If $L \ge n$ and $K \ge n$, the Vandermonde matrices have full column rank, so $\operatorname{rank} H_L(x) = n$ [11]. A real oscillatory mode corresponds to a complex conjugate pair of exponential atoms. Therefore, $r$ real oscillatory modes generate a Hankel rank of $2r$ in the noise-free case. This justifies the use of SSA/MSSA for estimating the number of active oscillatory sources.

For several channels the statement takes the following form.

Theorem 1 (rank of the block-Hankel matrix). Let every channel be a combination of the same $n$ exponential atoms, $x_j(t) = \sum_{i=1}^{n} c_{ji} z_i^{\,t}$, with pairwise distinct $z_i \in \mathbb{C} \setminus \{0\}$ and non-zero loading vectors $c_i = (c_{1i}, \dots, c_{mi})^{\top}$. If $L \ge n$ and $K \ge n$, then $\operatorname{rank} \mathcal{H}_L(X) = n$ for every number of channels $m$.

Proof. The entry of $\mathcal{H}_L(X)$ in channel block $j$, row $l$ and column $k$ is $x_j(l + k - 1) = \sum_i c_{ji}\, z_i^{\,l-1} \cdot z_i \cdot z_i^{\,k-1}$, hence

$$\mathcal{H}_L(X) = \sum_{i=1}^{n} z_i\, \bigl(c_i \otimes v_L(z_i)\bigr)\, v_K(z_i)^{\top}.$$

The rank is at most $n$. The vectors $c_i \otimes v_L(z_i)$ are linearly independent: if $\sum_i \beta_i\, c_i \otimes v_L(z_i) = 0$, then $\sum_i \beta_i c_{ji}\, v_L(z_i) = 0$ in every channel $j$; Vandermonde vectors with distinct nodes are independent for $L \ge n$, so $\beta_i c_{ji} = 0$ for all $i$ and $j$, and since every $c_i$ has a non-zero entry, $\beta_i = 0$. The vectors $v_K(z_i)$ are independent for $K \ge n$. $\square$

Corollary 1 (modes and sources). In model (1) without noise, let the sources have poles $z_k = e^{\lambda_k + i\omega_k}$ with $0 < \omega_k < \pi$. If the $r$ poles are pairwise distinct, $\operatorname{rank} \mathcal{H}_L(X) = 2r$. If several sources share a pole, they contribute a single conjugate pair of atoms between them, and the rank is twice the number of distinct poles (provided that the joint loading of every pole is non-zero).

The second part matters in practice. Two sources of the same frequency and damping but of different phase are linearly independent as time series — $\operatorname{rank} S = r$ holds — and are nevertheless one mode of the trajectory matrix. Conversely, a source whose waveform is not a single modulated sinusoid, such as a rhythm with harmonics, contributes one mode for each of its components. The quantity recovered by the first inverse problem is therefore the number of distinct oscillatory modes. It is the number of sources when every source is one mode and no two sources share a mode.

In the language of automata the rank is a number of states: a series whose Hankel matrix has rank $n$ is computed by a weighted automaton with $n$ states and by none with fewer [29]. Theorem 1 then says that observing the same modes through more channels adds no states. State counts do not always behave so mildly. For the deterministic finite automata that recognize regular languages the number of states can grow exponentially with the length of the regular expression, the exponential explosion problem treated in [8].

4. First inverse problem: estimating the number of oscillatory sources

The first inverse problem is

$$X \longmapsto r .$$

That is, given the observed multichannel signal, estimate the number of active latent oscillatory sources. This problem was previously addressed in Bernadotte’s Hankel-embedding framework for EEG/BCI source-number estimation [1, 2]. There, the effective rank or $\varepsilon$-rank of the time-series unfolding was used to infer the number of active oscillatory sources. Since a real oscillatory component contributes two Hankel dimensions, the number of oscillators is estimated as one half of the effective rank.

In MSSA form, define the intrinsic oscillatory dimensionality [6] as

$$\widehat{d}_{\mathrm{osc}}(X; L, \tau) = \bigl|\{ i : \sigma_i(\mathcal{H}_L(X)) > \tau \}\bigr|, \tag{2}$$

where $\tau$ separates structured components from noise. If significant components form conjugate oscillatory pairs, the estimated number of oscillatory sources is

$$\widehat{r} = \tfrac{1}{2}\, \widehat{d}_{\mathrm{osc}} .$$

By Corollary 1 this is the number of distinct oscillatory modes, and the number of sources when their poles are distinct.

The threshold decides the result, and two ways of setting it are compared in Section 11.8.

Calibrated random-matrix threshold. Let $\ell_i = \sigma_i^2 / K$ and $\beta = mL/K < 1$. For white noise of variance $\sigma_\varepsilon^2$ and independent entries, the $\ell_i$ of the noise would follow the Marchenko–Pastur law of ratio $\beta$ and the largest of them would converge to $\sigma_\varepsilon^2 (1 + \sqrt{\beta})^2$ [30, 31]. The noise variance is estimated by matching the median of the $\ell_i$ to the median of that law, which a few signal components do not disturb, and a component is declared significant when

$$\ell_i > \widehat{\sigma}_\varepsilon^{2}\, (1 + \sqrt{\beta})^{2}\, (1 + \kappa).$$

The correction $\kappa$ is necessary. The largest noise eigenvalue fluctuates about the edge [32], and the entries of a trajectory matrix are not independent, every noise sample being repeated along an anti-diagonal [33]; with $\kappa = 0$ the rule reports components in pure white noise in 16–80 % of trajectory matrices (Appendix A). We take for $\kappa$ the 95 % quantile, over simulated white-noise recordings of the same $m$, $N$ and $L$, of the relative excess of the largest eigenvalue over the edge; this takes under a second and fixes the false-positive rate at 5 %.

Surrogate threshold. A more robust form uses surrogate spectra. Let $X_{\mathrm{sur}}^{(b)}$ be surrogate signals, $b = 1, \dots, B$, obtained by permuting the samples of every channel independently, which destroys all temporal structure and keeps the distribution of amplitudes. For each surrogate, compute singular values $\sigma_i^{(b)} = \sigma_i(\mathcal{H}_L(X_{\mathrm{sur}}^{(b)}))$. Then

$$\widehat{d}_{\mathrm{osc}}(X) = \bigl|\{ i : \sigma_i(\mathcal{H}_L(X)) > Q_{1-\eta}(\sigma_i^{\mathrm{surrogate}}) \}\bigr| ,$$

where $Q_{1-\eta}$ is the quantile of level $1 - \eta$ over the surrogates. This estimates the number of statistically significant lagged oscillatory dimensions.

The two thresholds behave differently. The surrogate spectrum is that of white noise with the total variance of the recording, signal included. It is therefore conservative: with sources of comparable strength it returns the correct number in all trials from $-5$ dB upward, but it masks a weak source beside a strong one — at an amplitude ratio of ten the weak source is missed in every trial, while the calibrated threshold finds both in 96 %. Since attenuation makes distant sources weak, we use the calibrated threshold as the default and the surrogate threshold as a check. Both assume a white noise floor. A coloured background is itself a set of large eigenvalues and has to be whitened, or tested against a null that carries its colour [34], before the count is interpreted (Appendix A).

Thus, the first inverse problem is solved structurally:

$$X \xrightarrow{\ \mathrm{MSSA}\ } \widehat{d}_{\mathrm{osc}} \longrightarrow \widehat{r}.$$

5. Second inverse problem: recovering coordinates of attenuating oscillatory sources

We now consider the second inverse problem:

$$X \longmapsto Q = (q_1, \dots, q_r),$$

where $q_k$ are the coordinates of the latent oscillatory sources. The earlier work [1, 2, 7] estimated the number of sources and extracted their time courses; their positions were outside its scope.

The key assumption is that the source signal undergoes distance-dependent amplitude attenuation before reaching each sensor. Therefore, different sensors record the same temporal mode with different amplitudes. These amplitudes encode the geometry of the source-sensor configuration. Attenuation here means the decay of amplitude with distance; decay in time is the damping $\lambda_k$ of Section 3.1 and is a property of the source.

5.1 Forward model with attenuation

Let sensors have known coordinates $p_1, \dots, p_m \in \mathbb{R}^{d}$. Let the $k$-th source have unknown coordinate $q_k \in \Omega \subset \mathbb{R}^{d}$. Assume that the contribution of source $k$ to sensor $j$ is

$$a_{jk} = \alpha_k\, g(\lVert p_j - q_k \rVert), \tag{3}$$

where $\alpha_k > 0$ is the intrinsic amplitude of source $k$, $g(\rho)$ is a known attenuation function of the distance $\rho$, and $g$ is typically strictly decreasing. Then the observed signal is

$$x_j(t) = \sum_{k=1}^{r} \alpha_k\, g(\lVert p_j - q_k \rVert)\, s_k(t) + \varepsilon_j(t),$$

or $X = A(Q) S + E$ with $A_{jk}(Q) = \alpha_k\, g(\lVert p_j - q_k \rVert)$. For a position $q$ define the attenuation profile

$$g_q = \bigl(g(\lVert p_1 - q \rVert), \dots, g(\lVert p_m - q \rVert)\bigr)^{\top}.$$

The $k$-th column of $A(Q)$ is $a_k = \alpha_k\, g_{q_k}$. This vector is the amplitude fingerprint of the source coordinate $q_k$.

5.2 MSSA spatial loading as an amplitude fingerprint

After MSSA separates the latent modes, each recovered oscillatory source is associated with a spatial loading vector $\widehat{a}_k = (\widehat{a}_{1k}, \dots, \widehat{a}_{mk})^{\top}$. Under the attenuation model,

$$\widehat{a}_k \approx \alpha_k\, g_{q_k}.$$

Thus, the vector of spatial loadings recovered by MSSA is an empirical amplitude fingerprint of the source coordinate. This is the central bridge between MSSA and the coordinate inverse problem.

The loading has to be defined with some care, because the left singular vectors do not provide it directly. By Theorem 1 the leading singular vectors span the same subspace as the vectors $c_i \otimes v_L(z_i)$, but each of them is a combination of these vectors, and when two sources have similar strength the combination mixes them: this is the failure of strong separability known in singular spectrum analysis [35]. The subspace, however, determines the poles, and the poles determine the loadings.

Let $U \in \mathbb{R}^{mL \times 2r}$ hold the leading left singular vectors. Let $U^{\uparrow}$ be $U$ without the last row of every channel block and $U^{\downarrow}$ be $U$ without the first row of every channel block.

Theorem 2 (recovery of the spatial loadings). Let $X = AS$ with sources $s_k(t) = e^{\lambda_k t} \sin(\omega_k t + \phi_k)$ whose poles $z_k = e^{\lambda_k + i\omega_k}$, $0 < \omega_k < \pi$, are pairwise distinct, and let no column of $A$ vanish. If $L \ge 2r + 1$ and $K \ge 2r$, then

(i) the eigenvalues of $\Phi = (U^{\uparrow})^{+} U^{\downarrow}$ are the $2r$ poles $z_k$ and $\bar{z}_k$;

(ii) the amplitudes $C \in \mathbb{C}^{m \times 2r}$ in $x_j(t) = \sum_i C_{ji} z_i^{\,t}$ are the unique solution of the linear system $X^{\top} = V_N(z)\, C^{\top}$;

(iii) the column $c_k$ of $C$ that belongs to $z_k$ equals $\tfrac{1}{2} e^{i\theta_k} a_k$ with $\theta_k = \phi_k - \pi/2$, so that

$$a_k = \pm\, 2\, \mathrm{Re}\bigl(e^{-i\theta_k} c_k\bigr), \qquad \theta_k = \tfrac{1}{2} \arg\bigl(c_k^{\top} c_k\bigr) \pmod{\pi}. \tag{4}$$

The loading vector of every source is thus recovered exactly up to sign, together with its waveform $s_k(t) = \pm\mathrm{Re}(e^{i\theta_k} z_k^{\,t})$. For a positive attenuation law the sign is fixed by requiring positive entries.

Proof. (i) By the factorization in the proof of Theorem 1, the columns of $W = [\, c_i \otimes v_L(z_i) \,]_{i=1}^{2r}$ span the column space of $\mathcal{H}_L(X)$, so $U = WT$ with $T$ invertible. Deleting the last row of every block gives $W^{\uparrow} = [\, c_i \otimes v_{L-1}(z_i) \,]$, and deleting the first gives $W^{\downarrow} = W^{\uparrow} D$ with $D = \operatorname{diag}(z_i)$. The matrix $W^{\uparrow}$ has full column rank because $L - 1 \ge 2r$, by the argument of Theorem 1. Hence $U^{\downarrow} = U^{\uparrow}\, T^{-1} D T$, the solution $\Phi = T^{-1} D T$ is unique, and its eigenvalues are the $z_i$. (ii) $V_N(z)$ has full column rank since $N \ge 2r$. (iii) $s_k(t) = \mathrm{Im}(e^{i\phi_k} z_k^{\,t}) = (e^{i\phi_k} z_k^{\,t} - e^{-i\phi_k} \bar{z}_k^{\,t})/(2i)$, so the coefficient of $z_k^{\,t}$ in channel $j$ is $a_{jk}\, e^{i\phi_k}/(2i)$; then $c_k^{\top} c_k = \tfrac{1}{4} e^{2i\theta_k} \lVert a_k \rVert^2$. $\square$

With noise, the same three steps are applied to the estimated subspace: this is the ESPRIT estimator [18, 36] applied to the block-Hankel matrix, followed by linear least squares. Once the poles are fixed the loading is linear in the data, with

$$\operatorname{Var}(\widehat{a}_{jk}) \approx 2\, \sigma_\varepsilon^{2}\, v_k, \qquad v_k = \bigl[(V_N^{\mathsf{H}} V_N)^{-1}\bigr]_{kk},$$

and $v_k \approx 1/N$ for an undamped mode, which is the bound for the amplitude of a sinusoid of unknown frequency [37].

Two properties of this definition are used below. It does not require the modes to be orthogonal or of different strength. And it yields a signed vector: for a positive attenuation law all entries have one sign, while for the forward models of Section 5.4 the signs carry information.

5.3 Sources that share a mode

Proposition 1. Let sources $k$ and $l$ share a pole, and let no other source have it. The loading of their common mode is $c = \tfrac{1}{2}\,(e^{i\theta_k} a_k + e^{i\theta_l} a_l)$. If $a_k$ and $a_l$ are not collinear and $\theta_k \not\equiv \theta_l \pmod{\pi}$, the real and imaginary parts of $c$ are linearly independent and span the plane of $a_k$ and $a_l$. If $\theta_k \equiv \theta_l \pmod{\pi}$, then $c$ is proportional to the real vector $a_k \pm a_l$, and the two sources cannot be distinguished from one source with that loading.

Proof. $[\mathrm{Re}\, c,\ \mathrm{Im}\, c] = \tfrac{1}{2}\,[a_k,\ a_l]\, R$ with $R = \begin{bmatrix} \cos\theta_k & \sin\theta_k \\ \cos\theta_l & \sin\theta_l \end{bmatrix}$ and $\det R = \sin(\theta_l - \theta_k)$. $\square$

A single source has an in-phase loading: after removal of the common phase by (4), the quadrature part $b = 2\, \mathrm{Im}(e^{-i\theta} c)$ vanishes. Under noise, the statistic

$$T = \frac{\lVert \widehat{b} \rVert^{2}}{2\, v\, \widehat{\sigma}_\varepsilon^{2}}$$

is distributed approximately as $\chi^2$ with $m - 1$ degrees of freedom when one source underlies the mode, the phase having been fitted. This gives a test: a mode whose $T$ exceeds the quantile of level $1 - \eta$ is not a single attenuated source. Two sources in one mode can be recovered (Section 6.2). With three or more, only a two-dimensional subspace of the span of their loadings is observed, and they cannot be recovered individually. A propagation delay that differs between sensors would also produce a quadrature part; the model considered here has none, like the quasi-static approximation of electrophysiology.

5.4 Forward models with a gain matrix

The scalar law (3) is the simplest forward model. More generally, let

$$a_k = G(q_k)\, \vartheta_k, \tag{5}$$

where $G(q) \in \mathbb{R}^{m \times n_\vartheta}$ is a known gain matrix and $\vartheta_k \in \mathbb{R}^{n_\vartheta}$ collects the unknown linear parameters of the source. Model (3) is the case $n_\vartheta = 1$, $G(q) = g_q$, $\vartheta = \alpha$. A current dipole in a volume conductor is the case $n_\vartheta = 3$: $\vartheta$ is the dipole moment, and the columns of $G(q)$, the lead field, are the potentials or magnetic fields produced at the sensors by unit dipoles along three axes [38, 39]. In a homogeneous conductor, for instance, the row of $G(q)$ for sensor $j$ is proportional to $(p_j - q)^{\top} / \lVert p_j - q \rVert^{3}$. Such a pattern changes sign across the array, depends on the orientation of the source, and need not be largest at the nearest sensor. A distance-only law cannot represent it (Section 11.7), and (5) is the form in which the framework applies to EEG and MEG.

6. Coordinate recovery by attenuation-profile matching

The oscillatory modes are obtained as in Sections 2–5, by the Hankel-embedding approach of [1, 2, 7]. This section adds the step that turns their loadings into coordinates.

6.1 One source per mode

If the intrinsic amplitude $\alpha_k$ is unknown, we compare normalized profiles. The coordinate of source $k$ is estimated by

$$\widehat{q}_k = \arg\max_{q \in \Omega} \frac{\langle \widehat{a}_k, g_q \rangle}{\lVert \widehat{a}_k \rVert\, \lVert g_q \rVert}. \tag{6}$$

Equivalently,

$$\widehat{q}_k = \arg\min_{q \in \Omega,\ \alpha > 0}\ \sum_{j=1}^{m} \bigl(\widehat{a}_{jk} - \alpha\, g(\lVert p_j - q \rVert)\bigr)^{2}.$$

For fixed $q$, the optimal amplitude is $\widehat{\alpha}(q) = \langle \widehat{a}_k, g_q \rangle / \lVert g_q \rVert^{2}$. Substituting this back gives a pure optimization problem over coordinates,

$$\widehat{q}_k = \arg\min_{q \in \Omega}\ \lVert \widehat{a}_k - \widehat{\alpha}(q)\, g_q \rVert^{2}.$$

This is a finite-dimensional geometric inverse problem. We solve it by a search over a grid in $\Omega$ followed by least-squares refinement.

The value of the optimum is a diagnostic. The residual fraction $1 - \max_q \cos^2(\widehat{a}_k, g_q)$ is the part of the loading that no position explains. A value well above the level expected from noise indicates a wrong attenuation law, a source outside $\Omega$, or two in-phase sources in one mode (Proposition 1).

6.2 Two sources in one mode

For a mode that fails the in-phase test, let $W \in \mathbb{R}^{m \times 2}$ be an orthonormal basis of the span of $\widehat{a}$ and $\widehat{b}$. By Proposition 1 the fingerprints of both sources lie in this plane, so the sources are among the minima of

$$\mathcal{J}(q) = 1 - \frac{\lVert W^{\top} g_q \rVert^{2}}{\lVert g_q \rVert^{2}} .$$

This is a scan of the kind used in MUSIC [17], in a subspace defined by a temporal mode instead of a spatial covariance. The two deepest minima are refined by the joint fit $\min_{q, q'} \lVert (I - \Pi_{[g_q,\, g_{q'}]})\, [\widehat{a}, \widehat{b}] \rVert_F$, where $\Pi$ denotes the orthogonal projector onto the span of the two profiles. Uniqueness requires that no third position have its profile in the plane, which by a count of conditions needs $m \ge d + 3$ sensors in general position.

The same scan applies when the two sources are not sinusoids but narrow-band signals that occupy one band. The loadings of all modes in the band then lie in the plane of the two fingerprints, and $W$ is taken from them jointly, or from the two leading spatial singular vectors of the band-passed recording (Section 11.9). No assumption about the statistical relation between the sources is involved, beyond their not being proportional.

6.3 Gain-matrix models

For model (5) the linear parameters are eliminated in the same way as $\alpha$:

$$\widehat{q}_k = \arg\max_{q \in \Omega} \frac{\lVert \Pi_{G(q)}\, \widehat{a}_k \rVert^{2}}{\lVert \widehat{a}_k \rVert^{2}}, \qquad \widehat{\vartheta}_k = G(\widehat{q}_k)^{+}\, \widehat{a}_k ,$$

where $\Pi_{G(q)}$ is the orthogonal projector onto the range of $G(q)$. This is the fit of an equivalent dipole to a spatial map [24, 25], applied to the map of an oscillatory mode.

7. Identifiability

In the first inverse problem, identifiability rests on the finite rank of the unfolding [1, 2] (Theorem 1). In the second, it is a property of the attenuation map.

Proposition 2 (reduction to the attenuation map). Let the observed multichannel signal be generated by (1) and (3), and let the loading vectors $a_k$ of the $r$ sources be recovered, as they are without noise under the conditions of Theorem 2. If the normalized attenuation map

$$q \mapsto \frac{g_q}{\lVert g_q \rVert}$$

is injective on $\Omega$, then each source coordinate $q_k$ is uniquely recoverable from its MSSA spatial loading vector $a_k$, up to amplitude scaling and permutation of source labels.

Proof. For source $k$, $a_k = \alpha_k g_{q_k}$, so $a_k / \lVert a_k \rVert = g_{q_k} / \lVert g_{q_k} \rVert$. If another point $q' \in \Omega$ gave the same normalized attenuation profile, then $g_{q'} / \lVert g_{q'} \rVert = g_{q_k} / \lVert g_{q_k} \rVert$ and, by injectivity, $q' = q_k$. Since the ordering and scaling of separated sources are arbitrary, recovery is unique up to permutation of source labels and amplitude scaling. $\square$

The proposition reduces identifiability to a property of the array and of the attenuation law. For the two laws of Section 8.2 this property can be stated geometrically.

Theorem 3 (power-law attenuation). Let $g(\rho) = \rho^{-\gamma}$ with $\gamma > 0$, and let $q \neq q'$ be two points distinct from the sensors. The profiles $g_q$ and $g_{q'}$ are proportional if and only if all sensors lie

(a) on a sphere with respect to which $q$ and $q'$ are mutually inverse, or

(b) on the hyperplane with respect to which $q$ and $q'$ are mirror images.

Moreover, the differential of the normalized attenuation map is singular at $q$ if and only if all sensors lie on a sphere or a hyperplane that passes through $q$.

Proof. Proportionality means $\lVert p_j - q' \rVert = c\, \lVert p_j - q \rVert$ for all $j$ with one constant $c > 0$. For $c = 1$ the locus $\{x : \lVert x - q' \rVert = \lVert x - q \rVert\}$ is the bisecting hyperplane. For $c \neq 1$ it is the sphere of Apollonius, with centre $o = (q' - c^2 q)/(1 - c^2)$ and radius $R = c\, \lVert q - q' \rVert / |1 - c^2|$. Since $q - o = (q - q')/(1 - c^2)$ and $q' - o = c^2 (q - q')/(1 - c^2)$, the points $q$ and $q'$ lie on one ray from $o$ and $\lVert q - o \rVert\, \lVert q' - o \rVert = R^2$: they are inverse with respect to the sphere [40]. Conversely, if $q'$ is the inverse of $q$ in a sphere of centre $o$ and radius $R$, then $\lVert x - q' \rVert / \lVert x - q \rVert = R / \lVert q - o \rVert$ for every $x$ on the sphere. For the last statement, the differential annihilates a direction $v$ exactly when $(\partial g_q / \partial q)\, v$ is proportional to $g_q$, that is, when $(q - p_j)^{\top} v = c\, \lVert q - p_j \rVert^{2}$ for all $j$. For $c \neq 0$ this places the sensors on the sphere through $q$ with centre $q - v/(2c)$; for $c = 0$, on the hyperplane through $q$ orthogonal to $v$. $\square$

Corollary 2. For power-law attenuation and unknown source amplitude:

(i) if the sensors lie neither on a common sphere nor on a common hyperplane — which requires $m \ge d + 2$ and holds for $d + 2$ sensors in general position — the normalized attenuation map is injective on $\mathbb{R}^{d}$ without the sensors, and its differential is nowhere singular;

(ii) $m = d + 1$ affinely independent sensors always lie on their circumscribed sphere, and every source other than the centre of that sphere has exactly one twin, its inverse in the sphere;

(iii) sensors on a sphere $\Sigma$ that do not lie in a hyperplane confuse a point only with its inverse in $\Sigma$; the map is injective inside $\Sigma$ and outside it, so restricting $\Omega$ to the interior, when the sources are known to lie inside the array, removes the ambiguity;

(iv) sensors that span a hyperplane and do not lie on a sphere confuse a point only with its mirror image, and restricting $\Omega$ to one side removes the ambiguity. The mirror ambiguity of a planar array holds for every attenuation law, since all distances coincide.

Theorem 4 (exponential attenuation). Let $g(\rho) = e^{-\mu\rho}$ with $\mu > 0$, and $q \neq q'$. The profiles $g_q$ and $g_{q'}$ are proportional if and only if the difference of distances $\lVert p_j - q' \rVert - \lVert p_j - q \rVert$ is the same number $\Delta$ for all sensors; that is, if and only if all sensors lie on the bisecting hyperplane of $q$ and $q'$ ($\Delta = 0$), on one sheet of a hyperboloid of revolution with foci $q$ and $q'$ ($0 < |\Delta| < \lVert q - q' \rVert$), or on the line through $q$ and $q'$ beyond one of the two points ($|\Delta| = \lVert q - q' \rVert$).

Proof. $e^{-\mu \rho_j'} = c\, e^{-\mu \rho_j}$ for all $j$ is equivalent to $\rho_j' - \rho_j = -\mu^{-1} \log c$, and $|\Delta| \le \lVert q - q' \rVert$ by the triangle inequality. $\square$

The situation is that of hyperbolic positioning from time differences of arrival [41, 42], with logarithms of amplitude ratios in place of range differences. A second solution is described by $q'$ and $\Delta$, which are $d + 1$ numbers, against $m$ conditions: $d + 1$ sensors may admit one, and $d + 2$ sensors in general position do not. For a general strictly monotone law the same count applies. The normalized fingerprint has $m - 1$ degrees of freedom and the position has $d$, so $m \ge d + 1$ sensors are necessary for a locally unique solution, and one sensor more is what excludes isolated second solutions.

Proposition 3 (local identifiability and first-order accuracy). Let $J(q) = \partial g_q / \partial q \in \mathbb{R}^{m \times d}$, with rows $g'(\rho_j)\,(q - p_j)^{\top} / \rho_j$, $\rho_j = \lVert p_j - q \rVert$, and let $\Pi_q^{\perp} = I - g_q g_q^{\top} / \lVert g_q \rVert^{2}$.

(a) If $\Pi_q^{\perp} J(q)$ has rank $d$, the normalized attenuation map is locally injective at $q$. This requires $m \ge d + 1$.

(b) If $\widehat{a} = \alpha g_q + e$ with $\mathbb{E}\, e = 0$ and $\operatorname{Cov} e = \sigma_a^{2} I$, the estimate (6) satisfies, to first order in $e$,

$$\operatorname{Cov}(\widehat{q}) = \frac{\sigma_a^{2}}{\alpha^{2}}\, \bigl(J^{\top} \Pi_q^{\perp} J\bigr)^{-1}. \tag{7}$$

Proof. (a) The differential of $q \mapsto g_q / \lVert g_q \rVert$ is $\Pi_q^{\perp} J / \lVert g_q \rVert$. (b) The Jacobian of the residual $\widehat{a} - \alpha g_q$ with respect to $(q, \alpha)$ is $-[\alpha J,\ g_q]$. The covariance of the least-squares estimate is $\sigma_a^{2}$ times the inverse of its Gram matrix, whose block for $q$ is $(\alpha^{2} J^{\top} \Pi_q^{\perp} J)^{-1}$ by the Schur complement. $\square$

For an undamped mode in white noise, $\sigma_a^{2} = 2\sigma_\varepsilon^{2}/N$ (Section 5.2), and the root-mean-square position error is

$$\mathrm{RMS}(\widehat{q}) = \frac{\sigma_\varepsilon}{\alpha} \sqrt{\frac{2}{N}}\ \mathcal{D}(q), \qquad \mathcal{D}(q) = \Bigl[\operatorname{tr}\bigl(J^{\top} \Pi_q^{\perp} J\bigr)^{-1}\Bigr]^{1/2}.$$

The factor $\mathcal{D}(q)$ depends only on the array and on the law. It plays the part of the geometric dilution of precision of positioning systems [42], and it is infinite exactly where the differential of the normalized map is singular. Figure 2c maps it for a grid of nine sensors.

8. Conditions required for coordinate recovery

The inverse problem is solvable under the following conditions. They add to the condition under which the number of sources is read from the rank of the unfolding [1, 2]: that the sources be finite combinations of exponential atoms, observed against a noise floor from which they can be separated.

8.1 Known sensor coordinates

The sensor coordinates $p_1, \dots, p_m$ must be known. Without sensor geometry, the attenuation profile cannot be mapped to physical or abstract coordinates.

8.2 Known or calibrated attenuation law

The function $g(\rho)$ must be known or estimated from calibration data. Typical examples include

$$g(\rho) = \frac{1}{\rho^{\gamma}}, \quad \gamma > 0, \qquad \text{or} \qquad g(\rho) = e^{-\mu\rho}, \quad \mu > 0 .$$

More complex attenuation laws are possible, including anisotropic or medium-dependent models, and the gain-matrix models of Section 5.4. A parameter of the law can also be estimated from the fingerprint itself, together with the position, when the number of sensors allows one more unknown: Section 11.6 does this for the exponent $\gamma$.

8.3 Monotonicity of attenuation

The attenuation law must be monotone over the relevant distance range: $g'(\rho) < 0$. This ensures that different distances produce different amplitudes. If $g$ is not injective over the relevant range, distance cannot be inferred from amplitude. Monotonicity concerns one sensor. It does not by itself make the normalized attenuation map injective, which is a property of the law and of the array together.

8.4 Non-degenerate sensor geometry

The sensor configuration must not create unresolved symmetries. If the source amplitude were known, each sensor would give a distance, and $d + 1$ affinely independent sensors would determine the position; in $\mathbb{R}^{3}$ this means four non-coplanar sensors. When the source amplitude is unknown and only relative amplitudes are available — the case considered here — one more sensor is required. For power-law attenuation, $d + 1$ sensors always leave two solutions, and $d + 2$ sensors that lie neither on a sphere nor on a hyperplane give a unique one (Corollary 2); arrays on a sphere or in a plane are unambiguous once $\Omega$ is restricted to the interior of the sphere or to one side of the plane. For exponential attenuation the surfaces are hyperboloids: $d + 1$ sensors may leave a second solution, and $d + 2$ sensors in general position do not (Theorem 4). In practice, one uses $m > d + 2$ sensors and solves an overdetermined least-squares problem, whose accuracy is given by (7).

8.5 Injectivity of the normalized attenuation map

The essential condition is

$$q \neq q' \quad \Longrightarrow \quad \frac{g_q}{\lVert g_q \rVert} \neq \frac{g_{q'}}{\lVert g_{q'} \rVert}.$$

Equivalently, different coordinates must generate different normalized amplitude fingerprints. If two points produce proportional attenuation profiles, $g_q = c\, g_{q'}$, then they cannot be distinguished when source amplitude is unknown. Section 7 states when this happens.

8.6 Temporal separability of sources

MSSA must be able to separate the temporal source signals. The condition for this is that the poles of the sources be pairwise distinct: the sources must differ in frequency or in damping (Theorem 2). Linear independence of the source signals, $\operatorname{rank} S = r$, is not sufficient: two sources of one frequency and damping with different phases are independent and still form a single mode. If two sources have identical temporal dynamics up to a factor, $s_1(t) = c\, s_2(t)$, then they appear as a single effective source,

$$a_1 s_1(t) + a_2 s_2(t) = (c\, a_1 + a_2)\, s_2(t),$$

and their individual coordinates cannot be recovered. Between these cases lie two sources that share a pole and differ in phase; they are detected and localized through the plane spanned by their loadings (Sections 5.3 and 6.2).

8.7 Spatial distinguishability of sources

Sources with well-separated poles need not have distinct fingerprints: each mode is fitted on its own, and two such sources at one position are both localized there (Section 11.4). Distinct fingerprints,

$$\frac{g_{q_k}}{\lVert g_{q_k} \rVert} \neq \frac{g_{q_l}}{\lVert g_{q_l} \rVert}, \qquad k \neq l,$$

are required in two situations: for sources that share a mode, and for sources whose frequencies differ by less than about half the Fourier resolution, which are separated only because their loadings differ. For stable recovery in these situations, a stronger separation condition is required, $\angle(g_{q_k}, g_{q_l}) > \delta > 0$. If two such sources are too close, their spatial loadings become nearly collinear and the inverse problem becomes ill-conditioned.

8.8 Sufficient signal-to-noise ratio

In real data, $X = A(Q)S + E$. The recovery is stable only if the structured components are separated from noise in the MSSA spectrum. A practical condition is the existence of a singular-value gap,

$$\sigma_{\widehat{d}_{\mathrm{osc}}} \gg \sigma_{\widehat{d}_{\mathrm{osc}} + 1},$$

which the test of Section 4 turns into a decision with a stated error rate. More formally, perturbation results imply that the recovered subspace is stable when the noise norm is small relative to the spectral gap. Once the modes are detected, the position error is proportional to the noise level and to the dilution factor of the array, by (7).

9. Special cases of attenuation-based localization

In both cases below the loadings $a_{jk}$ are those of the oscillatory modes of Section 5.2, extracted by the Hankel-embedding approach of [1, 2, 7].

9.1 Power-law attenuation

Suppose $g(\rho) = \rho^{-\gamma}$. Then $a_{jk} = \alpha_k / \lVert p_j - q_k \rVert^{\gamma}$. Taking ratios eliminates $\alpha_k$:

$$\frac{a_{jk}}{a_{1k}} = \left( \frac{\lVert p_1 - q_k \rVert}{\lVert p_j - q_k \rVert} \right)^{\gamma}.$$

Define $\rho_j = (a_{1k} / a_{jk})^{1/\gamma}$. Then $\lVert p_j - q_k \rVert = \rho_j \lVert p_1 - q_k \rVert$. This gives a system of Apollonius sphere equations:

$$\lVert p_j - q_k \rVert^{2} = \rho_j^{2}\, \lVert p_1 - q_k \rVert^{2}.$$

The source coordinate is recovered from the intersection of these surfaces, or in the noisy case by least squares.

With $p_1$ taken as the origin and $w = \lVert q_k \rVert^{2}$, the equations are linear in $(q_k, w)$:

$$\lVert p_j \rVert^{2} - 2\, p_j^{\top} q_k + (1 - \rho_j^{2})\, w = 0, \qquad j = 2, \dots, m .$$

For $m \ge d + 2$ sensors in general position the $m - 1$ equations determine the $d + 1$ unknowns, and the solution satisfies $w = \lVert q_k \rVert^{2}$ by itself. This closed form is exact without noise and can start the least-squares refinement. For $m = d + 1$ the solutions form a line, and the constraint $w = \lVert q_k \rVert^{2}$ selects the two points of Corollary 2(ii).

9.2 Exponential attenuation

Suppose $g(\rho) = e^{-\mu\rho}$. Then $a_{jk} = \alpha_k e^{-\mu \lVert p_j - q_k \rVert}$. Taking logarithmic ratios,

$$\log \frac{a_{jk}}{a_{1k}} = -\mu\, \bigl( \lVert p_j - q_k \rVert - \lVert p_1 - q_k \rVert \bigr).$$

Thus,

$$\lVert p_j - q_k \rVert - \lVert p_1 - q_k \rVert = -\frac{1}{\mu} \log \frac{a_{jk}}{a_{1k}} =: \Delta_j .$$

This is a system of hyperboloid equations, analogous to time-difference localization, but based on amplitude attenuation rather than arrival time. With $p_1$ as the origin and $R = \lVert q_k \rVert$, squaring $\lVert p_j - q_k \rVert = R + \Delta_j$ gives the linear system

$$-2\, p_j^{\top} q_k - 2\, \Delta_j R = \Delta_j^{2} - \lVert p_j \rVert^{2}, \qquad j = 2, \dots, m,$$

the closed form used in hyperbolic positioning [41].

10. MSSA-assisted inverse localization algorithm

Steps 1–3 are the source-number algorithm of [1, 2] in multichannel form, and the extraction of oscillatory modes in Step 4 continues the pipeline of [7]. Steps 5–7 are new.

Input: multichannel signal $X \in \mathbb{R}^{m \times N}$; sensor coordinates $p_1, \dots, p_m$; attenuation function $g(\rho)$ or gain matrix $G(q)$; embedding length $L$; level $\eta$ of the tests; candidate source domain $\Omega$.

Step 1. Construct the block-Hankel MSSA matrix $\mathcal{H}_L(X)$ of Section 2.3.

Step 2. Compute the decomposition $\mathcal{H}_L(X) = U \Sigma V^{\top}$, through the SVD or through the eigenvectors of the lagged scatter.

Step 3. Estimate the number of active oscillatory dimensions $\widehat{d}_{\mathrm{osc}}$ with the calibrated threshold of Section 4, or with the surrogate threshold.

Step 4. Estimate the poles. With $U$ the first $\widehat{d}_{\mathrm{osc}}$ left singular vectors, compute $\Phi = (U^{\uparrow})^{+} U^{\downarrow}$ and its eigenvalues $\widehat{z}_i$. The modes are the poles with positive frequency; their number is $\widehat{r}$, equal to $\tfrac{1}{2} \widehat{d}_{\mathrm{osc}}$ if components are paired.

Step 5. Extract MSSA spatial loadings. Solve $X^{\top} = V_N(\widehat{z})\, C^{\top}$ by least squares. For each mode $k$ compute the common phase and the loading by (4),

$$\widehat{a}_k = (\widehat{a}_{1k}, \dots, \widehat{a}_{mk})^{\top},$$

its quadrature part $\widehat{b}_k$, and the statistic $T_k$ of Section 5.3.

Step 6. Match the loading to attenuation profiles. For a mode that passes the in-phase test, compute $g_q$ for candidate coordinates $q \in \Omega$ and estimate

$$\widehat{q}_k = \arg\max_{q \in \Omega} \frac{\langle \widehat{a}_k, g_q \rangle}{\lVert \widehat{a}_k \rVert\, \lVert g_q \rVert}.$$

For a mode that fails it, scan the plane of $\widehat{a}_k$ and $\widehat{b}_k$ for two sources (Section 6.2). For a gain-matrix model use the projector of Section 6.3.

Step 7. Estimate amplitude and accuracy. Once $\widehat{q}_k$ is obtained,

$$\widehat{\alpha}_k = \frac{\langle \widehat{a}_k, g_{\widehat{q}_k} \rangle}{\lVert g_{\widehat{q}_k} \rVert^{2}},$$

the residual fraction of Section 6.1 measures the quality of the fit, and (7) with $\sigma_a^{2} = 2 \widehat{\sigma}_\varepsilon^{2} v_k$ gives the covariance of the position.

Output: $\{\widehat{q}_k, \widehat{\alpha}_k, \widehat{s}_k(t)\}_{k=1}^{\widehat{r}}$ with $\widehat{s}_k(t) = \mathrm{Re}(e^{i\widehat{\theta}_k} \widehat{z}_k^{\,t})$, together with the frequency and damping of every source and the diagnostics of Steps 5 and 7.

11. Numerical experiments

The number of sources was evaluated on recordings of a non-invasive brain-computer interface in [1, 2], and the Hankel-embedded pipeline on EEG classification in [7]. The positions of sources are not known in such recordings, so the experiments of this section are simulations.

11.1 Setting

Unless stated otherwise, nine sensors form a $3 \times 3$ grid on the unit square ($d = 2$), so that lengths are in units of the side of the array. The sampling rate is 100 Hz, $N = 500$ (five seconds) and $L = 40$. Sources are undamped sinusoids of random phase at random positions in $[0.1, 0.9]^2$, at least 0.12 from every sensor. The attenuation law is $g(\rho) = \rho^{-1}$, or $e^{-3\rho}$ where indicated. Sensor noise is white, Gaussian and of equal variance $\sigma_\varepsilon^{2}$ in all channels. The signal-to-noise ratio is defined at the reference distance $\rho_0 = 0.5$,

$$\mathrm{SNR}_0 = \frac{\alpha^{2}\, g(\rho_0)^{2}}{2\, \sigma_\varepsilon^{2}},$$

the ratio that a source of amplitude $\alpha$ has at a sensor half an array side away. Position errors above 0.1 are called gross. Each point is the result of 300 trials unless noted. In Sections 11.3–11.7 the rank is given, so that the localization step is tested separately; in Section 11.8 it is estimated.

11.2 The exact statements

Without noise the statements of Sections 3–9 hold to machine precision.