{'fraction_non_alphanumeric': 0.05893923185755202, 'fraction_numerical': 0.01235377342213327, 'mean_word_length': 4.282927718530671, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Blind Direction-of-Arrival Estimation in Acoustic Vector-Sensor Arrays via Tensor Decomposition and Kullback-Leibler Divergence Covariance Fitting

A blind Direction-of-Arrivals (DOAs) estimate of narrowband signals for Acoustic Vector-Sensor (AVS) arrays is proposed. Building upon the special structure of the signal measured by an AVS, we show that the covariance matrix of all the received signals from the array admits a natural low-rank 4-way tensor representation. Thus, rather than estimating the DOAs directly from the raw data, our estimate arises from the unique parametric Canonical Polyadic Decomposition (CPD) of the observations' Second-Order Statistics (SOSs) tensor. By exploiting results from fundamental statistics and the recently re-emerging tensor theory, we derive a consistent blind CPD-based DOAs estimate without prior assumptions on the array configuration. We show that this estimate is a solution to an equivalent approximate joint diagonalization problem, and propose an ad-hoc iterative solution. Additionally, we derive the Cram\'er-Rao lower bound for Gaussian signals, and use it to derive the iterative Fisher scoring algorithm for the computation of the Maximum Likelihood Estimate (MLE) in this particular signal model. We then show that the MLE for the Gaussian model can in fact be used to obtain improved DOAs estimates for non-Gaussian signals as well (under mild conditions), which are optimal under the Kullback-Leibler divergence covariance fitting criterion, harnessing additional information encapsulated in the SOSs. Our analytical results are corroborated by simulation experiments in various scenarios, which also demonstrate the considerable improved accuracy w.r.t. a competing state-of-the-art blind estimate for AVS arrays, reducing the resulting root mean squared error by up to more than an order of magnitude.


I. INTRODUCTION
An Acoustic Vector-Sensor (AVS) is comprised of an omni-directional microphone and two/three particle velocity transducers aligned along the orthogonal coordinate axes in a two/three dimensional space [2]. It measures the acoustic pressure as well as the acoustic particle velocities in each one of the coordinate axes, thus providing a full description of the acoustic field at a given location in space [3], [4]. As one would expect, in the context of passive array-processing, AVS arrays enable enhanced capabilities relative to their equivalentaperture traditional (scalar) microphone arrays [5]- [7]. Consequently, and since these devices are already practically feasible [8], AVS arrays are the cornerstone in a wide variety of emerging applications such as aircraft acoustic detection, The author is with the Dept. of Computer Science and Applied Mathematics, Faculty of Mathematics and Computer Science, Weizmann Institute of Science, 234 Herzl Street, Rehovot 7610001, Israel, e-mail: amir.weiss@weizmann.ac.il. localization and tracking [9], [10], battlefield acoustics classification [11], and underwater acoustic communication [12] and source enumeration [13], to name but a few.
Within many of these applications, Direction-of-Arrival (DOA) estimation is a fundamental task, which plays a key role in the overall successful operation (e.g., localization, which reduces to DOA and range estimation). Of course, as one of the most prominent signal processing problems in general, and specifically in array processing, DOA estimation in the context of AVS arrays has already been extensively addressed in the literature during the past two and a half decades, as briefly reviewed in what follows.

A. Previous Work: DOA Estimation with AVSs
In their seminal paper [2], Nehorai and Paldi derived the Cramér-Rao Lower Bound (CRLB) on the Mean Squared Error (MSE) of any unbiased DOAs estimate, and proposed two algorithms-the Intensity-Based and the Velocity-Covariance-Based Algorithms-for DOA estimation, though for a single source with a single AVS only. For this (limited) scenario, the maximum steered response power and the Maximum Likelihood Estimate (MLE, as a special case of the former), were derived by Levin et al. in [14] and [15], resp. Multisource DOA estimation in a reverberant environment, still using a single AVS, has been recently addressed by Wu et al. [16], where low-reverberant-single-source points in the timefrequency domain are exploited.
For the extended Multiple-Sources Multiple-Sensors (MSMS) scenario, Hawkes and Nehorai considered in [17] both the conventional and the minimum-variance distortionless response beamforming DOAs estimates to demonstrate the improvement attained by using AVSs rather than traditional pressure sensors. Following this work, an abundance of methods have been proposed for various specific scenarios (whether for array-or signal-related properties), such as linear [18], circular [19]- [21], sparse [22]- [24] and nested [25]- [27] arrays, coherent signals [28]- [30], one-bit measurements [31] and a variety of others [32]- [34]. However, all these methods require perfect or (at least) partial prior knowledge of the array configuration, or, equivalently 1 , of the steering vectors parametric structure, in particular w.r.t. the sources' DOAs. Hence, these methods are typically sensitive to model inaccuracies/errors of this sort. This is exactly the point where our work comes into play w.r.t. to our novel contributions.
In this work, we consider blind DOAs estimation in AVS arrays within the framework of narrowband signals. Here, the term "blind" implies that no prior assumption on the array configuration (/geometry) is used, hence the steering vectors parametric structure is deemed unknown. In fact, blind DOA estimation in the context of AVSs has so far been only sparsely addressed in the literature. A direction-finding and blind interference rejection algorithm was proposed by Wong in [35], but only for up to three fast frequency-hop spread spectrum signals of unknown hop sequences, and using only a single AVS. Xiao et al. also proposed a blind DOA estimate for a single AVS only [36], based on the Joint Approximate Diagonalization of Eigen-matrices (JADE, [37]) algorithm, and therefore does not provide a solution for Gaussian signals.
For AVS arrays, Zhang et al. proposed the Trilinear Decomposition-based (TriD) blind DOAs estimate for incoherent signals [1], and a similar TriD approach in [38] for coherent signals. However, the trilinear decomposition considered both in [1] and [38] is exact only for the (less practical) noiseless signal model. Although the TriD methods perform quite well, the DOA estimates which stem from this approach were not shown analytically to be consistent or optimal. In contrast, on top of providing superior performance, our novel estimate is shown analytically (and demonstrated empirically) to be consistent regardless of the SNR conditions, and provides optimal performance for Gaussian signals for any SNR.

B. Blind DOAs Estimation: Motivation and Contributions
The motivation to address this blind scenario and develop a solution algorithm under this framework arises from several considerations. Firstly, by not assuming a specific array configuration (geometry/structure), the resulting solution algorithm could properly operate in different systems with arbitrary array configurations, and without the need for specific tuning prior to operation. Secondly, unknown inter-AVSs' gain and/or phase offsets (e.g., due to sensors mis-locations), which otherwise require a calibration procedure (e.g., [39]), are totally transparent to such a blind estimate. For example, many of the DOA estimation methods tailored to Uniform Linear Arrays (ULAs) exploit the special Toeplitz structure of the observations' spatial covariance matrix (e.g., [40] as one representative example). Clearly, the performance of such methods deteriorates rapidly in the presence of sensors error positioning, in contrast to the performance of a blind estimate, which remains (almost) indifferent to these errors. In other words, a blind algorithm is robust w.r.t. sensors error positioning. This robustness implies a significant practical advantage over other algorithms in terms of simplicity for the end-user, and may save time and resources, for example, when blind ("online") calibration is not possible and ad-hoc transmission of a calibrating source is required. Thirdly, most algorithms are developed for an underlying signal model based on some physical approximation, e.g., the "near-field" or "far-field" approximations (e.g., [41] in the context of AVSs). As a result, in these cases the steering vectors' parametric representation (in particular, w.r.t. the DOAs) is also only an approximation, which brings along with analytical convenience an inherent modeling error [42], [43]. This modeling error impairs the per-formance even for an optimal solution (in some well-defined sense) for this particular, approximated model. In contrast, a blind approach, in which the solution is not developed based on such (potential) modeling errors, would yield model-based errors free 2 DOA estimates. Lastly, a blind estimate can also successfully cope with faulty elements in the array [44], and maintains proper functionality for partially damaged arrays.
Motivated by the merits above, in this work we propose a novel blind DOA estimation algorithm for AVS arrays with arbitrary configurations, while making only a few a-priori assumptions on the signal model. Building upon the AVS measurement model, we exploit the special (block) structure of the observations' covariance matrix, which naturally lends itself to the recently flourishing tensor formulation in the signal processing literature (e.g., [45]- [52]). As a natural continuum thereof, we employ a statistical approach and invoke tensorcalculus-related results, leading to a consistent DOAs estimate for any SNR. We then show that this estimate can be improved by a second refinement phase via Kullback-Leibler Divergence (KLD) covariance fitting, which yields our proposed estimate. The main contributions of this paper are as follows: • Consistent blind DOAs estimation via tensor decomposition: Based on the observations' empirical covariance matrix and a consistent noise variance estimate, we show that joint estimation of all the DOAs and their associated steering vectors in a MSMS scenario is (asymptotically) equivalent to a parametric Canonical Polyadic Decomposition (CPD, see [53] and reference therein)-sometimes termed as "tensor rank decomposition" or "parallel factor model"-of a 4-mode tensor statistic. We show that the uniqueness theorem of quadrilinear decompositions of 4mode arrays due to Sidiropoulos and Bro (Theorem 2 in [54]) grants this CPD-based estimate its consistency. • Iterative solution algorithm of the CPD-based estimate: We show that computation the CPD-based estimate amounts to a partially-parametric Approximate Joint Diagonalization (AJD) problem. Accordingly, we proposed an iterative solution, which is a modified version of the Alternating Columns-Diagonal Centers (AC-DC) algorithm proposed by Yeredor in [55]. Our modified algorithm results in more accurate DOAs estimates than the ones a generic AJD algorithm would yield due to our tailored parametric adaptation. Further, the algorithm inherently yields estimates of the steering vectors, considered as nuisance parameters, which consequently also enable consistent blind separation of the latent sources. • Performance bounds and optimal estimation for Gaussian signals: For the particular case of Gaussian signals, we derive the CRLB on the MSE matrix of any unbiased estimate in joint estimation of all the unknown deterministic model parameters, namely the steering vectors, the DOAs and the noise variance. In addition, based on the Fisher Information Matrix (FIM), and using the (already obtained) CPD-based consistent estimates as initial solutions, we propose a Maximum Likelihood (ML) "refine-ment" phase, in which all the unknown parameters' MLEs are pursued via the Fisher Scoring Algorithm (FSA, [56]). As demonstrated in simulations, these refined estimates are asymptotically efficient, attaining the CRLB. • KLD covariance fitting enhancement: We show that the MLEs for the particular Gaussian signal model are also optimal under the KLD covariance fitting criterion regardless of the underlying signal model, and can therefore be used for non-Gaussian signals as well in order to achieve significant performance enhancement. Accordingly, our final proposed DOAs estimate enjoys higher accuracy and robustness to the underlying signal model. The rest of this paper is organized as follows. The following subsection contains an outline of our notations. In Section II we present the model under consideration and formulate our blind DOA estimation problem. The CPD-based (phase 1) estimates are presented in Section III, followed by the iterative solution algorithm for their actual computation in Section IV. We then consider in Section V the Gaussian signal model, and derive its respective CRLB, as well as the update equations of the FSA for the computation of the MLEs. Our proposed KLD-based (phase 2) estimates are presented in Section VI, followed by simulation results in Section VII, substantiating and demonstrating empirically our analytical results. Concluding remarks are given in Section VIII.

C. Notations and Preliminaries
We use x, x, X and X for a scalar, column vector, matrix and tensor, resp. The superscripts (·) T , (·) * , (·) † , (·) −1 and (·) + denote the transposition, complex conjugation, conjugate transposition, inverse and Moore-Penrose pseudo-inverse operators, resp. We use I K to denote the K × K identity matrix, and the pinning vector e k denotes the k-th column of the identity matrix with context-dependent dimension. Further, δ k e T k e denotes the Kronecker delta of k and . E[·] denotes expectation, the Diag(·) operator forms an M × M diagonal matrix from its M -dimensional vector argument, and 0 M ∈ R M ×1 is the all-zeros vector. The Kronecker, Khatri-Rao (column-wise Kronecker) and tensor outer products (e.g., [49]) are denoted by ⊗, and •, resp. We use  (a dotless j) to denote √ −1; The operators {·} and {·} denote the real and imaginary parts (resp.) of their complex-valued argument. rank(Q) denotes the rank of the matrix Q. The Frobenius and 2 norms are denoted by · F and · 2 , resp. Convergence in probability and in distribution are denoted by The function atan2(y, x) returns the principal value of the argument function applied to the complex number x + y. The gradient of a matrix function F (x) ∈ C M ×N w.r.t. its scalar argument x ∈ R is denoted by ∇ x F ∈ C M ×N . Conversely, the gradient of a scalar function x(F ) ∈ R w.r.t. its matrix argument F ∈ C M ×N is denoted by ∇ F x ∈ C M ×N . The vec(·) operator concatenates the columns of an M × N matrix into an M N × 1 column vector. The vec * (·) operator, defined only for Hermitian matrices, is the invertible transformation which concatenates the columns of its M × M Hermitian matrix argument into an (M (M + 1)/2) × 1 column vector, but takes each element (conjugately) "duplicated" by conjugate symmetry only once, on its first occurrence. cum(w, x, y, z) denotes the fourth-order joint cumulant of its four scalar random variable arguments. Finally, as we make use of the somewhat less familiar Kruskal Rank, for convenience, we bring its definition, given as follows.
Rank of Q, denoted by k Q , is r if and only if every r columns are linearly independent, and this fails for at least one set of r + 1 columns. It follows that k Q ≤ rank(Q) ≤ min(I, D).

II. PROBLEM FORMULATION
Consider an array of M AVSs, where each AVS consists of three elements, one pressure and two particle velocity transducers in two perpendicular directions. The configuration of the array, which is not constrained to a particular structure (e.g., uniform linear), is assumed as unknown, which dictates a "blind" setup in this respect. Further, consider the presence of D < M − 1 unknown narrowband sources, centered around some common carrier frequency with a wavelength λ, where we assume that the number of sources D is known. Assuming the received signals are down-converted, Low-Pass Filtered (LPF) 3 and sampled at least at the Nyquist rate, the vector of sampled baseband signals from all 3M sensors is given by is the vector of sources impinging on the array from unknown azimuth angles θ [θ 1 · · · θ D ] T ∈ [−π, π) D×1 , assumed as distinct from one another, i.e, ∀d = : [a(θ 1 ) · · · a(θ D )] ∈ C 3M ×D is the array manifold matrix, whose columns are the steering vectors in which {a(θ d ) ∈ C M ×1 } D d=1 are the unknown equivalent acoustic pressure sensor array steering vectors [2]; (iii) v[t] ∈ C 3M ×1 is an additive noise vector, spatially and temporally independent, identically distributed (i.i.d.) zero-mean circular Complex Normal (CN) [57] with a covariance matrix where R x is the covariance matrix of the noiseless signal x[t].
Thus, the problem at hand can be formulated as follows:

III. PHASE 1: THE CPD-BASED BLIND DOAS ESTIMATE
Let us begin with a general, bird's-eye view description of our strategy for the proposed solution, which stems from two fundamental observations. First, observe that since the array configuration is assumed unknown, the "core" (pressure) steering vectors {a(θ d )} D d=1 are unknown, and therefore {a(θ d )} D d=1 are as well. However, the structure of each a(θ d ) as a function of a(θ d ) given in (2), which is determined inherently by the AVS's basic structure, is not only known, but also encapsulates the dependence on the desired DOA θ d (via cos(θ d ) and sin(θ d )) regardless of the particular parametric structure of a(θ d ). Next, notice that while the measured signal (1), which is a function of the DOAs θ, is random for all t ∈ {1, . . . , T }, the covariance matrix R y , given in (3), is a deterministic function of θ. Therefore, we shall work towards writing the explicit dependence of R y in θ based on the special structure of an AVS steering vector, as prescribed by (2). Since the empirical covariance matrix " is a consistent estimate of R y (under mild conditions), we may exploit the aforementioned special dependence to derive our blind, consistent DOAs estimate based on " R y only. More specifically, notice first that due to (2), we may write the array manifold matrix A(θ) as where A(θ) [a(θ 1 ) · · · a(θ D )] ∈ C M ×D , cos(θ) and sin(θ) operate elementwise, and we have defined the auxiliary matrix which is completely determined by θ only. Since in our framework the "core" steering vectors' parametric structure {a(θ d )} D d=1 is assumed unknown, we denote for brevity hereafter a(θ d ) := a d for all d ∈ {1, . . . D}, and A := A(θ) accordingly. Note, however, that we intentionally keep the notation A(θ), as its dependence on θ via C(θ), regardless of {a d } D d=1 , is known and given by (4). Note further that although A(θ) has 3M D complex-valued elements, (4) implies that it is completely determined only by the unknowns {a d , θ d } D d=1 , namely 2M D + D free parameters (/ degrees of freedom). This special, economical structure will be exploited shortly.
Next, let us consider the covariance matrix R x of the noiseless signal x[t]. For this, observe first that both A and the sources powers, i.e., the diagonal elements of the (diagonal) matrix R s , are unknown. Thus, without loss of generality 5 (w.l.o.g.), we may assume that R s = I D . With this, R x reads where we have used the mixed product rule (e.g., [49], Section II) in moving from (7) to (8), and defined the rank-1 matrices for all i, j ∈ {1, 2, 3}. Indeed, since {A d } are rank-1 matrices, one can also write the (i, j, m, n)-th element of R x as, to conclude that R x is a 4-mode tensor of rank D.

Remarks:
i. Notice that while M , the number of AVSs in the array, may be very large, the number of sources D is typically smaller. Hence, since R x is a rank-D tenor, it admits a CPD of order D, and is therefore another compact, economical (and equivalent) representation of R x . ii. The dorsal slab R x (i, j, :, :) is the auto-/ cross-covariance matrix between similar / different types of sensors. Thus, R x (i, i, :, :) are the auto-covariance matrices of all types of sensors (e.g., R x (2, 2, :, :) is the auto-covariance matrix of the x-velocity sensors). Similarly, R x (i, j, :, :), for i = j, are the cross-covariance matrices between different types of sensors (e.g., R x (1, 3, :, :) is the cross-covariance matrix between the pressure and y-velocity sensors). iii. The tensor R x depends on (and is determined by) θ even without knowledge of the explicit parametric dependence of {a d } D d=1 in θ, i.e., when treating A as a general unknown complex-valued matrix with M D elements, which are 2M D free parameters (/ degrees of freedom).
A well-known appealing property of tensors is the uniqueness of their CPD. This powerful property holds under relatively mild conditions, as formulated in the following theorem for 4-mode tensors due to Sidiropoulos and Bro [54]: with Q ∈ C I×D , U ∈ C J×D , W ∈ C M ×D and Z ∈ C N ×D , and suppose that X (i, j, m, n) cannot be represented using fewer than D components as in (13). Then, given X , the factor matrices Q, U , W and Z are unique up to permutation and complex scaling of their columns provided that The proof of Theorem 1 is given in [54], Section 3. Going back to our problem, in order to avoid the immaterial permutation and scaling ambiguities, which have no effect whatsoever on the DOAs estimation, we assume w.l.o.g. that the first element of each steering vector is non-negative, i.e., A 1d ∈ R ≥0 for all d ∈ {1, . . . , D}, and that the DOAs are ordered in an ascending order, i.e., θ 1 < . . . < θ D . Furthermore, we now assume that rank(A) = D, and refer to this condition as the sensor array regularity condition. Recall that θ d = θ for all d = ∈ {1, . . . , D} by assumption, hence the sensor array regularity condition is quite mild, and typically holds for any reasonable array configuration / geometry. This regularity condition grants the tensor R x in (12) its uniqueness, as we show in the following theorem.

Theorem 2.
[CPD Uniqueness of R x ] Consider the 4-mode tensor R x as defined in (11). Assume D > 1, and that the sensor array regularity condition holds, i.e., rank(A) = D. Then, the CPD (12) of the rank-D tensor R x is unique.
At this point, based on the uniqueness of the CPD (12) of the covariance tensor R x , we are ready to present the blind CPD-based DOAs estimate, whose intuitive definition comes naturally from Theorem 2. Given the measurements {y[t]} T t=1 , 1) Compute the empirical covariance matrix " R y ; 2) Estimate the noise variance σ 2 v via the MLE for Gaussian signals [60] or the improved (less biased) estimate [61]. Denote it as σ 2 v (note that this is a consistent estimate); 3) Define the covariance matrix estimate of x[t], and the corresponding block partitioning (according to (10)) and using (15), construct the estimated 4-mode covariance tensor " 4) Given the statistic tensor " R x , jointly estimate the DOAs vector θ and the "core" steering vectors matrix A via where with θ ∈ [−π, π) D×1 and ‹ A [ a 1 · · · a D ] ∈ C M ×D , and we use to emphasize that, in general, the arguments may be different from the true unknown estimands θ, A. Obviously, due to Theorem 2, given the true covariance tensor R x , which admits the exact CPD (12), the sources' DOAs are readily extracted from the vectors Therefore, intuitively, if " R y and σ 2 v are "good" estimates, than so does " R x , and it makes sense to define the estimate (17): the best approximate CPD of " R x in the Least Squares (LS) sense yields approximate versions of {c(θ d )} D d=1 , from which the DOAs' estimates arise. Fortunately, this rationale may be rigorously justified, as we prove in the following theorem.

Theorem 3.
[Consistency of the CPD-based estimates θ CPD and " A CPD ] Let " R y and σ 2 v be consistent estimates of R y and σ 2 v , resp. Further, assume that all the elements of A are finite, such that ∃ρ ∈ R + : and that the DOAs vector θ belong to a compact set, such that Then, the estimates θ CPD and " A CPD defined in (17) are consistent, i.e., with any fixed SNR level, for a sample size T → ∞, Proof. See Appendix B.
Notice that the conditions Theorem 3 requires are quite mild, and are fulfilled in practice for the most part. For distributions with finite fourth order moments, " R y and σ 2 v are consistent estimates (where the noise variance is estimated as in [60] or [61]). Furthermore, in practice, the steering vectors {a d } D d=1 are always finite. Lastly, condition (21) is required for technical considerations in the proof presented in Appendix B, but is meaningless from a practical point of view for a sufficiently small . Therefore, we conclude that under these mild conditions, and without knowledge of the explicit dependence of A in θ, (17) are consistent. In particular, we derived θ CPD , consistent blind DOAs estimates, as desired.
We note in passing that our approach also yield, as a by product, the nuisance parameters' estimate " A CPD , allowing for consistent separation of the latent sources s[t], by multiplying Having provided the blind DOAs estimates (17) and the analytical guarantees for their consistency, we now turn to present an iterative algorithm for their actual computation.

IV. COMPUTATION OF THE CPD-BASED ESTIMATES VIA
THE MODIFIED AC-DC ALGORITHM Our goal in this section is to find a computationally feasible algorithm in order to obtain the proposed estimates θ CPD and " A CPD , given the estimated covariance tensor " R x . The roadmap towards this goal is the following. First, we show that the optimization problem (17) is in fact equivalent to an AJD problem, with underlying diagonal matrices which admit a particular parametric structure. Then, we resort to the AJD AC-DC algorithm [55], which due to its alternating mode of operation, allows for a local convenient modification, tailored ad-hoc to the particular parametric structure of the aforementioned diagonal matrices. By this, we obtain an iterative solution algorithm for the optimization problem (17).
For the first step, let us define the LS cost function Now, since we may write where F ij (θ) is elementwise (e.g., F 23 (θ) = 1 2 sin(2θ)) and D ij (θ) ∈ R D×D , using (10)- (11) and (16), observe that Therefore, the optimization problem in (17) is equivalent to AJD in the LS sense. More specifically, since A is not restricted to a particular structure (e.g., orthogonal matrix), (26) accounts for a non-orthogonal AJD problem.

A. Review of the AC-DC Algorithm
One viable solution approach for such a non-orthogonal AJD problem, is using the iterative AC-DC algorithm, proposed by Yeredor [55]. In a nutshell, given a set of K "target-matrices" {Q k ∈ C N ×N } K k=1 , the algorithm seeks a "diagonalizing matrix" B ∈ C N ×L and K associated diagonal is minimized, where {w k ∈ R + } K k=1 are some positive weights. The algorithm alternates between the two following minimization schemes: • The AC ("alternating columns") phase minimizes C AC-DC w.r.t. a single column of B while keeping its other columns, as well as {Λ k } K k=1 , fixed. This phase is sequentially repeated, for all columns of B, for a prespecified number of "sweeps". • The DC ("diagonal centers") phase minimizes C AC-DC w.r.t. the diagonal matrices {Λ k } K k=1 while keeping B fixed. It is readily seen that (26) admits the same formulation as (27), with the simple mapping 6 AC-DC DOAs via CPD 3}. Yet, although it is possible to apply the AC-DC algorithm in its original form in order to solve our specific CPD optimization problem, a significant enhancement can be achieved by exploiting a subtle difference in these two (almost) identical problem, as follows. While the AC-DC attempts to minimize (27) for some general diagonal matrices {Λ k }, namely with KL free parameters (corresponding to 6D parameters according to (28)) in C, in our problem, the diagonal matrices {D ij (θ)} are parametrized by the DOAs vector θ, namely by only D free parameters in [−π, π). Therefore, it would be desirable if certain modifications in the AC-DC algorithm could be made, such that the optimization w.r.t. the diagonal matrices would actually be only w.r.t. θ. The reason for this is twofold. Firstly, the optimization would be for less parameters-D rather than 6D. Secondly, the optimization for each θ d would be confined to the interval [−π, π), rather than searching on the whole complex plane. Consequently, not only the computational cost would be reduced, but the variance in the resulting diagonal elements' estimates would be reduced as well, due to their known parametric structure (24).
Fortunately, the AC-DC algorithm operates iteratively in an alternating manner between two phases: While in the AC phase the diagonalizing matrix B is optimized with the diagonal matrices {Λ k } kept fixed, the optimization in the DC phase is w.r.t. {Λ k } only with B kept fixed. Therefore, we shall now redesign the DC phase so as to adjust it to the specifics of our problem, in order to enjoy the aforementioned advantages.

B. The Modified DC Phase
As mentioned above, in terms of our problem, during this phase the current estimate of A is held fixed. Therefore, in order to emphasize that this is not the final estimate " A CPD , we denote this (intermediate) estimate as " A for brevity. Hence, the problem under consideration in this phase is as follows namely a trigonometrical scalar function of a real-valued scalar argument. As we show in the reminder of this subsection, (31) may be solved efficiently, thus leading to our desired goal in deriving a modified DC phase, tailored specifically to our primary parametric CPD optimization problem. Starting our derivation, in Appendix C we show that differentiating ‹ C LS w.r.t. θ d yields after algebraic simplifications where (3γ + α)τ 4 + 2βτ 3 + 2γτ 2 + (4δ + 2β)τ − (α + γ) = 0. (38) Evidently, (38) is a 4-th order polynomial equation in τ , which has at most four real-valued solutions, and may be solved efficiently using various methods (e.g., [62]). Thus, the global minimizer (31), denoted as θ * , is the solution that minimizes ‹ C LS out of these (maximum four) solutions, and is the updated estimate of θ d . We stress that this approach does not guarantee that the output of the modified AC-DC algorithm would be the global maximizer (17). However, it does guarantee that C LS is non-increasing w.r.t. the iterations, exactly as for the original AC-DC algorithm, thus retaining its original weak convergence property (see [55], Section IV). Summarizing the above, the modified DC phase is given in Algorithm 1.
Note that due to the required matrix multiplications, evaluating ‹ C LS (θ d ) at four points amounts to O M 3 operations. Since the overall computational load per iteration 7 of the AC-DC algorithm, in terms of the left hand side of (28), is O KN 3 (see [55], Section IV), it is O 6M 3 in terms of the parameters of our problem, according to the right hand side of (28). Therefore, we conclude that the overall computational load per iteration of the AC-DC algorithm with the modified DC phase remain O M 3 , namely unchanged.

C. Initialization via Exact Joint Diagonalization
Naturally, whenever an iterative algorithm is proposed the issue of initialization must be addressed. To this end, recall first that (17) is a solution to the equivalent AJD problem (26). Therefore, we propose to use Yeredor's Exact Joint Diagonalization (EJD) for AJD method [63], given in our case by the following steps: Find the two largest eigenvalues and associated eigenvectors p 1 and p 2 of M x ; 4) Construct the matrices P k unsvec (p k ), k = 1, 2; and 5) Compute the eigenvalue decomposition of P 1 P −1 2 , and denote " A EJD ∈ C M ×D as the matrix with the D eigenvectors corresponding to the D largest eigenvalues. For the definitions of the svec(·) and unsvec(·) operators, see [63], Subsection II-B. Having obtained the EJD-based solution " A EJD , we further carry out the following additional steps: , . . . , D}, such that the first element of each steering vector is real-valued, as required; 7) Compute the pair of estimated diagonal matrices 8) Compute θ EJD , the initial EJD-based DOA estimates via (19), where the estimates of c 2 (θ d ), c 3 (θ d ) are the associated diagonal elements of " D 12 (θ), " D 13 (θ), resp. Note that although this non-iterative solution is "merely" an initial solution to the proposed modified AC-DC algorithm, it is instrumental for a successful operation in practice. Indeed, while the estimate (17) is consistent, computing it boils down to finding the global minimizer of C LS , which, in general, is not a trivial task. However, and as we show via simulations in Section VII, when the proposed initialization is used, the iterative algorithm yields the desired estimates.
As an intermediate summary, by establishing Algorithm 1, we have thus presented a consistent blind DOAs estimate (17), along with an iterative solution algorithm-the AC-DC algorithm with the modified DC phase-equipped with a noniterative "educated" initial solution, tailored ad-hod to our specific blind AVS DOA estimation problem.
We now turn to the particular case of Gaussian signals, in which significant performance enhancement can be attained by further refining the CPD-based estimates θ CPD , " A CPD from (17). Moreover, and quite interestingly, we show that the results obtained within the Gaussian signal model framework are valid for other signal models as well, enabling this performance enhancement for a wider class of signals' distributions.

V. OPTIMAL ESTIMATION FOR GAUSSIAN SIGNALS
In this section only, we further assume that s[t] is circular CN. As a consequence, it follows that Since the model is now fully specified, it is first instructive to study the CRLB on the MSE of any unbiased estimate in joint estimation of all the unknown deterministic parameters, namely A, θ and σ 2 v (recall that R s = I D by assumption).

A. Cramér-Rao Lower Bound in Blind DOA Estimation
For brevity, we define the vector of all real-valued unknowns where I M [0 M −1 e 1 · · · e M ] ∈ R (M −1)×M , when multiplies from the left, preserves all the rows except for the first one of a matrix with M rows (recall A 1d ∈ R ≥0 ), and K ϕ 2M D + 1 is the total number of unknown parameters.
Since {y[t]} T t=1 are all CN and i.i.d., the FIM elements corresponding to A, θ and σ 2 v are given by 8 (see, e.g., [64]) where J (ϕ) denotes the FIM. In Appendix D we show that 8 We specifically use a different notation for the FIM's elements, with slight abuse in notation also in (47), which is more natural w.r.t. their definition, and therefore easier to comprehend in this context, for the sake of clarity.
In addition, by the Woodbury matrix identity [65], we have (46) Therefore, all the required expressions for the computation of the FIM are at hand. For example, by (44) and (45), we have The CRLB on the MSE in unbiased joint estimation of A, θ and σ 2 v is readily given by the inverse of the FIM J (ϕ), whose elements are prescribed in (41), using (42)- (46).
Next, we further utilize this result in order to derive an approximate iterative solution algorithm for the computation of the MLE

B. MLE Computation via the Fisher Scoring Algorithm
Given an initial estimate of ϕ, the FSA can be used in order to obtain a stationary point of the log-likelihood (if it converges). Moreover, if this initial estimate is "close" enough to the global maximizer of the log-likelihood, the algorithm is likely to converge to the MLE. Formally, the update equation of the FSA for the n-th iteration is given by where ϕ (n) is the estimate of ϕ in the n-th iteration, and is the log-likelihood function, where c is a constant independent of ϕ. Hence, in order to carry out the iterations (49) such that they will successfully converge to the MLE ϕ ML , three ingredients are required: a sufficiently "good" initial solution ϕ (0) , and closed-form expressions of the FIM and the score function, i.e., J (ϕ) and ∇ ϕ L, resp. Now, recall that " A CPD , θ CPD and σ 2 v , specified in (17) and (15), resp., are consistent estimates. Therefore, when rearranged in vector form according to (40), can serve as a "good" initial estimate of ϕ, which is presumably "close" to the global maximizer of (50). Furthermore, note that we have already obtained closed-form expressions for the elements of the FIM J (ϕ), given in (41) while using (42)- (46), which can be computed for any ϕ. Moreover, using the chain rule, we have for all i ∈ {1, . . . , K ϕ } where ∇ Ry L ∈ C 3M ×3M and ∇ ϕi R y ∈ C 3M ×3M . As we show in Appendix D, and with the already obtained expressions (42)-(45), we have obtained a closed-form expression of the score (52) w.r.t. any element of ϕ, such that ∇ ϕ L can be computed for any ϕ. By this, we now have the three required ingredients for a successful operation of the FSA for the computation of the MLE ϕ ML . Note further that as long as M ·D is not too "large" (in terms of matrix inversion), the computation of J −1 Ä ϕ (n) ä is not very costly w.r.t. computational load.

VI. PHASE 2: THE KLD-BASED BLIND DOAS ESTIMATE
Our goal in this section is to show that (49) can successfully operate not only for Gaussian signals. Moreover, we would like to achieve analytical arguments which explain and justify two aspects. First, the rationale of invoking results obtained under the Gaussian model for non-Gaussian signals. And second, the resulting higher accuracy attained with this approach for our primary goal, DOAs estimation. To this end, we begin by showing that, in general, it is theoretically possible to attain better estimates, in terms of lower MSEs, than the CPD-based estimates (17) with only mild assumptions on the signal model.

A. Suboptimality of the Equally-Weighted LS Estimates
In our general framework, as long as the fourth-order joint cumulants of the measurements {y[t]} T t=1 are finite, i.e., ∃ ∈ R + : ∀i, j, k, ∈ {1, . . . , M } : the estimate " R y , which is not necessarily the MLE of R y , is still consistent by virtue of the law of large numbers [66]. Moreover, denoting r y vec * (R y ) ∈ C Kr×1 with K r 3M (3M + 1)/2, by virtue of the central limit theorem 9 [66], where Γ ε ∈ C Kr×Kr and C ε ∈ C Kr×Kr are the covariance and pseudo-covariance matrices of r y , resp. Note that r y , Γ ε and C ε all depend on the unknown parameters vector ϕ. Now, since the asymptotic distribution of r y , or, equivalently, " R y , is given by (55), as we show in Appendix E, ML estimation of ϕ based on " R y asymptotically amounts to where ε r y − r y denotes the vector of estimation errors E " R y − R y in estimation of R y , i.e., ε = vec * (E), and Evidently, as seen from (56), asymptotically optimal estimation of the unknown parameters based only on " R y is via the Optimally-Weighted Non-Linear LS (OWNLLS) criterion. Note that the optimal weight matrix R −1 ε is a particular case of a general (not necessarily optimal) positive-definite weight matrix, denoted, say, by W . Indeed, with equal weights, namely W = I 2Kr , the criterion yields the CPD-based estimates 9 Note that the diagonal elements of Ry are in fact real-valued, and are therefore normal, rather than CN. However, this does not affect the following derivation, so we allow this slight abuse in notation. (17). However, as we show in Appendix E, the optimal weight matrix R −1 ε is generally not equal to a scaled identity matrix. Therefore, in general, the equally-weighted LS criterion (17) yields suboptimal estimates. Thus, we conclude that even for the general (non-Gaussian) signal model, pursuing improved, more accurate estimates than θ CPD and " A CPD is not in vain. We stress that the estimate (56) is generally not the MLE w.r.t. the raw data {y[t]} T t=1 , since " R y is not necessarily a sufficient statistic. Nevertheless, despite the possible strict statistical insufficiency of " R y , it still encapsulates valuable information which is not fully used via equally-weighted LS fitting (17), and can be further exploited, as we show next.

B. Enhancement via the KLD Covariance Fitting Criterion
Resorting temporarily to the framework of Gaussian signals, as explained in Subsection V-B, the FSA (49) asymptotically yield the MLE 10 ϕ ML . As can be easily seen from the loglikelihood (50), the sufficient statistic in the Gaussian model is the sample covariance matrix " R y . Therefore, ML estimation of ϕ based only on " R y yields the MLE of ϕ based on the raw data {y[t]} T t=1 . Moreover, as shown in Appendix E, ML estimation of ϕ based on " R y asymptotically amount to OWNLLS estimation of ϕ based on " R y (56). Consequently, it follows that for the CN signal model (39), asymptotically, Accordingly, this implies that the iterations (49), namely the FSA for the Gaussian model, are (approximately) an implicit computation of the OWNLLS estimate (56). In other words, maximizing the log-likelihood (50) is asymptotically equivalent to minimizing the OWNLLS objective in (56). Hence (49) is an iterative procedure for a solution of either of these equivalent objectives. Furthermore, note that in this model, the optimal weight matrix R −1 ε is computed under the Gaussian model (39) as well, and is thus fully specified by the covariances and pseudo-covariances (for all i, j, k, ∈ {1, . . . , M }) For the detailed derivation of (59), see Appendix E. At this point, observe that the log-likelihood (50) reads where D KL (·, ·) denotes the KLD [67], and c is a constant independent of ϕ. Clearly, it follows from (60) that maximizing the log-likelihood L(ϕ) is equivalent to minimizing the KLD R y , R y ä for shorthand. However, since the CN distribution is fully characterized by Second-Order Statistics (SOSs), note that the KLD of two zero-mean multivariate CN distributions serves in itself as a plausible criterion for consistent covariance matrix estimation. Indeed, from Gibbs' inequality, the KLD is always non-negative (e.g., [67]), and regardless of the true underlying distributions governing " R y and R y . Thus, asymptotically, we now have where the optimal weight matrix R −1 ε in (62), now computed for a general (not necessarily Gaussian) signal model, is fully specified by (see Appendix E) for all i, j, k, ∈ {1, . . . , M }. Notice that whenever κ y [i, j, , k] = 0, e.g., for CN signals, (63)-(64) coincide with (59). Otherwise, (59), which depend only on the elements of R y , can be regarded as a SOS-based approximation for (63)- (64). Therefore, it follows from (62) that in the general case, namely regardless of the distribution of y[t] (under the mild condition (54)), the iterations (49) serve as a quasi-Newton algorithm (e.g., [68]) for an implicit computation of our proposed KLD-based estimate, which is also the SOSbased approximate OWNLLS estimate. That is, the weight matrix R −1 ε in (62) for this approximate OWNLLS estimate is specified in (59), which depend on the SOS R y , and T . Of course, (62) also implies that ϕ KLD is a consistent estimate.
Summarizing the above, our proposed blind DOAs estimate θ KLD is extracted from ϕ KLD (as its D penultimate entries, according to (40)), which is computed via the iterations (49)the FSA for CN signals. With this, the instrumental role of the CPD-based estimate is now revealed. As seen from (62), our desired KLD-based DOAs estimate is the solution of a nonconvex high-dimensional optimization problem. In contrast, although the CPD-based estimate is statistically inferior, it is nevertheless consistent, and can be obtained (relatively) efficiently via AJD. Thus, the KLD-based estimate and its iterative solution (49) are of practical value only when provided with a sufficiently "good" initial estimate, which lies in the basin of attraction of the global minimizer of D CN the key role of θ CPD , " A CPD in our overall proposed solution.

VII. SIMULATION RESULTS
In this section, we consider three simulation experiments of blind DOAs estimation for different scenarios in order to corroborate our analytical derivations by empirical results. First, for a ULA, we demonstrate the asymptotic optimality of the proposed estimate for Gaussian signals. Second, we demonstrate its consistency for non-Gaussian signals, this time with an uncalibrated ULA, exemplifying the robustness of our blind approach w.r.t. the signals' distributions. Lastly, we demonstrate the robustness of our blind estimate w.r.t. the array geometry and reliability by considering a Uniform Circular Array (UCA) with faulty sensors. In all three simulation experiments we also compare our proposed method to Zhang et al.'s TriD blind estimate [1], thus demonstrating the considerable accuracy enhancement w.r.t. a state-of-the-art competitor.
Throughout, we consider model (1) with unit variance sources, and our proposed estimates were obtained as follows: A CPD via the modified AC-DC algorithm (Section IV), with the initial estimates θ EJD , " A EJD (Subsection IV-C); and 3) Phase 2: Compute θ KLD , " A KLD via (49) (the FSA for CN signals), with the initial estimates θ CPD , " A CPD .
All empirical results were obtained by 10 4 independent trials.
In the first two experiments which follow (Subsections VII-A and VII-B), we consider a ULA with M = 7 AVSs and half wavelength inter-element spacing (i.e., λ/2), and D = 3 sources impinging from azimuth angles θ = [−56

A. Optimal Performance for Gaussian Signals
In this experiment CN sources are considered. Therefore, the received signals' distribution is prescribed by (39), and it follows that θ KLD = θ ML (and also "

B. QPSK Sources, Laplace Noise and an Uncalibrated Array
In the second experiment we consider the same ULA as in the first one, but now with inter-AVSs' positioning errors and gain offsets [39]. Formally, each element of the uncalibrated steering vectors matrix A now reads where g m ∈ R + and ∆ (m) x , ∆ (m) y ∈ R are the gain offset and the xand y-axis positioning errors of the m-th AVS, resp., k 2π/λ is the sources' wavenumber, and A md exp (π(m − 1) cos(θ d )). In our simulation, the gain offsets {g m } and the positioning errors, {∆ (m) x } and {∆ (m) y }, were independently drawn (once, and then fixed) from the uniform distributions U (0.7, 1.3) and U (−1, 1), resp. In order to demonstrate the robustness of our method w.r.t. the signals' distributions, we consider Quadrature Phase-Shift Keying (QPSK) sources. Moreover, here the real and imaginary parts of the noise v[t] are independent Laplace distributed, intentionally not in accordance with our CN assumption, in order to demonstrate our proposed estimate's robustness to the noise distribution. Fig. 2 presents the DOAs RMSE vs. the sample size T for a fixed SNR level of 5[dB], and implicitly demonstrate that (62) holds, and hence the preceding    derivation, particularly (58). Indeed, here as well both blind estimates are robust w.r.t. the gain and phase offsets, though similarly to the previous scenario examined in Subsection VII-A, a considerable improvement is still obtained relative to the TriD method up to ∼ 7[dB]. Additionally, in this experiment only, we also present the performance of the (phase 1) CPD-based estimates, which play the key role of initial estimates in the computation of the proposed enhanced KLDbased estimates. Similar trends as in Fig. 1b were obtained in simulations for different SNR levels in this scenario as well.

C. Gaussian Mixtures and a UCA with Faulty Elements
In our last experiment, we consider a UCA with M = 5 equiangular spaced AVSs placed on its circumference [69]. In addition to changing the array geometry, we further investigate the performance of the proposed method in the presence of sensors failures (e.g., [70]). Specifically, we assume that the 2 nd and 4 th AVSs are faulty, such that all their (three) elements do not receive any external signal, and thus contain only noise. Assuming for simplicity that the radius of the UCA is λ/2, , ∀d ∈ {1, 2}. As evident from Figs. 3a and 3b, presenting the DOAs RMSEs vs. the sample size (for fixed SNR of 0[dB]) and the SNR (for fixed T = 500), resp., similar trends are obtained as in the previous two experiments. This not only demonstrates the robustness of our proposed approach w.r.t. the array geometry, sensors malfunctioning, and the sources' distribution, but also indicates that the substantial accuracy improvement, which in this case reaches up to more than an order of magnitude, is obtained in various different scenarios.
Lastly, we consider the Interference-to-Source Ratio (ISR), which measures the residual energy of the j-th source in the reconstruction of the i-th source, and is a common measure for the separation performance. Fig. 4, presenting the ISR vs. the sample size, shows that our estimate for the mixing matrix A also yields better source separation than the TriD estimate.

VIII. CONCLUSION
In the context of passive AVS arrays, we presented a novel blind DOAs estimate. Rather than estimating the DOAs directly from the raw data, our estimate, obtained in a two-phase procedure, exploits the special structure of the observations' covariance matrix. In the first phase, the CPD-based estimates are obtained via the modified AC-DC algorithm, based on the the unique quadrilinear decomposition of the SOSs covariance tensor. In the second phase, the proposed estimates, which were shown to be asymptotically equivalent the OWNLLS estimates, are obtained via KLD covariance fitting. Since minimization of the KLD is equivalent to maximization of the Gaussian ML objective, our proposed KLD-based estimates are optimal for Gaussian signals, and can be computed via the FSA for non-Gaussian signals as well. Our analytical results were supported by simulation experiments in various scenarios, which demonstrated our estimate's robustness, as well as its superiority over another leading blind DOAs estimate.

APPENDIX A PROOF OF THEOREM 2
Proof. In order to prove Theorem 2, we shall invoke Theorem 1. To this end, note that R x , given explicitly in (12), may also be written in the form (13) of Theorem 1 as with C(θ), C(θ), A and A * as the four factor matrices of R x . Then, according to Theorem 1, the factor matrices, which are in this case C(θ) and A, are unique up to permutation and complex scaling of columns provided that From the sensor array regularity condition, we have rank(A) = D, which also implies rank(A * ) = rank(A) = D.
Since A (and therefore A * ) is full column rank, it follows that k A = k A * = D. In addition, since θ d = θ for all d = by assumption, it follows that C(θ) is also full column rank, which means that k C(θ) = rank(C(θ)) = min{3, D}. Therefore, in our case, the uniqueness condition (67) reads For D ≥ 3, (68) obviously hold. For the complementary case D < 3, (68) becomes 2D ≥ 3. Therefore, we conclude that for D > 1, which holds by assumption, the CPD (66) is unique up to permutation and complex scaling of columns. Finally, since the first row of C(θ) is the all-ones D-dimensional vector and {A 1d ∈ R ≥0 }, thus eliminating the complex scaling ambiguity, and since θ 1 < . . . < θ D , thus eliminating the permutation ambiguity, we conclude that (66) is unique.

APPENDIX B PROOF OF THEOREM 3
Proof. To prove Theorem 3, we shall use the basic consistency theorem for extremum estimators [71]. To this end, for shorthand, let us first define the vector of real-valued unknowns ×M preserves all the rows except for the first one of the matrix it is left-multiplied with, and further define " . Then, by virtue of Theorem 2.1 in [71] due to Newey and McFadden, if there exist a function Q 0 ( ϑ) such that: (i) Q 0 ( ϑ) is uniquely maximized at ϑ; (ii) ϑ ∈ Θ, where Θ is a compact set; (iii) Q 0 ( ϑ) is continuous; and (iv) " Q T ( ϑ) converges uniformly in probability to Q 0 ( ϑ), then namely (17) are consistent estimates.
We shall now show that the four conditions above (i)- First, condition (i) is an immediate consequence of Theorem 2, the uniqueness of the CPD of R x , proved in Appendix A. Second, conditions (20) and (21) stated in the Theorem imply that both θ and ‹ A belong to compact sets. The union of these compact sets, denoted Θ henceforth, is a compact set as well, hence ϑ belong to the compact set Θ, and condition (ii) is satisfied. Third, notice that R Ä θ, ‹ A ä is a continuous function w.r.t. θ and ‹ A by its definition (18), in the sense that each element of R Ä θ, ‹ A ä is a continuous function w.r.t.
each of the elements of θ and of ‹ A. Additionally, note that , the function Q 0 ( ϑ) is a composition of continuous functions, and is therefore continuous, thus (iii) holds. To show that condition (iv) holds, i.e., uniform convergence in probability [72], we first present the notion of stochastic equicontinuity ([72], Section 2): " Q T (ϑ) is said to be stochastically equicontinuous if for every , η > 0 there exist a sequence of random variables " ∆ T ( , η) and a samples size T 0 ( , η) such that for T ≥ T 0 ( , η), Pr " ∆ T ( , η) > < η, and for each ϑ there is an open set U ϑ containing ϑ with Now, according to Lemma 2.8 in [71], if ϑ belong to a compact set and Q 0 ( ϑ) is continuous, then condition (iv) holds if and only if As we have already shown above, ϑ belong to the compact set Θ and Q 0 ( ϑ) is continuous. Thus, if (a) and (b) hold, condition (iv) is fulfilled, and the proof is completed. By the assumption in Theorem 3, " R y and σ 2 v are consistent. Hence " R x defined in (15) is consistent, and also " R x . Further, note that " Q T ( ϑ) is continuous w.r.t. " R x . Therefore, by the Continuous Mapping Theorem (CMT) [73], for all ϑ ∈ Θ, thus (a) holds. We now turn to the final phase, proving (b).

APPENDIX C DERIVATION OF THE DOAS ALTERNATING LS EQUATIONS
Let us begin by rewriting the cost function (23) as where R x depends on θ and A as prescribed in (8). Now, differentiating (73) w.r.t. θ d gives where have have used R x = R † x as well as " Starting with the first term in (74), substituting R x with (8), where we have defined cos(2θ k ).
(76) Moving to the second term in (74), observe first that (78) Substituting (78) into (77), and using (15), we now have where we have defined Finally, substituting (75) and (79) into (74), we obtain the LS equations for θ d , For the alternating LS equations, solving for θ d while fixing all the other parameters (as in (30)), we substitute " A and { θ k } k =d with A and {θ k } k =d , resp., everywhere in (82).

APPENDIX D COMPUTATION OF THE SCORE AND THE CRLB
In order to obtain closed-form expressions of the score w.r.t. each element of the vector of unknowns ϕ, using the chain rule, we may alternatively use (52). Hence, our goal now is to compute the gradient of R y w.r.t. each element, and the gradient the log-likelihood L(ϕ) w.r.t. R y .
Starting with the gradient of R y w.r.t. the real and imaginary parts of {A md }, using (3) and (8) we have for all m ∈ {1, . . . , M }, and in the same fashion Lastly, the gradient of R y w.r.t. σ 2 v reads Having obtained ∇ ϕi R y for all i ∈ {1, . . . , K ϕ }, namely the gradient w.r.t. all the unknowns, we move proceed to the gradient of log-likelihood w.r.t. R y . Recall from (50) that where c is independent of ϕ and R y . Hence, based on wellknown matrix functions derivatives (e.g., [74]), using we obtain Substituting (86)-(89) and (92) into (52), we obtain closedform expressions for the score. Further, substituting (86)-(89) into (41), we obtain closed-form expressions for the FIM elements, which, upon inversion, gives the CRLB.

APPENDIX E ASYMPTOTIC ML ESTIMATION BASED ON "
R y As mentioned in Subsection VI-A, by virtue of the central limit theorem, Our goal here is to show that ML estimation of ϕ based only on r y , asymptotically amounts to OWNLLS (56). To this end, we first derive closed-form expressions of the elements of the covariance and pseudo-covariance matrices Γ ε and C ε , resp. By definition, Γ ε , C ε are also the covariance and pseudocovariance matrices of the vector of errors ε = vec (E) in estimation of R y . Therefore, we may equivalently compute the covariance and pseudo-covariance of E, which explicitly determine Γ ε and C ε by the injective vec * (·) mapping. Thus, Using the fact that y[t] is proper, we have hence we may write the summand in (95) as for all i, j, k, ∈ {1, . . . , M }. Note that in the particular case of CN sources, it follows that y[t] is CN, and therefore κ y [i, j, , k] vanishes for all i, j, k, ∈ {1, . . . , M }, which gives (59). Thus, (98) and (99) are the closed-form expressions for all the associated elements (in compliance with the vec * (·) mapping) of Γ ε and C ε , resp., as required. Next, we turn to the asymptotic ML estimation of ϕ based on " R y , whose asymptotic distribution is prescribed in (94). First, notice that, according to (98)-(99), all the covariances and pseudo-covariances are a multiplication of 1 T by a factor independent of T . Accordingly, we define such that Γ ε and ‹ C ε depend only on the elements of R y , and more specifically, are independent of T . With these notations, we further define and the auxiliary matrices (to be used shortly) which implies that ‹ R ε and ‹ P ε are also independent of T . Now, using the asymptotic distribution (94), the probability density function of r y reads pr y ( r y ; ϕ) Thus, by definition, the MLE ϕ ML of ϕ based on r y , is asymptotically given by which is the OWNLLS estimate, and we have used the fact that log is an increasing monotonic function in (i), omitted irrelevant constants w.r.t. ϕ in (ii), as well as in (iii) due to and the approximation (iv) holds for a sufficiently large sam-