# The Maximally-Coherent Reference technique and its application to sound source extraction without synchronous measurements

Muhammad N. Albezzawy, Jérôme Antoni, Quentin Leclère

► To cite this version:

Muhammad N. Albezzawy, Jérôme Antoni, Quentin Leclère. The Maximally-Coherent Reference technique and its application to sound source extraction without synchronous measurements. Journal of Sound and Vibration, 2024, 599, 10.1016/j.jsv.2024.118896. hal-04886260

HAL Id: hal-04886260

https://hal.science/hal-04886260v1

Submitted on 14 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Maximally-Coherent Reference technique and its application to sound source extraction without synchronous measurements

Muhammad N. Albezzawy<sup>ID*</sup>, Jérôme Antoni, Quentin Leclère<sup>ID</sup>

INSA Lyon, LVA, UR677, Villeurbanne, 69621, France

<figure><img src="image_3.png" /></figure>

Check for
updates

## ARTICLE INFO

**Keywords:**

Source extraction

Resynchronization

Reference techniques

Coherence methods

Microphone array

Electric engines

## ABSTRACT

In source identification, it is often necessary to perform source extraction, and in cases involving sequential measurements, to also perform resynchronization. Coherence techniques, which are based on the use of references (i.e., fixed sensors), are widely used to solve these two equivalent problems. However, when the number of references surpasses the number of sources, the cross-spectral matrix becomes ill-conditioned, invalidating the popular least squares (LS) solution. Although the truncated singular value decomposition (TSVD) was successfully applied in the literature to solve this problem, its validity is limited to the case of scalar noise on the references. It is also difficult to apply, when the singular values are gradually decreasing. This paper proposes a solution based on a set of virtual references that is maximally correlated with the measurements, named the Maximally-Coherent Reference (MCR) technique, accompanied with a technique for estimating the number of sources. The method is validated using both numerical and physical laboratory experiments, and by using real acoustical data from an e-motor. It is shown to return better results than LS and TSVD when employed for the same purpose.

## 1. Introduction

In applications involving the identification and ranking of sources originating from different physical origins through remote sensing (e.g. in acoustic source reconstruction [1]), the extraction/separation of mutually incoherent sources is crucial [2–4]. In many scenarios, the sources tend to overlap in both the space and frequency domains, which poses a significant challenge for their extraction. Another similar scenario is sequential measurements [5–7], where the loss of the phase requires some sort of resynchronization of all the measurements [1]. The equivalence between the two problems (source extraction and resynchronization) lies in the fact that the two problems boils down to extracting or restoring the pure contribution of the sources to the sensors in the presence of additive noise. As a solution for the two problems, the reference techniques (also known as coherence techniques) [8] were introduced in the 70’s [9]. They use so-called reference sensors to extract the coherent spectra (due to sources) from noisy measurements, by restricting the rank of the cross-power spectral matrix (CSM) of the measurements to the number of references [10], through the least squares (LS) solution between the estimated source contributions and the measurements.

However this solution came with restrictive assumptions on the references, which are hard to meet in practice: (1) the purity of references (i.e. high signal to noise ratio (SNR)), (2) each reference is perfectly correlated with one source and (3) totally uncorrelated with other sources, (4) the number of references is not less than the number of sources. However, in practice, when the number of references exceeds the number of sources, the references are certainly inter-correlated. In addition, noise always pollutes them,

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: muhammad.albezzawy@insa-lyon.fr (M.N. Albezzawy).</footnote>
leading to a full-rank CSM, with rank equal to the number of reference sensors and not to the number of sources. Hence, a rank reduction strategy needs to be implemented on the CSM of the references.

Two main rank reduction techniques were proposed in the literature: the partial coherence technique [11–13], and the virtual coherence technique [14]. The first one requires an ordering of the references, which is essentially problematic, as a model with arbitrary N inputs will theoretically require analyzing N! possible orderings. To avoid this, one has to rely on a priori knowledge on the cause-and-effect relations between every pair of inputs, which is practically difficult [9]. Besides, in the case of systems with large number of references, the probability of existence of perfectly correlated (linearly dependent) inputs arises, which would lead to numerical problems due to singularity [14,15]. Eventually, if there is noise on the references, the technique is conceptually no longer valid. On the other hand, the virtual coherence technique, utilizing principal component analysis (PCA), through a truncated singular value decomposition (TSVD) of the CSM of references, is valid only when dealing with spatially white and decorrelated noise (i.e. having its covariance matrix proportional to identity), hereafter coined as “scalar noise”. In cases where the noise on the references deviates from this characteristic, the resulting virtual spectra (i.e. the principal spectra obtained from TSVD) get biased. Practically, this means that the technique fails in the cases of correlated noise on the references, and/or noise of different variances (e.g. as encountered with reference sensors of different types) [16]. Recent attempts were reported in the literature to tackle this problem, by averaging the reference auto-spectra over several scans [17,18], or through regularizing the inverse of the CSM [19,20].

Another strong challenge, which is nearly not discussed in the literature, facing the rank reduction techniques as a whole, arises when the noise levels on the references result in a well-conditioned CSM (i.e. a full rank matrix with gradually decreasing singular values). In this situation, there is no clear threshold that divides principal spectra of the sources from those of the noise.

This paper presents the Maximally-Coherent Reference (MCR) technique, recently introduced by the authors [16], as a solution for the source extraction problem, by finding a set of virtual spectra, which are linear combinations of the measured references, that are maximally correlated with the measurements and so with the actual sources. MCR, once tested by simulations [16], showed its immunity against the limitations of both the LS and the TSVD solutions, which makes it a robust basis for source extraction. This work also presents a novel technique for selecting the number of sources, using a new likelihood ratio (LR) formulation based on testing against a saturated model. A bootstrap strategy is adopted to make the technique applicable to the cases of short records. The validity of the method, and its supremacy over LS and TSVD, are proved by numerical simulations and by real acoustical experiments. The method is also applied for solving the resynchronization problem of sequential measurements on numerical simulations, on real acoustical experiments, and on data from an e-motor industrial test bench.

The paper is structured as follows. Mathematical formulation of the problem is detailed in Section 2. Numerical validation is undertaken in Section 3. Validation by the laboratory experiment is presented in Section 4. Application to the e-motor data is reported in Section 5. Finally, the work is concluded in Section 6.

## 2. Mathematical formulation of the problem

### 2.1. Problem statement

Let $Y \in \mathbb{C}^{M \times I}$ denote a set of complex-valued measurements (e.g. Fourier coefficients) returned by M sensors and recorded for I independent snapshots, where $I \ge M$. The data Y are supposed to be produced by S sources, say $S \in \mathbb{C}^{S \times I}$, whose contributions are noted $X \in \mathbb{C}^{M \times I}$, such that

$$X = HS, \qquad (1)$$

for some linear but unknown operator $H \in \mathbb{C}^{M \times S}$, and corrupted by additive disturbances $N \in \mathbb{C}^{M \times I}$ uncorrelated with X, such that

$$Y = X + N. \qquad (2)$$

$H$ is assumed of full rank (i.e. rank($H$) = $S$) without loss of generality (otherwise there would be the effect of less than $S$ sources in $X$). It is also assumed that the number of output sensors $M$ exceeds or is equal to the number of sources, i.e. $M \ge S$.

The data Y comes with a set of R references, $R \in \mathbb{C}^{R \times I}$, which are supposed to be perfectly correlated with the S sources, in the sense that their exists a linear but unknown operator $L \in \mathbb{C}^{R \times S}$ such that

$$R = LS, \qquad (3)$$

where $L$ is assumed of full rank (i.e. rank($L$) = $S$) without loss of generality (otherwise there would be the effect of less than $S$ sources in $R$). In practice, a small amount of additive noise may be present in the references, such that

$$R = LS + \epsilon E, \qquad (4)$$

with $\epsilon \ge 0$ a deterministic scaling factor and $E$ a matrix of complex random variables such that the noise covariance matrix, $\Sigma_E = E[\epsilon \epsilon^H]$ is identical for any column $e$ of $E$, but is not necessarily proportional to the identity.

It is further assumed that the number of references exceeds or is equal to the number of sources, i.e. $R \ge S$.

The aim is to predict the contributions $X$ from the references $R$, i.e. to find an operator $G \in \mathbb{C}^{M \times R}$ such that

$$\hat{X} = GR \qquad (5)$$

is an estimate of $X$.
## 2.2. The proposed Maximally-Coherent Reference technique

A simple solution to the above problem that is often used in practice is provided by the least squares estimate

$$G = \operatorname{Arg min}_{A} \|Y - AR\|^2 = S_{YR}S_{RR}^{-1} \qquad (6)$$

with $S_{RR} = RR^H/I$ the CSM of the references and $S_{YR} = YR^H/I$ the CSM between the references and the measurements. Unfortunately, $S_{RR}$ might be badly conditioned when $R > S$ and $\epsilon$ is small; it tends to be singular when $\epsilon \to 0$ and $R > S$.

A popular way to deal with the aforementioned situation is to replace the set of R references R by their first S principal components, as returned by TSVD. However, the so-selected S components do not need to be the best ones to predict the contributions X, i.e., the most correlated with the sources S.

In light of this discussion, a better choice is to select S linear combinations of references R that best predict the data Y. This amounts to finding a vector $b_1$ such that $z_1 = b_1^H R$ is maximally correlated with Y, then to remove the contribution of $z_1$ from the data and to find the next linear combination of references defined by vector $b_2$ such that $z_2 = b_2^H R$ is maximally correlated with the residual, etc. The proposed algorithm reads

### Algorithm 1

* Set $Y_0 = Y$

* FOR $s = 1$ to $S$

- Find $b_s = \operatorname{argmax}_b \operatorname{corr}(Y_{s-1}, b^H R)$

- Set $z_s = b_s^H R$

- Find $c_s = \operatorname{argmin}_c \|Y_{s-1} - cz_s\|^2$

- Set $Y_s = Y_{s-1} - c_s z_s$

* $s \leftarrow s + 1$

* END if $s > S$

Let denote $B = [b_1 \cdots b_S]$ and $C = [c_1 \cdots c_S]$ the values obtained from the above algorithm. Then, an estimate of the contribution X in the data Y is returned by

$$\hat{X} = \sum_{s=1}^{S} c_s b_s^H R = CB^H R. \qquad (7)$$

This algorithm is equivalent to the generalized eigenvalue decomposition (GEVD) of the pair of matrices $(S_{RY} \Gamma S_{YR}, S_{RR})$, for some positive definite (weight) matrix $\Gamma \in \mathbb{C}^{M \times M}$, where only the eigenvectors associated to the S largest eigenvalues $\lambda_1 \ge \lambda_2 \ge \cdots \lambda_S \ge 0$ are conserved (see Appendix A).

It is clear by looking at the first matrix $(S_{RY} \Gamma S_{YR})$ that, if $M < R$, then its rank is M, and there is no hope to search for more than M sources. This means that one must have $M \ge S$ if the effects of all sources are to be recovered, as initially assumed.

The matrix B then will consist of the first S generalized eigenvectors, and $C = S_{YR}B$. Eq. (7) then becomes

$$\hat{X} = CB^H R = S_{YR}BB^H R. \qquad (8)$$

As compared to Eq. (6), the above result replaces the (possibly unstable or non-existing) inverse operator $S_{RR}^{-1}$ by the projector $BB^H$.

Different strategies can be implemented to define the weight matrix $\Gamma$, that enters into the definition of the correlation $\operatorname{corr}(Y_{s-1}, b^H R)$ introduced in Alg. 1:

* A simple choice is $\Gamma = I$.

* Another choice is $\Gamma = S_{YY}^{-1}$ with $S_{YY} = YY^H/I$ the CSM of the measurements, which entails some useful properties. In particular, it then holds that $0 \le \lambda_s \le 1$, which can be interpreted as the partial coherence between the data Y and the sth virtual reference $z_s$. Hence, $\|Y - \hat{X}\|^2/I = \sum_{s=S+1}^{R} \lambda_s$ measures the overall prediction error.

* From the Gauss-Markov theorem, the minimum prediction error is returned by $\Gamma = \Omega^{-1}$, where $\Omega = E\{nn^H\}$ is the noise covariance matrix, with n a column of N. Since $\Omega$ will rarely be known beforehand, it can be estimated iteratively from the data by estimating $\hat{X}$ from the GEVD of $(S_{RY} \Gamma S_{YR}, S_{RR})$, and $\Gamma$ from $((Y - \hat{X})(Y - \hat{X})^H/I)^{-1}$, until convergence. Interestingly, it can be proven that using $\Gamma = \Omega^{-1}$ is equivalent to $\Gamma = S_{YY}^{-1}$, (see Appendix B).

## 2.3. Failure of partial coherence and TSVD in the case of spatially-colored and/or correlated noises on the references

The partial coherence technique is not a good candidate to solve the problem under consideration. Firstly, when the references are in greater number than the sources, thus partially correlated, an ambiguous ordering operation is required. This means there exists a triangular matrix L such that $R = LS$. Second, it requires the availability of pure references, without additive noise. The
TSVD technique is more general as it requires **R** = **LS** + **E**, with **L** full-rank, but **E** with a scalar covariance matrix. To explain the latter condition, let us consider the energetic form of Eq. (4),

$$S_{RR} = \tilde{S}_{RR} + S_{EE}, \qquad (9)$$

where $\tilde{S}_{RR}$ represents the noise-free part, and $S_{EE}$ is the CSM of the noises. Replacing $\tilde{S}_{RR}$ with $LS_{SS}L^H$, where $S_{SS}$ is the CSM of the sources, one gets

$$S_{RR} = LS_{SS}L^H + S_{EE}. \qquad (10)$$

Keeping in mind that the sources are uncorrelated, $S_{SS}$ is diagonal. Moreover, even though the diagonal elements would be of different values, they can be absorbed into the factor **L**; hence, and without any loss of generality, $S_{SS}$ can be taken as the identity matrix. Then,

$$S_{RR} = LL^H + S_{EE}. \qquad (11)$$

Now let us consider the very special case of spatially white (scalar) noise, when $S_{EE}$ is proportional to the identity, i.e. $S_{EE} = \epsilon^2 I$, then the previous equation becomes

$$S_{RR} = LL^H + \epsilon^2 I. \qquad (12)$$

Replacing **L** with its singular value decomposition $U\Sigma V^H$, where $UU^H = I$, $V^H V = I$, and $\Sigma$ is diagonal, and, replacing **I** in the second term with $UU^H$, one gets after some algebra

$$S_{RR} = U(\Sigma^2 + \epsilon^2 I)U^H. \qquad (13)$$

This is recognized as the eigenvalue decomposition of $S_{RR}$, where the variation among the eigenvalues is coming only from the sources. That is why TSVD exclusively works fine in the very special situation of $S_{EE} = \epsilon^2 I$. However, when noises are inter-correlated, $S_{EE} \neq \epsilon^2 I$, and the variations among the eigenvalues is no longer coming only from the sources, but also from the noise.

Therefore, TSVD requires two conditions to be met:

* all the off-diagonal noise elements $[S_{EE}]_{i,j} = 0$ for all $i \neq j$,

* a constant noise power on all the references, i.e. $[S_{EE}]_{i,j} = \epsilon^2$ for all $i$.

The first condition is theoretically met if and only if the noises on all the references are completely uncorrelated with each other, and that an infinite amount of averaging is conducted (i.e. $I \to \infty$). The second condition is even harder to meet. The situation actually gets worse as the noise level increases. Hence, MCR provides a good alternative to TSVD by finding an alternative reference set, that maximally correlates to the sources.

## 2.4. Expected supremacy of MCR in the case of spatially-colored and/or correlated noises on the references

The model presented in Section 2.1 is rewritten, for simplicity, considering one column **s** of **S**, **x** of **X**, **n** of **N**, **y** of **Y**, and **r** of **R**. Correlation among noise elements in $\epsilon E$ is modeled by the linear relationship $\epsilon E = Fv$, with $v$ composed of $Q \le R$ independent noise sources on the references. Then, the model reads

$$x = Hs, \qquad (14)$$

$$y = Hs + n, \qquad (15)$$

$$r = Ls + Fv, \qquad (16)$$

where, $H \in \mathbb{C}^{M \times S}$, $M \ge S$, $L \in \mathbb{C}^{R \times S}$, $R \ge S$, and $F \in \mathbb{C}^{R \times Q}$, $R \ge Q$.

The processes **s**, **n** and **v** are assumed mutually uncorrelated. Without loss of generality, **s** and **v** have covariances $E\{ss^H\} = I$ and $E\{vv^H\} = I$, respectively. **H** and **L** are assumed of full rank (i.e. rank($H$) = rank($L$) = $S$) without loss of generality (otherwise there would be the effect of less than $S$ sources in **y** and **r**). Given $Q$ independent noise sources, **F** is assumed of full rank (i.e. rank($F$) = $Q$) without loss of generality (otherwise there would be the effect of less than $Q$ noise sources in **r**). The proposed MCR method finds the matrix **B** $\in \mathbb{C}^{R \times S}$ which minimizes

$$ \begin{aligned} e &= \langle \|y - S_{YR}BB^Hr\|^2 \rangle \\ &= \langle \|Hs + n - S_{YR}BB^H(Ls + Fv)\|^2 \rangle \\ &\approx \underbrace{\langle \| (H - S_{YR}BB^H L)s \|^2 \rangle}_{(1)} + \underbrace{\langle \|n\|^2 \rangle}_{(2)} + \underbrace{\langle \|S_{YR}BB^H Fv\|^2 \rangle}_{(3)} \end{aligned} \qquad (17) $$
as $\langle sn^* \rangle \to 0$, $\langle sv^* \rangle \to 0$ and $\langle nv^* \rangle \to 0$, with $\langle \cdot \cdot \cdot \rangle$ the average over $I$ snapshots. The problem requires the joint minimization of terms (1) and (3), while term (2) is constant. With $S_{YR} = E\{yr^H\} = HL^H$, (1) $= \langle \|H(I - L^H BB^H L)s\|^2 \rangle$ and (3) $= \langle \|HL^H BB^H Fv\|^2 \rangle$. Conditions to zeroe (1) and (3) are:

$$ (i) \quad L^H BB^H L = I_S \Leftrightarrow B^H L \text{ is a square root of } I_S, $$

$$ (ii) \quad B^H F = O_{S \times Q}. $$

Note that if the two conditions are satisfied, then $x = Hs$ is exactly recovered by MCR, as

$$ \begin{aligned} \hat{x} &= S_{YR} BB^H r = HL^H BB^H (Ls + Fv) \\ &= H \underbrace{L^H BB^H L}_{I} s + HL^H B \underbrace{B^H F}_{O} v \\ &= Hs = x. \end{aligned} \qquad (18) $$

It is proved in Appendix C, that $L^H F = O_{S \times Q}$ is the only condition to satisfy to meet (i) and (ii). Since the orthogonal subspace of $L$ has dimension $dim = R - S$, one must have $Q \le R - S \Leftrightarrow R \ge Q + S$ for this condition to be met. Hence, the exact recovery is guaranteed in the case of correlated noise on the references, given that a sufficient number of references is used. In practice, one neither know the true number of sources $S$ nor the true number of noise sources on the references $Q$, so it is advised to use as many references as possible to approach exact recovery.

## 2.5. Selection of the number of sources using bootstrapped likelihood

Prior knowledge of the true number of sources ($S$) is needed for any rank reduction technique. In practice, the number of sources is rarely known beforehand, hence, it has to be estimated from the data. The likelihood ratio test (LRT) is usually used for solving this problem. It compares the fit of nested models based on Wilks' theorem [21]. However, the conditions for Wilks' theorem, which states that the LRT statistic asymptotically follows a $\chi^2$-distribution with degrees of freedom equal to the difference in the number of free parameters between the two nested models under comparison, showed often not to hold for models with latent variables as the one considered here [22]. Hence, a bootstrapped likelihood ratio (LR) is proposed which enables successively testing models with $K = 0, 1, 2, \dots, R - 1$ tentative sources against the saturated model $K = R$. The estimated source contributions for a model of $K \le R$ sources is

$$ \hat{X}(K) = C(K)B(K)^H R, \qquad (19) $$

as returned by Eq. (8), with $C(K) = [c_1 \dots c_K]$ and $B(K) = [b_1 \dots b_K]$, consisting of the first $K$ generalized eigenvectors of $(S_{RY} \Gamma S_{YR}, S_{RR})$. The corresponding covariance matrix is estimated as

$$ \hat{\Omega}(K) = (Y - \hat{X}(K))(Y - \hat{X}(K))^H / I. \qquad (20) $$

The data $Y$ is assumed to follow a multivariate complex Gaussian, (e.g. a justified assumption in the case of Fourier coefficients according to the central limit theorem [23]) with mean $\hat{X}(K)$ (resp. $\hat{X}(R)$), with $K$ (resp. $R$) components involved, and covariance matrix $\hat{\Omega}(K)$ (resp. $\hat{\Omega}(R)$) under the $K^{th}$ (resp. saturated) model. Then the log-likelihood ratio (LLR) of the saturated model to the $K^{th}$ model reads

$$ \Lambda(K) = -\|Y - \hat{X}(R)\|_{\Omega(R)}^2 + \|Y - \hat{X}(K)\|_{\hat{\Omega}(K)}^2 - I \ln |\Omega(R)| + I \ln |\hat{\Omega}(K)|, \qquad (21) $$

where $\|A\|_{\Omega}^2 = \text{tr}\{A^H \Omega^{-1} A\}$. By substituting for $\Omega(R)$ and $\hat{\Omega}(K)$ with their estimates, from Eq. (20), respectively, the expression in Eq. (21) reduces to

$$ \Lambda(K) = I(\ln |\hat{\Omega}(K)| - \ln |\hat{\Omega}(R)|). \qquad (22) $$

Even though this LLR asymptotically (with $I$) follows a $\chi^2$-distribution under milder conditions than the one used in Ref. [16], it is preferred to use a bootstrapping technique, in order to avoid the need for long records. Hence, parametric bootstrap is used. The LLR defined in Eq. (22) has a reference distribution, say $D(K)$, that one can bootstrap under the restricted model. First, one estimates $\hat{\Omega}(K)$ as in Eq. (20). Next, computes the singular value decomposition of $\hat{\Omega}(K) = U \Sigma V^H$, and then draws new noise matrices as $N_{new} = U \Sigma^{1/2} N_1$, where the elements of $N_1$ follow a complex Gaussian distribution with zero mean and a unit variance. A new set of references $R_{new}$ is also synthesized in the same way using the $S_{RR}$ matrix. Finally the new output data matrix is computed as $Y_{new} = C(K)B(K)^H R_{new} + N_{new}$. Hence, one rejects all the models until finding the first model with $K$ sources that satisfies $\Lambda(K) < D(K)_{1-\alpha}$, with $D(K)_{1-\alpha}$ the quantile of the bootstrapped distribution associated with probability $1 - \alpha$. If so, $K$ sources are selected; if not, the next model with $K + 1$ is tested.

The algorithm for the selection of the number of sources reads:
## Algorithm 2

* Inputs: $S_{YY}$, $S_{RR}$ and $S_{YR}$

* Estimate $C(R)$ and $B(R)$ from $S_{YY}$, $S_{RR}$ and $S_{YR}$ assuming $R$ sources

* Estimate $\hat{X}(R)$ and $\hat{\Omega}(R)$ from Equations (19) and (20) with $K = R$

* Set $K = 0$

* WHILE $K < R$

- Construct $C(K)$ and $B(K)$ from $C(R)$ and $B(R)$ assuming $K$ sources

- Estimate $\dot{X}(K)$ and $\dot{\Omega}(K)$ from Equations (19) and (20)

- Compute $\Lambda(K)$ from Eq. (22)

- Set $i = 1$

- WHILE $i \le N$ (Number of Iteration)

* Draw $\mathbf{R}_{\text{new}}$ using $S_{RR}$

* Draw $\mathbf{N}_{\text{new}}$ using $\hat{\Omega}(K)$

* Compute $\mathbf{Y}_{\text{new}} = C(K)\mathbf{B}(K)^H\mathbf{R}_{\text{new}} + \mathbf{N}_{\text{new}}$

* Estimate $C(R)_{\text{new}}$ and $B(R)_{\text{new}}$ from $S_{Y_{\text{new}}Y_{\text{new}}}$, $S_{R_{\text{new}}R_{\text{new}}}$ and $S_{Y_{\text{new}}R_{\text{new}}}$ assuming $R$ sources

* Construct $C(K)_{\text{new}}$ and $B(K)_{\text{new}}$ from $C(R)_{\text{new}}$ and $B(R)_{\text{new}}$ assuming $K$ sources

* Estimate $\dot{X}(R)_{\text{new}}$ and $\dot{\Omega}(R)_{\text{new}}$ from Equations (19) and (20) applied to $C(R)_{\text{new}}$, $B(R)_{\text{new}}$ and $\mathbf{R}_{\text{new}}$, with $K = R$

* Estimate $\dot{X}(K)_{\text{new}}$ and $\dot{\Omega}(K)_{\text{new}}$ from Equations (19) and (20) applied to $C(K)_{\text{new}}$, $B(K)_{\text{new}}$ and $\mathbf{R}_{\text{new}}$

* Compute $\Lambda(K)_{\text{new}}$ from Eq. (22) applied to $\hat{\Omega}(R)_{\text{new}}$ and $\hat{\Omega}(K)_{\text{new}}$

* Set $D(K)^{(i)} = \Lambda(K)_{\text{new}}$

* $i \leftarrow i + 1$

- IF $\Lambda(K) < D(K)_{1-a}$ STOP where $D(K)_{1-a}$ is quantile of $\{D(K)^{(i)}; i = 1, ..., N\}$ associated with probability $1 - \alpha$

- ELSE $K \leftarrow K + 1$

* IF $K = R$, $S \leftarrow R$,

* ELSE $S \leftarrow K$

## 3. Numerical validation

### 3.1. Simulation A: A comparative sensitivity study

Simulations were performed to compare the performance of MCR, LS and TSVD. Two sources $S = 2$ were simulated as independent and identically distributed (iid) complex Gaussian random processes with zero mean and unit variance, which contribute to the measurements $Y$ and the references $R$ through complex-valued transfer matrices $H \in \mathbb{C}^{M \times S}$ and $L \in \mathbb{C}^{R \times S}$, whose elements were randomly drawn from iid complex Gaussian random processes with zero mean and unit variance (see Eqs. (1) and (3)). The number of sensors was set to $M = 50$, and the number of snapshots to $I = 1k$. Simulations were performed while varying the number of references ($R = 5, 10, 25$), the signal-to-noise ratio (SNR) of the measurements, also referred to as outputs, ($SNR_{\text{out}} = -20, -19, ..., 10$ dB), and the SNR of the references ($SNR_{\text{ref}} = -20, 0, 20$ dB). The behavior of the methods is judged by comparing the estimated $SNR_{\text{out}}$ to its actual value. The estimated $SNR_{\text{out}}$ for the $i^{th}$ microphone is

$$ \hat{\mathbf{SNR}}_{\text{out}} = \frac{[S_{Y_R} \mathbf{BB}^H S_{R_Y}]_{ii}}{[S_{Y_Y}]_{ii} - [S_{Y_R} \mathbf{BB}^H S_{R_Y}]_{ii}} \quad (23) $$

where the subscript $(ii)$ indicates the diagonal elements of a matrix. Two simulations were performed. In the first, uncorrelated noises on the output and the reference sensors were simulated as random complex Gaussian variables with zero mean, and their variances were specified so as to give the aforementioned SNRs. Whereas, in the second, an 0.9 cross-correlation factor was imposed on the noises of the references. The number of the preserved eigenvalues in MCR and TSVD was set to the number of simulated sources. The results of the two simulations are shown in Figs. 1 and 2, where $\hat{\mathbf{SNR}}_{\text{out}}$ is displayed as a function of $SNR_{\text{out}}$.

In Fig. 1, it is clear that when using clean references (i.e. with high $SNR_{\text{ref}}$) in a number close to the number of sources, all the methods return excellent results. However once deviating from these two conditions, two types of deterioration are observed. First, increasing the number of references deteriorates the results of LS, as it tends to overestimate $SNR_{\text{out}}$. Second, by increasing the noise level on the references, the result of TSVD deteriorates, as it tends to underestimates $SNR_{\text{out}}$. The proposed MCR always mediates the over/under estimation of the two other solutions, giving the closest solution to the ground truth, especially when the number of references gets large. In this simulation case of uncorrelated noise on the references, all the three methods perform more or less equivalently and fairly well, except for the case of very low $SNR_{\text{ref}} = -20$ dB, where all the methods were badly impacted by the high levels of noise on the references, especially TSVD.
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Results of simulation A, with uncorrelated noise on the references: $\hat{\text{SNR}}_{\text{out}}$ (estimated value) as a function of $\text{SNR}_{\text{out}}$ (actual value), where the ground truth is the black diagonal, LS in magenta, TSVD in red, and MCR in dotted blue. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

The results of the second simulation are depicted in Fig. 2, where it is noticed that introducing correlation between the noises on the references somehow improved the results of both LS and MCR, except from the persistent overestimation of the former. This can be explained by the possibility that the correlation between the noises on the references compensates for their low SNR. In other words, given identical levels, the effect of *R* correlated noises is less than *R* uncorrelated noises. TSVD gives worse results (significant underestimation) than in the first simulation, in coherence with the discussion of Section 2.3; the case of correlated noises. In conclusion, MCR gives the best results, as it does not underestimate $\text{SNR}_{\text{out}}$ as TSVD does, while avoiding the overestimation of LS, and this is even more prominent in the cases of correlated and/or spatially colored noises on the references, especially when the noise levels are high.

## 3.2. Simulation B: The resynchronization problem

The necessity of resynchronizing sequential measurements always arises when simultaneous measurements are not achievable due to hardware limitations (such as the lack of sensors or channels), or other technical constraints (e.g. array dimension). Given the availability of fixed-in-space sensors (references), the phase lost between sequential measurements can be retrieved using the reference/coherence techniques. Numerical simulations were conducted to assess the efficacy of MCR in addressing the resynchronization problem, as compared to both LS and TSVD. The number of preserved eigenvalues in MCR and TSVD was set to match the number of simulated sources, in this case $S = 1$. The number of references was $R = 10$, the number of sensors $M = 25$, and the number of snapshots $I = 70$.

The sources were simulated as iid complex Gaussian random processes with zero mean and unit variance, and were mapped to the sensors and references using complex-valued transfer matrices $\mathbf{H} \in \mathbb{C}^{M \times S}$ and $\mathbf{L} \in \mathbb{C}^{R \times S}$, whose elements were randomly drawn from iid complex Gaussian random processes with zero mean and unit variance. To introduce spatial characteristics into the simulation, the sensors were arranged in the configuration of a $5 \times 5$ array. To model a field emitted by a single point source facing exactly the middle of the array, the sensor located in the center of the array (sensor no. 13) had its power increased by a factor of 100.
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Results of simulation A, with correlated noise on the references: $\hat{\text{SNR}}_{\text{out}}$ (estimated value) as a function of $\text{SNR}_{\text{out}}$ (actual value), where the ground truth is the black diagonal, LS in magenta, TSVD in red, and MCR in dotted blue. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Uncorrelated noises on the sensors were simulated as random complex Gaussian variables with zero mean, and their variances were adjusted such that $\text{SNR}_{\text{out}} = -10$ dB. Noises on the references were simulated as random complex Gaussian variables with zero mean, and their variances were adjusted such that $\text{SNR}_{\text{ref}} = 0$ dB. It is important to note that using the same SNR for all the sensors led to different noise power levels since the signals themselves had varying power levels. Furthermore, a cross-correlation factor of 0.9 was introduced on the noises of the references. As seen in Section 2.3, these conditions represent a challenge for TSVD.

**Fig. 3** displays radiation maps, which represent interpolated figures of $5 \times 5$ simulated complex-valued measurements. The first row shows the true coherent part in the measurements (a single snapshot x of $X$), while the second row depicts the noisy measurements (a single snapshot y of $Y$). The third, fourth and fifth rows illustrate the estimated coherent parts in the measurements using MCR, TSVD, and LS, respectively. This presentation allows for a concise comparison through visual inspection of the magnitude, the phase, and the real and imaginary parts.

The second row of Fig. 3 demonstrates significant amplitude and phase distortion resulting from the addition of noise. However, MCR shows excellent performance, accurately retrieving the coherent part in terms of both amplitude and phase. On the other hand, TSVD fails to recover either the amplitude or the phase, leading to substantial underestimation. Regarding LS, its behavior is in line with its mathematical definition, tending to return the full noisy measurements, i.e. overestimating $\text{SNR}_{\text{out}}$. It also performs poorly in restoring the phase of the coherent part. In summary, the simulations confirm that MCR better handles the challenges faced by LS and TSVD, including low SNR, correlated references with correlated noise, and heterogeneous noise powers.

## 4. Laboratory experiment

To validate the MCR-based resynchronization technique on real data, an experiment was conducted in the coupled testing rooms of the Laboratory of Vibration and Acoustics (LVA), INSA de Lyon. One of the two rooms is semi-anechoic and the other one is reverberant.
<figure><img src="image_6.png" /><figcaption>Fig. 3. Results of Simulation B: the radiation maps of the (a) magnitude, (b) phase, (c) real, and (d) imaginary parts of the true coherent part of the measurements x (first row), the noisy measurements y (second row), the MCR estimate (third row), the TSVD estimate (fourth row), and the LS estimate (fifth row).</figcaption></figure>

## 4.1. Experimental apparatus

The experimental apparatus is schematically illustrated in Fig. 4(a). A microphone array was placed in the semi-anechoic room with 18 microphones, as shown in Fig. 4(b), in order to capture the sound radiated by a vibrating metallic plate which was fixed in the window between the two rooms and was excited by a shaker. It is depicted in Fig. 4(c) that two shakers were attached to the plate, yet only the one on the right was fed with a white Gaussian noise, with an upper cutoff frequency of 20 kHz, during the experiment to simulate a single source, i.e. $S = 1$. Two accelerometers were fixed to the vibrating plate, as shown in Fig. 4(c), to be used as references. On the reverberant side, three reference microphones were installed near the vibrating plate, as shown in Fig. 4(d). Hence, in this experiment, there were a total of five references $R = 5$, consisting of three microphones and two accelerometers. Finally, the data was acquired with a sampling frequency of 51.2 kHz. The sensitivities of all the sensors were set to 1, so as to keep the same unit (V) for signals coming from different types of sensors.

## 4.2. Spectral and principal component analysis of the experimental data

The signal length was 1,335,652 samples, which is equivalent to a record length of about 26 s. Spectral analysis was done using Welch's averaged periodogram method. A fine frequency resolution of 0.7813 Hz was achieved using a hanning window of length $L = 2^{16}$ samples (1.28 s), with 2/3 overlap, which produced an effective number of snapshots $I_{\text{eff}} = 41$ (dimension of an equivalent set of uncorrelated snapshots [24]) defined as

$$I_{\text{eff}} = \frac{J}{1 + 2 \sum_{j=1}^{J-1} \frac{J-j}{J} \rho(j)}, \qquad (24)$$

where $J > I_{\text{eff}}$ is the actual number of segments (or windows), and

$$\rho(j) = \left[ \sum_{n=0}^{L-1} W(n)W(n + jD_{\text{sh}}) \right]^2 \left[ \sum_{n=0}^{L-1} W^2(n) \right]^2, \qquad (25)$$

is the correlation function, with $W(n)$ the time window and $D_{\text{sh}}$ the shift between consecutive segments. Hence, in this case, $I_{\text{eff}}$ replaces $I$ in Eq. (22).
<figure><img src="image_7.png" /><figcaption>Fig. 4. The LVA experimental apparatus: (a) schematic illustration of the experimental setup in the coupled rooms, (b) the microphone array in the semi-anechoic room, (c) a close-up on the shaker and the accelerometers, and (d) the reference microphones in the reverberant room.</figcaption></figure>

The power spectral densities (PSDs) of the array signals and the reference signals are shown in Figs. 5(a) and 5(b), respectively. In Fig. 5(a), it is seen that some of the microphones contain high levels of noise. It is also obvious that there are two regions of high amplitude (below 3 kHz and [6–10] kHz), which also appear on the reference microphones in Fig. 5(b). This could be due to either an amplification of the system or some external sources. The answer to this question is delayed until discussing the estimation of the number of sources. It is also clear in Fig. 5(b) that the powers of the two reference acceleration signals are higher than those of the reference microphones, which may pose a problem to TSVD. Figs. 5(c) and 5(d) display the principal spectra (eigenvalues of the CSM) of the array signals and of the reference signals, respectively. The principal spectra of the array microphones show a single high principal component, however the high levels on some of the microphones caused the appearance of other principal components high enough for not being discarded. Thus, the number of dominant principal spectra of the array microphones alone does not reflect the actual number of the uncorrelated sources in this experiment, i.e. *S* = 1. Fig. 5(d) shows the principal spectra of the reference set of sensors. Although it was hoped that the latter might have given a clearer idea about the number of sources, they turned out to have a gradual dispersion, where it is not possible to distinguish dominant principal components from the other ones. The authors believe this is due to the difference in power levels between the acceleration and the pressure references (see Fig. 5(b)).

### 4.3. MCR analysis of the experimental data

MCR was performed through the computation of the GEVD ($S_{RY}\Gamma S_{YR}\cdot S_{RR}$), where the weight matrix was chosen to be the identity, $\Gamma = I_M$. The MCR-principal spectra are shown in Fig. 6(a), where the existence of one dominant MCR-principal component is obvious, which corresponds to the true stochastic dimension of the field as there is one source in the experiment. The dimensionality of the field is now easily assessable by visual inspection. However, it is always preferable to have an algorithm that automatizes the selection of the important principal components. Hence, the bootstrap algorithm (Algorithm 2) was used, with a risk $\alpha = 0.01$, in order to automatically estimate the number of sources. The estimation results are shown in Fig. 6(b), where one can see that $\hat{S} = S = 1$ all along the frequency axis, except for few frequency lines within two ranges (below 3 kHz and in [6–10] kHz), where it gave $\hat{S} = 2$ to 3. These ranges correspond to the energetic ranges in the PSDs of the microphone signals. This answers the question posed in the previous section about the origin of this power amplification: the rise in $\hat{S}$ indicates the existence of one to two extra sources external to the system. One possibility is that they come from the data acquisition system. Therefore, the estimated number of sources is taken as $\hat{S} = 1$ for all frequencies up to 20 kHz. For the estimated number of sources above 20 kHz, as it is seen in Fig. 6(b), the algorithm returned $\hat{S} = 0$, what further proves the effectiveness of the proposed algorithm, as that was the upper cutoff frequency of the excitation signal fed to the shaker. Hence, for the resynchronization analysis in the next subsection only frequencies up to 20 kHz are considered.
<figure><img src="image_8.png" /><figcaption>Fig. 5. Spectral analysis of the LVA experimental data: (a) PSDs of the array microphones, (b) PSDs of the references (3 microphones + 2 accelerometers), (c) principal spectra of the CSM of the array microphones, and (d) of the references.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. MCR analysis of the LVA experimental data: (a) MCR-based principal spectra using the five references, and (b) estimated number of sources $\hat{S}$.</figcaption></figure>

### 4.4. Application to resynchronization

To test the ability of MCR to resynchronize sequential measurements, the signals of the array (Y) were divided into two sets, $Y_1$ and $Y_2$, each of which consists of half the number of signals of the original set, i.e. $(18/2 = 9)$ signals. Then, the phase between the two sets was distorted by shifting the signals of $Y_1$ in the time domain by the signal length, thus the two subsets were no longer synchronous, however $Y_2$ was still synchronous to the original data set Y. Two reference sets were used, the original measured set $R$, which is synchronous to Y and $Y_2$, and its time-shifted version $R_1$, synchronous to $Y_1$. In this manner, successive asynchronous measurements were created, to mimick the scenario where, for instance, an array is moved sequentially in front of a radiating object [25]. Analyzing such artificially phase-distorted data has the merit of providing the original synchronous data, needed to assess the performance of the resynchronization. Spectral analysis was performed, using the same signal length and spectral analysis parameters as before.

Fig. 7 shows all the data sets estimated in the resynchronization analysis, where asynchronous data are shown in different colors. First, the true source contributions $X$ in the original set $Y$ is recovered as $\hat{X} = GR$, where $G$ is estimated from MCR applied to (Y, R) (i.e. $\hat{X} = \hat{X}_{MCR}$). Next, the source contributions in each subset are recovered. $\hat{X}_1 = G_1R_1$, is the estimate of the source contribution in $Y_1$, where $G_1$ is estimated from $(Y_1, R_1)$ using MCR (i.e. $\hat{X}_1 = \hat{X}_{1MCR}$), whereas, $\hat{X}_2 = G_2R$, is the estimate of the source contribution in $Y_2$, where $G_2$ is estimated from $(Y_2, R)$ using LS, TSVD and MCR, thus providing $\hat{X}_{2LS}$, $\hat{X}_{2TSVD}$ and $\hat{X}_{2MCR}$. The concatenated matrix

$$ \hat{X} = [\hat{X}_{1MCR}^T \quad \hat{X}_{2MCR}^T]^T \qquad (26) $$
<figure><img src="image_10.png" /><figcaption>Fig. 7. Different matrices used in the resynchronization of the out-of-phase two data subsets, where different colors means out-of-phase, and the shaded parts are the parts used for calculating the adjusted resynchronization error. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

represents the asynchronous solution. Finally, a resynchronized estimate $\tilde{X}_1 = G_1 R$ of the source contributions in $Y_1$ is estimated, where $G_1$ is obtained from $(Y_1, R_1)$ using LS, TSVD and MCR, leading to $\tilde{X}_{1LS}, \tilde{X}_{1TSVD}$ and $\tilde{X}_{1MCR}$. Three resynchronized solutions are then obtained,

$$ \tilde{X}_{LS} = [\tilde{X}_{1LS} \quad \tilde{X}_{2LS}^T]^T, \qquad (27) $$

$$ \tilde{X}_{TSVD} = [\tilde{X}_{1_{TSVD}}^T \quad \tilde{X}_{2_{TSVD}}^T]^T, \qquad (28) $$

$$ \tilde{X}_{MCR} = [\tilde{X}_{1_{MCR}}^T \quad \tilde{X}_{2_{MCR}}^T]^T. \qquad (29) $$

The relative error $e$ between each of the resynchronized solutions, $\tilde{X}$, and the target $\hat{X}$ is calculated as $e = \|\tilde{X} - \hat{X}\|_F / \|\hat{X}\|_F$. The relative error values of the three methods are plotted against frequency in Fig. 8(a), along with the relative error of the asynchronous solution $\hat{X}$. It is clear that MCR returns the lowest resynchronization error, for nearly all frequencies. However, it could be argued that the error $e$ is essentially biased towards the MCR solution, as the submatrix $\hat{X}^{(2)}$, which is the second-half signals of the target $\hat{X}$ (see Fig. 7), is theoretically equal to the matrix $\hat{X}_{2MCR}$, as both are MCR estimates of the same set of signals. To avoid this bias, an adjusted error formula $e_{ad} = \|\hat{X}_1 - \hat{X}^{(1)}\|_F / \|\hat{X}^{(1)}\|_F$ is used to evaluate the resynchronization error of the compared techniques, which involves only the resynchronized subset $\hat{X}_1$ and its counterpart subset of $\hat{X}$ (i.e., $\hat{X}^{(1)}$). The values of the adjusted error against frequency for the three compared methods, as well as the asynchronous solution, are shown in Fig. 8(b), where the supremacy of MCR is again evident, as it returns the lowest error for nearly all frequencies. The average relative errors over all frequencies, are shown in Table 1, from which the total improvement in resynchronization obtained by MCR is quantitatively validated.

It is also useful to compare the resynchronization effectiveness among the methods using radiation maps, as done in the previous simulation. Figs. 9 and 10 show the resynchronization results at an arbitrary chosen frequency of $f = 7$ kHz, at snapshots $i = 33$ and $i = 56$, respectively. In these two figures, the first row shows a single snapshot $\hat{x}$ of the target, the second row the phase-distorted (asynchronous) solution $\hat{x}$, the third row the MCR solution, the fourth row the TSVD solution, and the fifth row the LS solution.

It is again observed that MCR is the only one able to accurately retrieve the phase of the target. The phase maps of TSVD and LS are distorted, while MCR retrieves nearly exactly the correct phase. The most remarkable result is the large amplitude underestimation returned by TSVD (see Fig. 10). Besides, LS appears to give noisy results. It is still clear that MCR gives the best results, while TSVD completely fails.

## 5. Application to electric motor data

Data from the Siemens SimRod electric motor (e-motor) test bench were used to validate the proposed method, where the noise emitted by the motor was measured using an array consisting of 50 microphones, located 2 cm away from the motor. For this
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Resynchronization relative error (a) and its adjusted version (b) for the three solutions under comparison: asynchronous data (black), LS (orange), TSVD (purple), and MCR (green), for the LVA experimental data. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<caption>**Table 1**</caption>

Mean relative error $e$, and adjusted mean relative error $e_{ad}$, over all frequencies, as returned by each compared resynchronization method, as well as by the asynchronous solution $\dot{X}$, for the LVA experimental data.

<table><thead><tr><th></th><th>X&#775;</th><th>LS</th><th>TSVD</th><th>MCR</th></tr></thead><tbody><tr><td>e</td><td>1.0252</td><td>0.2264</td><td>0.1622</td><td>0.1316</td></tr><tr><td>e<sub>ad</sub></td><td>1.4478</td><td>0.2524</td><td>0.1915</td><td>0.1765</td></tr></tbody></table>

analysis, the array was positioned as depicted by the blue rectangle in Fig. 11. Five fixed microphones were used as references. The microphone signals were sampled at a frequency of 40960 Hz. Besides, data from two tri-axial accelerometers were used as additional references, resulting in a total of 11 references. The acceleration tracks were sampled at 10240 Hz, which is only one-fourth of the microphone signals sampling frequency, necessitating an up-sampling by a factor of 4. In order to avoid numerical problems due to the empty acceleration spectra, beyond their initial Nyquist frequency of 5120 Hz, white Gaussian noises, with an SNR of 50 dB, were added to the acceleration signals. The time records had a signal length of approximately 14.6 s. However, for the sake of performing resynchronization analysis by splitting the array signals into two asynchronous sets, in the same manner as in Section 4.4, all the time records were split into two equal portions, leaving 7.3 seconds/300k samples-long signals.

## 5.1. Spectral, principal component, and MCR analyses

For a better conditioning of the CSM of the references, the Fourier spectra of all signals were normalized by dividing them by their corresponding RMS values, allowing the incorporation of acceleration signals alongside microphone (pressure) signals. The PSDs were computed using Welch's averaged periodogram method with a hanning window of $2^{13}$ samples, and using 2/3 overlap, which yielded a frequency resolution of 5 Hz and an effective number of snapshots of 74. The analysis was restricted only up to 5 kHz (the initial Nyquist frequency of the acceleration signals), as above this frequency, there was no information in the acceleration signals.

The PSDs of the array signals and the reference signals are shown in Figs. 12(a) and 12(b), respectively, whereas, Figs. 12(c) and 12(d) display the principal spectra (calculated as the eigenvalues of the CSM) of each set, respectively. The gradual decrease
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Resynchronization of two subsets at *f* = 7 kHz and *i* = 33: radiation maps of the (a) magnitude, (b) phase, (c) real, and (d) imaginary parts, of the MCR estimate of the full array $\bar{x}$ (first row), the phase-distorted data $\dot{x}$ (second row), the MCR-based resynchronization (third row), the TSVD-based resynchronization (fourth row), and the LS-based resynchronization (fifth row), for the LVA experimental data.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Resynchronization of two subsets at *f* = 7 kHz and *i* = 56: radiation maps of the (a) magnitude, (b) phase, (c) real, and (d) imaginary parts, of the MCR estimate of the full array $\bar{x}$ (first row), the phase-distorted data $\dot{x}$ (second row), the MCR-based resynchronization (third row), the TSVD-based resynchronization (fourth row), and the LS-based resynchronization (fifth row), for the LVA experimental data.</figcaption></figure>

in the principal spectra, of both sets, hinders drawing any conclusions about the underlying number of sources. Fig. 13(a) shows the MCR-principal spectra, which can be recognized, with careful observation, as a coherence-based mix of the principal spectra of the array and the reference sets. In this specific case, there is no clear threshold that distinguishes dominant principal components due to sources from negligible ones due to noise. Thus, the bootstrap algorithm (Algorithm 2) was used for the estimation of the number of sources. The estimated number of sources $\hat{S}$ is shown in Fig. 13(b), which is about 5 to 6 on average.

## 5.2. Application to resynchronization

A resynchronization test similar to the one in the previous section was performed, where successive measurements were created out of a single array position by dividing the signals of the array into two subsets, $Y_1$ and $Y_2$, each of which consisting of half the number of signals of the original set, i.e. (50/2 = 25) signals. Then, the phase between the two sets was distorted (by shifting the signals of $Y_1$ in the time domain), to make them asynchronous. A reference set $R_1$ was created from the original reference set $R$, so as to be synchronous to $Y_1$, consisting of the full number of signals of the original set (11 signals). Spectral analysis was performed,
<figure><img src="image_14.png" /><figcaption>Fig. 11. Siemens e-motor test bench, the array position is marked by the blue frame. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Spectral analysis of the Siemens e-motor data: (a) PSDs of the array microphones, (b) PSDs of the references (5 microphones + 6 accelerometers), (c) principal spectra of the CSM of the array microphones, and (d) of the CSM of the references.</figcaption></figure>

using the same signal length and spectral analysis parameters as before, and using the number of sources estimated in the previous subsection. The relative error *e* and its adjusted version *e*<sub>ad</sub> are compared for LS, TSVD and MCR, as shown in Fig. 14. It is seen that MCR returns the lowest error for nearly all frequencies. The average errors, over all frequencies, are displayed in Table 2, from which the improvement in resynchronization obtained by MCR is quantitatively validated.

A radiation-map-based comparison at a frequency of *f* = 1 kHz is also shown in Figs. 15 and 16, for the snapshots *i* = 95 and *i* = 100, respectively. It is seen that MCR is the only technique able to accurately retrieve the phase of the target.

## 6. Conclusions

The difficulties that arise when dealing with source extraction or resynchronization in remote sensing primarily stem from the conditioning of the CSM of the references. An ill-conditioned CSM leads to the failure of LS, while a well-conditioned one leads to
<figure><img src="image_16.png" /><figcaption>Fig. 13. MCR analysis of the Siemens e-motor data: (a) MCR-based principal spectra using the eleven references, and (b) estimated number of sources $\hat{S}$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Resynchronization relative error (a) and its adjusted version (b) for the three solutions under comparison: asynchronous data (black), LS (orange), TSVD (purple), and MCR (green), for the Siemens e-motor data. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<caption>Table 2</caption>

Mean relative error $e$, and adjusted mean relative error $e_{ad}$, over all frequencies, as returned by each compared resynchronization method, as well as by the asynchronous solution $\hat{X}$, for the Siemens e-motor data.

<table><thead><tr><th></th><th>X&#x0307;</th><th>LS</th><th>TSVD</th><th>MCR</th></tr></thead><tbody><tr><td>e</td><td>1.0097</td><td>0.1986</td><td>0.2220</td><td>0.1834</td></tr><tr><td>e<sub>ad</sub></td><td>1.3919</td><td>0.2554</td><td>0.2651</td><td>0.2408</td></tr></tbody></table>

the failure of TSVD. Therefore, the motivation behind this work was to find an optimal solution that combines the strengths of both approaches and mitigates their weaknesses as much as possible. The MCR technique was developed by searching a solution that maximizes the coherent spectra. This involves identifying a set of virtual references, which are essentially linear combinations of the actual references, that exhibit maximum correlation with the measurements and, consequently, with the true sources. When tested
<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Resynchronization of two subsets at *f* = 1 kHz and *i* = 95: radiation maps of the (a) magnitude, (b) phase, (c) real, and (d) imaginary parts, of the MCR estimate of the full array *k* (first row), the phase-distorted data *k* (second row), the MCR-based resynchronization (third row), the TSVD-based resynchronization (fourth row), and the LS-based resynchronization (fifth row), for the Siemens e-motor data.</figcaption></figure>

in simulations and on real experimental data, MCR yielded excellent results in comparison to both LS and TSVD, demonstrating significant improvements in terms of both source extraction and resynchronization.

CRediT authorship contribution statement

**Muhammad N. Albezzawy:** Writing – original draft, Validation, Methodology, Conceptualization. **Jérôme Antoni:** Writing – review & editing, Supervision, Methodology, Conceptualization. **Quentin Leclère:** Writing – review & editing, Supervision, Methodology, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

We gratefully acknowledge the European Commission for its support of the Marie Skłodowska-Curie program through the H2020 ETN ECO-DRIVE project (GA 858018). We also would like to express our special thanks to SIEMENS for providing us with the e-motor data. This work was performed within the framework of the Labex CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ANR-11-IDEX-0007).
<figure><img src="image_19.png" /><figcaption>Fig. 16. Resynchronization of two subsets at f = 1 kHz and i = 100: radiation maps of the (a) magnitude, (b) phase, (c) real, and (d) imaginary parts, of the MCR estimate of the full array $\hat{x}$ (first row), the phase-distorted data $\hat{x}$ (second row), the MCR-based resynchronization (third row), the TSVD-based resynchronization (fourth row), and the LS-based resynchronization (fifth row), for the Siemens e-motor data.</figcaption></figure>

## Appendix A. Generalized eigenvalue decomposition

We now prove that the vectors $b_s$, $s = 1, ..., S$ are directly obtained by the GEVD of the pair of matrices $(S_{RY}, \Gamma S_{YR}, S_{RR})$ for some weight matrix $\Gamma$, with $S_{YR} = YR^H / I$ and $S_{RY} = S_{YR}^H$.

It is first noted that the set of virtual references $\{z_1, ..., z_S\}$ is defined up to any scaling (changing the norm of any $z_s$ will be directly compensated by the scaling $c_s$ in Eq. (7)). Therefore, one can equivalently place an arbitrary scaling on the vectors $b_s$. Here, the following choice is made:

$$ \mathbf{b}_s^H \mathbf{S}_{RR} \mathbf{b}_s = 1. \qquad (A.1) $$

Now, assume $\mathbf{b}_1$ is known and find $c_1$ under constraint (A.1). Define the quadratic norm of a matrix $\mathbf{A}$ as $\|\mathbf{A}\|_{\Gamma}^2 = \text{tr}\{\mathbf{A}^H \mathbf{\Gamma} \mathbf{A}\}$ for some positive definite matrix $\mathbf{\Gamma} \in \mathbb{C}^{M \times M}$. Therefore,

$$ c_1 = \arg \min_{\mathbf{c}} \|\mathbf{Y} - \mathbf{c}_1 \mathbf{b}_1^H \mathbf{R} \|_{\Gamma}^2 = \mathbf{S}_{YR} \mathbf{b}_1. \qquad (A.2) $$

The next step is to substitute $c_1$ with this expression and now minimize $\|\mathbf{Y} - \mathbf{c}_1 \mathbf{b}_1^H \mathbf{R} \|_{\Gamma}^2 = \|\mathbf{Y} - \mathbf{S}_{YR} \mathbf{b}_1 \mathbf{b}_1^H \mathbf{R} \|_{\Gamma}^2$, with respect to $\mathbf{b}_1$, still under the constraint (A.1). The solution is given by the equation

$$ \mathbf{S}_{RY} \mathbf{\Gamma} \mathbf{S}_{YR} \mathbf{b}_1 = \lambda_1 \mathbf{S}_{RR} \mathbf{b}_1 \qquad (A.3) $$

for the largest possible scalar $\lambda_1$. This is the eigenvector associated with the largest eigenvalue, $\lambda_1$, of the pair of matrices $(S_{RY}, \Gamma S_{YR}, S_{RR})$.
The process is repeated for finding $c_2$ and $b_2$, where $Y$ is now replaced by the residual $Y_1 = Y - c_1 z_1$, where the contribution of $z_1 = b_1^H R$ has been removed. By construction, $z_2 = b_2^H R$ will be orthogonal to $z_1$, which implies that $b_2$ will be orthogonal to $b_1$ with respect to $S_{RR}$, i.e.

$$b_2^H S_{RR} b_1 = 0. \qquad (A.4)$$

Overall, the complete process is recognized as the GEVD of the pair of matrices ($S_{RY}\Gamma S_{YR}, S_{RR}$), where only the eigenvectors associated to the $S$ largest eigenvalues $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_S \ge 0$ are conserved.

Eq. (7) then becomes

$$\hat{X} = CB^H R = S_{YR} BB^H R.$$

## Appendix B. Equivalence between $\Gamma = S_{YY}^{-1}$ and $\Gamma = \hat{\Omega}^{-1}$

An interesting result is that the GEVD

$$S_{RY}\Gamma S_{YR}\mathbf{B} = S_{RR}\mathbf{B}\mathbf{\Lambda}, \quad \text{with diagonal } \mathbf{\Lambda} \qquad (B.1)$$

returns the same solution $\mathbf{B}$ whether $\Gamma = S_{YY}^{-1}$ or $\Gamma = \hat{\Omega}^{-1}$. To prove it, let suppose that $\mathbf{B}$ satisfies (B.1) with $\Gamma = S_{YY}^{-1}$, and check that it then also satisfies it for $\Gamma = \hat{\Omega}^{-1}$ for some (generally different) diagonal matrix $\mathbf{\Lambda}$.

First note that

$$\begin{aligned} \hat{\Omega} &= (Y - CB^H R)(Y - CB^H R)^H / I \\ &= S_{YY} - S_{YR} BB^H S_{RY}, \end{aligned} \qquad (B.2)$$

where $C = S_{YR}\mathbf{B}$ and the constraint $\mathbf{B}^H S_{RR}\mathbf{B} = \mathbf{I}$ have been used. Using the matrix inversion lemma,

$$\hat{\Omega}^{-1} = S_{YY}^{-1} - S_{YY}^{-1} S_{YR} \mathbf{B} (\mathbf{B}^H S_{RY} S_{YR} \mathbf{B} - \mathbf{I})^{-1} \mathbf{B}^H S_{RY} S_{YY}^{-1}. \qquad (B.3)$$

Therefore,

$$S_{RY}\hat{\Omega}^{-1}S_{YR}\mathbf{B} = S_{RY}S_{YY}^{-1}S_{YR}\mathbf{B} - S_{RY}S_{YY}^{-1}S_{YR}\mathbf{B}(\mathbf{B}^H S_{RY}S_{YR}\mathbf{B} - \mathbf{I})^{-1}\mathbf{B}^H S_{RY}S_{YY}^{-1}S_{YR}\mathbf{B}. \qquad (B.4)$$

Using the equality (B.1) with $\Gamma = S_{YY}^{-1}$, this becomes

$$\begin{aligned} S_{RY}\hat{\Omega}^{-1}S_{YR}\mathbf{B} &= S_{RR}\mathbf{B}\mathbf{\Lambda} - S_{RR}\mathbf{B}(\mathbf{B}^H S_{RY}S_{YR}\mathbf{B}\mathbf{\Lambda} - \mathbf{I})^{-1}\mathbf{\Lambda} \\ &= S_{RR}\mathbf{B}(\mathbf{I} - (\mathbf{B}^H S_{RY}S_{YR}\mathbf{B}\mathbf{\Lambda} - \mathbf{I})^{-1})\mathbf{\Lambda}. \end{aligned} \qquad (B.5)$$

Since the equality (B.1) is verified with $\Gamma = S_{YY}^{-1}$, it comes (after pre-multiplying it with $\mathbf{B}^H$) that

$$\mathbf{B}^H S_{RY} S_{YR} \mathbf{B} = \mathbf{B}^H S_{RR} \mathbf{B} \mathbf{\Lambda} = \mathbf{\Lambda}. \qquad (B.6)$$

Therefore,

$$S_{RY}\hat{\Omega}^{-1}S_{YR}\mathbf{B} = S_{RR}\mathbf{B}\mathbf{\Lambda}. \qquad (B.7)$$

with $\mathbf{\Lambda} = \mathbf{I} - (\mathbf{\Lambda} - \mathbf{I})^{-1}\mathbf{\Lambda}$ a diagonal matrix. This proves that $\mathbf{B}$ solves the same GEVD when $\Gamma = S_{YY}^{-1}$ or $\Gamma = \hat{\Omega}^{-1}$.

As a consequence, the weights $\Gamma = S_{YY}^{-1}$ and $\Gamma = \hat{\Omega}^{-1}$ also return the same quadratic error. The quadratic error reads

$$e = \|Y - CBR\|^2 / I = \text{tr}\{\hat{\Omega}\} = \text{tr}\{S_{YY} - S_{RY} BB^H S_{YR}\}, \qquad (B.8)$$

where $C = S_{YR}\mathbf{B}^H$ and the constraint $\mathbf{B}^H S_{RR}\mathbf{B} = \mathbf{I}$ have been used. Since it does only depend on $\mathbf{B}$, it is the same for $\Gamma = S_{YY}^{-1}$ and $\Gamma = \hat{\Omega}^{-1}$.

## Appendix C. The condition for exact recovery of $x$

Analysis of the conditions (i) and (ii) for exact recovery of $x$ reveals the following constraints:

(i) $L^H BB^H L = I_S \Leftrightarrow \mathbf{B}^H L$ is a square root of $I_S$: this imposes $\frac{S(S^{-1})}{2} \times 2 + S = S^2$ constraints. $\mathbf{B}$ has $R \times S$ unknowns.

(ii) $\mathbf{B}^H \mathbf{F} = O_{S \times Q}$: the orthogonal subspace of $\mathbf{F}$ has $\text{dim } R - Q$ (i.e., there exist $R - Q$ independent vectors $\mu_k$ such that $\mu_k^H \mathbf{F} = O_{(R-Q) \times Q}$). Since $\mathbf{B} \in \mathbb{C}^{R \times S}$ has $S$ columns, this condition can only be met if $R - Q \ge S \Leftrightarrow R \ge Q + S$.

(iii) There is an additional constraint, existing in the original assumptions of the MCR technique, $\mathbf{B}^H S_{RR}\mathbf{B} = I_S$, which removes $S^2$ degrees of freedom.

* Condition (i) is met if and only if $\mathbf{B}^H = \mathbf{W}\mathbf{L}^+$, $\mathbf{W} \in \mathbb{C}^{S \times S}$: $\mathbf{W}\mathbf{W}^H = \mathbf{I}$, and $\mathbf{L}^+ \in \mathbb{C}^{S \times R}$: $\mathbf{L}^+\mathbf{L} = \mathbf{I}$, hence, $\mathbf{L}^H \mathbf{B}\mathbf{B}^H \mathbf{L} = \mathbf{L}^H \mathbf{L}^{+H} \mathbf{W}^H \mathbf{W}\mathbf{L}^{+}\mathbf{L} = \mathbf{W}^H \mathbf{W} = \mathbf{I}$.
• Then with $S_{RR} = LL^H + FF^H$, the additional constraint (iii) is satisfied as follows,

$$ 
\begin{aligned}
B^H S_{RR} B &= B^H LL^H B + B^H FF^H B \\
&= (WL^+L)(WL^+L)^H + B^H FF^H B \\
&= WW^H + B^H FF^H B \\
&= I + O,
\end{aligned}
 $$

if $B^H F = O$, which is condition (ii).

Therefore, with $B^H F = O$ and $B^H = WL^+$, all constraints are met, but this requires that

$$ 
\begin{aligned}
WL^+F = O &\Leftrightarrow L^+F = O \\
&\Leftrightarrow (L^H L)^{-1} L^H F = O \\
&\Leftrightarrow L^H F = O.
\end{aligned}
 $$

Hence, a necessary condition to achieve exact recovery of $x$ is that $L_{S \times R}^H F_{R \times Q} = O_{S \times Q}$. Since the orthogonal subspace of $L$ has $dim = R - S$, one must have $Q \le R - S \Leftrightarrow R \ge Q + S$.

## Data availability

The authors do not have permission to share data.

## References

[1] J. Antoni, Y. Liang, Q. Leclère, Reconstruction of sound quadratic properties from non-synchronous measurements with insufficient or without references: Proof of concept, J. Sound Vib. 349 (2015) 123–149, MAG ID: 2043944389.

[2] B. Dong, J. Antoni, E. Zhang, Blind separation of sound sources from the principle of least spatial entropy, J. Sound Vib. 333 (9) (2014) 2643–2668.

[3] B. Dong, J. Antoni, A. Pereira, W. Kellermann, Blind separation of incoherent and spatially disjoint sound sources, J. Sound Vib. 383 (2016) 414–445.

[4] Hongjie Hou, Fangli Ning, Qingbo Zhai, Jing Li, Keqiang Yao, Wenxun Li, Juan Wei, Cross-spectral matrix denoising for beamforming based on Schatten-p norm, Appl. Acoust. 197 (2022) 108938, MAG ID: 4288057960.

[5] L. Yu, Liang Yu, L. Yu, J. Antoni, Q. Leclère, Spectral matrix completion by Cyclic Projection and application to sound source reconstruction from non-synchronous measurements, J. Sound Vib. 372 (372) (2016) 31–49, MAG ID: 2292417508.

[6] L. Yu, J. Antoni, Q. Leclère, W. Jiang, Acoustical source reconstruction from non-synchronous sequential measurements by Fast Iterative Shrinkage Thresholding Algorithm, J. Sound Vib. 408 (2017) 351–367, MAG ID: 2739534060.

[7] Dingyu Hu, Jinjin Ding, Han Zhao, Liang Yu, Achieving high-resolution 3D acoustic imaging in a large-aspect-ratio cabin by the non-synchronous measurements, Appl. Acoust. 197 (2022) 108953, MAG ID: 4289526712.

[8] Di Wang, Jun Guo, X. Xiao, X. Sheng, CSA-based acoustic beamforming for the contribution analysis of air-borne tyre noise, Mech. Syst. Signal Process. 166 (2022) 108409, MAG ID: 3199973575.

[9] J. Bendat, A. Piersol, Random data: Analysis and measurement procedures, in: Wiley Series in Probability and Statistics, Wiley, 2011, URL: https://books.google.fr/books?id=qYsviFRNMlwC.

[10] S.H. Yoon, P.A. Nelson, Method for the efficient construction of acoustic pressure cross-spectral matrices, J. Sound Vib. 233 (2000) 897–920.

[11] J.S. Bendat, Solutions for the multiple input/output problem, J. Sound Vib. 44 (3) (1976) 311–325.

[12] C. Dodds, J. Robson, Partial coherence in multivariate random processes, J. Sound Vib. 42 (2) (1975) 243–249, MAG ID: 2095035951.

[13] J. Bendat, A. Piersol, Engineering applications of correlation and spectral analysis, in: A Wiley-interscience publication, Wiley, 1993, URL: https://books.google.fr/books?id=DI1RAAAAMAAJ.

[14] S. Price, R. Bernhard, Virtual coherence: A digital signal processing technique for incoherent source identification, in: Proceedings of IMAC, Vol. 4, 1986, pp. 3–6.

[15] D. Otte, Development and Evaluation of Singular Value Analysis Methodologies for Studying Multivariate Noise and Vibration Problems (Ph.D. thesis), KUL, 1994.

[16] M.N. Albezzawy, J. Antoni, Q. Leclère, The Maximally-Coherent Reference Technique: a solution to deal with large sets of references, in: Proceedings of ISMA2022, Leuven, Belgium, 2022, pp. 1683–1697, URL: https://hal.science/hal-04063060.

[17] H.-S. Kwon, Y.-J. Kim, J.S. Bolton, Compensation for source nonstationarity in multireference, scan-based near-field acoustical holography, J. Acoust. Soc. Am. 113 (1) (2003) 360–368.

[18] Q. Leclère, Multi-channel spectral analysis of multi-pass acquisition measurements, Mech. Syst. Signal Process. 23 (5) (2009) 1415–1422.

[19] M. Lee, J.S. Bolton, Scan-based near-field acoustical holography and partial field decomposition in the presence of noise and source level variation, J. Acoust. Soc. Am. 119 (1) (2006) 382–393.

[20] Q. Leclère, A. Dinsenmeyer, J. Antoni, E. Julliard, A. Pintado-Peño, Thresholded multiple coherence as a tool for source separation and denoising: Theory and aeroacoustic applications, Appl. Acoust. 178 (2021) 108021.

[21] S.S. Wilks, The large-sample distribution of the likelihood ratio for testing composite hypotheses, Ann. Math. Stat. 9 (1) (1956) 60–62.

[22] Y. Chen, I. Moustaki, H. Zhang, A note on likelihood ratio tests for models with latent variables, 2020, URL: https://arxiv.org/abs/2008.03971.

[23] M. Peligrad, W.B. Wu, Central limit theorem for fourier transforms of stationary processes, Ann. Probab. 38 (5) (2010) 2009–2022.

[24] P. Welch, The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms, IEEE Trans. Audio Electroacoust. 15 (2) (1967) 70–73.

[25] Q. Leclère, B. Laulagnet, L. Polac, Application of an innovative acoustic imaging technique to assess acoustic power maps of a gasoline engine, in: Forum Acusticum 2005, Hungary, 2005, URL: https://hal.science/hal-00414817.