# Source enumeration in reduced-rank reference techniques for source extraction

Muhammad N. Albezzawy, Jérôme Antoni, Quentin Leclere

► To cite this version:

Muhammad N. Albezzawy, Jérôme Antoni, Quentin Leclere. Source enumeration in reduced-rank reference techniques for source extraction. Journal of Sound and Vibration, In press, 10.1016/j.jsv.2025.119167. hal-05071476

HAL Id: hal-05071476

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

Submitted on 16 May 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
Journal Pre-proof

Source enumeration in reduced-rank reference techniques for source
extraction

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

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

PII: S0022-460X(25)00241-X

DOI: https://doi.org/10.1016/j.jsp.2025.119167

Reference: YJSVI 119167

To appear in: *Journal of Sound and Vibration*

Received date : 25 August 2024

Revised date : 24 April 2025

Accepted date : 29 April 2025

Please cite this article as: M.N. Albezzawy, J. Antoni and Q. Leclère, Source enumeration in reduced-rank reference techniques for source extraction, *Journal of Sound and Vibration* (2025), doi: https://doi.org/10.1016/j.jsp.2025.119167.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2025 Published by Elsevier Ltd.
# Source enumeration in reduced-rank reference techniques for source extraction

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

<sup>a</sup>INSA Lyon, LVA, UR677, Villeurbanne, 69621, France.

## Abstract

In source extraction, source enumeration is essential when the true number of sources is unknown. This is the case in reduced-rank reference (coherence) techniques, where a number of references higher than the number of sources is used. The estimation of the number of sources is critical for accurate source extraction. However, this ill-posed inverse problem has not been sufficiently addressed in the literature within the framework of reference techniques. In this paper, after providing a unified formalism for all reference techniques in the literature, three alternative source enumeration methods applicable to all reference techniques are presented: a direct likelihood ratio test (LRT) against the saturated model, a parametric bootstrap technique, and a cross-validation approach. A comparative study is conducted among the three methods based on simulated numerical data, real sound experimental data, and real industrial data from an electric motor. The results reveal two important findings. First, the number of snapshots (spectral windows) used in spectral analysis significantly affects the performance of the three methods, and they behave differently for the same number of snapshots. Second, parametric bootstrapping proves to be the best method in terms of both estimation accuracy and robustness concerning the number of snapshots used.

**Keywords:** Source enumeration, Reference techniques, Sound source separation, Likelihood ratio test, Parametric bootstrap, Cross-validation

<footnote>*Corresponding author.</footnote>

<footnote>E-mail address: muhammad.albezzawy@insa-lyon.fr</footnote>
# 1. Introduction

In remote sensing of physical fields, where the contributions from mutually incoherent sources overlap in both the spatial and the frequency domains, no definite single phase value can be assigned at any sensor, and hence sources can be hardly identified, in a phenomenon called spatial noise [1]. Hence, a pre-processing source separation step is needed before retro-propagating the measured field from the measurement plane to the source plane. The existence of incoherent output noise is an example of this problem, where the measurements are corrupted by additive noise, and one is interested in breaking down the measured signal into two parts; the contribution of the sources of interest (SOI) and noise. In this case, source separation can be called source extraction. The sound source extraction/separation problem is a typical example where the source field is remotely sensed using a microphone array [2]. This problem could be solved either in a supervised or a blind manner, using the coherence/reference techniques [3] or blind source separation [4], respectively.

In source extraction using remotely-sensed multi-sensor data, where the signals are modeled as a linear combination of the contributions of some uncorrelated sources corrupted by additive noise [5], the accurate estimation of the source contributions is conditioned by the accurate source enumeration, i.e., the accurate estimation of the number of incoherent sources (the stochastic dimension of the field) [6]. Whether a supervised or a blind source separation approach is pursued, it is always essential the estimation of the number of sources for the separation results to be correct. In supervised source separation techniques, the data are available in the form of two sets of signals, namely the array signals (the targets) and the so-called reference signals. The reference signals are believed to carry more information about the sources than the array signals. This could be due to a better position (being closer to the sources), due to being of higher signal-to-noise ratio (SNR), or simply due to being of lower number than the array signal, which makes their cross-spectral matrix (CSM) better conditioned than the array signals. Assuming that the noise on these reference signals are uncorrelated to the noise on the array, the coherent part between the two sets is then coming from the sources. Hence, the coherence/reference signal processing techniques could be used (using the reference signals as inputs and the array signals as outputs) to quantify the part of the array signals coherent to the sources (which represents the source contributions to the measurements, later denoted as $X$) [7]. This is conventionally done using a least squares (LS) solution. Considering the necessity of using reference sensors of a number equal to or larger than the number of sources, in
order to ensure capturing the stochastic dimension of the field [8], a reduced-rank
solution is needed, otherwise the above technique become numerically invalid due
to a badly-conditioned reference CSM. Techniques in the literature (e.g., the par-
tial coherence [9, 10], the virtual coherence [11, 12] and the maximally-coherent
reference (MCR) techniques [13, 14]) attempted to solve this problem, providing
a reduced-rank solution. The rank-reduction step, involved in all these solutions,
requires either a priori knowledge or the estimation of the true number of sources
(source enumeration). This source enumeration problem, in the current context,
is an ill-posed inverse problem, which is insufficiently investigated in the context
of coherence/reference techniques, what motivates, in this paper, the presentation
of three alternative solutions for it. For these three source enumeration solutions
to be universal, in the sense that, they are applicable to all the reduced-rank ref-
erence techniques, a unified formalism for all the reference techniques is provided.

The source enumeration problem, in the context of the reference techniques,
is poorly addressed in the literature, compared to that in the context of other tech-
niques, such as blind source separation, which is based on a similar, but not iden-
tical model (a linear mixture model) [2, 4, 15–19]. For nearly all the literature
timeline of the reduced-rank reference techniques, the cut-off number of partial
(or virtual) references was solely found based on the value of multiple coherence
of the preceding partial (or virtual) references approaching unity. This method
showed to fail, for example, in the case of aeroacoustic sources [20], or in the case
of a CSM with gradually decreasing eigenvalues. To solve this, Leclere et al. [21]
established a threshold using a statistical hypothesis test based on the probability
law of the estimated coherence between uncorrelated Gaussian signals given by
Carter and Nuttall [22]. Solving the source enumeration problem using statistical
inference is a good idea, as the problem is equivalent to a model order selection
problem. Hence, to choose which statistical techniques to use, the model order
selection (model identification) literature is reviewed hereafter.

The literature on model order selection (source enumeration) is vast and dates
back to the emergence of the communication, the detection and the information
theories during the first half of the twentieth century. The solutions provided dur-
ing that period were mainly based on the maximum likelihood estimation (MLE).
This MLE-based solutions can be further categorized under two main branches,
namely; the likelihood ratio testing method [23–25] and the information criteria
(IC) method [26–28]. The LRT method consists of a sequence of hypothesis tests
(i.e., LRTs) based on the likelihood ratios between each two consecutive mod-
els while replacing unknown parameters by their maximum likelihood estimates.
The IC method is based on selecting the model which yields the minimum of an information criterion which is in turn a function of the log-likelihood of the maximum-likelihood estimator of the parameters of the model. The IC method is then an extension of the maximum likelihood principle [26], that was advised as an alternative to implementing LRTs [29], to avoid the struggle of setting a threshold in the cases where the statistical specification of the data is not clear, and hence, a subjective judgment is needed to set a threshold [30, 31]. The IC method represents a good alternative to LRTs also in the cases when no sufficient data samples are available [32]. However, when sufficiently long data instances are available, and when the statistical specification of the data is known (e.g., its statistical model), the inference step becomes relatively easy, in terms of finding a threshold, especially if the data is describable with a parametric model, in this case the only remaining difficulty is to count the degrees of freedom (DOF) (i.e., number of free parameters) for each specific model under consideration. In these cases, the LRTs are simple, independent and sufficient, and opting to the IC method is computationally and statistically complicated due to the consisting in statistically correlated tests associated with minimizing the IC [33, 34]. The source enumeration problem in the rank-reduction framework is described by a nested model, hence, based on Wilks' theorem [23, 35], 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. However, this theorem showed often not to hold for factor models, which are models describing the data in terms of a number of latent factors smaller than the number of measurements (such as the one considered in this work), due to violating certain regularity conditions (such as, parameters identifiability with respect to the likelihood function and non-singularity of the Fisher information matrix). Such models possess irregularity points in their parameter space, i.e., specific parameter values or conditions where the regularity assumptions, required for the LRT, do not hold: e.g., parameters taking values on the boundary of its parameter space, leading to a possibly singular Fisher information matrix, which invalidates the use of LRT. Specifically, such issues arise when comparing models with different number of factors rather than comparing a factor model against the saturated model (i.e., a model with a number of factors equal to the number of measurements) [36].

For cases where the likelihood function is not tractable, or where the data samples are insufficiently few for an asymptotic framework (i.e., very short records in the signal processing context), one can resort to another main stream of solutions (rather say techniques) which was evolved with the advancement in the computational power, during the second half of the twentieth century. These techniques
are aggregately called the computer-intensive techniques (also called resampling techniques), and were based on the idea of resampling (the reuse of the samples). The literature on this main stream is also vast and contain many techniques, that range from the non-parametric (cross-validation, jackknife, and non-parametric bootstrap) to the parametric ones (parametric bootstrap), which evolved from one another. They all started from the idea of cross-validation [37], passing through Quenouille-Tukey's jackknife [38, 39], and ending by the emergence of the different bootstrapping techniques [40].

The cross-validation, in its simplest forms, consists in splitting the data into two sets (an inference set and a validation set), the estimation step is first done using the inference set, then the validation is done by applying the estimated statistics (in the current context, the estimated number of sources) to the validation set and then evaluating some error. The earliest reported application of this technique was to evaluate the shrinkage of the coefficient of multiple correlation (the R coefficient) [37] leading to the formulation of an adjusted R coefficient [41]. After years of distinct applications, the term “Cross-Validation” had started to become loose, therefore it had to be further well-defined and distinguished from similar statistical validation techniques, what was done in the work of Mosier [42]. An important corner stone in the evolution journey of the resampling techniques is the development of the Quenouille-Tukey's jackknife, initially presented by Quenouille [38] and then extended by Tukey [39]. This technique is done through systematically leaving out one observation at a time from the dataset to create multiple subsets (jackknife samples), the statistics are then estimated for each jackknife sample. The estimates from these jackknife samples are then used to calculate the jackknife estimate of the parameter of interest. Efron [40] generalized the jackknife method into a more robust method, which he coined “Bootstrap”. It is done by repeatedly generating, a large number of times (typically thousands of times), multiple samples (bootstrap samples) with replacement from the observed data. The statistics of interest are then calculated for each bootstrap sample. In comparison to jackknife, in bootstrap, the same data point can be selected multiple times in a single bootstrap sample, and some data points may not be selected at all. Using a comparative study, Efron [40] showed the supremacy of the bootstrap method over both the cross-validation and the jackknife techniques.

All the resampling techniques mentioned above (cross-validation, jackknife, and bootstrap) despite of being efficient, are non-parametric, and it is known that a good parametric analysis, when parametric assumptions are easy to justify, can be far more efficient than its non parametric counterpart. Hence, it was the development of the parametric bootstrap very early after the (non-parametric) original
bootstrap invention [43]. It involves resampling by generating new datasets from an assumed parametric distribution, using the covariance structure estimated from the measured dataset. Hence, if the data was proved to follow a specific parametric model, it is more appropriate to use a parametric bootstrap strategy.

Hence, this work provides three distinct solutions for the source enumeration problems, in the context of reference techniques, which span the reviewed literature, ranging from parametric to non-parametric, and from asymptotic to exact (computer-intensive) solutions:

**An MLE-based solution (asymptotic: i.e., based on an approximation of the true Probability Density Function (PDF)):**

1. A modified LRT (parametric).

**and two computer-intensive solutions (exact: i.e., based on draws from the true PDF):**

2. Cross-validation (non-parametric).

3. Parametric bootstrap (parametric).

Thanks to the unified formalism of the reference techniques provided in this paper, the three presented source enumeration solutions are equally applicable to all the reduced-rank reference techniques in the literature. However, a single reference technique has to be selected as a datum for the analysis and comparison between the presented source enumeration methods, which is the MCR technique, recently provided by the authors, which was formulated as an optimal solution for the source extraction problem, and it succeeded to overcome the inherent shortcomings of the conventional coherence-based solutions (namely, the partial and the virtual coherence) [13, 14]. The performance of the presented source enumeration techniques is assessed and compared using both numerical simulations, laboratory experiments and electric motor (e-motor) data.

The structure of the rest of the paper is as follows. Section 2 defines the problem. The unified formalism for the reference techniques is presented in Section 3. The presentation of the different source enumeration algorithms, and assessing and comparing their performances using simulations are in Section 4. The laboratory experiments and its results are presented in Section 5. Comparison of the techniques based on e-motor data is in Section 6. Section 7 concludes the paper. Detailed mathematical derivations are presented in the appendices.
## 2. The source extraction problem using references

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

$$ \mathbf{X} = \mathbf{HS}, \qquad (1) $$

for some linear but unknown operator $\mathbf{H} \in \mathbb{C}^{M \times S}$, which is assumed of full rank (i.e., $\text{rank}(\mathbf{H}) = S$) without loss of generality (otherwise there would be the effect of less than $S$ sources in $\mathbf{X}$), and corrupted by additive disturbances $\mathbf{N} \in \mathbb{C}^{M \times I}$ uncorrelated with $\mathbf{X}$, such that

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

It is assumed that the number of output sensors $M$ exceeds or is equal to the number of sources, i.e., $M \ge S$.

The data $\mathbf{Y}$ comes with a set of $R$ references, $\mathbf{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 $\mathbf{L}$ such that

$$ \mathbf{R} = \mathbf{LS}, \qquad (3) $$

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

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

with $\epsilon \ge 0$. Given that, $\mathbf{e}$ is a single realization/snapshot (i.e., one column vector) of $\mathbf{E}$, the noise covariance matrix, $\Sigma_{\mathbf{E}} = \mathbb{E}\{\mathbf{e}\mathbf{e}^{\mathrm{T}}\}$, is not necessarily proportional to the identity. It is also 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 $\mathbf{X}$ from the references $\mathbf{R}$, i.e., to find an operator $\mathbf{G} \in \mathbb{C}^{M \times R}$ such that

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

is an estimate of $\mathbf{X}$. The simplest solution to the above problem is the LS solution

$$ \mathbf{G} = \underset{\mathbf{A}}{\operatorname{Arg min}} \|\mathbf{Y} - \mathbf{A}\mathbf{R}\|^2 = \mathbf{S}_{\mathrm{YR}}\mathbf{S}_{\mathrm{RR}}^{-1}, \qquad (6) $$
with $S_{RR} = RR^H/I$ and $S_{YR} = YR^H/I$.

Unfortunately, the matrix $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$, hence, its inversion is numerically problematic. The solutions in the literature for this problem are: a subset LS [44], the truncated singular value decomposition (TSVD), also known as virtual coherence (i.e., LS on the principal components (PC) of $S_{RR}$) [11], the partial coherence method [9, 45], and the MCR method [13, 14]. All these solutions can be summarized by avoiding the problematic inversion of $S_{RR}$ by finding a reduced-rank factorization of the operator $G$ in Eq. (5). For the clarification of this concept, a unified formalism that unites all the solutions in the literature (including the LS solution) is presented in the next section. This unified formalism allows for the design of a general statistical test, which is necessary for the reduced-rank solution, applicable to all the solutions.

### 3. A unified formalism for all the reduced-rank reference techniques in the literature

All the reduced-rank reference techniques in the literature can be generalized by factorizing the operator $G$, in Eq. (5), as $G = CP^H$, with a weighting matrix $C \in \mathbb{C}^{M \times R}$, and a projection matrix $P \in \mathbb{C}^{R \times R}$ for the LS solution (or any other solution evaluated for the (so-called saturated) model of full-rank equals $R$, i.e., keeping all the $R$ components), and with $C \in \mathbb{C}^{M \times S}$, and $P \in \mathbb{C}^{R \times S}$, for the reduced-rank solution, where $S < R$ is the number of sources (i.e., discarding $R-S$ noise-related components). Hence,

$$ \hat{X} = CP^H R. \qquad (7) $$

For identifiability of the solution, a scaling is placed on the matrix $P$ such that,

$$ P^H MP = I, \qquad (8) $$

for any positive definite matrix $M_{R \times R}$. Given $P$, and by defining the quadratic norm of a matrix $Q$ as $\|Q\|_I^2 = \text{tr}\{Q^H \Gamma Q\}$ for some positive definite matrix $\Gamma \in \mathbb{C}^{M \times M}$, $C$ can be found as

$$ C = \arg \min_{A} \|Y - AP^H R\|_I^2 = S_{YR} P (P^H S_{RR} P)^{-1}. \qquad (9) $$

Hence, the operator $G$ becomes

$$ G = S_{YR} P (P^H S_{RR} P)^{-1} P^H, \qquad (10) $$

which is interpreted, for the reduced-rank solution case, as using the stable projector $P(P^H S_{RR} P)^{-1} P^H$ to avoid the problematic inversion of the potentially unstable $S_{RR}$.
### 3.1. Solutions in the literature as particular cases

The solutions provided in the literature (including the LS solution) actually are special cases of the presented general solution in Eqs. (7) and (8), in essence, each solution has different choice for the matrices $\mathbf{P}^H$ and $\mathbf{M}$, as follows.

1) LS: $\mathbf{P}_{R \times R}^H = \mathbf{I}_{R \times R}$ and $\mathbf{M} = \mathbf{I}$.

2) Subset LS: $\mathbf{P}_{S \times R}^H = \tilde{\mathbf{I}} \mathbf{I}_{R \times R}$, where $\tilde{\mathbf{I}} = [\mathbf{I}_{S \times S} \quad \mathbf{O}_{S \times (R-S)}]$ is a selection matrix, and $\mathbf{M} = \mathbf{I}$.

3) TSVD: $\mathbf{P}_{S \times R}^H = \tilde{\mathbf{I}} \mathbf{U}_{R \times R}^H$, where $\mathbf{U}$ is the matrix containing the eigenvectors associated to the eigenvalues, sorted in descending order, of $\mathbf{S}_{RR}$, and $\mathbf{M} = \mathbf{I}$.

4) Partial coherence: $\mathbf{P}_{S \times R}^H = \tilde{\mathbf{I}} \mathbf{T}^{-1}$, where $\mathbf{T}_{R \times R}$ is the Cholesky factor of $I \cdot \mathbf{S}_{RR}$, with $I$ the number of snapshots, and $\mathbf{M} = I \cdot \mathbf{S}_{RR} = \mathbf{R}\mathbf{R}^H$.

To briefly explain where this came from, Let $\mathbf{R} = \mathbf{T}_{R \times R}\mathbf{Q}_{R \times I}$ be the LU decomposition of $\mathbf{R}$, with $\mathbf{T}$ lower triangular and $\mathbf{Q}$ unitary, then, $\mathbf{T}^{-1}\mathbf{R} = \mathbf{Q}$. Hence, the matrix $\mathbf{T}^{-1}$ orthogonalizes $\mathbf{R}$. It also follows that, $\mathbf{R}\mathbf{R}^H = I \cdot \mathbf{S}_{RR} = \mathbf{T}\mathbf{T}^H$ is the Cholesky decomposition of $I \cdot \mathbf{S}_{RR}$.

5) MCR: $\mathbf{P}_{S \times R}^H = \mathbf{B}^H$, where $\mathbf{B} = [\mathbf{b}_1 \quad \dots \quad \mathbf{b}_S]$ is the matrix containing the generalized eigenvectors associated to the $S$ largest generalized eigenvalues, sorted in a descending order, of the pair of matrices $(\mathbf{S}_{RY}\mathbf{\Gamma}\mathbf{S}_{YR}, \mathbf{S}_{RR})$, with $\mathbf{\Gamma}$ is any positive definite matrix (usually taken as the identity $\mathbf{I}$ or the inverse of the CSM of the outputs $\mathbf{S}_{YY}^{-1}$), and $\mathbf{M} = \mathbf{S}_{RR}$ [46].

It is clear that, the solutions (3) to (5) are different orthogonalizations of the references $\mathbf{R}$, which result in a new set of orthogonal (uncorrelated) references $\tilde{\mathbf{R}} = \mathbf{P}^H\mathbf{R}$, so that, $\tilde{\mathbf{R}}\tilde{\mathbf{R}}^H$ is diagonal.

It is also clear that, for all the reduced-rank solutions from (2) to (5) to be correct, it is essential the knowledge of the true number of sources $S$, however in the absence of this information, these solutions can be evaluated for any model of an arbitrary assumed rank ($K$), where $K = 0, 1, 2, \dots, R$, then using the data, the model fitness can be judged using any of the proposed source enumeration techniques, i.e., $S$ can be estimated.

### 3.2. The degrees of freedom in a reduced-rank model

When evaluating the fitness of a reduced-rank model, the model DOF, which is the number of free parameters (i.e., parameters whose values are allowed to
vary), changes with the model order (rank), affecting the LRT results, and hence,
the correct estimation of S. Thus, it has to be accurately counted. For the general
solution in Eq. (7), and regardless of the particular cases, the DOF of the reduced-
rank model of rank S (i.e., of which $dim(C) = M \times S$ and $dim(P) = R \times S$) is
counted as follows (note that, every restriction imposed by a constraint will be
assigned a negative sign, i.e., will be counted as -1 free parameter).

* $CP^H$ contains a total of $(M + R)S$ complex elements, which contributes by $2(M + R)S$ free parameters.

* These elements are not totally free, as they are subjected to the constraint $P^HMP = I_{S \times S}$, which imposes $-2(S^2 - S)$ free parameters, from the complex (off-diagonal) elements, and $-S$ from the real (diagonal) ones. However, due to the hermitian property of $P^HMP$, there are only $-(S^2 - S)$ free parameter from the complex elements, along with the $-S$ from the real ones, thus, we now have a total of $-S^2 + S - S = -S^2$ parameter imposed by the constraint.

* Recognizing the fact that the same constraint affects the matrix $C$ as $C = S_{YR}P(P^H S_{RR}P)^{-1}$, i.e., the number and values of the elements in $P$ predetermines the number and values of the free parameters in $C$, then the correct total number of free parameter imposed by the constraint is $-2S^2$, not only $-S^2$.

Adding the positive and negative free parameters from the previous discussion,
the total DOF of the reduced model reads DOF=$2(M + R)S - 2S^2$, this is when
the saturated model (the model of a full-rank equals $R$) is reduced to a model
of a rank equal to the true number of sources $S$ (the correct model). Yet, this
formula similarly applies to any model with an arbitrary rank ($K$), where $K =
0, 1, 2, \dots, R$, just by replacing $S$ with $K$.

**4. Selection of the number of sources**

As it can be seen from the previous section that, a priori knowledge of the true
number of sources (S) is needed, or it needs to be estimated. In practice, the num-
ber of sources is rarely known and must be estimated from the data, which is a hard
ill-posed inverse problem. Hence, three source enumeration techniques, namely, a
modified LRT, a parametric bootstrapping technique and a cross-validation tech-
nique are presented, validated and compared by numerical simulations in the fol-
lowing subsections. Note that, even if all the algorithms are formulated to be
equally applicable to all the reduced-rank reference techniques, however, for the sake of implementation, all the following results (simulations and real data) are obtained based on the MCR technique.

### 4.1. The likelihood ratio test

An asymptotic parametric approach that is usually used for source enumeration is the LRT, which compares the fit of nested models based on Wilks' theorem [23]. However, Wilks' theorem, which states that the LRT statistic (i.e., the Log-Likelihood Ratio (LLR)) 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, shows often not to hold for factor models, because of the existence of irregularity points in the models parameter space that in some cases invalidate some regularity conditions of the theorem (as discussed in the introduction). However, it was noted by Chen et al. [36] that, comparing each sub-model against the saturated model, rather than comparing each two successive sub-model, could solve this issue, as the LLR in this case follows a $\chi^2$ distribution under milder conditions. Hence, a modified LRT is proposed here based on comparing against the saturated model. In the following subsections, the conventional and the modified LRTs are detailed, their respective statistics are derived for the case of complex data, and their performances are evaluated and compared by numerical simulations.

#### 4.1.1. The conventional LRT between two successive models involving complex data

One assumes the presence of $K$ sources and wants to test whether a model with $K + 1$ is more likely using the LRT, based on the LLR between the compared models

$$ \Lambda(K) = \ln \frac{P(Y|\hat{\theta}_{K+1})}{P(Y|\hat{\theta}_K)}, \qquad (11) $$

where $\hat{\theta}_i$ is the maximum-likelihood estimates of the free parameters in the model with $i$ number of sources.

The estimate of source contributions for the model of $K$ sources is obtained from any of the reduced-rank reference techniques, truncating after the $k^{th}$ component, as follows.

$$ \hat{\mathbf{X}}(K) = \mathbf{C}(K)\mathbf{P}(K)^H\mathbf{R}, \qquad (12) $$

as returned by Eq. (7), but using $\mathbf{C}(K) = [\mathbf{c}_1 \cdots \mathbf{c}_K]$, instead of the full $\mathbf{C}_{M \times R}$ (denoted hereafter as $\mathbf{C}(R)$), and $\mathbf{P}(K) = [\mathbf{p}_1 \cdots \mathbf{p}_K]$, instead of the full $\mathbf{P}_{R \times R}$
(denoted hereafter as $\mathbf{P}(R)$). The covariance matrix of the $K$ model is then estimated as

$$ \hat{\Omega}(K) = (\mathbf{Y} - \hat{\mathbf{X}}(K))(\mathbf{Y} - \hat{\mathbf{X}}(K))^{\mathrm{H}}/I. \quad (13) $$

For the model of $K + 1$ sources, the estimate of source contributions is

$$ \hat{\mathbf{X}}(K + 1) = \mathbf{C}(K + 1)\mathbf{P}(K + 1)^{\mathrm{H}}\mathbf{R}, \quad (14) $$

as returned by Eq. (7), but using $\mathbf{C}(K + 1) = [\mathbf{c}_1 \cdots \mathbf{c}_{K+1}]$, instead of the full $\mathbf{C}(R)$, and $\mathbf{P}(K + 1) = [\mathbf{p}_1 \cdots \mathbf{p}_{K+1}]$, instead of the full $\mathbf{P}(R)$. Then the covariance matrix of the $K + 1$ model is estimated as

$$ \hat{\Omega}(K + 1) = (\mathbf{Y} - \hat{\mathbf{X}}(K + 1))(\mathbf{Y} - \hat{\mathbf{X}}(K + 1))^{\mathrm{H}}/I. \quad (15) $$

Given that the data are complex, which resembles Fourier coefficients, then $\mathbf{Y}$ is assumed to follow a multivariate complex Gaussian, according to the central limit theorem [47], with mean $\hat{\mathbf{X}}(K)$ (resp. $\hat{\mathbf{X}}(K + 1)$), and covariance matrix $\hat{\Omega}(K)$ (resp. $\hat{\Omega}(K + 1)$) under the first (resp. second) model.

Then, the likelihood expression for the first model reads

$$ P(\mathbf{Y}|\hat{\theta}_K) = \frac{1}{\pi^{MI}|\hat{\Omega}(K)|^I} \exp(-\mathrm{tr}\{(\mathbf{Y} - \hat{\mathbf{X}}(K))^{\mathrm{H}}\hat{\Omega}(K)^{-1}(\mathbf{Y} - \hat{\mathbf{X}}(K))\}), \quad (16) $$

and the likelihood expression for the second model reads

$$ P(\mathbf{Y}|\hat{\theta}_{K+1}) = \frac{1}{\pi^{MI}|\hat{\Omega}(K+1)|^I} \exp(-\mathrm{tr}\{(\mathbf{Y} - \hat{\mathbf{X}}(K+1))^{\mathrm{H}}\hat{\Omega}(K+1)^{-1}(\mathbf{Y} - \hat{\mathbf{X}}(K+1))\}), \quad (17) $$

substituting in LLR expression (11), then the LLR of the second model (of $K + 1$ sources) to the first model (of $K$ sources) reads

$$ \Lambda(K) = -\|\mathbf{Y} - \hat{\mathbf{X}}(K+1)\|_{\hat{\Omega}(K+1)}^2 + \|\mathbf{Y} - \hat{\mathbf{X}}(K)\|_{\hat{\Omega}(K)}^2 - I \ln |\hat{\Omega}(K+1)| + I \ln |\hat{\Omega}(K)| \quad (18) $$

where $\|\mathbf{A}\|_{\hat{\Omega}}^2 = \mathrm{tr}\{\mathbf{A}^{\mathrm{H}}\hat{\Omega}^{-1}\mathbf{A}\}$. Note that the $I$ snapshots used in the previous equation have to be independent (uncorrelated), however during the estimation of the Power Spectral Density (PSD), weighting and overlapping between the time windows are performed to avoid spectral leakage, this leads to having correlated snapshots. In this situation Eq. (18) becomes

$$ \Lambda(K) = -(I_{\mathrm{eff}}/I)\|\mathbf{Y} - \hat{\mathbf{X}}(K+1)\|_{\hat{\Omega}(K+1)}^2 + (I_{\mathrm{eff}}/I)\|\mathbf{Y} - \hat{\mathbf{X}}(K)\|_{\hat{\Omega}(K)}^2 \\ - I_{\mathrm{eff}} \ln |\hat{\Omega}(K+1)| + I_{\mathrm{eff}} \ln |\hat{\Omega}(K)| \quad (19) $$
with $I_{\text{eff}}$ the effective number of equivalent uncorrelated snapshots, as calculated by Welch [48]:

$$I_{\text{eff}} = I/\{1 + 2 \sum_{i=1}^{I-1} \frac{I-i}{I} \rho(i)\}, \qquad (20)$$

where $\rho(i)$ is the correlation function defined, in the case of a Gaussian process, as

$$\rho(i) = \left[ \sum_{n=0}^{L-1} W(n)W(n + iD_{sh}) \right]^2 \left[ \sum_{n=0}^{L-1} W^2(n) \right]^2, \qquad (21)$$

where $W(n)$ is the weighting time window, $L$ is the window length, and $D_{sh} = L - \text{overlap}$ is the shift between the overlapping segments.

For the mathematical basis, Appendix A works out the equivalence between a set of correlated random variables and a reduced set of uncorrelated random variables.

The LLR expression in Eq. (19) could be simplified, assuming a full covariance matrix. Hence, the covariance estimates $\hat{\Omega}(K)$ and $\hat{\Omega}(K + 1)$ in this equation are taken as the estimates in Eqs. (13) and (15), respectively. Then, it follows that

$$\Lambda(K) = I_{\text{eff}}(\ln|\hat{\Omega}(K)| - \ln|\hat{\Omega}(K + 1)|). \qquad (22)$$

The above expression might still be difficult to evaluate numerically, because the computation of the logarithm of the determinant of large matrices tends to be unstable. However further simplification could be achieved under the MCR solution, by noting that

$$\hat{\Omega}(K + 1) = S_{YY} - \sum_{k=1}^{K+1} S_{YR} b_k b_k^H S_{RY} = \hat{\Omega}(K) - S_{YR} b_{K+1} b_{K+1}^H S_{RY}. \qquad (23)$$

Setting $\mathbf{t} = S_{YR}b_{K+1}$ and using the determinant lemma,

$$|\hat{\Omega}(K + 1)| = |\hat{\Omega}(K)| \cdot (1 - \mathbf{t}^H \hat{\Omega}(K)^{-1} \mathbf{t}) \qquad (24)$$

and so

$$\Lambda(K) = -I_{\text{eff}} \ln(1 - \mathbf{t}^H \hat{\Omega}(K)^{-1} \mathbf{t}). \qquad (25)$$

Wilk's LRT was originally defined for real random variables. One of the novelties of the present work is presenting the mathematical derivation in the case
of vectors of complex random variables, which is up to the best of the authors'
knowledge never introduced before in the literature. According to the mathe-
matical findings in Appendix B.1, the above LLR is asymptotically (with $I_{\text{eff}}$) dis-
tributed as $k \cdot \chi_{2\nu}^2$, where $k = \frac{1}{2}$, and $\chi_{2\nu}^2$ is a Chi-squared with $2\nu = 2(M+R)-4K-2$
degrees of freedom, as found in Appendix B.2. Therefore, one accepts the new
model with $K+1$ sources if $\Lambda(K) > \chi_{2\nu,1-\alpha}^2/2$, with $\chi_{2\nu,1-\alpha}^2$ the quantile of the
Chi-squared associated with probability $1-\alpha$. If so, the next model with $K+2$
is tested; if not, $K$ sources are selected. Algorithm 1 gives the procedure of the
selection of the number of sources.

For confirming the mathematical findings, $S = 4$ independent sources were
simulated, which contributes to $M = 50$ output sensors and $R = 11$ references
through a random transfer matrix, whose elements were randomly drawn from an
i.i.d. complex Gaussian random process with zero mean and unit variance.

**Algorithm 1** The conventional LRT

• Inputs: $\mathbf{Y}, \mathbf{R}, \mathbf{S}_{YY}, \mathbf{S}_{RR}$ and $\mathbf{S}_{YR}$

• Estimate $\mathbf{C}(\mathbf{R})$ and $\mathbf{P}(\mathbf{R})$ from $\mathbf{S}_{YY}, \mathbf{S}_{RR}$ and $\mathbf{S}_{YR}$ using a chosen reduced-rank reference technique, and keeping all the $R$ components.

• Set $K = 0$

• WHILE $K < R$

- Construct $\mathbf{C}(K)$ and $\mathbf{P}(K)$ from $\mathbf{C}(\mathbf{R})$ and $\mathbf{P}(\mathbf{R})$ assuming $K$ sources

- Estimate $\hat{\mathbf{X}}(K)$ and $\hat{\mathbf{\Omega}}(K)$ from Eqs. (12) and (13)

- Construct $\mathbf{C}(K+1)$ and $\mathbf{P}(K+1)$ from $\mathbf{C}(\mathbf{R})$ and $\mathbf{P}(\mathbf{R})$ assuming $K+1$ sources

- Estimate $\hat{\mathbf{X}}(K+1)$ and $\hat{\mathbf{\Omega}}(K+1)$ from Eqs. (14) and (15)

- Compute $\Lambda(K)$ from Eq. (22) for all the reduced-rank reference techniques (or alternatively from Eq. (25) only for the MCR technique, in this case the two previous steps are skipped).

- IF $\Lambda(K) < \chi_{2\nu,1-\alpha}^2/2$ STOP

- ELSE $K \leftarrow K + 1$

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

• ELSE $S \leftarrow K$

Uncorrelated complex Gaussian noise, with zero mean, was added to the out-
puts with random heterogeneous variance of $\sigma_M^2 = 10^{-2}$ order of magnitude.
Scalar complex Gaussian noise (i.e., $\Omega = \sigma_R^2 I$) was added to the references with variance $\sigma_R^2 = 10^{-6}$. The simulation was repeated for varying number of snapshots $M < I < 10^5$. The data was simulated 1k times. After applying the MCR technique, $\Lambda(S)$, i.e., the LLR between a model of $S + 1$ and the model of $S$ sources, was computed each time, ending up having a 1k sample of LLRs.

As one can see from Fig. 1(a), that the average value of a 1k sample of simulated LLR/$k$ (where $k$ is the constant in front of the $\chi_{2\nu}^2$) nearly coincided with the DOF computed as $2\nu = 2(M + R) - 4S - 2$, and this seems to be unaffected by the used number of snapshots $I$. Fig. 1(b) shows that the simulated value of the constant $k$, calculated as $k = \mathbb{E}\{\Lambda(K)\}/(2 \times \text{Var}\{\Lambda(K)\})$, converges to 0.5 given that a sufficient number of snapshots are used in the analysis (in this example, approximately $I > 1$k snapshot, when the outputs are 50 microphones). Figures 1(c), 1(d) and 1(e) show the histograms of the simulated LLR superimposed on the histogram of an equivalent sample drawn from $\frac{1}{2}\chi_{2\nu}^2$, at $I = 100$, 1k, 10k respectively. Although the mode (note that: for the $0.5\chi^2$ distribution, the mean = $0.5 \times \text{DOF}$, and the mode = $0.5 \times (\text{DOF}-2)$) of the simulated LLR nearly converged to that of the $\frac{1}{2}\chi_{2\nu}^2$ for relatively low number of snapshots ($I = 1$k), the histograms show shape nonconformity between the simulated LLR distribution and the $\frac{1}{2}\chi_{2\nu}^2$, even for the case of using a very high number of snapshots ($I = 10$k). What conforms with the findings of Chen et al. [36] that Wilk's theorem fails for factor models. The shape mismatch, especially in the upper tail, leads to accepting wrong hypotheses, thus returning wrong estimate of the number of sources.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Conventional LRT: (a) the average simulated value of LLR/k, where k is the constant in $k \cdot \chi_{2\nu}^2$, (in dashed red), as a function of the number of snapshots *I*, superimposed on the DOF as computed by the analytically-found formula (in black), (b) the constant *k*, as a function of the number of snapshots *I*, the simulated value (in dashed blue) and the theoretical value (in black), (c) the histogram of the simulated LLR (the hollow red) superimposed on $\frac{1}{2}\chi_{2\nu}^2$ (the solid violet) for data with *I* = 100, (d) *I* = 1k, and (e) *I* = 10k snapshots.</figcaption></figure>
### 4.1.2. The proposed modified LRT

As was pointed out by Chen et al. [36] that the LRT for comparing a K-sources model against the saturated model does follow a $\chi^2$-distribution under mild conditions. A modified LRT is proposed which successively tests models with $K = 0, 1, 2, ..., R-1$ sources against the saturated model with an unrestricted covariance matrix (i.e., with all the R components). In this case, the LLR becomes

$$ \Lambda(K) = \ln \frac{P(Y|\hat{\theta}_R)}{P(Y|\hat{\theta}_K)}. \qquad (26) $$

The estimate of source contributions for the model of R sources is obtained by

$$ \hat{\mathbf{X}}(R) = \mathbf{C}(R)\mathbf{P}(R)^H\mathbf{R}, \qquad (27) $$

as returned by Eq. (7), but using the full $\mathbf{C}(R) = \mathbf{C}_{M \times R}$, and the full $\mathbf{P}(R) = \mathbf{P}_{R \times R}$. The covariance matrix of the R model is then estimated as

$$ \hat{\Omega}(R) = (\mathbf{Y} - \hat{\mathbf{X}}(R))(\mathbf{Y} - \hat{\mathbf{X}}(R))^H / I. \qquad (28) $$

Following the same lines as in the previous subsection, and assuming having correlated snapshots, the LLR becomes

$$ \Lambda(K) = -(I_{\text{eff}}/I)\|\mathbf{Y} - \hat{\mathbf{X}}(R)\|_{\hat{\Omega}(R)}^2 + (I_{\text{eff}}/I)\|\mathbf{Y} - \hat{\mathbf{X}}(K)\|_{\hat{\Omega}(K)}^2 - I_{\text{eff}} \ln |\hat{\Omega}(R)| + I_{\text{eff}} \ln |\hat{\Omega}(K)|, \quad (29) $$

assuming a full covariance matrix, and substituting with the covariance estimates in Eqs. (13) and (28), it follows that

$$ \Lambda(K) = I_{\text{eff}}(\ln |\hat{\Omega}(K)| - \ln |\hat{\Omega}(R)|). \qquad (30) $$

To avoid numerical problems, when the $\hat{\Omega}$ matrix is large, the quantity $\ln |\hat{\Omega}|$ can be computed as $\sum_{m=1}^{M} \ln(\sigma_m)$, where $\sigma_m$ are the singular values of $\hat{\Omega}$. The quantity in Eq. (30) is asymptotically (with $I_{\text{eff}}$) distributed as $k \chi_{2\nu}^2$, where the constant $k = \frac{1}{2}$ as before, but now the formula for the DOF of the Chi-squared reads $2\nu = 2(R-K)(M-K)$, as analytically found in Appendix B.3. Therefore, one rejects all the models until finding the first model with K sources that satisfies $\Lambda(K) < \frac{1}{2}\chi_{2\nu,1-\alpha}^2$, with $\chi_{2\nu,1-\alpha}$ the percentile of the Chi-squared associated with probability $1-\alpha$. If so, K sources are selected. If not, the next model with $K+1$ is tested. Algorithm 2 explains the procedure of the selection of the number of sources.

Numerical simulations, with the same simulation parameters described in Section 4.1.1, were carried out to check the correctness of the DOF formula. As one
can see from Fig. 2(a), the average of a simulated 1k sample of LLR/k, for a reduced model of S sources obtained with the MCR technique, nearly coincides with the DOF computed as $2\nu = 2(R-S)(M-S)$, and this seems to be unaffected by the used number of snapshots I.

<caption>**Algorithm 2** The modified LRT (testing against the saturated model)</caption>

* **Inputs:** $\mathbf{Y}$, $\mathbf{R}$, $\mathbf{S}_{YY}$, $\mathbf{S}_{RR}$ and $\mathbf{S}_{YR}$
* Estimate $\mathbf{C}(R)$ and $\mathbf{P}(R)$ from $\mathbf{S}_{YY}$, $\mathbf{S}_{RR}$ and $\mathbf{S}_{YR}$ using a chosen reduced-rank reference technique, and assuming $R$ sources.
* Estimate $\hat{\mathbf{X}}(R)$ and $\hat{\mathbf{\Omega}}(R)$ from Eqs. (27) and (28)
* Set $K = 0$
* WHILE $K < R$
    * Construct $\mathbf{C}(K)$ and $\mathbf{P}(K)$ from $\mathbf{C}(R)$ and $\mathbf{P}(R)$ assuming $K$ sources
    * Estimate $\hat{\mathbf{X}}(K)$ and $\hat{\mathbf{\Omega}}(K)$ from Eqs. (12) and (13)
    * Compute $\Lambda(K)$ from Eq. (30).
    * IF $\Lambda(K) < \chi_{2\nu,1-\alpha}^2/2$ STOP
    * ELSE $K \leftarrow K + 1$
* IF $K = R$, $S \leftarrow R$,
* ELSE $S \leftarrow K$

It is clear from Fig. 2(b), that the constant $k$ in front of the $\chi_{2\nu}^2$ converges to 0.5 given that a sufficient number of snapshots is used in the analysis (in this example, for approximately $I > 1000$ snapshot, when the outputs are 50 microphones). This also can be seen from the histograms, in Figs. 2(c), 2(d) and 2(e), of the simulated LLR superimposed on the histogram of an equivalent sample drawn from $\frac{1}{2}\chi_{2\nu}^2$. Besides the convergence of mode of the simulated LLR distribution to that of the $\frac{1}{2}\chi_{2\nu}^2$, these histograms also show good shape agreement between the simulated LLR histogram and the histogram of the equivalent sample drawn from the $\frac{1}{2}\chi_{2\nu}^2$ (given sufficient number of snapshots $I > 1k$ is used). Hence, the fitness of the proposed $\frac{1}{2}\chi_{2\nu}^2$ model for the LLR is assured (analytically and by numerical simulations), given a sufficient number of snapshots is available, and so is the effectiveness of the proposed saturated LRT.
<figure><img src="image_4.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e)</figcaption></figure>

**Fig. 2.** The proposed modified LRT (against a saturated model): (a) the average simulated value of LLR/$k$, where $k$ is the constant in $k \cdot \chi_{2\nu}^2$, (in dashed red), as a function of the number of snapshots $I$, superimposed on the DOF as computed by the analytically-found formula (in black), (b) the constant $k$, as a function of the number of snapshots $I$, the simulated value (in dashed blue) and the theoretical value (in black), (c) the histogram of the simulated LLR (the hollow red) superimposed on $\frac{1}{2}\chi_{2\nu}^2$ (the solid violet) for data with $I = 100$, (d) $I = 1k$, and (e) $I = 10k$ snapshots.
4.1.3. Comparing the proposed modified LRT to the conventional LRT

A simulated experiment was performed to compare the performance of the
modified LRT to the conventional LRT, where the same simulation parameters as
in Section 4.1.1 were used, but with varied number of outputs, M = 1 to 100, in
order to compare the capability of the two tests to accommodate a large number
of outputs. Algorithms 1 and 2 were used for the conventional and the modified
LRTs, respectively, both using the MCR technique.

Figures 3(a), 3(c) and 3(e) show the estimated number of sources $\hat{S}$ (in dashed red) using the conventional LRT, as a function of the number of outputs, along with the number of references $R$ (in black) and the true number of sources $S$ (in dash-dot blue), for numbers of snapshots $I = 500, 1500$ and $3000$ respectively. Figures 3(b), 3(d) and 3(f) show the corresponding estimated number of sources $\hat{S}$ using the modified LRT. The relative errors in the estimated coherent output, defined as $e(\hat{S}) = \|\hat{X}(\hat{S}) - X\|_F^2/\|X\|_F^2$, when keeping only $\hat{S}$ components in the MCR solution, as returned by the conventional and the modified tests, for only the case of $I = 3000$, are displayed in Figs. 3(g) and 3(h) respectively, along with the error $e(R) = \|\hat{X}(R) - X\|_F^2/\|X\|_F^2$ returned by keeping all the MCR components (R), and the error $e(S) = \|\hat{X}(S) - X\|_F^2/\|X\|_F^2$ obtained by taking the true number of sources (S) in MCR. The perfection of the results of the proposed modified LRT, given the availability of a sufficiently large number of snapshots, is clear in Figs. 3(f) and 3(h). Whereas the results of the conventional LRT overestimates the number of sources, even for the case of a few number of outputs $M$, and no improvement is witnessed by increasing the number of snapshots, even when $I = 3000$ (in Figs. 3(e) and 3(g)), which further proves that the statistic of the conventional LRT, in Eq. (18), does not asymptotically follow the Chi-squared distribution defined in Section 4.1.1. This failure of fit is caused by the failure of Wilk's theorem in the case of factor analysis models. The proposed modified statistic in Eq. (29) proved to asymptotically follow the Chi-squared distribution defined in Section 4.1.2. However, for allowing the proposed statistical test to work well while accommodating a large number of outputs, a very large number of snapshots has to be used.

This simulation was repeated four times for four different references' noise
characteristics, with increasing complexity. The first time was with a scalar noise
on the references, the second time with uncorrelated noise with heterogeneous
variance, the third time with correlated noise (full noise covariance matrix) with
homogeneous variance, and the last time with correlated noise with heterogeneous
variance.
<figure><img src="image_9.png" /><figcaption>**Fig. 3.** Comparing the simulation results of the modified LRT to those of the conventional LRT: (a), (c), and (e) show the estimated number of sources $\hat{S}$ (in dashed red) by the conventional LRT, as a function of the number of outputs ($M$), along with the number of references $R$ (in black) and the true number of sources $S$ (in dash-dot blue), for number of snapshots of $I = 500, 1500$ and $3000$, respectively, (b), (d), and (f) show the corresponding estimated number of sources $\hat{S}$ by the modified LRT, and (g) and (h) show, for the case of $I = 3000$, the corresponding relative errors in the estimated coherent output, as returned by the conventional and the modified tests, respectively, when keeping $\hat{S}$ components (in dashed red), along with the errors returned when keeping $R$ (in black) and $S$ (in dash-dot blue) components.</figcaption></figure>
The four input noise scenarios yielded the same results as in Fig. 3, which indicates the robustness of the proposed modified LRT, given the availability of a sufficiently high number of snapshots (a long record). In practical situations where only length-limited records are available, the frequency resolution has to be sacrificed in order to divide the signal into a sufficient number of snapshots for the asymptotic LRT to work, or one could resort to computer-intensive solutions (e.g., bootstrapping and cross-validation), however with higher cost on time, which are discussed in the following subsections.

## 4.2. The parametric bootstrap

An alternative technique for source enumeration is the parametric bootstrap, where a large number of data are bootstrapped under the restricted (reduced) model, and the corresponding LRT statistics is computed, yielding a distribution for the LRT statistic. The parametric bootstrap has the merit of better approximation for finite sample data than the asymptotic distribution. This would provide a good alternative when only short records are available or when a high frequency resolution, that can not be achieved using the direct LRT, is required. Hence, two parametric bootstrap algorithms, one based on the conventional LLR, and the other one based on the modified LLR, are detailed and validated by simulations in the following subsections.

### 4.2.1. Parametric bootstrap of the conventional LRT

Bootstrapping the distribution of the conventional LLR, defined in Eq. (19), can be considered as a fix for the problematic conventional LRT. The quantity in Eq. (22) has a reference distribution, which is denoted hereafter as $D(K)$, that could be computed using parametric bootstrap under the restricted (reduced) model (i.e., the one with $K$ sources), through estimating the covariance $\hat{\Omega}(K)$ using Eq. (13), and then synthesizing $I_{\text{eff}}$ independent vectors of $M$ complex random variables (i.e., bootstrapped output signals), with zero means and imposed covariance matrix $\hat{\Omega}(K)$. This is done through computing the singular value decomposition (SVD) of $\hat{\Omega}(K) = U\Sigma V^H$, then computing the new noise signals as $N_{\text{new}} = U\Sigma^{1/2}N_1$, where the elements of $N_1$ are drawn from a complex Gaussian distribution with zero mean and unit variance. A bootstrapped set of references $R_{\text{new}}$ is synthesized in the same manner, starting from the SVD of the $S_{RR}$ matrix. Finally the bootstrapped output data matrix is computed as $Y_{\text{new}} = C(K)P(K)^H R_{\text{new}} + N_{\text{new}}$. Hence, one accepts the new model with $K + 1$ sources if $\Lambda(K) > D(K)_{1-\alpha}$, with $D(K)_{1-\alpha}$ the quantile of the bootstrapped distribution associated with probability $1 - \alpha$. If so, the next model with $K + 2$ is
tested; if not, *K* sources are selected. Algorithm 3 gives in detail the steps for the selection of the number of sources.

**Algorithm 3 Bootstrapping the conventional LLR**

* Inputs: $\mathbf{Y}, \mathbf{R}, \mathbf{S}_{YY}, \mathbf{S}_{RR}$ and $\mathbf{S}_{YR}$

* Estimate $\mathbf{C}(\mathbf{R})$ and $\mathbf{P}(\mathbf{R})$ from $\mathbf{S}_{YY}, \mathbf{S}_{RR}$ and $\mathbf{S}_{YR}$ using a chosen reduced-rank reference technique, and keeping all the $\mathbf{R}$ components

* Set $K = 0$

* WHILE $K < R$

- Construct $\mathbf{C}(\mathbf{K})$ and $\mathbf{P}(\mathbf{K})$ from $\mathbf{C}(\mathbf{R})$ and $\mathbf{P}(\mathbf{R})$, assuming $\mathbf{K}$ sources

- Estimate $\hat{\mathbf{X}}(\mathbf{K})$ and $\hat{\Omega}(\mathbf{K})$ from Eqs. (12) and (13)

- Construct $\mathbf{C}(\mathbf{K} + 1)$ and $\mathbf{P}(\mathbf{K} + 1)$ from $\mathbf{C}(\mathbf{R})$ and $\mathbf{P}(\mathbf{R})$, assuming $\mathbf{K} + 1$ sources

- Estimate $\hat{\mathbf{X}}(\mathbf{K} + 1)$ and $\hat{\Omega}(\mathbf{K} + 1)$ from Eqs. (14) and (15)

- Compute $\Lambda(\mathbf{K})$ from Eq. (22) for all the reduced-rank reference techniques (or alternatively from Eq. (25) only for the MCR technique, in this case the two previous steps are skipped).

- Set $i = 1$

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

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

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

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

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

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

* Estimate $\hat{\mathbf{X}}(\mathbf{K})_{\text{new}}$ and $\hat{\Omega}(\mathbf{K})_{\text{new}}$ from Eqs. (12) and (13) applied to $\mathbf{C}(\mathbf{K})_{\text{new}}, \mathbf{P}(\mathbf{K})_{\text{new}}$ and $\mathbf{R}_{\text{new}}$

* Construct $\mathbf{C}(\mathbf{K} + 1)_{\text{new}}$ and $\mathbf{P}(\mathbf{K} + 1)_{\text{new}}$ from $\mathbf{C}(\mathbf{R})_{\text{new}}$ and $\mathbf{P}(\mathbf{R})_{\text{new}}$, assuming $\mathbf{K} + 1$ sources

* Estimate $\hat{\mathbf{X}}(\mathbf{K} + 1)_{\text{new}}$ and $\hat{\Omega}(\mathbf{K} + 1)_{\text{new}}$ from Eqs. (14) and (15) applied to $\mathbf{C}(\mathbf{K} + 1)_{\text{new}}, \mathbf{P}(\mathbf{K} + 1)_{\text{new}}$ and $\mathbf{R}_{\text{new}}$

* Compute $\Lambda(\mathbf{K})_{\text{new}}$ from Eq. (22) applied to $\hat{\Omega}(\mathbf{K})_{\text{new}}$ and $\hat{\Omega}(\mathbf{K} + 1)_{\text{new}}$, for all the reduced-rank reference techniques (or alternatively from Eq. (25) applied to $\hat{\Omega}(\mathbf{K})_{\text{new}}$, only for the MCR technique, in this case the two previous steps are skipped).

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

* $i \leftarrow i + 1$

- IF $\Lambda(\mathbf{K}) < D(\mathbf{K})_{1-\alpha}$ STOP where $D(\mathbf{K})_{1-\alpha}$ is quantile of $\{D(\mathbf{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$

In order to verify the algorithm, a numerical simulation was performed, using
a number of sources, $S = 4$, a number of output sensors, $M = 50$, a number of references, $R = 11$, and a number of snapshots, $I = 74$. The sources and the transfer functions were randomly drawn from i.i.d. complex Gaussian random processes with zero mean and unit variance. The noises were simulated as random complex Gaussian variables with zero mean, and their variances were specified so as to give output sensors' SNR$_{out}$ = 0 dB, on each output signal, and to give references' SNR$_{ref}$ = 30 dB, on each reference signal. The simulation was repeated twice, the first time using uncorrelated noises on the references, and the second time using 90% cross-correlated noises (i.e., with 0.9 correlation coefficient) on the references. A 99% confidence level was used in the statistical test, i.e., $\alpha = 0.01$.

The results obtained while using uncorrelated and correlated noises on the references turned out to be very similar; hence, only one of them is shown in Fig. 4. The figure consists of eleven sub-figures, each one representing a statistical test, for $K = 0, ..., 10$. For example, Fig. 4(a) displays the bootstrapped distribution $D(K = 0)$ of the LLR ($\Lambda(0)$) obtained under the restricted model with $K = 0$ source; hence, it can be used for comparing a model with $K + 1 = 1$ source to a model with no sources $K = 0$, i.e., a model that describes the data as pure noise. This is indicated in the figure subcaptions by $D(0) : Model(1)/Model(0)$. Similarly, Fig. 4(b) displays the distribution $D(1)$ and can be used to compare a model with $K = 2$ sources, *Model(2)*, to a model with only $K = 1$ source, *Model(1)*, ...etc. The red dotted vertical line is the current value of $\Lambda(K)$ to be tested, such that, one accepts the new model with $K+1$ sources if $\Lambda(K) > D(K)_{1-\alpha}$, and then the next model with $K + 2$ is tested. It is clear that, $\Lambda(K) > D(K)_{1-\alpha}$ up to Fig. 4(d), which shows the results of testing the model with 4 sources (the correct model) against the one with only 3 sources. However, as seen in Fig. 4(e), the value of $\Lambda(4)$ is approximately in the middle of the distribution $D(4)$, i.e., $\Lambda(4) << D(4)_{1-\alpha}$, hence, the model with 5 sources is very unlikely. Based on the previous discussion the model with 5 sources is rejected, thus the estimated number of sources $\hat{S} = 4$. Hence, the algorithm shows effectiveness in correctly estimating the true number of sources, as it returns $\hat{S} = S = 4$. It also shows an excellent accuracy, as it can be seen from the figure that, one could have got the same correct results using a lower confidence level, say 95%. Bearing in mind that, these excellent results were obtained using a relatively small number of snapshots (74 snapshots) compared to the number of sensors (50 output sensors). Any asymptotic-distribution-based test, such as the LRT, usually fails with such a low number of snapshots. However, the only drawback of bootstrapping is that it is more time consuming than asymptotic tests.
<figure><img src="image_10.png" /><figcaption>**Fig. 4.** Simulation results of bootstrapping the conventional LRT; comparison of $\Lambda(K)$ (red dotted line) with the bootstrapped distribution $D(K)$ (blue bar-chart) for the 11 virtual references (as returned by the MCR technique). The designation $Model(K + 1)/Model(K)$ indicates the two models in comparison: (a) $D(0) : Model(1)/Model(0)$, (b) $D(1) : Model(2)/Model(1)$, (c) $D(2) : Model(3)/Model(2)$, (d) $D(3) : Model(4)/Model(3)$, (e) $D(4) : Model(5)/Model(4)$, (f) $D(5) : Model(6)/Model(5)$, (g) $D(6) : Model(7)/Model(6)$, (h) $D(7) : Model(8)/Model(7)$, (i) $D(8) : Model(9)/Model(8)$, (j) $D(9) : Model(10)/Model(9)$, and (k) $D(10) : Model(11)/Model(10)$.</figcaption></figure>
4.2.2. Bootstrapping the proposed modified statistic

Another possibility worth investigating is to bootstrap the modified LLR, defined in Eq. (29), which proved to be more effective in likelihood ratio testing of nested factor models than the conventional statistic defined in Eq. (19). Therefore a bootstrap algorithm similar to Algorithm 3, with the modified LLR replacing the conventional one, is presented. The quantity defined in Eq. (30) has a reference distribution, $D(K)$, that one can obtain using parametric bootstrap under the restricted model (i.e., the one with $K$ sources), in the same way described in the previous subsection. 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. Algorithm 4 states the steps for the selection of the number of sources.
**Algorithm 4** Bootstrapping the modified LLR

* Inputs: $\mathbf{Y}, \mathbf{R}, \mathbf{S}_{YY}, \mathbf{S}_{RR}$ and $\mathbf{S}_{YR}$

* Estimate $\mathbf{C}(\mathbf{R})$ and $\mathbf{P}(\mathbf{R})$ from $\mathbf{S}_{YY}, \mathbf{S}_{RR}$ and $\mathbf{S}_{YR}$ assuming $\mathbf{R}$ sources

* Estimate $\hat{\mathbf{X}}(\mathbf{R})$ and $\hat{\mathbf{\Omega}}(\mathbf{R})$ from Eqs. (27) and (28)

* Set $K = 0$

* WHILE $K < R$

- Construct $\mathbf{C}(\mathbf{K})$ and $\mathbf{P}(\mathbf{K})$ from $\mathbf{C}(\mathbf{R})$ and $\mathbf{P}(\mathbf{R})$ assuming $\mathbf{K}$ sources

- Estimate $\hat{\mathbf{X}}(\mathbf{K})$ and $\hat{\mathbf{\Omega}}(\mathbf{K})$ from Eqs. (12) and (13)

- Compute $\Lambda(\mathbf{K})$ from Eq. (30)

- Set $i = 1$

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

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

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

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

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

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

* Estimate $\hat{\mathbf{X}}(\mathbf{R})_{\text{new}}$ and $\hat{\mathbf{\Omega}}(\mathbf{R})_{\text{new}}$ from Eqs. (27) and (28) applied to $\mathbf{C}(\mathbf{R})_{\text{new}}, \mathbf{P}(\mathbf{R})_{\text{new}}$ and $\mathbf{R}_{\text{new}}$

* Estimate $\hat{\mathbf{X}}(\mathbf{K})_{\text{new}}$ and $\hat{\mathbf{\Omega}}(\mathbf{K})_{\text{new}}$ from Eqs. (12) and (13) applied to $\mathbf{C}(\mathbf{K})_{\text{new}}, \mathbf{P}(\mathbf{K})_{\text{new}}$ and $\mathbf{R}_{\text{new}}$

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

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

* $i \leftarrow i + 1$

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

- ELSE $K \leftarrow K + 1$

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

* ELSE $S \leftarrow K$

The same simulations as in Section 4.2.1 was performed to show the boot-
strapping results when the modified statistic is used, and the results are shown in
Fig. 5. The results clearly show that the saturated model (with $\mathbf{R}$ sources) does not
provide any further information about the data than the model with $K = \mathbf{S} = 4$,
which is the true number of sources, see Fig. 5(e). Hence, bootstrapping the mod-
ified LRT excellently works in estimating the number of sources that makes up
the data.
<figure><img src="image_11.png" /><figcaption>**Fig. 5.** Simulation results of bootstrapping the modified LRT; comparison of $\Lambda(K)$ (red dotted line) with the bootstrapped distribution $D(K)$ (blue bar-chart) for the 11 virtual references (as returned by the MCR technique). The designation *Model(R)/Model(K)* indicates the two models in comparison: (a) *D(0) : Model(11)/Model(0)*, (b) *D(1) : Model(11)/Model(1)*, (c) *D(2) : Model(11)/Model(2)*, (d) *D(3) : Model(11)/Model(3)*, (e) *D(4) : Model(11)/Model(4)*, (f) *D(5) : Model(11)/Model(5)*, (g) *D(6) : Model(11)/Model(6)*, (h) *D(7) : Model(11)/Model(7)*, (i) *D(8) : Model(11)/Model(8)*, (j) *D(9) : Model(11)/Model(9)*, and (k) *D(10) : Model(11)/Model(10)*, note: here *Model(R)* is the saturated model with *R* = 11 (the full set of references) for all test steps.</figcaption></figure>
### 4.3. The cross-validation-based approach

The last presented alternative solution for the source enumeration problem is cross-validation, which requires the existence of two data sets (one for the estimation step and the other for the evaluation step) generated by the same stochastic processes, in the current case ($\mathbf{Y}_1, \mathbf{R}_1$) and ($\mathbf{Y}_2, \mathbf{R}_2$). Hence, starting from one of the two data sets, say ($\mathbf{Y}_1, \mathbf{R}_1$), one estimates $\mathbf{G}_1(R) = \mathbf{C}_1(R)\mathbf{P}_1(R)^H$ for the saturated model, then constructs $\mathbf{G}_1(K) = \sum_{s=1}^K \mathbf{c}_{1s}\mathbf{p}_1^H$ using the first $K$ columns in $\mathbf{C}_1(R)$ and $\mathbf{P}_1(R)$ for each sub-model of $K$ sources, with $K = 0, 1, ..., R-1$. One also estimates $\hat{\mathbf{\Omega}}_1(R) = (\mathbf{Y}_1 - \mathbf{G}_1(R)\mathbf{R}_1)(\mathbf{Y}_1 - \mathbf{G}_1(R)\mathbf{R}_1)^H/I^{est}$ for the saturated model, and $\hat{\mathbf{\Omega}}_1(K) = (\mathbf{Y}_1 - \mathbf{G}_1(K)\mathbf{R}_1)(\mathbf{Y}_1 - \mathbf{G}_1(K)\mathbf{R}_1)^H/I^{est}$ for each sub-model, with $I^{est}$ the number of snapshots (windows) in the estimation subset. A metric is needed for the evaluation step of each sub-model, the one used here is defined by the following LLR:

$$ \Lambda(K) = \ln \frac{P(\mathbf{Y}_2|\mathbf{R}_2, \hat{\theta}_{1_R})}{P(\mathbf{Y}_2|\mathbf{R}_2, \hat{\theta}_{1_K})} \quad (31) $$

where $\hat{\theta}_{1i}$ is the maximum likelihood estimates of the free parameters of the model with $i$ number of sources, as estimated from ($\mathbf{Y}_1, \mathbf{R}_1$). Taking the same assumptions about the data as before (in Section 4.1.1), and substituting with the maximum likelihood estimates one gets

$$ \Lambda(K) = -(I_{\text{eff}}^{\text{evl}}/I^{\text{evl}})\|\mathbf{Y}_2 - \mathbf{G}_1(R)\mathbf{R}_2\|_{\hat{\mathbf{\Omega}}_1(R)}^2 + (I_{\text{eff}}^{\text{evl}}/I^{\text{evl}})\|\mathbf{Y}_2 - \mathbf{G}_1(K)\mathbf{R}_2\|_{\hat{\mathbf{\Omega}}_1(K)}^2 \\ - I_{\text{eff}}^{\text{evl}} \ln |\hat{\mathbf{\Omega}}_1(R)| + I_{\text{eff}}^{\text{evl}} \ln |\hat{\mathbf{\Omega}}_1(K)| \quad (32) $$

where $I_{\text{eff}}^{\text{evl}}$ and $I_{\text{eff}}^{\text{evl}}$ are the number of snapshots (windows) and the effective number of snapshots in the evaluation set, respectively. The choice of evaluating the sub-models against the saturated model is again used here, because of the effectiveness, which it showed earlier, in evaluating sub-models. Hence, $\Lambda(K)$ is calculated, from Eq. (32), for each sub-model, i.e., for $K = 0, 1, ..., R-1$, and for the saturated model itself $\Lambda(K = R) = 0$, by definition. Finally, the model with the least $\Lambda(K)$, with number of sources $\hat{S} = K$, is selected. The model with the correct number of sources will return the least $\Lambda(K)$ value because the denominator of Eq. (31) ($P(\mathbf{Y}_2|\mathbf{R}_2, \hat{\theta}_{1_K})$) is maximized for the correct model, i.e., when $K = S$ (true number of sources). As the likelihood of the data is maximized, given $\hat{\theta}_{1_K} = \hat{\theta}_{1_S}$ of the correct model. While under-fitting (when $K < S$) produces relatively large error, over-fitting (when $K > S$) produces relatively smaller error, given that the data is of a reasonable SNR, i.e., the noise is not dominant. That is why it is expected a fast rate of decrease in error for ($K < S$) and a slow rate
increase of error for ($K > S$). The details for implementing cross-validation are described in Algorithm 5. There are several possible cross-validation strategies for a given data set. One possibility, which is used in this work, is to apply a 2-fold cross validation strategy by splitting the data into two equal subsets. The number of snapshots of each subset is $I^{sub} = I^{est} = I^{evl} = I/2$, with $I$ total number of snapshots. Then, subset #1 is used for estimation, and subset #2 is used for sub-model evaluation, resulting in $\Lambda(K)_1$. After that, the process is repeated with the two subsets changing roles, resulting in $\Lambda(K)_2$. Finally, the metric used for evaluating a given sub-model is the mean value of the two LLRs. In this manner, a good balance between estimation accuracy and processing time is achieved.

**Algorithm 5 Cross-validation**

* Inputs: $\mathbf{Y}_1, \mathbf{R}_1, \mathbf{Y}_2, \mathbf{R}_2$ and $\mathbf{S}_{Y_1 Y_1}, \mathbf{S}_{R_1 R_1}, \mathbf{S}_{Y_1 R_1}$, using $(\mathbf{Y}_1, \mathbf{R}_1)$

* Estimate $\mathbf{C}_1(R)$ and $\mathbf{P}_1(R)$ from $\mathbf{S}_{Y_1 Y_1}, \mathbf{S}_{R_1 R_1}, \mathbf{S}_{Y_1 R_1}$ assuming $R$ sources

* Estimate $\mathbf{G}_1(R) = \mathbf{C}_1(R)\mathbf{P}_1(R)^H$ and $\hat{\mathbf{\Omega}}_1(R) = (\mathbf{Y}_1 - \mathbf{G}_1(R)\mathbf{R}_1)(\mathbf{Y}_1 - \mathbf{G}_1(R)\mathbf{R}_1)^H/I^{est}$

* FOR $K = 0$ to $R$

- IF $K = R$

$$ \ast \Lambda(K) = 0 $$

- ELSE

* Construct $\mathbf{G}_1(K) = \sum_{s=1}^K \mathbf{c}_{1s}\mathbf{p}_1^H$ from the first $K$ columns of $\mathbf{C}_1(R)$ and $\mathbf{P}_1(R)$

* Calculate $\hat{\mathbf{\Omega}}_1(K) = (\mathbf{Y}_1 - \mathbf{G}_1(K)\mathbf{R}_1)(\mathbf{Y}_1 - \mathbf{G}_1(K)\mathbf{R}_1)^H/I^{est}$

* Calculate $\Lambda(K)$ from Eq. (32)

* $\hat{S} = \text{argmin}_K \Lambda(K)$

Simulations similar to the one in section 4.2.1 were performed to validate the proposed algorithm, where $S = 4$ sources, $M = 50$ outputs, and $R = 11$ references were simulated with 90% correlated noise. However, this time, to validate the cross-validation technique, two data sets were randomly generated using two independently simulated source matrices, $\mathbf{S}_1 \in \mathbb{C}^{S \times I/2}$ and $\mathbf{S}_2 \in \mathbb{C}^{S \times I/2}$, where $I$ is the total number of snapshots, two independently simulated output noises, $\mathbf{N}_1$ and $\mathbf{N}_2$, and two independently simulated reference noises, $\mathbf{E}_1$ and $\mathbf{E}_2$, but using the same transfer functions, $\mathbf{H}$ and $\mathbf{L}$ for the two instances. The simulation was repeated four times with different numbers of snapshots $I = 150, 300, 600$ and $1200$. The entire simulation was run $100$ times and the average values of the LLR
$(\Lambda(K))$ are shown in Fig. 6 for $K = 0, 1, ..., 11$, for different number of snapshots. It is clear from Fig. 6(a), that a number of snapshots $I = 150$ is not enough for the cross-validation technique to work, as it returned a minimum LLR value of $\Lambda(0)$ that does not correspond to the true number of sources, i.e., at $K = 0 \neq S$. Whereas, only 74 snapshots were enough for bootstrapping to work (in Sections 4.2.1 and 4.2.2). However, it can be seen from Fig. 6(b), that 300 snapshots are enough for cross-validation to correctly estimate the true number of sources, as it gives a minimum LLR of $\Lambda(4)$, which corresponds to the true number of sources, $K = S = 4$. Further increase of the number of snapshots, in Figs. 6(c) and 6(d), proves the stability of the solution, and gives the characteristic behavior expected from a cross-validation error value.

<figure><img src="image_12.png" /><figcaption>**Fig. 6.** The simulation results of the cross-validation technique using (a) $I = 150$, (b) $I = 300$, (c) $I = 600$ and (d) $I = 1200$ number of snapshots.</figcaption></figure>

## 4.4. A comparative simulation

Comparing the results of the previous simulations reveals important observations regarding the sensitivity of the three compared source enumeration method to the number of snapshots. Bootstrapping perfectly works with a few snapshots, roughly as few as 150% of the number of output sensors (e.g., for $M = 50, I = 74$). Cross-validation perfectly works for a moderate number of snapshots, roughly 6
times the number of output sensors (e.g., for *M* = 50, *I* = 300). The direct imple-
mentation of the asymptotic LRT requires a huge number of snapshots to work,
roughly 12 fold of the number of output sensors (e.g., for *M* = 50, *I* = 600), still
with potential slight over-estimation. Although these observations are clear, they
are seen from independent simulations. Therefore, it is more accurate to perform
a comparison based on the same simulation instance.

A data set with statistical characteristics and SNRs as in Section 4.2.1 was sim-
ulated to compare the proposed three algorithms (the modified LRT, bootstrapping
and 2-fold cross-validation). The data were simulated with *S* = 4 sources, *M* = 50
outputs, and *R* = 20 references with 90% correlated noise. A higher number of
references used here, compared to previous simulations, in order to test if it has
effect on the estimation of the number of sources. The data were simulated in the
form of two subsets of the same size (i.e., of same number of snapshots *I*<sup>sub</sup> = *I*/2),
in the same manner as in the previous subsection, for the purpose of applying the
2-fold cross-validation. While the full data of *I* number of snapshots were fed
to the other two techniques. The simulation was repeated with different numbers
of snapshots that quasi-quadratically span a range of *I* = 50 to 10000. The data
were processed using the three compared methods, namely, the modified LRT,
bootstrapping the modified LLR (both with 1 − *α* = 0.99), and cross-validation,
all based on MCR as the decomposition tool. The entire simulation was repeated
20 times to test repeatability.

The results of the comparison are depicted in Fig. 7(a), in terms of range
(shaded area) and mean value (line), which shows the convergence of the results
of the three compared methods to the true number of sources (Ś = *S* = 4), at a
distinct number of snapshot for each. It also shows a great stability, accuracy and
precision of the three techniques, for numbers of snapshots beyond convergence,
except from a potential for a slight overestimation in the case of LRT. In order to
compare the convergence points of the three methods, a zoomed view is depicted
in Fig. 7(b). Note that, there are no cross-validation results below *I* = 100, as
this point corresponds to cross-validation subsets of *I*<sup>sub</sup> = *I*/2 = 50 snapshots,
which is the minimum required for a non-singular output cross-spectral matrix (as
*M* = 50). It is seen that, in this simulation instance, bootstrapping needs only
*I* = 100 snapshots (i.e., *I* = 2 × *M*) to work. Cross-validation needs *I* = 300 (i.e.,
*I* = 6 × *M*), to work. While the asymptotic modified LRT needs *I* = 620 snapshots
(i.e., as high as *I* ≈ 12 × *M*) to work. These results confirms the observations seen
from the previous stand-alone simulations. The methods also showed robustness
against the change in the number of references, as the results are nearly similar
(up to random variations) between this simulation of *R/S* = 5 and the previous
simulations of R/S = 11/4 = 2.75.

<figure><img src="image_13.png" /><figcaption>**Fig. 7.** Comparison results of the three source enumeration methods (a), and a zoomed version (b): the estimated number of sources $\hat{S}$ is plotted, in terms of mean (line) and range (shaded area), against the number of snapshots $I$, the true number of sources is $S = 4$, the number of output sensors is $M = 50$ (represented by the dotted light gray vertical line), the results of the modified LRT (in blue) converge at $I = 620$, of bootstrap (in dashed red) converge at $I = 100$, and of cross-validation (in dash-dot green) converge at $I = 300$.</figcaption></figure>

### 4.5. Remarks on the simulations

Finally, it can be concluded from all the previous simulations that bootstrap-ping needs a ratio of the number of independent snapshots to the number of output sensors $I/M > 2$ to work. Cross-validation needs $I/M > 6$ to work. The modified LRT requires $I/M > 12$ to work, still with potential slight over-estimation. Noting
that, *I*, in this ratio, is replaced by *I*<sub>eff</sub> in the case of correlated snapshots. It was also shown that the conventional sequential LRT (with the conventional LLR between two consecutive models) is not valid for the current factor analysis model, even when using extremely high number of snapshots (e.g., *I* = 3000). It was also shown that the proposed modified LRT (with the modified LLR of a specific model against the saturated model) is valid for the current factor analysis model and works perfectly fine for a relatively fewer snapshots (compared to the conventional LRT). Hence, from now on, only the modified LLR defined in Eq. (29) is considered in the upcoming sections, where laboratory experimental and real industrial data are used for the comparison among the LRT, parametric bootstrap and cross-validation.

**5. Laboratory Experiments**

This section describes a laboratory experiment conducted in the semi-anechoic
room of the Laboratory of Vibration and Acoustics (LVA), INSA Lyon. The data
are used for the comparison of the three methods, namely, the direct LRT based
on the modified LLR (Algorithm 2), the bootstrapping of the modified LLR (Al-
gorithm 4), and cross-validation (Algorithm 5), and using MCR as a datum for
comparison. The experiment contained two sound sources, i.e., *S* = 2. The exper-
imental apparatus is depicted in Fig. 8(a) and 8(b). An array of 20 microphones
was placed in the semi-anechoic room, in order to capture the sound radiated by
two speakers. Four reference microphones were installed near to the speakers.
Finally, an acquisition system was used to collect the data with a sampling fre-
quency of 51.2 kHz, and to feed the speakers with white Gaussian noise with an
upper cut-off frequency of 20 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.
The spectral analysis was done using Welch method, with frequency resolution
of 0.8 Hz, which corresponds to a hanning window of length *L* = 2<sup>16</sup> samples
(1.28 seconds), with overlap of 2/3. To test the compared methods under different
numbers of snapshots, the record length was varied from the sufficiently long to
the very short. The record lengths used in the comparison, and the corresponding
number of snapshots are stated in Table 1. As a validation set is needed for the
cross-validation technique, the data ware divided into two equivalent subsets (i.e.,
of the same numbers of snapshots), and a 2-fold cross-validation strategy was ap-
plied. Note that, for a fair comparison, the cross-validation results at *I*<sup>sub</sup> = *n*
subset's snapshots (which counts to a total number of snapshots *I* = 2*n*) compares
to the LRT (and bootstrapping) results at *I* = 2*n*. Note that, for the cases (i) and (j)
in Table 1, the cross-validation-related columns are empty, because the values of $I_{\text{eff}}^{\text{sub}}$ in these cases are less than the number of sensors $M = 20$, causing numerical problems to the spectral estimation, hence, cross-validation is not applied to these cases. The results of the comparison, for the two experiments, are shown in the following subsections.

<figure><img src="image_14.png" /><figcaption>**Fig. 8.** Apparatus of the laboratory experiment: (a) Schematic of experiment, (b) Apparatus of experiment.</figcaption></figure>

<table><caption>**Table 1.** Different signal lengths used in the comparison.</caption><thead><tr><th></th><th>Total signal length (in samples)</th><th>Total signal length (in seconds)</th><th>Total number of snapshots (I)</th><th>Total effective number of snapshots (I<sub>eff</sub>)</th><th>Number of snapshots in the cross-validation subset (I<sup>sub</sup>)</th><th>Effective number of snapshots in the cross-validation subset (I<sub>eff</sub><sup>sub</sup>)</th></tr></thead><tbody><tr><td>(a)</td><td>41779200</td><td>816</td><td>1910</td><td>1321</td><td>954</td><td>660</td></tr><tr><td>(b)</td><td>20889600</td><td>408</td><td>954</td><td>660</td><td>476</td><td>329</td></tr><tr><td>(c)</td><td>10444800</td><td>204</td><td>476</td><td>329</td><td>237</td><td>164</td></tr><tr><td>(d)</td><td>5222400</td><td>102</td><td>237</td><td>164</td><td>117</td><td>81</td></tr><tr><td>(e)</td><td>2611200</td><td>51</td><td>117</td><td>81</td><td>57</td><td>39</td></tr><tr><td>(f)</td><td>1899054</td><td>37.1</td><td>84</td><td>58</td><td>41</td><td>28</td></tr><tr><td>(g)</td><td>1740800</td><td>34</td><td>77</td><td>53</td><td>37</td><td>25</td></tr><tr><td>(h)</td><td>1606892</td><td>31.4</td><td>71</td><td>49</td><td>34</td><td>23</td></tr><tr><td>(i)</td><td>1160533</td><td>22.66</td><td>51</td><td>35</td><td>—</td><td>—</td></tr><tr><td>(j)</td><td>949527</td><td>18.54</td><td>41</td><td>28</td><td>—</td><td>—</td></tr></tbody></table>

## 5.1. The effect of the number of snapshots on the MCR-based principal spectra

Applying the MCR technique to the data, one ends up with the MCR-principal components, which are the generalized eigenvalues of $(S_{RY}\Gamma S_{YR}, S_{RR})$. Figure 9 shows these principal spectra. The results are shown for an effective number of
$I_{\text{eff}}$ = 660, and for $I_{\text{eff}}$ = 28. First, and before analyzing the results, one should keep in mind that the excitation was only up to 20 kHz. In general, it is expected that a number of principal spectra equal to the number of sources ($S$ = 2) is protrude from the remaining principal spectra. However, the results are seen to be highly affected by the used number of snapshots. As seen in Fig. 9(a), for $I_{\text{eff}}$ = 660, the $S$ principal components are more than 10 dB higher than the noise floor, whereas for $I_{\text{eff}}$ = 28, in Fig. 9(b), the difference is less than 1dB. This is because of the fact that using fewer snapshots, in Welch's method, results in a higher variance in the PSD estimate (i.e., higher random error). It is noticed also in the case of relatively high number of snapshots ($I_{\text{eff}}$ = 660), the reduction of the random error revealed some band-limited components (in the range from 3 kHz to 13 kHz), which were buried under the random error in the case of $I_{\text{eff}}$ = 28. These band-limited high amplitudes could be originating from weak source/sources extraneous to the experiment (e.g., from the measurement system). The dependence of the behavior of the MCR-principal spectra on the used number of snapshots certainly, in turn, affects the estimation of the number of sources, which is shown in the following subsection for the three compared methods.

<figure><img src="image_15.png" /><figcaption>Fig. 9. The MCR-based principle components, for the LVA experiment at (a) $I_{\text{eff}}$ = 660 and (b) $I_{\text{eff}}$ = 28.</figcaption></figure>

## 5.2. *The estimated number of sources*

Figures 10, 11 and 12 show the results of the experiment, for the direct LRT, parametric bootstrap and cross-validation, respectively. Each figure shows the results for the different values of $I_{\text{eff}}$ stated in Table 1. Observations similar to the ones drawn from the simulations (stated in section 4.5) are witnessed in the results of the laboratory experiment. It is obvious that the performance of the estimation process is affected by the used number of snapshots, and it is also obvious that
this effect varies from an estimation method to the other. As seen from the results of the direct LRT in Fig. 10, there exists a frequent overestimation of the number of sources (i.e., $\hat{S}$) for the whole frequency range due to the random error accompanied with the PSD estimation, which could not be observed in the performed (frequency-based) simulations, where data were directly simulated in the frequency domain. The results get better with increasing the number of snapshots. The best results, at Fig. 10(a), are hence achieved when using the longest available record, still with a frequent overestimation by 1 source (i.e., $\hat{S} = 3$ instead of 2). A local rise of 1 to 2 sources in $\hat{S}$ is also observed in the range from 3 kHz to 13 kHz, which correspond to the rise of the MCR-principal spectra witnessed in Fig. 9(a), originating from sources extraneous to the experiment (potentially from the measurement chain). By decreasing the number of snapshots the increased random error ruin the results for the whole frequency range causing an increased overestimation of $S$. The susceptibility of the LRT to the random error is understood, as the whole method is asymptotic with the number of snapshots $I_{\text{eff}}$, hence, a small $I_{\text{eff}}$ means a departure of the data distribution from the $\frac{1}{2}\chi_{2\nu}^2$ distribution, resulting in overestimation.

The results of bootstrapping, in Fig. 11 show great immunity against the random error, as the method returns $\hat{S} = 2$ for few snapshots (i.e., for $I_{\text{eff}} = 35$ to 58). It even returns very good results (in Fig. 11(j)) at $I_{\text{eff}} = 28$, which is as few as $1.4 \times M$ (where $M = 20$ array microphone), which represents excellent results at such low number of snapshots. Considering the witnessed immunity of bootstrapping against random error, and by looking at the bootstrapping results at relatively high numbers of snapshots (i.e., $I_{\text{eff}} = 1321, 660$ and 329), in Figs. 11(a), 11(b) and 11(c), it is seen both the global rise of $\hat{S}$ by 1 source, and the band-limited rise (in the range from 3 kHz to 13 kHz), which also exist in the results of the LRT (in Figs. 10(a), 10(b) and 10(c)). This suggests that the rise of 1 source over the whole frequency range is not due to the random error, but due to a weak broadband source extraneous to the experiment, that reveals when the random error lessens due to the use of a high number of snapshots.

The results of cross-validation, in Fig. 12, show underestimation of the number of sources at low values of $I_{\text{eff}}$, as seen in Figs. 12(f), 12(g) and 12(h). Regarding the sufficient $I_{\text{eff}}$ needed for the method to work, it gives good results starting from a moderate total number of snapshots, $I_{\text{eff}} = 81$, i.e., $I_{\text{eff}}/M = 4$ (in Fig. 12(e)), compared to a lower ratio of $I_{\text{eff}}/M = 1.4$, at which bootstrapping gives good results (in Fig. 11(j)), which conforms with what was found by simulations. In general, it returns results of the same level of goodness (on average over all frequencies) as the other two techniques at high values of $I_{\text{eff}}$, but it suffers from
underestimation at low $I_{\text{eff}}$.

Overall, the parametric bootstrap seems to be the most robust technique among the three, as it returns the best on-average results over all values of $I_{\text{eff}}$.

It should be noted also that the upper cut-off frequency of the excitation is 20 kHz, which expects $\hat{S} = 2$ below and $\hat{S} = 0$ above this frequency. This is a good test to check the capability of the methods to returns an estimate of $S$ that depends on frequency. By visual inspection, it is clear that all the three methods return $\hat{S}$ that is different on average above and below 20 kHz.

### 5.3. Remarks on the laboratory experiments

Aside from the observations that were initially witnessed from simulation, and then confirmed by the two experiments, there is a practical advice that could be concluded from the results. The best way to find out the true number of sources of a real data set is to process the data with different source enumeration method, and by using various signal lengths (i.e., various $I_{\text{eff}}$). In this way, and by comparison, one could arrive at the true number of sources. However, this needs long records, which are not always available when dealing with real industrial data. Hence, when only short records are available, bootstrapping is the optimum choice. This will be also confirmed by the analysis of the real industrial case in the next section, where only short records are available.

<figure><img src="image_16.png" /><figcaption>**Fig. 10.** Experiment 2: the estimated number of sources using the direct LRT using (a) $I_{\text{eff}} = 1321$, (b) $I_{\text{eff}} = 660$, (c) $I_{\text{eff}} = 329$, (d) $I_{\text{eff}} = 164$, (e) $I_{\text{eff}} = 81$, (f) $I_{\text{eff}} = 58$, (g) $I_{\text{eff}} = 53$, (h) $I_{\text{eff}} = 49$, (i) $I_{\text{eff}} = 35$ and (j) $I_{\text{eff}} = 28$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 11.** Experiment 2: the estimated number of sources using parametric bootstrap using (a) *I*<sub>eff</sub> = 1321, (b) *I*<sub>eff</sub> = 660, (c) *I*<sub>eff</sub> = 329, (d) *I*<sub>eff</sub> = 164, (e) *I*<sub>eff</sub> = 81, (f) *I*<sub>eff</sub> = 58, (g) *I*<sub>eff</sub> = 53, (h) *I*<sub>eff</sub> = 49, (i) *I*<sub>eff</sub> = 35 and (j) *I*<sub>eff</sub> = 28.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 12.** Experiment 2: the estimated number of sources using cross-validation using (a) *I*<sub>eff</sub> = 1321, (b) *I*<sub>eff</sub> = 660, (c) *I*<sub>eff</sub> = 329, (d) *I*<sub>eff</sub> = 164, (e) *I*<sub>eff</sub> = 81, (f) *I*<sub>eff</sub> = 58, (g) *I*<sub>eff</sub> = 53 and (h) *I*<sub>eff</sub> = 49.</figcaption></figure>
# 6. Application to industrial data

Data measured on a SimRod e-motor test bench from SIEMENS were used for comparing the performance of the three methods proposed for the estimation of the number of sources.

## 6.1. *Description of SIEMENS' test bench and the data*

The sound emitted from the electric motor was measured using an array consisting of 50 microphones, located 2 cm away from the motor, which was running at 2500 rpm. The considered position for analysis is illustrated by the blue rectangle in Fig. 13(a) and by a top view of the array in-place in Fig. 13(b). Five fixed microphones were used as references. The sampling frequency of the microphone signals was 40960 Hz. Data from two tri-axial accelerometers (6 tracks) were also used as references; hence, there exists a total of 11 references. The sampling frequency of the acceleration tracks was 10240 Hz, which is only one-fourth the sampling frequency of the microphone signals, hence, they were upsampled by a factor of 4. White Gaussian noise was added to every acceleration signal with SNR of 50 dB, to keep a full-rank CSM of the references. The time records had a length of approximately 14.6 seconds, which corresponds to 600,000 samples. In cross-validation, the total record was divided into two subsets of the same length (i.e., 7.3 seconds). For a better conditioning of the CSM of the references, the Fourier spectra of all signals (pressure and acceleration) were normalized by dividing them by their corresponding (root mean square) RMS values, allowing the incorporation of acceleration signals alongside microphone (pressure) signals. A frequency resolution of 5 Hz was obtained using a hanning window of length $L = 2^{13}$ samples (0.2 second), and 2/3 overlap, when computing the power spectra with the Welch's method, which corresponds to an effective number of snapshots of $I_{\text{eff}} = 150$.

The PSDs of the array and reference signals are shown in Figs. 14(a) and 14(b), respectively. Where it is clearly seen the loss of the acceleration signals above 5 kHz, which corresponds to the initial (i.e., before upsampling) Nyquist frequency of the acceleration signals. This loss will have an impact on the estimation of the number of sources.
<figure><img src="image_19.png" /><figcaption>Fig. 13. SIEMENS e-motor test bench: (a) the analyzed array position indicated by the blue rectangle, and (b) the array situated at the position for measurement.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 14. PSDs: (a) of the array signals, and (b) of the reference signals.</figcaption></figure>
## 6.2. The results

The three methods proposed for estimating the number of sources was applied to the SIEMENS e-motor data. The MCR-based principal spectra, and the estimated number of sources $\hat{S}$, as returned by each of the three methods are displayed in Fig. 15. In this industrial example, and contrary to all the previous numerical and laboratory experimental examples, the true number of sources is unknown. However, a known fact is that above 5 kHz (i.e., the initial Nyquist frequency of the acceleration tracks) the acceleration signals contain pure noise uncorrelated with the true sources. Hence, $\hat{S}$ can not exceed the number of physical references, which is 11 in the frequency range below 5 kHz, and can not exceed 5 in the frequency range above 5 kHz. This information can be used as a criterion for judging the performance of the methods under comparison. From the first sight on Fig. 15(a), one notices, above 5 kHz, the separation between the five source-related principal components and the six noise-related principal components. The direct application of the LRT, in Fig. 15(b), overestimates the number of sources, as it gives $\hat{S} > 5$ for frequencies above 5 kHz. This overestimation is due to the small number of snapshots used in the LRT relatively to the number of microphones. While cross-validation, in Fig. 15(d), suffers from an underestimation of $S$, as it gives $\hat{S} = 0$ at many frequencies. This is also due to the small number of snapshots, as already observed in the previous experiment (in Figs. 12(f), 12(g) and 12(h)). It could be argued that, the total $I_{\text{eff}} = 150$ in the current case (of the e-motor data) is more than the total $I_{\text{eff}} = 58$ in the experimental case (in Fig. 12(f)), however one should keep in mind the fact that in the current case the ratio $I_{\text{eff}}/M = 150/50 = 3$ (where $M = 50$ array microphone) nearly equals the ratio $I_{\text{eff}}/M = 58/20 = 2.9$, in the experimental cases (where $M = 20$ microphone only), which is lower than the minimal value of 6, required for cross-validation to work, found by simulations. Finally, bootstrapping returns, as expected, the only reliable results, in Fig. 15(c). It does not overestimate $S$, as it gives $\hat{S} < 11$ below 5 kHz, and $\hat{S} \le 5$ above 5 kHz. It does not underestimate $S$, as it gives $\hat{S} > 0$ for all frequencies. It seems to converge around a specific number of sources for all frequencies, $\hat{S} = 9$ below 5 kHz, and $\hat{S} = 5$ above 5 kHz. This means that the true number of sources is most likely to be $S = 9$ for all frequencies, and $\hat{S}$ dropped to 5 after 5 kHz just due to the loss of the acceleration signals.

To visualize the effect of losing the acceleration signals on the source extraction results, acoustic imaging results just before and after 5 kHz are presented in Fig. 16. The figure depicts two radiation maps ($5 \times 10$ interpolated plots) of the magnitudes of Fourier coefficients of the extracted source contributions at the 50 array microphones, for a single snapshot, at 4 kHz and 6 kHz. The source ex-
<figure><img src="image_21.png" /><figcaption>**Fig. 15.** Estimated number of sources $\hat{S}$ of SIEMENS data: (a) MCR-based principal spectra, (b) the direct LRT results, (c) bootstrapping results, and (d) cross-validation results.</figcaption></figure>

traction results at 4 kHz (Fig. 16(a)) were obtained using an estimated 8 sources, whereas the results at 6 kHz (Fig. 16(b)) were obtained using an estimated 4 sources, as returned by bootstrapping. As shown in Fig. 16(b), the drop in source energy is evident due to the use of an underestimated number of sources, resulting from the loss of the acceleration signals. This highlights the importance of correctly estimating the number of sources and, prior to that, ensuring the use of a number of references at least equal to the number of sources.

Another important point is that the sources are not necessarily spatially disjoint (i.e., they may overlap in space), particularly in the case of a compact machine. In fact, the assumption of spatially disjoint sources is rarely satisfied in practice, especially on the measurement plane (the array). As a result, a single high-energy spot on the acoustic image may correspond to multiple sources. Conversely, multiple high-energy spots could be reflections of the same source, leading to the appearance of two or more spots for a single source. Therefore, relying solely on an acoustic image to count the number of sources is unreliable. This highlights the value of the statistical source enumeration methods presented in this paper.
<figure><img src="image_22.png" /><figcaption>Fig. 16. Acoustic imaging results of the analyzed array position: (a) at 4 kHz, and (b) at 6 kHz.</figcaption></figure>

7. Conclusions

Source enumeration is crucial for reference-based source extraction techniques while the number of references exceeds the number of sources. The estimation of number of sources is essential for reducing the rank of the signal's subspace in order to isolate the sources contribution from the noise while avoiding any numerical problems. Source enumeration, in this context, is a hard ill-posed inverse problem, which has not tackled enough in the literature within the framework of the reference techniques. What motivated the presentation of three alternative techniques for source enumeration. One of them is an asymptotic solution based on an LRT against the saturated model, and the other two are exact-pdf computer-intensive solutions built upon the same LLR from the proposed saturated LRT. The first exact-pdf solution is parametric bootstrap, which is a parametric solution, while the second is non-parametric which is cross-validation. For permitting these techniques to be applicable to all reference based techniques, a unified formalism for all the reference techniques was provided. After performing a comparative study based on numerical simulations, laboratory experiments and real industrial data, the number of snapshots (hence, the signal length) turned out to have a great effect on the performance of the three methods. A minimal value of the ratio of the number of independent (otherwise, effective) snapshots to the number of output sensors is necessary for a method to work, which differs from a method to the other. For bootstrapping it is 2, for cross-validation it is 6, and for the modified
LRT it is 12. Hence, parametric bootstrap turned out to be the best, especially when only short records are available.

Acknowledgements

We gratefully acknowledge the European commission for its support of the Marie Sklodowska 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).

Appendix A. The equivalence between a set of correlated random variables and a reduced set of uncorrelated random variables

This appendix works out the equivalence between a set of correlated random variables and a reduced set of uncorrelated random variables. Let {$Y_i$; $i = 1, ..., I$} be a set of $I$ zero-mean, correlated, Gaussian random variables. The aim is to replace it by a set of $I_{eff} < I$ zero-mean, uncorrelated, Gaussian random variables, {$Z_i$; $i = 1, ..., I_{eff}$}, with same mean statistical properties on the second-order. The $Z_i$'s are such that $V\{Y_i\} = V\{Z_i\} = \sigma^2$. Therefore, one imposes

$$ \frac{1}{I} \sum_{i=1}^{I} |Y_i|^2 = \frac{1}{I_{eff}} \sum_{i=1}^{I_{eff}} |Z_i|^2. \qquad (A.1) $$

Note that, taking the expected value of the above equality, one correctly recovers $\sum_{i=1}^{I} \sigma^2/I = \sum_{i=1}^{I_{eff}} \sigma^2/I_{eff}$. The sufficient statistic $\sum_{i=1}^{I_{eff}} |Z_i|^2$ can then be replaced by $(I_{eff}/I) \sum_{i=1}^{I} |Y_i|^2$.

For instance, Eq. (18) applied to correlated random variables becomes

$$ \Lambda(K) = -(I_{eff}/I)\|\mathbf{Y} - \hat{\mathbf{X}}(K+1)\|^2_{\hat{\Omega}(K+1)} + (I_{eff}/I)\|\mathbf{Y} - \hat{\mathbf{X}}(K)\|^2_{\hat{\Omega}(K)} \\ - I_{eff} \ln |\hat{\Omega}(K+1)| + I_{eff} \ln |\hat{\Omega}(K)|. \qquad (A.2) $$
# Appendix B. Derivation of LRT for data of complex random vectors

## Appendix B.1. Evaluating the asymptotic distribution

$\theta \in \mathbb{C}^p$ is the vector that lists the free values in $\mathbf{C}$ and $\mathbf{P}$. Recognizing the property of the MLE to asymptotically follow a Gaussian law, i.e.,

$$ \hat{\theta} \sim CN(\theta, \mathcal{I}(\theta)^{-1}), \qquad (B.1) $$

where $\mathcal{I}(\theta)$ is the Fisher information matrix whose elements are generally defined as $\mathcal{I}_{kl} = [\nabla_{\phi} \nabla_{\psi} l(\cdot)]_{kl} = -\mathbb{E}\{\frac{\partial^2}{\partial\phi_k \partial\psi_l} l(\cdot)\}$, for any two vectors $\phi$ and $\psi$, and where $l(\cdot)$ is the log-likelihood function which is assumed to be twice differentiable with respect to $\phi$ and $\psi$. Then, $\mathcal{I}(\theta) = -\mathbb{E}\{\nabla_{\theta} \nabla_{\theta'} l(\theta)\}$, assuming that the likelihood function asymptotically follows a circular complex Gaussian law, then

$$ \begin{aligned} l(\theta) &= \ln P(Y|\theta) = \ln \left( \frac{\exp(-\|Y - \mathbf{X}(\theta)\|_{\Omega(\theta)}^2)}{\pi \Omega(\theta)^l} \right) \\ &= -\|\mathbf{Y} - \mathbf{X}(\theta)\|_{\Omega(\theta)}^2 - I \ln |\Omega(\theta)| - I \ln \pi, \end{aligned} \qquad (B.2) $$

which simplifies to

$$ l(\theta) = \ln \left( \frac{\exp(-\|\theta - \hat{\theta}\|_{\mathcal{I}(\theta)^{-1}}^2)}{\pi \mathcal{I}(\theta)^{-1}} \right) = -(\theta - \hat{\theta})^H \mathcal{I}(\theta) (\theta - \hat{\theta}) + \text{constant}, \qquad (B.3) $$

hence, and recognizing that $\mathcal{I}(\theta)$ is a constant matrix,

$$ \nabla_{\theta} l(\theta) = -(\theta - \hat{\theta})^H \mathcal{I}(\theta) \qquad (B.4) $$

$$ \nabla_{\theta} \nabla_{\theta'} l(\theta) = \mathbf{0}, \qquad (B.5) $$

and similarly,

$$ \nabla_{\theta'} l(\theta) = -\mathcal{I}(\theta)(\theta - \hat{\theta}) \qquad (B.6) $$

$$ \nabla_{\theta'} \nabla_{\theta} l(\theta) = \mathbf{0}, \qquad (B.7) $$

and

$$ \nabla_{\theta} \nabla_{\theta'} l(\theta) = \nabla_{\theta'} \nabla_{\theta} l(\theta) = -\mathcal{I}(\theta). \qquad (B.8) $$

From Eq. (B.1), $\mathcal{I}(\theta)^{1/2} (\hat{\theta} - \theta) \sim CN(\mathbf{0}, \mathbf{I})$ and

$$ (\hat{\theta} - \theta)^H \mathcal{I}(\theta) (\hat{\theta} - \theta) \sim \frac{1}{2} \chi_{2p}^2, \qquad (B.9) $$

because if $X \sim CN(\mathbf{0}, 1)$, then $|X|^2 \sim \frac{1}{2} \chi_2^2$.
The test hypotheses are as
$$
\begin{cases}
H_0: \theta_s, \theta_r = 0 \\
H_1: \theta_s, \theta_r \neq 0
\end{cases}
$$

where $\theta_s$ lists the values of $\hat{C}_{M \times K}$ and $\hat{P}_{K \times R}$, and $\theta_r$ lists the values in the $K+1$ column (in the $R-K$ columns) of $\hat{C}_{M \times R}$ and $\hat{P}_{R \times R}$ in the case of the traditional test (in the case of the saturated test), respectively.

The logarithm of the Generalized Likelihood Ratio (GLR) is

$$
\ln \frac{P(\mathbf{Y}|\hat{\theta}_1)}{P(\mathbf{Y}|\hat{\theta}_0)} = \ln P(\mathbf{Y}|\hat{\theta}_1) - \ln P(\mathbf{Y}|\hat{\theta}_0) = \ell(\hat{\theta}_1) - \ell(\hat{\theta}_0), \quad (B.10)
$$

where $\hat{\theta}_i$ is the maximum-likelihood estimate of $(\theta_s, \theta_r)$ under $H_i$, $i = 0, 1$, i.e., such that $\nabla_{\theta^*} l(\theta) = 0$ at $\theta = \hat{\theta}_i$.

A general result, using a Taylor’s expansion on the second order around $\hat{\theta}_i$

$$
\begin{equation}
\begin{split}
\ell(\theta) \approx \ell(\hat{\theta}_i) &+ (\theta - \hat{\theta}_i) \nabla_{\theta^*} \ell(\hat{\theta}_i) + (\theta - \hat{\theta}_i)^T \nabla_{\theta^*} \ell(\hat{\theta}_i) \\
& + \frac{1}{2} \begin{bmatrix} (\theta - \hat{\theta}_i) \\ (\theta - \hat{\theta}_i)^* \end{bmatrix}^H \begin{bmatrix} \nabla_{\theta^*} \nabla_{\theta^*} \ell(\hat{\theta}_i) & \nabla_{\theta^*} \nabla_{\theta} \ell(\hat{\theta}_i) \\ \nabla_{\theta^*} \nabla_{\theta^*} \ell(\hat{\theta}_i) & \nabla_{\theta^*} \nabla_{\theta} \ell(\hat{\theta}_i) \end{bmatrix}^H \begin{bmatrix} (\theta - \hat{\theta}_i) \\ (\theta - \hat{\theta}_i)^* \end{bmatrix}
\end{split}
\tag{B.11}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
\ell(\theta) \approx {}& \ell(\hat{\theta}_i) + (\theta - \hat{\theta}_i) \nabla_{\theta^*} \ell(\hat{\theta}_i) + (\theta - \hat{\theta}_i)^T \nabla_{\theta^*} \ell(\hat{\theta}_i) \\
& + \frac{1}{2} (\theta - \hat{\theta}_i)^H \nabla_{\theta^*} \nabla_{\theta^*} \ell(\hat{\theta}_i) (\theta - \hat{\theta}_i) + \frac{1}{2} (\theta - \hat{\theta}_i)^T \nabla_{\theta^*} \nabla_{\theta^*} \ell(\hat{\theta}_i) (\theta - \hat{\theta}_i)^*
\end{split}
\tag{B.12}
\end{equation}
$$

Recognizing that the last term in the previous equation is scalar; hence, taking the
conjugate of the last term do not alter the equation

$$
\ell(\theta) \approx \ell(\hat{\theta}_i) + \frac{1}{2}(\theta - \hat{\theta}_i)^H \nabla_{\theta^*} \nabla_{\theta^*} \ell(\hat{\theta}_i) (\theta - \hat{\theta}_i) + \frac{1}{2}(\theta - \hat{\theta}_i)^H \nabla_{\theta^*} \nabla_{\theta^*} \ell(\hat{\theta}_i) (\theta - \hat{\theta}_i) \quad (B.13)
$$

hence, because $(\theta - \hat{\theta}_i)^H \nabla_{\theta^*} \nabla_{\theta^*} \ell(\hat{\theta}_i) (\theta - \hat{\theta}_i)$ is a real scalar,

$$
\ell(\theta) \approx \ell(\hat{\theta}_i) + (\theta - \hat{\theta}_i)^H \nabla_{\theta^*} \nabla_{\theta^*} \ell(\hat{\theta}_i) (\theta - \hat{\theta}_i) \quad (B.14)
$$

Replacing $\hat{\theta}_i$ by $\hat{\theta}_1$ and $\theta$ by $\hat{\theta}_0$

$$
\ell(\hat{\theta}_0) - \ell(\hat{\theta}_1) \approx (\hat{\theta}_0 - \hat{\theta}_1)^H \nabla_{\theta^*} \nabla_{\theta^*} \ell(\hat{\theta}_1) (\hat{\theta}_0 - \hat{\theta}_1) \quad (B.15)
$$
hence, the log-likelihood ratio defined in Eq. (B.10) is

$$ \ell(\hat{\theta}_1) - \ell(\hat{\theta}_0) \approx (\hat{\theta}_0 - \hat{\theta}_1)^H (-\nabla_{\theta} \nabla_{\theta^*} \ell(\hat{\theta}_1))(\hat{\theta}_0 - \hat{\theta}_1) \quad (B.16) $$

Recognizing that $\mathcal{I}(\hat{\theta}_1) = -\nabla_{\theta} \nabla_{\theta^*} \ell(\hat{\theta}_1)$, one gets

$$ \ell(\hat{\theta}_1) - \ell(\hat{\theta}_0) \approx (\hat{\theta}_0 - \hat{\theta}_1)^H \mathcal{I}(\hat{\theta}_1)(\hat{\theta}_0 - \hat{\theta}_1) \quad (B.17) $$

where

$$ \hat{\theta}_0 - \hat{\theta}_1 = \begin{bmatrix} \mathbf{0} \\ \hat{\theta}_{s_0} \end{bmatrix} - \begin{bmatrix} \hat{\theta}_r \\ \hat{\theta}_{s_1} \end{bmatrix}. \quad (B.18) $$

We now show that $\hat{\theta}_{s_0} - \hat{\theta}_{s_1}$ can be expanded by means of $\hat{\theta}_r$ using Eq. (B.14), with $\hat{\theta}_i = \hat{\theta}_1$, but now with $\theta = \begin{bmatrix} \mathbf{0} \\ \theta_{s_0} \end{bmatrix}$, it comes

$$ \ell(\begin{bmatrix} \mathbf{0} \\ \theta_{s_0} \end{bmatrix}) \approx \ell(\hat{\theta}_1) + \begin{bmatrix} -\hat{\theta}_r \\ \theta_{s_0} - \hat{\theta}_{s_1} \end{bmatrix}^H (-\mathcal{I}(\hat{\theta}_1)) \begin{bmatrix} -\hat{\theta}_r \\ \theta_{s_0} - \hat{\theta}_{s_1} \end{bmatrix}. \quad (B.19) $$

Recognizing $\mathcal{I}(\hat{\theta}_1)$ as a matrix of four block, where the diagonal block $\mathcal{I}(\hat{\theta}_1)_{11}$ corresponds to the element $-\hat{\theta}_r$, whereas the diagonal block $\mathcal{I}(\hat{\theta}_1)_{22}$ corresponds to the element $\theta_{s_0} - \hat{\theta}_{s_1}$,

$$ \ell(\begin{bmatrix} \mathbf{0} \\ \theta_{s_0} \end{bmatrix}) \approx \ell(\hat{\theta}_1) - \begin{bmatrix} -\hat{\theta}_r \\ \theta_{s_0} - \hat{\theta}_{s_1} \end{bmatrix}^H \begin{bmatrix} \mathcal{I}(\hat{\theta}_1)_{11} & \mathcal{I}(\hat{\theta}_1)_{12} \\ \mathcal{I}(\hat{\theta}_1)_{21} & \mathcal{I}(\hat{\theta}_1)_{22} \end{bmatrix} \begin{bmatrix} -\hat{\theta}_r \\ \theta_{s_0} - \hat{\theta}_{s_1} \end{bmatrix}. \quad (B.20) $$

The MLE $\hat{\theta}_{s_0}$ is the value that maximizes $\ell(\begin{bmatrix} \mathbf{0} \\ \theta_{s_0} \end{bmatrix})$, hence, taking the gradient of the above equation w.r.t. $\theta_{s_0}^H$ and setting it to zero gives

$$ \nabla_{\theta_{s_0}^H} \ell(\begin{bmatrix} \mathbf{0} \\ \theta_{s_0} \end{bmatrix}) \approx \nabla_{\theta_{s_0}^H} \ell(\hat{\theta}_1) - \nabla_{\theta_{s_0}^H} (\hat{\theta}_r^H \mathcal{I}(\hat{\theta}_1)_{11} \hat{\theta}_r + (\theta_{s_0} - \hat{\theta}_{s_1})^H \mathcal{I}(\hat{\theta}_1)_{22} (\theta_{s_0} - \hat{\theta}_{s_1}) \\ - \hat{\theta}_r^H \mathcal{I}(\hat{\theta}_1)_{12} (\theta_{s_0} - \hat{\theta}_{s_1}) - (\theta_{s_0} - \hat{\theta}_{s_1})^H \mathcal{I}(\hat{\theta}_1)_{21} \hat{\theta}_r) \quad (B.21) $$

$$ 0 \approx 0 - (0 + \mathcal{I}(\hat{\theta}_1)_{22}(\theta_{s_0} - \hat{\theta}_{s_1}) - 0 - \mathcal{I}(\hat{\theta}_1)_{21}\hat{\theta}_r) \quad (B.22) $$

$$ \hat{\theta}_{s_0} - \hat{\theta}_{s_1} \approx \mathcal{I}(\hat{\theta}_1)_{22}^{-1} \mathcal{I}(\hat{\theta}_1)_{21} \hat{\theta}_r. \quad (B.23) $$
Therefore, plugging Eq. (B.23) into Eq. (B.18), one gets

$$ \hat{\theta}_0 - \hat{\theta}_1 = \left[ \mathcal{I}(\hat{\theta}_1)_{22}^{-1} \mathcal{I}(\hat{\theta}_1)_{21} \right] \hat{\theta}_r \quad (B.24) $$

and then plugging this into the log-likelihood ratio in Eq. (B.17)

$$ \begin{aligned} \ell(\hat{\theta}_1) - \ell(\hat{\theta}_0) &\approx \hat{\theta}_r^H [-\mathbf{I} \quad \mathcal{I}(\hat{\theta}_1)_{21}^H \mathcal{I}(\hat{\theta}_1)_{22}^{-1}] \left[ \begin{array}{cc} \mathcal{I}(\hat{\theta}_1)_{11} & \mathcal{I}(\hat{\theta}_1)_{12} \\ \mathcal{I}(\hat{\theta}_1)_{21} & \mathcal{I}(\hat{\theta}_1)_{22} \end{array} \right] \left[ \begin{array}{c} -\mathbf{I} \\ \mathcal{I}(\hat{\theta}_1)_{22}^{-1} \mathcal{I}(\hat{\theta}_1)_{21} \end{array} \right] \hat{\theta}_r \\ &\approx \hat{\theta}_r^H (\mathcal{I}(\hat{\theta}_1)_{11} - \mathcal{I}(\hat{\theta}_1)_{21}^H \mathcal{I}(\hat{\theta}_1)_{22}^{-1} \mathcal{I}(\hat{\theta}_1)_{21}) \hat{\theta}_r \end{aligned} \quad (B.25) $$

where the quantity inside the brackets is recognized as the inverse of $\mathcal{I}(\hat{\theta}_1)^{11}$ which is the first diagonal block in $\mathcal{I}(\hat{\theta}_1)^{-1}$ (i.e., $\mathrm{Cov}(\hat{\theta}_1)$), then the above equation becomes

$$ \ell(\hat{\theta}_1) - \ell(\hat{\theta}_0) \approx \hat{\theta}_r^H [\mathcal{I}(\hat{\theta}_1)^{11}]^{-1} \hat{\theta}_r \quad (B.26) $$

But from Eq. (B.1),

$$ \hat{\theta}_r \sim CN(\theta_r, \mathcal{I}(\hat{\theta}_1)^{11}), \quad (B.27) $$

and therefore from Eq. (B.9),

$$ \hat{\theta}_r^H [\mathcal{I}(\hat{\theta}_1)^{11}]^{-1} \hat{\theta}_r \sim \frac{1}{2} \chi_{2\nu}^2 \quad (B.28) $$

with $\nu = \dim(\theta_r)$. Hence

$$ \ell(\hat{\theta}_1) - \ell(\hat{\theta}_0) \sim \frac{1}{2} \chi_{2\nu}^2. \quad (B.29) $$

## Appendix B.2. Counting the degrees of freedom for the conventional LRT

Recall that, the estimated quantity in a reduced-rank reference solution is $CP^H$, where $C \in \mathbb{C}^{M \times K}$ and $P \in \mathbb{C}^{R \times K}$, which contains a total of $2MK + 2RK$ complex elements with only $2(M + R)$ elements corresponding to $2 \times \dim(\theta_r)$, in the case of the conventional LRT. The value of these elements are not totally free, as they are subjected to the constraint $P^HMP = I_{K \times K}$, which imposes $-2(K^2 - K)$ real and imaginary free parameters, relating to the complex off-diagonal elements, and $-K$ real ones, relating to the diagonal elements. However, due to the hermitian property of the matrix $P^HMP$, there remain only $-(K^2 - K)$ free parameters, relating to the complex of diagonals, along with the $-K$ real ones, thus there is a total
of $-(K^2 - K) - K = -K^2$ constraint-related free parameters. Recognizing the fact that, the same constraint applies to the matrix $C$ as $C = S_{YR}P(P^H S_{RR}P)^{-1}$, then the correct total number of constraint-related free parameters is $-2K^2$, where $K$ represents here the number of preserved virtual sources in the reduced model. Similarly, there exist $-2(K + 1)^2$ constraint-related free parameters in the tested model of $K + 1$ preserved sources. Therefore when comparing a model with $K$ sources to a model with $K + 1$ sources, the degree of freedom ($2\nu$) of the $\chi_{2\nu}^2$ is as follows

$$ \begin{aligned} 2\nu &= 2(M + R) - 2(K + 1)^2 - (-2K^2) \\ &= 2(M + R) - 4K - 2. \end{aligned} \qquad (B.30) $$

### Appendix B.3. Counting the degrees of freedom for the modified LRT

Recall that, the estimated quantity in a reduced-rank reference solution is $CP^H$, where $C \in \mathbb{C}^{M \times K}$ and $P \in \mathbb{C}^{R \times K}$, which contains a total of $2MK + 2RK$ complex elements with only $2M(R-K)+2R(R-K)$ elements corresponding to $2 \times \dim(\theta_r)$, in the case of the proposed modified LRT against the saturated model. The value of these elements are not totally free, as they are subjected to the constraint $P^HMP = I_{K \times K}$, which imposes $-2(K^2 - K)$ real and imaginary free parameters, relating to the complex off-diagonal elements, and $-K$ real ones, relating to the diagonal elements. However, due to the hermitian property of the matrix $P^HMP$, there remain only $-(K^2 - K)$ free parameters, relating to the complex off-diagonals, along with the $-K$ real ones, thus there is a total of $-(K^2-K)-K = -K^2$ constraint-related free parameters. Recognizing the fact that, the same constraint applies to the matrix $C$ as $C = S_{YR}P(P^H S_{RR}P)^{-1}$, then the correct total number of constraint-related free parameters is $-2K^2$, where $K$ represents here the number of preserved virtual sources in each of the tested models. Similarly, there are $-2K^2 = -2R^2$ constraint-related free parameters for the saturated model, of $R$ preserved virtual sources. Therefore when comparing a sub-model with $K$ sources to the saturated model (i.e., with $R$ sources), the degree of freedom ($2\nu$) of the $\chi_{2\nu}^2$ is as follows

$$ \begin{aligned} 2\nu &= 2M(R - K) + 2R(R - K) - 2R^2 - (-2K^2) \\ &= 2(R - K)(M - K). \end{aligned} \qquad (B.31) $$

### References

[1] D. Halman, J. Bolton, Multi-reference nearfield acoustical holography, in: PrOC, of Inter-Noise, volume 92, 1992, pp. 1165–1170.
[2] R. Wang, Y. Zhang, L. Yu, J. Antoni, Q. Leclère, W. Jiang, A probability model with Variational Bayesian Inference for the complex interference suppression in the acoustic array measurement, Mech. Syst. Signal Process. 191 (2023) 110181. DOI: 10.1016/j.ymssp.2023.110181

[3] 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. DOI: https://doi.org/10.1016/j.apacoust.2021.108021

[4] B. Dong, J. Antoni, E. Zhang, Blind separation of sound sources from the principle of least spatial entropy, J. Sound Vibration 333 (2014) 2643–2668. DOI: 10.1016/j.jsv.2013.12.011

[5] Q. Leclère, Multi-channel spectral analysis of multi-pass acquisition measurements, Mech. Syst. Signal Process. 23 (2009) 1415–1422. DOI: https://doi.org/10.1016/j.ymssp.2008.12.002

[6] 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 Vibration 349 (2015) 123–149. DOI: 10.1016/j.jsv.2015.03.008

[7] J. Bendat, A. Piersol, Random Data: Analysis and Measurement Procedures, Wiley Series in Probability and Statistics, Wiley, 2011. URL: https://books.google.fr/books?id=qYSViFRNM1wC.

[8] M. Lee, J. S. Bolton, Scan-based near-field acoustical holography and partial field decomposition in the presence of noise and source level variation, The J. Acoust. Soc. America 119 (2006) 382–393. DOI: 10.1121/1.2133717

[9] C. Dodds, J. Robson, Partial coherence in multivariate random processes, J. Sound Vibration 42 (1975) 243–249. DOI: https://doi.org/10.1016/0022-460X(75)90219-9

[10] R. Fan, Z. Su, G. Meng, C. He, Application of sound intensity and partial coherence to identify interior noise sources on the high speed train, Mech. Syst. Signal Process. 46 (2014) 481–493. DOI: 10.1016/j.ymssp.2013.11.014
[11] S. Price, R. Bernhard, Virtual coherence: A digital signal processing technique for incoherent source identification, in: Proceedings of IMAC, volume 4, 1986, pp. 3–6.

[12] X. Pan, W. Jiang, The separation and identification of distributed noise sources with local coherence by a subspace equivalent-source method, Mech. Syst. Signal Process. 168 (2022) 108601. DOI: 10.1016/j.ymssp.2021.108601

[13] M. N. Albezzawy, J. Antoni, Q. Leclere, The Maximally-Coherent Reference Technique: a solution to deal with large sets of references, in: Proceedings of ISMA2022, Leuven, Belgium, 2022. URL: https://hal.science/hal-04063060.

[14] M. N. Albezzawy, J. Antoni, Q. Leclere, Resynchronization of sequential measurements using the Maximally-Coherent Reference technique, in: Surveillance, Vibrations, Shock and Noise, Institut Supérieur de l'Aéronautique et de l'Espace [ISAE-SUPAERO], Toulouse, France, 2023. URL: https://hal.science/hal-04165667.

[15] R. Balan, Estimator for Number of Sources Using Minimum Description Length Criterion for Blind Sparse Source Mixtures, in: M. E. Davies, C. J. James, S. A. Abdallah, M. D. Plumbley (Eds.), Independent Component Analysis and Signal Separation, volume 4666, Springer Berlin Heidelberg, Berlin, Heidelberg, 2007, pp. 333–340. doi:10.1007/978-3-540-74494-8-42, series Title: Lecture Notes in Computer Science.

[16] M. Chiani, M. Z. Win, Estimating the number of signals observed by multiple sensors, in: 2010 2nd International Workshop on Cognitive Information Processing, IEEE, Elba Island, Italy, 2010, pp. 156–161. URL: http://ieeexplore.ieee.org/document/5604227/. doi:10.1109/CIP.2010.5604227.

[17] L. Sun, Q. Cheng, Indoor sound source localization and number estimation using infinite Gaussian mixture models, in: 2014 48th Asilomar Conference on Signals, Systems and Computers, IEEE, Pacific Grove, CA, USA, 2014, pp. 1189–1193. URL: http://ieeexplore.ieee.org/document/7094646/. doi:10.1109/ACSSC.2014.7094646.
[18] Y. Kawakami, T. Ishibashi, Estimation for the Number of Sound Sources Using Directional Histogram of Observed Mixture Signals, in: Proceedings of The 7th IIAE International Conference on Industrial Application Engineering 2019, The Institute of Industrial Applications Engineers, 2019, p. 225. doi:10.12792/iciae2019.043.

[19] L. Yu, J. Antoni, J. Deng, C. Li, W. Jiang, Low-rank Gaussian mixture modeling of space-snapshot representation of microphone array measurements for acoustic imaging in a complex noisy environment, Mech. Syst. Signal Process. 165 (2022) 108294. DOI: 10.1016/j.ymssp.2021.108294

[20] M. Lee, J. S. Bolton, L. Mongeau, Application of cylindrical near-field acoustical holography to the visualization of aeroacoustic sources, The J. Acoust. Soc. America 114 (2003) 842–858. DOI: 10.1121/1.1587735

[21] Q. Leclere, 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. MAG ID: 3136681030 S2ID: 0686e7587e76b19fccec5d88191e3972c726c3c1. DOI: 10.1016/j.apacoust.2021.108021

[22] G. Carter, A. Nuttall, Statistics of the estimate of coherence, Proc. IEEE 60 (1972) 465–466. DOI: 10.1109/PROC.1972.8671

[23] S. S. Wilks, The large-sample distribution of the likelihood ratio for testing composite hypotheses, The Ann. Math. Stat. 9 (1956) 60–62. DOI: 10.1214/aoms/1177732360

[24] M. S. Bartlett, A note on the multiplying factors for various $\chi^2$ approximations, J. Royal Stat. Soc. Ser. B (Methodological) 16 (1954) 296–298. Publisher: John Wiley & Sons, Ltd. DOI: 10.1111/j.2517-6161.1954.tb00174.x

[25] D. N. Lawley, Tests of significance for the latent roots of covariance and correlation matrices, Biometrika 43 (1956) 128–136. MAG ID: 1996340537. DOI: 10.1093/biomet/43.1-2.128

[26] H. Akaike, Information theory and an extension of the maximum likelihood principle, in: Proceeding of the second international Symposium on Information Theory, 1973, pp. 199–213. doi:10.1007/978-1-4612-0919-5-38.
[27] J. Rissanen, Modeling by shortest data description, Automatica 14 (1978) 465-471. DOI: 10.1016/0005-1098(78)90005-5

[28] G. Schwarz, Estimating the dimension of a model, Ann. Stat. 6 (1978) 461-464. MAG ID: 2168175751. DOI: 10.1214/aos/1176344136

[29] P. Stoica, Y. Selen, Jian Li, J. Li, On information criteria and the generalized likelihood ratio test of model order selection, IEEE Signal Process. Lett. 11 (2004) 794-797. MAG ID: 2144106004 S2ID: e2a45652457879e54c967362eb2091eb7fcb5c0d. DOI: 10.1109/lsp.2004.835468

[30] K. M. Wong, Q. -T. Zhang, J. P. Reilly, P. C. Yip, On information theoretic criteria for determining the number of signals in high resolution array processing, IEEE Trans. on Acoust. Speech, Signal Process. 38 (1990) 1959-1971. DOI: 10.1109/29.103097

[31] R. R. Nadakuditi, A. Edelman, Sample eigenvalue based detection of high-dimensional signals in white noise using relatively few samples, IEEE Trans. on Signal Process. 56 (2008) 2625-2638. DOI: 10.1109/TSP.2008.917356

[32] S. Kritchman, B. Nadler, Determining the number of components in a factor model from limited noisy data, Chemom. Intell. Lab. Syst. 94 (2008) 19-32. DOI: 10.1016/j.chemolab.2008.06.002

[33] G. Huaijin, Estimating the number of signals and signal resolution, IEEE Trans. on Signal Process. 46 (1998) 2267-2270. DOI: 10.1109/78.705462

[34] Q. T. Zhang, K. M. Wong, P. C. Yip, J. P. Reilly, Statistical analysis of the performance of information theoretic criteria in the detection of the number of signals in array processing, IEEE Trans. on Acoust. Speech, Signal Process. 37 (1989) 1557-1567. DOI: 10.1109/29.35394

[35] S. S. Wilks, On the independence of k sets of normally distributed statistical variables, Econometrica 3 (1935) 309. MAG ID: 2314024971. DOI: 10.2307/1905324

[36] 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. doi:10.48550/ARXIV.2008.03971.

[37] S. C. Larson, The shrinkage of the coefficient of multiple correlation., J. Educ. Psychol. 22 (1931) 45–55. Place: US Publisher: Warwick & York. DOI: 10.1037/h0072400

[38] M. H. Quenouille, Approximate tests of correlation in time-series, J. Royal Stat. Soc. Ser. B (Methodological) 11 (1949) 68–84. Publisher: [Royal Statistical Society, Wiley].

[39] J. W. Tukey, Bias and confidence in not quite large samples (abstract), The Ann. Math. Stat. 29 (1958) 614–623. DOI: 10.1214/aoms/1177706647

[40] B. Efron, Bootstrap methods: Another look at the jackknife, The Ann. Stat. 7 (1979) 1–26. Publisher: Institute of Mathematical Statistics.

[41] R. J. Wherry, A new formula for predicting the shrinkage of the coefficient of multiple correlation, The Ann. Math. Stat. 2 (1931) 440–457. Publisher: Institute of Mathematical Statistics.

[42] C. I. Mosier, I. problems and designs of cross-validation 1, Educ. Psychol. Meas. 11 (1951) 5–11. Publisher: SAGE Publications Inc. DOI: 10.1177/001316445101100101

[43] B. Efron, The Jackknife, the Bootstrap and Other Resampling Plans, Society for Industrial and Applied Mathematics, 1982. doi:10.1137/1.9781611970319, eprint: https://epubs.siam.org/doi/pdf/10.1137/1.9781611970319.

[44] S. YOON, P. NELSON, A METHOD FOR THE EFFICIENT CONSTRUCTION OF ACOUSTIC PRESSURE CROSS-SPECTRAL MATRICES, J. Sound Vibration 233 (2000) 897–920. DOI: 10.1006/jsvi.1999.2888

[45] J. Bendat, Solutions for the multiple input/output problem, J. Sound Vibration 44 (1976) 311–325. DOI: https://doi.org/10.1016/0022-460X(76)90504-6

[46] M. N. Albezzawy, J. Antoni, Q. Leclère, The maximally-coherent reference technique and its application to sound source extraction without synchronous measurements, J. Sound Vibration 599 (2025) 118896. DOI: https://doi.org/10.1016/j.jsv.2024.118896
[47] M. Peligrad, W. B. Wu, Central limit theorem for fourier transforms of stationary processes, The Ann. Probab. 38 (2010) 2009–2022.

[48] 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. on Audio Electroacoustics 15 (1967) 70–73.
DOI: 10.1109/TAU.1967.1161901
# Highlights

## Source enumeration in reduced-rank reference techniques for source extraction

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

* A unified formalism for all the reference techniques for source extraction.

* Wilk's theorem is reformulated for vectors of complex random variables.

* A robust likelihood ratio test against the saturated model for source enumeration.

* A parametric bootstrap technique as an alternative in the case of short records.

* A cross-validation alternative, and the performance of the techniques is compared.
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.