{'fraction_non_alphanumeric': 0.049296046287367404, 'fraction_numerical': 0.017357762777242044, 'mean_word_length': 4.139968279143537, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sparse Bayesian Learning for DOA Estimation in Heteroscedastic Noise

The paper considers direction of arrival (DOA) estimation from long-term observations in a noisy environment. In such an environment the noise source might evolve, causing the stationary models to fail. Therefore a heteroscedastic Gaussian noise model is introduced where the variance can vary across observations and sensors. The source amplitudes are assumed independent zero-mean complex Gaussian distributed with unknown variances (i.e. the source powers), inspiring stochastic maximum likelihood DOA estimation. The DOAs of plane waves are estimated from multi-snapshot sensor array data using sparse Bayesian learning (SBL) where the noise is estimated across both sensors and snapshots. This SBL approach is more flexible and performs better than high-resolution methods since they cannot estimate the heteroscedastic noise process. An alternative to SBL is simple data normalization, whereby only the phase across the array is utilized. Simulations demonstrate that taking the heteroscedastic noise into account improves DOA estimation.


I. INTRODUCTION
With long observation times weak signals can be extracted in a noisy environment. Most analytic treatments analyze these cases assuming Gaussian noise with constant variance. For long observation times the noise process though is likely to change with time causing the noise variance to evolve. This is called a heteroscedastic Gaussian process, meaning that the noise variance is evolving. While the noise variance is a nuisance parameter that we are not interested in, it still needs to be estimated or included in the processing in order to obtain an accurate estimate of the weaker signals.
Accounting for the noise variation is certainly important for machine learning [1], [2] and related to robust statistics [3], [4]. Heteroscedastic noise models have been used in e.g. finance [5] and image processing [6]. In statistical signal processing, the noise has been assumed to vary spatially [7], [8], [9], but spatiotemporally varying noise as considered here has not been studied. The proposed processing could be applied to spatial coherence loss [10], [11], [12] or to wavefront decorrelation, where turbulence causes the wave front to be incoherent for certain observations (thus more noisy). This has lead to so-called lucky imaging in astronomy [13] or lucky ranging in ocean acoustics [14], where only the measurements giving good results are used. As a result an involved hypothesis testing is needed to determine the measurements to be used. In contrast, we propose to use all measurements.
In applications, a simple way to account for noise power variations is to normalize the data magnitude to only contain phase information as demonstrated for beamforming in seismology [15], [16], noise cross correlation in seismology [17], [18], [19], [20], [21] and acoustics [22], source deconvolution in ocean acoustics [23], [24] and speaker localization [25], [26]. High-resolution beamformers such as MUSIC [27] relying on a sample covariance matrix are not likely to perform well for heteroscedastic noise as the loud noise might dominate the sample covariance matrix. More advanced methods than an eigenvalue decomposition are needed to separate the signal and noise subspaces. We demonstrate that for well-separated sources normalizing the data to only contain phase information works well.
When the sources are closely spaced, more advanced parametric methods are needed for DOA estimation when the noise power is varying in space and time and the sources are weak. We derive and demonstrate this for the application of multiple measurement vector (MMV, or multi-snapshot) compressive beamforming [28], [29], [30], [31]. We solve the MMV problem using the sparse Bayesian learning (SBL) framework [30], [32], [33] and use the maximum-a-posteriori (MAP) estimate for DOA reconstruction. We assume the source signals to jointly follow a zero-mean multivariate complex normal distribution with unknown power levels. The noise across sensors and snapshots also follows a zero-mean multivariate normal distribution with unknown variances. These assumptions lead to a Gaussian likelihood function.

A. Heteroscedastic noise observation model
For the lth observation snapshot, we assume the linear model where the dictionary A ∈ C N ×M is constant and known, and the source vector x l ∈ C M contains the physical information of interest. Further, n l ∈ C N is additive zero-mean circularly symmetric complex Gaussian noise, which is generated from a heteroscedastic Gaussian process n l ∼ CN (n l ; 0, Σ n l ). Due to the circular symmetry of the noise the phase is uniformly distributed. We specialize to diagonal covariance matrices, parameterized as where J n = diag(e n ) = e n e T n with e n the nth standard basis vector. Note that the covariance matrices Σ n l are varying over the snapshot index l = 1, . . . , L and we introduce the matrix of all noise standard deviations as where R 0+ denotes non-negative real numbers. We consider three special cases for the a priori knowledge on the noise covariance model (2) which are expressed as constraints on V N . Case I We assume wide-sense stationarity of the noise in space and time: σ 2 nl =σ 2 . The model is homoscedastic, Case II We assume wide-sense stationarity of the noise in space only, i.e., the noise variance for all sensor elements is equal across the array, σ 2 nl = σ 2 0l and it varies over snapshots. The noise variance is heteroscedastic in time (across snapshots), Case III No additional constraints other than (3). The noise variance is heteroscedastic across both time and space (sensors and snapshots.). In this case V N ∈ V III = R N ×L 0+ . The relation between these noise cases is V I ⊂V II ⊂V III = R N ×L 0+ . From the sets V d (d is I, II, or III) both V N in (3) and Σ n l in (2) can be constructed.

B. Array model
Let X = [x 1 , . . . , x L ] ∈ C M ×L be the complex source amplitudes, x ml = [X] m,l = [x l ] m with m ∈ {1, · · · , M } and l ∈ {1, · · · , L}, at M DOAs (e.g., θ m = −90 • + m−1 M 180 • ) and L snapshots for a frequency ω. We observe narrowband waves on N sensors for L snapshots Y = [y 1 , . . . , y L ] ∈ C N ×L . A linear regression model relates the array data Y to the source amplitudes X as The dictionary A=[a 1 ,...,a M ]∈C N ×M contains the array steering vectors for all hypothetical DOAs as columns, with the (n, m)th element given by e −j ωdn c sin θm (d n is the distance to the reference element and c the sound speed).
We assume M >N and thus (6) is underdetermined. In the presence of only few stationary sources, the source vector x l is K-sparse with K M . We define the lth active set and assume M l =M={m 1 ,...,m K } is constant across all snapshots l. Also, we define A M ∈C N ×K which contains only the K "active" columns of A. In the following, · p denotes the vector p-norm and · F the matrix Frobenius norm. Similar to other DOA estimators, K can be estimated by model order selection criteria or by examining the angular spectrum. The parameter K is required only for the noise power in the SBL algorithm. An inaccurate estimate influences the algorithm's convergence.

C. Prior on the sources
We assume that the complex source amplitudes x ml are independent both across snapshots and across DOAs and follow a zero-mean circularly symmetric complex Gaussian distribution with DOA-dependent variance γ m , m = 1, . . . , M , i.e., the source vector x l at each snapshot l∈{1,···,L} is multivariate Gaussian with potentially singular covariance matrix, as rank(Γ)=card(M)=K≤M (typically K M ). Note that the diagonal elements of Γ, denoted as γ γ γ γ γ, represent source powers and thus γ γ γ γ γ ≥ 0. When the variance γ m =0, then x ml =0 with probability 1. The sparsity of the model is thus controlled with the parameter γ γ γ γ γ, and the active set M is equivalently The SBL algorithm ultimately estimates γ γ γ γ γ rather than the complex source amplitudes X. This amounts to a significant reduction of the degrees of freedom resulting in a low variance of the DOA estimates.

D. Stochastic likelihood
We here derive the well-known stochastic maximum likelihood function [44], [45], [46]. Given the linear model (6) with Gaussian source (9) and noise (2) the array data Y is Gaussian with for each snapshot l the covariance Σ y l given by The probability density function of Y is thus given by The L-snapshot log-likelihood for estimating γ γ γ γ γ and V N is This likelihood function is identical to the Type II likelihood function (evidence) in standard SBL [36], [35], [32] which is obtained by treating γ γ γ γ γ as a hyperparameter. The Type II likelihood is obtained by integrating the likelihood function over the complex source amplitudes, cf. (29) in [32]. The stochastic maximum likelihood approach is used here as it is more direct.
The parameter estimatesγ γ γ γ γ and V N are obtained by maximizing the likelihood, leading to where V d is the feasible set of noise variances V N in (3) corresponding to the noise cases (d = I, II, or III, see Sec. I-A).The likelihood function (14) is similar to the ones derived for SBL and LIKES [34]. If γ γ γ γ γ and Σ n l are known, then the MAP estimate is the posterior meanx MAP l and covariance Σ x l [47], [32], The diagonal elements of Γ, i.e., γ γ γ γ γ, control the row-sparsity ofx MAP l as for γ m = 0 the corresponding mth row ofx MAP l becomes 0 T .

E. LASSSO versus SBL
Both LASSO and SBL use the linear model (6) with complex zero-mean Gaussian random noise but they differ in the modeling of the source matrix X.
The LASSO approach assumes a priori X random with uniformly i.i.d. distributed phase and Laplace-like prior amplitudes, Thus only the summed amplitudes (19) are Laplacian. The elements in X are unknown and must be estimated. The LASSO approach uses the conditional likelihood (Type I) for p(Y|X) and applies Bayes rule with the prior p(X) giving the MAP estimate The LASSO approach (20) estimates the realization of x l for each snapshot l. The number of parameters to be estimated for the LASSO approach grows linearly with the number of snapshots.
The SBL approach on the other hand, assumes that every column of X is random with the same complex zeromean Gaussian a priori distribution CN (x l ; 0, Γ). The SBL approach uses the likelihood p(Y|γ γ γ γ γ, Σ N ) in (15) and estimates γ γ γ γ γ but not the realization of X. Therefore, the number of parameters to be estimated for the SBL approach is independent of the number of snapshots. Thus the quality of the estimate improves faster than with LASSO.
For the heteroscedastic noise model, a different covariance matrix could be included in the data fit of (20). However, since the number of parameters to be estimated grows with the number of snapshots, this seems less attractive than using SBL.

F. Pre-whitening
The purpose of this section is to motivate the empirical evidence [15]- [26] that phase-only processing might provide improved estimates over using also amplitudes. This is most likely to happen when sources are not closely spaced and at low SNR as demonstrated in the examples, Sec. V.
Let us introduce the factorization Σ −1 n l = W H l W l , where W l is a square and non-singular matrix. For our particular setup, we have The matrix W l is useful for pre-whitening the sensor data. The corresponding whitened sensor data and dictionary matrix are For known diagonal noise covariance Σ n l the above means we have to normalize each row in (1) with σ nl as then the noise satisfies n l ∼ CN ( n l ; 0, I), and thus all entries are identically distributed.
If the noise covariance matrices are not known, they can be estimated as we will show later on. More specifically, at low SNR, the Noise Cases I, II, and III, lead to the noise variance estimates (43), (44). The corresponding pre-whitening matrices can then be computed as leading to Empirically, it has been found that applying pre-whitening to the data y l only (the dictionary A is non-whitened) effective in finding the strongest DOA.
For Noise Case I, pre-whitening does not play a role and the conventional beamformer is formulated using the power spectrum at DOA θ m where the sample covariance matrix is For Noise Cases II and III on the other hand, we will work with the pre-whitened data. In those cases, , the conventional beamformer leads to the powerspectra P CBF2 and P CBF-Phase , respectively, which can be formulated as where the whitened sample covariance matrix is The weighting W l is given by (22). For Case III, only phase is used as can be observed from (23), which results in phaseonly processing. As demonstrated in Section V-B this simple pre-whitening can improve the DOA performance at low SNR.

II. SOURCE POWER ESTIMATION
Let us now focus on the SBL algorithm solving (15). The algorithm iterates between the source power estimateŝ γ γ γ γ γ derived in this section and the noise variance estimates V N computed in Sec. III. After detailing these estimation procedures, the full algorithm is summarized in Sec. IV.
We impose the diagonal structure Γ=diag(γ γ γ γ γ), in agreement with (9), and form derivatives of (14) with respect to the diagonal elements γ m , cf. [44]. Using the derivative of (14) is formulated as For the solution to (15), we impose the necessary condition Assuming γ old m and Σ y l given (from previous iterations or initialization) and forcing (32) to zero, we obtain the following fixed point iteration [48] for the γ m : We use b = 0.5, but have not carefully tested for optimal values and this value depends on many factors such as which noise estimate is used and closeness of the DOAs. A value of b = 1 gives the update equation used in [47], [30], [33] and b = 0.5 gives the update equation used in [32].

III. NOISE VARIANCE ESTIMATION
While there has been more focus on estimating the DOAs or γ γ γ γ γ, noise is an important part of the physical system and a correct estimate is needed for good convergence properties. In SBL, the noise variance controls the sharpness of the peaks in the γ γ γ γ γ spectrum, with higher noise levels giving broader peaks. Thus, as we optimize the DOAs we expect the noise levels to decrease and the γ γ γ γ γ spectrum to become sharper. In this section we estimate the noise variance for the three noise cases in Sec. I-A. In Secs. III-A-III-C, we will assume the support of γ γ γ γ γ is known.

A. Noise estimate, Case I
Under Noise Case I, where Σ n l = σ 2 I N with I N the identity matrix of size N , stochastic maximum likelihood [38], [41], [43] can provide an asymptotically efficient estimate of σ 2 if the set of active DOAs M is known.
Let Γ M =diag(γ γ γ γ γ new M ) be the covariance matrix of the K active sources obtained above with corresponding active steering matrix A M which maximizes (14). The corresponding data covariance matrix is Note that for Noise Case I, the data covariance matrices (12) and (34) are identical. Following [42], we continue from (30), for all active sources (m ∈ M). Since range(Σ −1 To obtain Jaffer's condition below, we impose the following additional constrains a H m (S y − Σ y l ) a p = 0, ∀m, p ∈ M, m = p.
Together, the conditions (37a), (37b) give Jaffer's condition ( [42]:Eq.(6)), i.e., which we enforce at the optimal solution (Γ M , σ 2 ). Jaffer's condition follows from allowing arbitrary correlations among the source signals, i.e. when the Γ matrix is not restricted to be diagonal. Substituting (34) into (38) gives Let us then define the projection matrix onto the subspace spanned by the active steering vectors Evaluating the trace, using tr PP H =K and tr PS y P H = tr(PS y ), gives where we have defined =tr[Σ y l −S y ].
The above approximation motivates the noise power estimate for Noise Case I (43), which is error-free if tr[Σ y ]= tr[S y ], unbiased because E[ ] = 0, consistent since also its variance tends to zero for L→∞ [49], and asymptotically efficient as it approaches the CRLB for L→∞ [50]. Note that, the Noise Case I estimate (43) is valid for any number of snapshots, even for just one snapshot.

B. Noise estimate, Case II
For Noise Case II, where Σ l = σ 2 l I N , we apply (43) for each snapshot l individually, leading tô Several alternative estimators for the noise variance are proposed based on EM [36], [30], [37], [47], [51]. For a comparative illustration in Sec. V-D1 we will use the iterative noise estimate from [36], given by where the posterior mean x MAP l and covariance Σ x l are given in (16) and (17). Empirically, EM noise estimates such as (45) significantly underestimate the true noise variance in our applications.

C. Noise estimate, Case III
Let us start from the definition of the noise covariance This motivates plugging-in the single-observation signal estimate x M,l = A + M y l ∈ C K for the active (non-zero) entries in x l . This estimate is based on the single observation y l and the projection matrix (40), giving the rank-1 estimate Σ n l = (I − P)y l y H l (I − P).
Since the signal estimate x M,l maximizes the estimated signal power, this noise covariance estimate is biased and the noise level is likely underestimated.
Since we assume the noise independent across sensors, all off-diagonal elements of Σ n l are known to be zero. With this constraint in mind, we modify (47) as The estimate (48) is demanding as for all the N × L complexvalued observations in Y, we obtain N × L estimates of the noise variance. Note that the estimate Σ n l in (47) is not invertible whereas the diagonal constraint in (48) leads to a non-singular estimate of Σ n l with high probability (it is singular only if an element of y l is 0). As a result, the expression for Σ y l that is used for estimating γ γ γ γ γ in (33) is likely invertible.
On the other hand, an overestimate of the noise is easily obtained by assuming x l = 0 which is equivalent to setting P = 0 in (48), resulting in Σ n l = diag diag y l y H l (49) or σ nl = |y nl |.
This can be shown be the maximum likelihood estimate for no sources (M = 0) or very low power sources.

IV. SBL ALGORITHM
The complete SBL algorithm is summarized in Table I. The same algorithm is valid for Noise Cases I, II, and III (at high SNR). Given the observed Y, we iteratively update Σ y l (12) by using the current γ γ γ γ γ and Σ n l . The Σ −1 y l is computed directly as the numerical inverse of Σ y l . For updating γ m , m = 1, . . . , M we use (33). For the initialization of γ γ γ γ γ we use the conventional beamformer (CBF) power Based on the corresponding noise case, (43), (44), (45), or (48) is used to estimate Σ n l . The noise is initialized using (43), (44), (45), or (48) with P = 0, K = 0, which provides an over estimate of the noise variance.
The convergence rate measures the relative change in estimated total source power, The algorithm stops when ≤ min and the output is the active set M (see (11)) from which all source parameters are computed. In both cases the standard deviation has mean 1.

V. EXAMPLES
When the noise process is stationary and only one source is active, the peak of the CBF spectrum provides the optimal estimate of the source location. With heteroscedastic noise this may not hold true. We propose to use a robust CBF (Sec.V-B) or find the DOA with SBL, which takes the heteroscedastic noise into account (Sec. V-C). With multiple sources and heteroscedastic noise we use SBL (Sec. V-D).

A. Simulation setup
In the analysis of seismic data the noise for each snapshot was observed to be log-normal distributed [52]. Noise has also been modeled with extreme distributions [53]. In the simulations here, the noise follows a normal-uniform hierarchical model. The noise is complex zero-mean Gaussian with the noise standard deviation uniformly distributed over two decades, i.e., log 10 σ nl ∼U(−1,1), where U is the uniform distribution. Three noise cases are simulated: (a) Case I: constant noise standard deviation over snapshots and sensors, (b) Case II: standard deviation changes across snapshots with log 10 σ l ∼ U(−1,1), and (c) Case III: standard deviation changes across both snapshots and sensors with log 10 σ nl ∼U (−1,1).
A realization of the noise standard deviation is shown for Noise Cases II (Fig. 1a) and III (Fig. 1b) for N =20 sensors and L=50 snapshots. The corresponding histogram is presented in Fig. 1c.

B. Single DOA with CBF and MUSIC
The single source is located at −3 • . The complex source amplitude is stochastic and there is additive heteroscedastic Gaussian noise with SNR variation from −40 to 0 dB. The sensor array has N =20 elements with half wavelength spacing. We process L=50 snapshots. The angle space [−90, 90] • is divided into a 0.5 • grid (M =360). The single-snapshot array signal-to-noise ratio (SNR) is We first compute the beam spectra using CBF (24), CBF2 and CBF-Phase both (26) but with different weighting (23). When the noise is homoscedastic (constant standard deviation), the beam spectra for the three processors behave similarly (first row in Fig. 2). For heteroscedastic noise CBF2 and CBF-Phase give much better discrimination between signal to noise, see Fig. 2 middle (Noise Case II) and bottom row (Noise Case III).
The root mean squared error (RMSE) of the DOA estimates over N sim =100 runs with random noise realizations is used for evaluating the algorithm where θ 0 i is the true DOA andθ n i is the estimated DOA for the ith source when K sources are present.  The SNR curves (Fig. 3) demonstrate increased robustness of CBF2 and CBF-Phase, failing 20 dB later. Due to the heteroscedastic noise, MUSIC performs worse than CBF for a single source. CBF-Phase diverges at an SNR 15-20 dB later than CBF for Noise Cases II and III.

C. Single DOA estimation with SBL
We use the following SBL methods, with γ γ γ γ γ update (33) and SBL: Solves Case I using standard SBL, with σ from (43). SBL2: Solves Case II, with σ l for each snapshot from (44). SBL3: Solves Case III, with σ nl from (48). In addition to these methods we also use basis pursuit (BP) as implemented in [31].
For Noise Case I, all the methods fail near the same SNR of −12.5 dB, Fig. 4a. When the noise is heteroscedastic across snapshots, CBF, BP, and SBL fail early. Since both SBL2 and SBL3 correctly model the noise, they perform the best for Noise Case II, Fig. 4b. CBF-Phase also performs well. For heteroscedastic Noise Case III, SBL3 fails the last and as before CBF-Phase also performs well. SBL3 fails 15 dB later than any other method. This demonstrates the usefulness of accurate noise models in DOA processing. The histograms (Fig. 5) of the DOA (location of peak in power spectrum) at SNR −25 dB for Noise Case III demonstrate that when the heteroscedastic noise is accounted for, the histograms concentrate near the true DOA (θ = −3 • ).
An example statistic of the heteroscedastic noise standard deviation is shown in Fig. 6. The standard deviation for each sensor is either 0 or √ 2. SBL3 estimates the standard deviation from (48). Average noise in Fig. 6b resembles well the true noise (Fig. 6a) whereas the sample estimate (Fig. 6c) has high variance. Fig. 6d plots the mean across simulations and snapshots of the estimated standard deviation. On average, the noise estimate is close to the true noise.

D. Three DOA estimation with SBL
Now consider three sources located at [−3, 2, 50] • with power [10,22,20] dB, see Fig. 7. Relative to the single source case, the CBF-Phase performs significantly worse than SBL3, as the sources at −3 • and 2 • are both in the CBF main lobe.
The localization ability of the algorithms can also be gauged from the histogram of the top three peaks, see Fig. 8. Since SBL3 accounts for the heteroscedastic Noise Case III, its histogram is concentrated near the true DOAs.
1) Noise estimate convergence: The performance of the SBL methods is strongly related to the accuracy of the noise estimates. For the lth snapshot, the true noise power σ 2 l,T (Noise Case II) is The estimated (σ 2 l ) new (44) deviates from σ 2 l,T (55) randomly. We also consider the noise estimate (45) based on the EM method. Fig. 9 compares the two noise estimates for noise generated with log σ l ∼ U(−1, 1). The evolution of the histograms of the relative noise variance σ 2 Est /σ 2 T with iterations is in Fig. 9b. The mean of the ratio of σ 2 Est /σ 2 T is close to 1 for SBL2 but is much lower for the EM noise estimate (45), this likely causes the SBL2 (EM noise, (45)) to fails 5 dB earlier. 2) Number of snapshots: The RMSE versus number of snapshots for Noise Case III (Fig. 10) shows that SBL3 is most accurate with CBF-Phase following.
3) Noise distribution: For Noise Case III, we are just using one observation to estimate the standard deviation for SBL3 (48). Thus the estimates are not good, see Fig. 11. The true noise standard deviation is generated from log 10 σ T,nl ∼ U(−1,1), see Fig. 11a. The distribution of the deviation from the true standard deviation σ Est −σ T (Fig. 11c) is well-centered (mean 0.007).

VI. CONCLUSION
Likelihood based methods for single and multiple DOA estimation from long-term observations corrupted by nonstationary additive noise are discussed. In such a setting, the DOA estimators for a stationary noise model perform poorly. Therefore a heteroscedastic Gaussian noise model is introduced where the noise variance varies across sensors and snapshots. We develop Sparse Bayesian Learning (SBL) approaches to estimate source powers, source DOAs, and the noise variance parameters. Simulations show that normalizing the array data magnitude (such that only the phase information is retained) is simple and useful for estimating a single source in heteroscedastic Gaussian noise. For the estimation of (multiple) closely spaced sources, a problem specific SBL approach gives a much lower RMSE in the DOA estimates.