A Degenerate Spatial ARMA Process of External Sources Impinging Upon an Azimuth-Only ULA and the Estimation of Doas and Noise Variances

Array processing is an interdisciplinary topic of both physics and signal processing. Physical basis of array processing is the orderly-setup of sensors in space that induces regular change of the phase of external sources incident upon sensor array. In the current paper the direction-of-arrival (DoA) estimation of an azimuth-only uniform linear array (ULA) is investigated in theory and method from the fused perspective of signal and physical properties of the incident sources. By organically fusing the stationary assumption on source signals and noises with difference operation on the phases, a distinctive system of linear equations satisfied by the incident sources is theoretically derived for the azimuth-only ULA with the Hankel-block-matrix of signal correlations as coefficient matrix and the elementary power-sum symmetric functions of the propagators of incident sources as the unknowns. Based on the derived system of linear equations, signal model of the incident sources is first proved as a degenerate spatial ARMA process subject to the identical autoregressive and moving average parameters and simultaneously obeying the dimensional homogeneity principle (DHP) in physics. The explicit root-finding polynomial is proposed with the unknowns of the system of linear equations as polynomial coefficients and the propagators as the roots. No extraneous roots and conjugate symmetry constraint on polynomial coefficients are involved. The DoAs and noise variances can be separately estimated under the backgrounds of spatially white and colored noises, which are numerically analyzed with the different coherent lengths of noises. A simple sound experiment is designed and performed to verify the proposed DoA estimation method. It is promising to investigate the DoA estimation of the ULA model, particularly of the ULAs of the multi-dimensional array from the fused perspective of physical and signal properties of the incident sources.

calibration. The phase is proportional to the corresponding phase-difference. In contrast, for the multi-dimensional (M-D) array or near-field source, function form of the phase becomes highly complicated due to the introduction of azimuth, elevation angles and propagation distance.
Meanwhile, array processing is signal model-based, which is reflected by the statistical property of source signals and sensor noises. In practice, the well-known concepts such as non-stationary/stationary, deterministic/stachastic, independent/correlated/coherent, narrow-band and wide-band and so forth have been frequently introduced to reflect signal properties of source and noise and impacted on array processing including the DoA estimation in theory and method.
Physical and signal properties have formed a pair to fully describe the incident sources. In this sense, array processing can be taken as a problem of both physics and signal processing. It is beneficial and meaningful to study array processing from the fused perspective of physical and signal properties of the incident sources.
In past, various popular DoA estimation methods had been proposed from signal perspective of the incident sources and noises and mainly classified into beamforming technique [1]- [3], linear prediction (LP) [4]- [6], the subspace-based methods such as the conventional MUSIC (Multiple signal classification) and ESPRIT (Estimation of signal parameters via rotational invariance techniques) [7]- [12], Maximum likelihood (ML) [13] and their variants and combinations. The reviews could be found in [14]- [16].
It was noted that a well-known concept was the so-called subspace in the DoA estimation, in particularly the signal subspace that stemmed from the partition of data covariance matrix. Signal subspace was proved to be spanned both by the implicit steering vectors of the ULA representation and by signal eigenvectors after eigenvalue decomposition (EVD) on data covariance matrix, which showed signal subspace in effect had played a role of bridge connecting steering vector and the ULA data. After considering the one-to-one map from steering vector to azimuth angle, azimuth angles could be estimated from the ULA data via signal subspace in theory. The DoA estimation of the azimuth-only ULA had been elegantly accomplished by the MUSIC and ESPRIT methods, which also resulted in the popularity of the subspace-based methods in array processing.
In the recent years, the sparsity of the incident source signals in the spatial domain has drawn growing attention, spare signal processing techniques based on CS (Compressed sensing) theory [17] have been extensively introduced and exploited for the DoA estimation of array processing. Under the CS frame, the DoA estimation had been formulated as a sparse signal representation problem after discretizing angular area and formulating an overcomplete dictionary and then performed by the convex optimization problems of the l 1 -norm penalty [18] or l p -norm penalty with p ≤ 1 [19], [20] and the sparse Bayesian Learning (SBL) [21]- [23]. The CS-based DoA estimation methods exhibited some advantages including the increased resolution, improved robustness to noise and limitations in data quantity, as well as not requiring an accurate initialization. On the other hand, the accuracy and performance of the above methods relied on the exact-alignment requirement of the true DoAs of the incident sources with the discretized grids and the tuning of regularization parameter balancing the sparsity of the spectrum and the residual norm, which have been the focuses in study and also resulted in the promising off-grid models based on the Taylor series expansion method and the linear interpolation method.
Physical properties of the incident sources induced by the orderly-setup of sensors has been more-or-less exploited in some estimation methods. The noted one was the rotationalinvariance structure of array manifold developed in the ESPRIT. The propagators of all the sources were wholly extracted out from array manifolds of the displaced or partitioned subarrays and explicitly formed a diagonal unitary matrix. Owing to avoiding array calibration and reducing computational complexity, the ESPRIT had been the other representative subspace method, the successive methods such as MPM (Matrix pencil method) [24], [25], PM (Propagator method) [26]- [28], and SUMWE (The subspace-based method without eigen-decomposition) method [29] as well as CODE (The cross-correlation based 2-D DOA estimation) method [30] had borrowed ESPRIT idea. The rotational transformation was a whole operation acting on the displaced or partitioned subarrays. Analogous to the MUSIC, the estimation of propagators still required to make the connection with the ULA data, which was completed by exploiting the rotational invariance of the underlying signal subspace.
Array processing is an interdisciplinary problem of both physics and signal processing. As a physical parameter containing the unknown DoAs, phase-difference can be treated as a result of difference operation on phases of the same source impinging upon the adjacent sensors. Once such a connection is made, except for the rotational-invariance structure of array manifold, by making use of difference operation on phases, phase-difference also can be individually extracted out from the ULA representation and explicitly expressed for the subsequent estimation. By fusing difference operation on phases with the appropriate signal assumption on source and noise, here we take the azimuth-only ULA as an object and probe the DoA estimation in theory and method.
Indeed, the more deep and significant insight into the DoA estimation benefits from the fused perspective. It is from the fused perspective that the MYW system of linear equations and the root-finding polynomial are first derived for the azimuth-only ULA. Signal model of the incident sources is revealed as a degenerate spatial ARMA process (Autoregressive moving-average process), neither AR (Autoregressive) nor MA (Moving-average) models. Furthermore, the ARMA process is not arbitrary, but explicit in function forms of model coefficients. The degenerate ARMA process constrains the DoA estimation to possess its own distinctive method.
By contrast, in past signal model had been uncertain in principle, which led various AR processes or LP equations to be introduced into the DoA estimation [4]- [6], [8].
The paper is structured as follows: The extraction manners of phase-difference are described in Section II. The systems of linear equations and the root-finding polynomial are derived for the azimuth-only ULA, respectively, with the forward and backward difference operations (FDO and BDO) in Section III. A degenerate spatial ARMA process of the incident sources is developed and further analyzed from the standpoints of physics and signal processing in Section IV. The properties of the MYW system of linear equations are discussed in Section V. The implementation of the DoA estimation is described in Section VI. Numerical analysis and experimental verification are performed in Section VII and VIII.

A. THE AZIMUTH-ONLY ULA MODEL
Consider a typical one-dimensional azimuth-only ULA comprising the M identical and omnidirectional sensors, the K far-field narrow-band sources {s k (t)} K k=1 with azimuth angles {θ k } K k=1 are impinging upon the array. The azimuth angle θ k is defined as the one between the incident direction of the k-th source and the normal of the ULA. The representation of the azimuth-only ULA is expressed in matrix notation as [14] where the column vectors of sensor data, source signals and sensor noises are X(t) = [x 1 (t), x 2 (t), · · · , x M (t)] T , S(t) = [s 1 (t), s 2 (t), · · · , s K (t)] T and W (t) = [w 1 (t), w 2 (t), · · · , w M (t)] T with t as time variable, (·) T denotes transpose operator. The statistical property of source signal is reflected via s k (t). A(θ ) stands for array manifold or array steering matrix with θ = [θ 1 , θ 2 , · · · , θ K ] T . Array manifold A(θ ) of the azimuth-only ULA is a Vandermonde matrix given by where steering vector a(θ k ) = [e jϕ 1 (θ k ) , · · · , e jϕ m (θ k ) , · · · , e jϕ M (θ k ) ] T with j as imaginary unit. Steering vector a(θ k ) of array manifold A(θ ) in (2) contains the two physical parameters, one is the phase ϕ(θ ), the other is phase-difference ϕ(θ ). Suppose ϕ m (θ k ) in a(θ k ) is the phase of the k-th source impinging upon the m-th sensor after the first sensor is calibrated as the zero-phase sensor, then where ϕ(θ k ) is the corresponding k-th phase-difference.
For the far-field narrow-band source, the k-th phasedifference ϕ(θ k ) is simplified as In (4) the k-th phase-difference ϕ(θ k ) is a function of azimuth angle θ k , sensor interspacing d, carrier frequency f and propagation velocity c. The corresponding propagator v k is defined as It is shown in (3) that the phase ϕ m (θ ) is the product of phase-difference ϕ(θ ) and the index (m − 1) after sensor calibration (To simplify the notation, here we omit the index k of ϕ m (θ k ) and ϕ(θ k ); we will reinstate the index k later on, when needed.). For the azimuth-only ULA, only one phasedifference ϕ(θ ) appearing in the phase ϕ m (θ) guarantees the one-to-one map from the propagator v to azimuth angle θ, while the index (m − 1) in (3) describes sensor position and is linearly-varied. The function form of the phase ϕ m (θ ) is induced by the artificial and orderly setup of sensors.
In contrast to the phase ϕ m (θ ), owing to independent of position index of sensor and also containing the DoAs, phasedifference ϕ(θ ) is the more essential parameter in physics. The DoA estimation of the azimuth-only ULA in effect has been viewed as the estimation of propagator or phasedifference from the noise-contaminated ULA data.

B. THE EXTRACTION OF PHASE-DIFFERENCE WITH DIFFERENCE OPERATION
Except for the rotational-invariance structure of array manifold proposed by the ESPRIT which was used to wholly extract the propagators of the incident sources, in accordance with the above definition, here we provide the other manner of explicit extraction of phase-difference.
The proportional relation between the phase and phasedifference in (3) yields and also For comparison, with positive phase ϕ m (θ ), phasedifference ϕ(θ) is extracted by using the forward difference operation (FDO) in (6), while with negative phase −ϕ m (θ ), ϕ(θ) obtained by using the backward difference operation (BDO) in (7). The BDO can be implemented after conjugate preprocessing of the ULA data, which will be discussed next.
Here we probe the DoA estimation of the azimuth-only ULA in theory and method via implementing difference operation on phases and explicitly extracting phase-difference. Apparently the implementation of difference operation requires the appropriate condition on source signals and noises since the phases {ϕ m (θ )} M m=1 and the signals {s k (t)} K k=1 and the noises {w m (t)} M m=1 coexist into the representation of the azimuth-only ULA in (1).

A. THE CONSTRUCTION OF CORRELATION SEQUENCE OF THE ULA DATA
In order to study the DoA estimation via performing the FDO, a raw sequence S f 0 composed of correlation functions VOLUME 8, 2020 of sensor data {x m (t)} M m=1 is first constructed and given by S where the m-th element R x m x n = E{x m (t)x H n (t)}. E{·} and (·) H denote the expectation operation and conjugate transpose, respectively. For the sake of convenience, the data x n (t) used in R x m x n is called reference signal and can be chosen from sensor data as needed.
After the first sensor is calibrated as the zero-phase sensor, the m-th element R x m x n is expanded as where R s k x n = E{s k (t)x H n (t)} denotes the correlation function between the k-th source signal s k (t) and x n (t), R w m x n = E{w m (t)x H n (t)} is the one between the m-th noise w m (t) and x n (t).
The property of R s k x n and R w m x n in R x m x n are usually unknown. However, under the stationary assumption on source signals and noises, R s k x n and R w m x n become the timeindependent constants.

For the first element
Similarly, for the second element R x 2 x n , Dividing the right-side and left-side hands of (11) by the corresponding sides of (10), we preform the FDO on the phases of s 1 (t) incident on the first and second sensors. To guarantee the implementation of the FDO, the involved correlation function R s 1 x n is required to be a time-independent constant, which implies that source signals and noises are stationary. After the FDO, we get R 1 x 1 x n = R 1 s 2 x n + · · · + R 1 s k x n + · · · + R 1 It is shown in (12) that the unknown correlation function R s 1 x n appearing in (10) and (11) is eliminated by the FDO under the stationary assumption on source signals and noises.
For the m-th and (m + 1)-th sensors, it follows (12) that For the convenience next, R 1 x m x n (m = 1, · · · , M −1) in (13) is named as the compound correlation function in respect to the index m after the first FDO.
After the first FDO, it is noted that the original sequence S f 0 in (8) changes as S x m x n in (13) remains the similar function structure as R x m x n in (9), which shows that the FDO can be performed on the phases of R 1 s 2 x n in R 1 x m−1 x n and R 1 x m x n . After the k-th FDO, by derivation, the correlation functions R s 1 x n , R s 2 x n , · · · , R s k x n are shown to be orderly eliminated under the stationary assumption on source signals and sensor noises, the compound correlation function R k x m x n is derived as where the involved recurrence formulas satisfy The sequence after the k-th FDO is given by , which similarly shows that the FDO can continue to be performed for the elements in S f k . After the K -th FDO, all the K correlation functions {R s k x n } K k=1 are eliminated in derivation, the compound correlation function R K x m x n in respect to the index m contains no signal correlations {R s k x n } K k=1 and satisfies The sequence after the K -th FDO is Employing recurrence formulas in (15), (16), and (17), Equation (18) in respect to the index m is expanded as where the unknown coefficients {b k } K k=1 denote the elementary power-sum symmetric functions of the propagators and have the following distinctive forms Equation (19) is a difference equation satisfied by the propagators of the K sources impinging upon the azimuthonly ULA under the stationary assumption on source signals and noises.
Writing (19) for the indices m = 1, · · · , M − K gives the following system of linear equations in matrix notation where H f , b and r f denote coefficient matrix, the vectors of the unknowns and constant terms and are given by where In order to solve the K unknowns {b k } K k=1 , in accordance with the concrete form of (22), the relation between M and K is required to satisfy K ≤ M 2 . As K = M 2 , the system of the K equations in K unknowns is obtained and all the sensors of the ULA are employed for estimation.
Although the azimuth-only ULA has the highly complicated representation and array manifold A(θ) contains no explicit propagator and phase-difference, by fusing the FDO with the stationary assumption on source signals and noises, it is proven by theoretical derivation that the azimuthonly ULA model has a system of linear equations with the block-Hankel matrix as coefficient matrix and the elementary power-sum symmetric functions of the propagators as the unknowns.
Matrix equation (21) and the corresponding derivation manner are based on the equal-interspacing azimuth-only ULA model. Since any ULA is always one-dimensional, the linearly-combined relation between the phase and its component phase-difference still holds, no matter how complicated function forms of the phase and phase-difference are, therefore the derivation manner for the azimuth-only ULA including difference operation can be extended to any ULA of the M-D array and further used to investigate the complicated 2-D DoA estimation problem.
If {h fm } M m=1 are known, the unknowns {b k } K k=1 can be estimated and explicitly visualized by solving the system of linear equations in (21).

C. THE SYSTEM OF LINEAR EQUATIONS WITH THE BDO
As shown in (7), phase-differences { ϕ(θ k )} K k=1 also can be explicitly extracted out by making using of the BDO on the negative phase. To generate the negative phase, here correlation function R x n x m between reference signal x n (t) and where the m-th correlation function R x n x m = E{x n (t)x H m (t)}. The correlation function R x n x m is expanded as where R x n s k = E{x n (t)s H k (t)} (k = 1, · · · , K ) and R x n w m = E{x n (t)w H m (t)}. In (24) the phase ϕ m (θ k ) of the k-th signal s k (t) incident on the m-th sensor is negative and equal to −(m − 1) ϕ(θ k ).
Owing to the introduction of negative phase, the BDO can be performed on the elements of the sequence S b 0 . Under the stationary assumption on source signals and noises, by making use of the similar derivation manner as used for the system of linear equations in (21), the difference equation in respect to m (m = M , · · · , K + 1) is given by Equation (25) is the other difference equation satisfied by the propagator of the incident source with the BDO on negative phase. Similarly, writing (25) for the indices m = M , · · · , K + 1 gives the other system of linear equations where the vector

D. THE ROOT-FINDING POLYNOMIAL OF THE AZIMUTH-ONLY ULA
After {b k } K k=1 are estimated and known by solving (21) or (26), there still requires one step to obtain the propagators {v k } K k=1 from the known {b k } K k=1 . In accordance with the distinctive function forms of {b k } K k=1 in (20), we construct the following K -degree monic polynomial in propagator variable v as It is not difficult to show that f (v) is factored as follows: Thus the propagators {v k } K k=1 are proven as the roots of the polynomial f (v). With the known propagators {v k } K k=1 , the individual azimuth angle θ k is straightforward calculated via Hence the DoA estimation of the azimuth-only ULA is taken as a root-finding problem of the higher-order monic polynomial with {b k } K k=1 as the explicit polynomial coefficients and the propagators {v k } K k=1 as the roots. As a simple example, let the number of sources K = 2, from (20), the coefficients b 1 and b 2 have the following simple and explicit forms as After b 1 and b 2 are estimated from the ULA data, the 2-degree monic polynomial f (v) is constructed as The closed-form solutions of two propagators are .

IV. THE DEGENERATE SPATIAL ARMA PROCESS OF THE INCIDENT SOURCES
A. DIMENSIONAL ANALYSIS ON DIFFERENCE EQUATIONS (19) AND (25) Difference equations (19) and (25) can be analyzed from the standpoint of physics. An evident difference in function forms of (19) and (25) (25) are arranged in the ascending order from m − K + 1 to m. The reason for the above difference is simple to explain from the dimensional analysis on two difference equations.
In effect, after noise correlation R w m x n is subtracted from R x m x n , the m-th element h fm is given by Then h fm is contributed only by source correlations (19) and (25) should physically obey the dimensional homogeneity principle (DHP) in which the propagator in each term of the same equation remains identical in dimension.
The dimension of propagator is reflected by the corresponding power. For example, the dimension of b k (k = 1, · · · , K ) is equal to k. Similarly, for h fm in (31), without considering x n (t) appearing in all the K terms, the dimension of propagator in each term is equal to (m − 1), while for h bm , is equal to −(m − 1). Now we perform the dimensional analysis on difference equation (19). The dimension of propagator in the left-hand term h f (m+K ) is equal to m + K − 1, each of the right-hand terms is the product of b i and h f (m+K −i) (i = 1, · · · , K ), then the dimension of propagator in each of right-hand terms is the sum of i and K + m − i − 1 and equal to m + K − 1. Hence difference equation (19) physically obeys the DHP with the arrangement of {h fi }

B. SIGNAL MODEL OF THE INCIDENT SOURCES OF THE AZIMUTH-ONLY ULA
Difference equations (19) and (25) also can be analyzed from the standpoint of signal processing.
Let c m+K = R x m+K x n and e m+K = R w m+K x n , c m+K −k = R x m+K −k x n and e m+K −k = R w m+K −k x n , then (19) is rewritten as Similarly, difference equation (25) is rewritten as As shown in (32) and (33), under the stationary assumption, via the FDO or BDO, signal model of the external sources impinging upon the azimuth-only ULA is definitely revealed as the so-called degenerate spatial ARMA(K , K ) process subject to the identical autoregressive and moving average parameters. The definition of degenerate ARMA(K , K ) process can be found in [31].
As the stationary assumption is narrowed as the conventional one in which source signals and noises are uncorrelated and the noises are supposed as the GWN (Gaussian white noise) with zero mean and identical variance σ 2 , then (32) is further simplified by reasonably choosing reference signal. For example, let reference index n = m + K , then only e m+K = σ 2 , other e m+K −1 , · · · , e m are equal to zero, from (32) we get The degenerate ARMA(K , K ) process is shown as the AR(K ) process with the following power spectral density In contrast to (35), indeed, the MUSIC was the spatial pseudospectral method.
In general, signal model of the incident sources is the degenerate spatial ARMA process under stationary assumption on source signals and noises. Meanwhile, it should be stressed that such an ARMA process is not a common one. Model coefficients {b k } K k=1 are required to remain the explicit and distinctive function forms and satisfy the DHP.
For convenience next, Equations (21) is called the modified Yule-Walker (MYW) system of linear equations with the FDO.

C. THE DISCUSSION OF THE ROOT-FINDING POLYNOMIALS IN THE DoA ESTIMATION
The various root-finding polynomials for the DoA estimation had been developed in the Root-Music [11], Min-norm method [8], IQML (Iterative quadratic maximum likelihood) algorithm [32], [33], Root-WSF (Weighted subspace fitting) and MODE (Method of direction estimation) [13], [34]. The Root-WSF algorithm had been recognized as 'a strong candidate for the best method for ULAs' [14].
It was an ordinary operation to prespecify the propagators as the roots of the polynomial. However, a key problem was how to construct the reasonable root-finding polynomial. There mainly existed two constructing manners.
In the Root-MUSIC and Min-norm method, the rootfinding polynomials stemmed from noise eigenvectors, the propagators and noise zeroes became the roots of the polynomial in accordance with the orthogonality between steering vectors and noise eigenvectors. Therefore the polynomials in the Root-MUSIC and Min-norm method were not the pure one in propagator. In estimation, the roots distributed in and on the unit circle, it was required to separate the propagators from noise zeroes. Otherwise, when noise zeroes got closer to the unit circle, the spurious DoA estimates arose.
In the IQML, Root-WSF and MODE, the ML criterions in terms of angles were parameterised in terms of another parameter vector whose elements were set as polynomial coefficients, which was the other constructing manner. To guarantee the roots on the unit circle, polynomial coefficients were imposed to satisfy the conjugate symmetry constraint.
The diverse constructing manners, nonuniform polynomial forms and the attached constraints all showed that the rootfinding polynomial still deserved further study, even for the simple azimuth-only ULA. In addition, from the standpoint of physics, the propagator is a physical parameter, one source corresponds to one propagator. Once the K propagators are prespecified as the roots of polynomial, the dimensions of propagator in all the terms of the root-finding polynomial are required to be equal in accordance with the DHP, which is a sound way of checking the validity of the root-finding polynomial.
Different from the root-finding polynomials mentioned above, the proposed polynomial shown in (23) is derived from the fused perspective of physics and signal processing. From the signal perspective, polynomial coefficients are explicit and distinctive in form. The propagators are clearly proved as the roots of the proposed polynomial. No other extraneous roots are involved and no conjugate symmetry constraint is imposed on polynomial coefficients. Meanwhile, the propagators are both the poles and zeroes of the so-called degenerate spatial ARMA(K , K ) process that are only distributed on the unit circle. From the standpoint of physics, it is shown that the proposed polynomial obeys the DHP in which the dimensions of all the terms of f (v) are equal to the number K of sources.
In addition, for the simple azimuth-only ULA modle, the reason on the different function forms of the root-finding polynomials induced by difference operation and ML parameterization is an interesting problem.

V. THE PROPERTIES OF THE MYW SYSTEM OF LINEAR EQUATIONS A. THE DoA ESTIMATION UNDER MULTIPATH SCENARIO
Multipath propagation due to various reflections was frequently encountered in the DoA estimation. Similar to the unknown DoAs, the uncorrelated/coherent properties of source signals reflecting multipath propagation are also unknown to the processor.
Rank deficiency of source covariance matrix induced by the presence of multipath propagation was main reason for performance degradation and the failure of the subspacebased methods. Considerable efforts had been spent to handle the DoA estimation of the highly correlated and coherent sources.
The spatial smoothing operation had been recognized as an effective preprocessing technique in which the total array was partitioned into the subarrays and generated the average of covariance matrices of subarray data, however, which still depended on the uncorrelated condition between source signals and noises.
By contrast, the proposed difference operation on phases can protect the applicability of the DoA estimation to multipath propagation under the more looser stationary condition on source signals and noises.
As an illustration, suppose the k-th source signal s k (t) is rewritten as s k (t) = α ki s i (t), where the coefficient α ki denotes the complex-valued attenuation of the k-th s k (t) with respect to the i-th signal s i (t), (i, k = 1, · · · , K ), then correlation function R s k x n changes as R s k x n = α ki E{s i (t)x H n (t)}, under the stationary assumption on source signals and noises, R s k x n containing α ki is still a time-independent constant and can be eliminated by difference operation as we use above. Hence Multipath propagation do not pose a problem and is unnecessary to be considered in estimation.

B. THE MAXIMUM NUMBER OF DETECTABLE SOURCES WITH THE COMBINED SYSTEM OF LINEAR EQUATIONS
To solve the MYW system of linear equations in (21), coefficient matrix H f and the vector r f are required to be known. However, owing to the presence of the unknown noise correlation R w m x n (m = 1, · · · , M ) in the corresponding element h fm , H f and r f comprising {h fm } M m=1 in effect are unknown, the MYW system of linear equations in effect is unsolvable.
The more concrete signal property of sources and noises is desirable or needed. When the stationary source signals and noises are further assumed to be uncorrelated and suppose the noises are the spatially GWNs with zero mean and unequal variances {σ 2 m } M m=1 , the DoA estimation is readily implemented by constructing and solving the eligible MYW system of linear equations without the presence of noise variance (PNV).

1) THE MAXIMUM NUMBER OF DETECTABLE SOURCES WITH THE MYW SYSTEM OF LINEAR EQUATIONS
Under the above assumption, the construction of the eligible MYW system of linear equations depends on the chosen VOLUME 8, 2020 reference signal. When the first sensor is chosen as reference one, namely n = 1, the unknown R w 1 x 1 is equal to the nonzero σ 2 1 . Except for the first element h f 1 containing R w 1 x 1 , other M − 1 elements {h fm } M m=2 without PNV can be employed to construct the eligible MYW system of linear equations. Similarly, when n = M , the elements without PNV are {h fm } M −1 m=1 . For other sensors chosen as reference sensor, the number of the elements without PNV is less than M − 1.
By analysis, the number n c of the elements without PNV, the number n e of the equations without PNV and the number K satisfy n e = n c − K (36) As reference index n = 1 or n = M , n c = M − 1.
With the single FDO, as n = 1, owing to noise variance σ 2 1 only appearing into h f 1 , in the system of linear equations in (21), only the first equation containing h f 1 is not eligible. However, the first equation can be replaced by the counterpart with the M -th sensor chosen as reference sensor. Thus there have a total of n e + 1 equations without PNV. As n e + 1 = K and n c = M − 1, the maximum number K max of detectable sources is reached, then After n c = M − 1 is substituted into (37), the maximum number K max of detectable sources is obtained and equal to int(M )/2. The same K max can be reached with the BDO.

2) THE MAXIMUM NUMBER OF DETECTABLE SOURCES WITH THE COMBINED SYSTEM OF LINEAR EQUATIONS
Each equation of the systems of linear equations in (21) and (26) has the different dimension of propagator, then all the equations without PNV are independent and readily combined to increase the number of detectable sources. The eligible combined system of linear equations is given by where H and r stand for the combined coefficient matrix and the vector of constant terms without PNV.
As an example, the first sensor is chosen as reference sensor to investigate the maximum number of detectable sources with the combined system of linear equations. The number n c of the elements without PNV is equal to 2(M − 1). In accordance to the construction of (21) and (26), the 2(M − 1) − 2K equations without PNV are constructed by the 2(M − 1) eligible elements, namely, n e = 2(M − 1) − 2K . After including the two counterpart equations with n = M to replace the first equations of (21) and (26) with n = 1, respectively, the number of the equations without PNV becomes n e + 2.
To estimate the K sources requires n e + 2 ≥ K , then The maximum number K max of detectable sources is given by It is shown that the maximum number K max is equal to that by making use of the subspace-based method based on the FSBB, however, is completed by the portion of the elements of data covariance matrix in the proposed method.

C. THE CALCULATION OF NOISE VARIANCES
Signal model of the external sources impinging upon the azimuth-only ULA is revealed as the degenerate spatial ARMA(K , K ) process. As source signals and noises are assumed to be uncorrelated and suppose the noises are the spatially GWNs with unequal variances {σ 2 m } M m=1 , by reasonably choosing reference sensors, the DoA estimation and the calculation of noise variances are performed separately. Noise variances of all the sensors are calculated from the ULA data with the known {b k } K k=1 . When the first sensor x 1 (t) is chosen as reference one, the unknown noise variance σ 2 1 only appears in the first equation of the MYW system of linear equations in (21) and given by Once {b k } K k=1 are determined, the variance σ 2 1 is calculated as Similarly, as the M -th sensor is chosen as reference one, the variance σ 2 M can be extracted from the last equation and is given by With the similar manner, we choose other sensors as reference sensor, respectively, and calculate the corresponding noise variance from the corresponding equation. Each of noise variances has the different calculation formula.
It is a distinct advantage of the proposed method to extract more parameters including the DoAs and noise variances from the ULA data.

VI. THE IMPLEMENTATION OF THE DoA ESTIMATION OF THE AZIMUTH-ONLY ULA A. THE TLS ESTIMATION FOR THE COMBINED SYSTEM OF LINEAR EQUATIONS
In practice, the true data correlations {R x m x n } M m=1 are usually unavailable and replaced by the finite sampling ones where P is the number of snapshots, (·) * denotes conjugate operation. x m (p) and x n (p) are the sampling data with p as the discreted time index. After noise correlation is removed from the corresponding data correlation, the combined system of linear equations in (38) is still subject to various unknown 'residual' errors. Hence the determination of {b k } K k=1 is a classical total least squares (TLS) problem. Central to the TLS is stated as perturbing both H and r in a minimal fashion to obtain a consistent matrix equation.
where E and δ are the minimal TLS perturbations. Equation (48) is rewritten as The TLS aims to solve the following constrained minimization problem where · F denotes the Frobenius norm. The TLS property does not depend on the distribution of 'residual' errors. The TLS estimates can be obtained from the SVD of augmented matrix, which was discussed in [31], [35]- [37].

B. THE IMPLEMENTATION OF THE DoA ESTIMATION AND THE CALCULATION OF NOISE VARIANCES
Here the uncorrelated assumption on source signals and noises is considered and the noises are preset as the GWNs with zero mean and unequal variances {σ 2 m } M m=1 . As the degenerate spatial ARMA(K , K ) process, the construction of the combined system of linear equations without PNV is readily realized by flexibly choosing reference signals, especially for the ULA with the large number of sensors.
To implement the DoA estimation, it is first assumed that the number K of incident sources is known or estimated by the existing enumeration techniques in advance.
The step-by-step description of the proposed method is summarized as follows: Step 1. Construct and solve the sampling combined system of linear equations without PNV to obtain the estimates of {b k } K k=1 .
Step 3. Construct the root-finding polynomial with the known {b k } K k=1 , obtain the estimates of {v k } K k=1 and {θ k } K k=1 via the formulas of roots of polynomial or the numerical rootfinding algorithm.
If only the DoAs need to be estimated, Step 2 can be cancelled.

VII. NUMERICAL ANALYSIS
The validity and performance of the proposed method are first evaluated by several numerical examples. We consider four far-field narrow-band sources (K = 4) impinging upon an azimuth-only ULA comprising seven sensors (M = 7) with sensor separation of half the wavelength of signal. Azimuth angles [θ 1 , θ 2 , θ 3 , θ 4 ] are preset as [−60 0 , −40 0 , 20 0 , 35 0 ]. The ULA are corrupted by the so-called spatially colored noise. For simplicity, here the sources are supposed to have zero mean and unitary power, while sensor noises have the unequal powers {σ 2 m } 7 m=1 . In simulation, source signals and noises all are generated by the complex-valued zero-mean Gaussian white random process. The number of snapshots is equal to 1000.
For the seven-sensor ULA, only three DoAs can be estimated by the single system of linear equations with the FDO or BDO. However, in accordance with (40), the number of the detectable DoAs can be up to four with the combined system of linear equations.
In simulation, following Step1, by means of the first sensor as the reference sensor, a total of six correlation elements {h 1 fm } 7 m=2 and the corresponding conjugate counterparts {h 1 bm } 7 m=2 without PNV in data covariance matrix are chosen to construct the four independent equations in which the superscript 1 is used to describe the first sensor data as reference data. The four independent equations are enough to form the eligible combined system of linear equations to estimate the unknowns {b k } 4 k=1 . So the equations with the seven-th sensor as reference sensor are not considered.
The concrete form of the used combined system of linear equations in simulation is constructed as where h 1 fm = R x m x 1 − R w m x 1 and h 1 bm = R x 1 x m − R x 1 w m with the first sensor data x 1 (t) as reference data. According to (20), the unknowns {b k } 4 k=1 is given by Under the uncorrelated assumption between source signals and noises and further suppose the noises are the GWNs, theoretically, R w m x 1 and R x 1 w m are all equal to zero as m = 1, then the used elements   Carlo trial, with the known {b k } 4 k=1 , the root-finding polynomial is given by Here numerical approach in which the root-finding problem of the complex-valued propagator is converted into the real-valued angle-searching operation is proposed. We construct the cost function h(v) = 1/ log 10 [|f (v)|] and plot its amplitude for points versus θ in the interval of [−π/2, π/2] and pick peak positions as the estimated angles. In anglesearching process, in order to obtain more accurate estimates,  the effective off-gird searching algorithms are highly recommended [38], [39].

A. THE UNCORRELATED CASES OF SENSOR NOISES
Example 1: We first assume that sensor noises are uncorrelated with one another and have an identical noise variance, which is termed the so-called spatially white. In simulation, noise variance is preset via the unitary signal power and SNR = 2dB. Since the FBSS (Forward and backward spatial smoothing) had been recognized as an effective preprocessing technique to combat multipath scenario, the comparison is focused on the proposed method and the MUSIC based on the FBSS (FBSS-MUSIC). The DoAs of the uncorrelated and coherent sources can be well estimated by two methods in Fig.1(a) and (b). Besides the estimated DoAs, {σ 2 m } 7 m=1 also can be directly calculated by the proposed method and shown in Fig.1(c) and (d).
Example 2: Sensor noises are still assumed to be uncorrelated but with unequal variances {σ 2 m } M m=1 , which is the simple spatially-colored noise background. In simulation, noise variances are preset via the unitary power of signal and the SNR vector randomly generated by computer program within a prespecified range [−5dB, 5dB]. The azimuth-searching results of the uncorrelated and coherent sources are plotted in Fig.2(a) and (b), the calculated noise variances with the proposed method are in Fig.2(c) and (d).
As shown in Fig.1 and 2, the estimates of azimuth angles, regardless of under the equal or unequal SNRs, all can be obtained by the proposed method and FBSS-MUSIC, which shows that the spatially homogeneous or colored noises in effect have the less impact on the DoA estimation when the noises are uncorrelated with one another. The extraction of noise variances from the ULA data is an advantage of the proposed method.

B. THE COHERENT CASES OF SENSOR NOISES
Example 3: The proposed method can be applied to the background in which the noises are correlated/coherent with a 'correlated/coherent length' [40], [41] when source signals and noises are still uncorrelated. In simulation, sensor noises are assumed to have an identical variance under equal SNR with two coherent lengths. Owing to employing the different elements and amounts of data, the RMSEs ofθ are separately plotted to show the behaviors of two methods. As shown in Fig.4(a) and (c), as the coherent length increases from 3 to 4, estimation performance gradually degrades, which spreads from the low to moderate SNR. By contrast, the proposed method performs more stable.
Example 4: The azimuth-searching results with the coherent length of six noises under equal SNR and unequal SNRs are obtained and plotted in Fig.4. As source signals are uncorrelated or coherent, the FBSS-MUSIC all breaks down, regardless of equal or unequal SNRs. It is shown that the coherent length of noises has a significant impact on the FBSS-MUSIC. By contrast, the proposed method still provides the correct DoA estimation even with one independent noises. Hence the proposed method also can be applicable to the so-called partly calibrated arrays in spatially correlated noise fields in which the noise in the calibrated sensors is uncorrelated with the noise in the other sensors [42].

VIII. EXPERIMENTAL VERIFICATION OF THE PROPOSED DoA ESTIMATION METHOD
The proposed DoA estimation method is verified by a simple sound experiment. Fig.5 shows the schematic diagram of the experimental set-up comprising three microphones and one sound source. Cellphone 2 is taken as one source emitting a sinusoidal waveform. The Cellphone 1 is added to align with  three microphones for time synchronization. The frequency of sinusoidal waveform f = 1.0kHz and microphone interspacing d = 0.17m. With time synchronization, the DoA is determined by the proposed method after the received waveforms of the microphones are transformed into the complex data with the Hilbert transform. Meanwhile, the differences between time-of-arrivals (DToAs) of three received waveforms also can be readily marked to calculate the DoA. The two types of the DoA results are compared to analyze and verify the proposed method. Fig.6 shows the signals received by three microphones in one experiment in which the signals emitted by Cellphone 1 is employed for time synchronization. After time synchronization, the used waveforms and the corresponding results of three experiments are displayed in Fig.7. The source has the different incident angle in each experiment. For the first experiment, the starting waveforms are intentionally chosen and plotted in Fig.7(a). For other experiments, the non-starting waveforms are plotted in Fig.7(c) and (e), respectively.
We take Fig.7(a) for an example to show the DToAs method, the time t i (i = 1, 2, 3) of the first peak of the waveform received by the i-th sensor is marked as the corresponding time-of-arrival, so t im = t i − t m (i > m) denotes time-difference between the i-th and m-th sensors, then the angle θ t im is determined via t im from t i and t m . Each concrete θ t im is calculated as described above. Here the averaging valuē θ t of the corresponding θ t 32 , θ t 21 and θ t 31 is used to describe the DoA estimate with the DToAs.
In the proposed method, data covariance matrix is first constructed by windowing three received waveform. the lengths of the windowed waveform increase from 200 points to 1000 points with the increment of 20 points and the zero starting point. Employing (21) and (22), each windowing operation provides one result for each concrete angle θ a im with the following expression where x h i (p) denotes the Hilbert transform of the i-th sampling data x i (p), the correlation functionR x h i x h n is obtained by making use of (44) with x h n (p) as the reference signal and n = i, n = m. Re{.} stands for real part.
The results of θ a 32 , θ a 21 and θ a 31 of each incident angle versus the windowing times are plotted in Fig.7(b), (d) and (f), respectively. Because the non-signals before the received waveforms in Fig.7(a) are included in the windowing operations and result in the inaccurate and unstable results in calculation, only as the windowing times is larger than 12, the impact of the non-signals disappears in Fig.7(b). It is shown that all the results of each concrete angle θ a im are almost equal and nearly form the straight line with the averaging valueθ a im . Similarly, the total averaging valueθ a of the correspondingθ a 32 ,θ a 21 andθ a 31 is taken as the DoA estimate with the proposed method. The angle results with the DToAs and the proposed method are listed in TABLE 1. It is shown that VOLUME 8, 2020 each concreteθ a is well consistent with the correspondingθ t , which verifies the effectiveness of the proposed method.
As shown in TABLE 1, an interesting phenomenon is observed for Source 2, the large value difference arises between θ t 32 andθ a 32 , and also between θ t 21 andθ a 21 . However, qualitatively and quantitatively, the values ofθ a 32 andθ a 21 in Fig.7(d) are observed to be more consistent with the corresponding time-of-arrivals of the received waveforms shown in Fig.7(c). The main reason lies in that the time-of-arrival of the waveform of sensor 2 is slightly varied in propagation in the practical experiment and become more closer to that of sensor 3 and away from that of sensor 1 in Fig.7(c), which is well reflected byθ a 32 andθ a 21 with the proposed method.

IX. CONCLUSION
The DoA estimation of the azimuth-only ULA is investigated in theory and method from the fused perspective of physical and signal properties of the incident sources induced by the orderly-setup of array sensors. Under the mild stationary assumption on source signals and noises, it is theoretically proved that signal model of the external sources incident on the azimuth-only ULA is a degenerate spatial ARMA process subject to the identical autoregressive and moving average parameters. The azimuthonly ULA has its own distinctive MYW system of linear equations and the root-finding polynomial. The DoAs of the azimuth-only ULA are estimated directly from the ULA data via the chain comprising the MYW system of linear equations and the root-finding polynomial. The EVD is not involved.
The proposed method can be applied to multipath scenario. The DoAs estimation and the calculation of noise variances are separately performed under the backgrounds of spatially white and colored noises. The maximum number of the detectable sources is equal to floor[ 2M 3 ] by making use of the combined system of linear equations with the FDO and BDO, where M is the number of sensors.
The MYW system of linear equations and the root-finding polynomial are analyzed to obey the DHP, which manifests that the conclusion of the degenerate spatial ARMA process is self-consistent both in physics and signal processing.