Robust estimation of DOA from array data at low SNR

Signal Processing 166 (2020) 107262

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Robust estimation of DOA from array data at low SNRR Christoph F. Mecklenbräuker a,∗, Peter Gerstoft b, Erich Zöchmann a, Herbert Groll a a b

Institute of Telecommunications, TU Wien, Vienna 1040, Austria University of California San Diego, La Jolla, CA 92093-0238, USA

a r t i c l e

i n f o

Article history: Received 14 April 2019 Revised 14 August 2019 Accepted 19 August 2019 Available online 28 August 2019 Keywords: Array processing DOA Estimation Least absolute deviation LAD Laplace-like noise Heteroscedastic noise Phase-only processing

a b s t r a c t We consider direction of arrival (DOA) estimation for a plane wave hidden in additive circularly symmetric noise at low signal to noise ratio. Starting point is the maximum-likelihood DOA estimator for a deterministic signal carried by a plane wave in noise with a Laplace-like distribution. This leads to the formulation of a DOA estimator based on the Least Absolute Deviation (LAD) criterion. The phase-only beamformer (which ignores the magnitude of the observed array data) turns out to be an approximation to the LAD-based DOA estimator. We show that the phase-only beamformer is a well performing DOA estimator at low SNR for additive homoscedastic and heteroscedastic Gaussian noise, as well as Laplacelike noise. We compare the root mean squared error of several different DOA estimators versus SNR in a simulation study: the conventional beamformer, the phase-only beamformer, and the weighted phaseonly beamformer. The simulations indicate that the phase-only DOA estimator has desirable properties when the additive noise deviates from the Laplace-like assumption. The qualitative robustness of these DOA estimators is investigated by comparing the empirical influence functions. Finally, the estimators are applied to passive sonar measurements acquired with a horizontal array in the Baltic Sea. © 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license. (http://creativecommons.org/licenses/by-nc-nd/4.0/)

1. Introduction In a series of contributions, Böhme and students [1–4] developed several types of approximate maximum likelihood estimators for locating wideband sources in the presence of partly unknown noise fields from array data. In order to define approximate likelihood functions, the asymptotic normality of finite Fourier transformed array data is exploited [5]. Asymptotic normality holds when the so-called mixing conditions are fulfilled by the array data [5]. However, the effects of deviations from normality on approximate maximum likelihood estimates have been little investigated In this work, we investigate the robustness of several Direction of Arrival (DOA) estimators which are strongly related to Least Absolute Deviation (LAD) regression. It is well-known that LADregression is resistant to outliers in the data and becomes maximum likelihood when the additive noise follows a Laplace distribution [6]. Laplace distributions are often used to model impulsive noise in seismic, acoustic, and radio environments [7–9]. Multivari-

R ∗

Johann Böhme gewidmet zum 80. Geburtstag. Corresponding author. E-mail addresses: [email protected], [email protected] (C.F. Mecklenbräuker). URL: http://noiselab.ucsd.edu (P. Gerstoft)

ate Laplace distributions and their generalizations are discussed in Refs. [10,11]. In the context of compressive beamforming, a complex-valued Laplace-like probability density has previously been used to model the prior amplitude distribution of plane waves hidden in additive white Gaussian noise in [12–14]. The complex-valued Laplace-like distribution is similar to the multivariate Laplace distribution, but they should not be confused. In [15], the unknown complex amplitudes of multiple monochromatic plane waves and their unknown DOAs are modeled as unknown parameters and a Laplace-like probability density is assumed for the additive noise. We model the unknown complex amplitudes at multiple frequencies and its DOA of a single wideband plane wave arrival as unknown parameters and assume a Laplace-like probability density for the additive noise. Application of the maximum-likelihood principle leads to complex-valued least absolute deviation (LAD) regression for the unknown complex amplitude. The influence function [16] for LAD regression is universally bounded which guarantees that outliers in the observation data have a limited influence on the estimate of the unknown complex wave amplitude. This shows that LAD regression features desirable qualitative robustness [6]. In this work, we use LAD regression for the unknown complex wave amplitude (which we regard as a nuissance parameter) to formulate a DOA

https://doi.org/10.1016/j.sigpro.2019.107262 0165-1684/© 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license. (http://creativecommons.org/licenses/by-nc-nd/4.0/)

2

C.F. Mecklenbräuker, P. Gerstoft and E. Zöchmann et al. / Signal Processing 166 (2020) 107262

estimator based on multiple measurement vectors and multiple frequencies. This work is structured as follows: Section 2 introduces the array data model and derives several DOA estimators. Next, the relation between those DOA estimators and approximate prewhitening is discussed. Section 4 discusses simulations with synthetic data and various noise models. The empirical influence function is used in Section 5 to assess the qualitative robustness of the DOA estimators. Finally, the DOA estimators are applied to realworld data with and without additional contamination.

negative single-frequency log-likelihoods where the summation is taken over all available frequencies ω ∈  = {ω1 , . . . , ωJ },

L (θ , x , λ ) =

y  ( ω ) = a ( θ , ω ) x  ( ω ) + n ( ω ) ,

(  = 1, . . . , L )

θˆLAD =

ω dn ej c sin θ

and the nth element of a(θ , ω) is given by an (θ , ω ) = (dn is the distance to the reference element and c the phase speed). All elements of the steering vector a(θ , ω) have magnitude one and n is a complex random vector with independent identically distributed (i.i.d.) zero-mean circularly-symmetric Laplace-like noise. Its probability density function (pdf) is

p(n (ω )) =

N  k=1

=



λ (ω ) √ 2π

2

e−λ(ω )|nk (ω )|

(λ(ω ))2N −λ(ω )n (ω )1 e ( 2π )N

(2)

where λ(ω ) ∈ R+ is the unknown inverse scale parameter and the  1 -norm of a vector x ∈ CN is defined as x1 = N k=1 |xk |. We assume the noise to be independent across all measurement vectors ( = 1, . . . , L ). Then the pdf of the array data matrix Y (ω ) = [y1 (ω ), . . . , yL (ω )] with typical element yn (ω) is

p(Y (ω ); θ , x(ω )) =

L  N  =1 n=1

=



λ (ω ) √ 2π

2

e−λ(ω )|yn (ω )−an (θ ,ω )x (ω )|

L 

λ2N (ω ) −λ(ω )y (ω )−a(θ ,ω ) x (ω )1 e . ( 2π )N =1

Lω (θ , x(ω ), λ(ω )) = − ln p(Y (ω ); θ , x(ω )) = λ (ω )

L  =1

y (ω ) − a(θ , ω )x (ω )1 − LN ln

=−

λ (ω ) 2

aH (θ , ω ) W ,ω (θ ) [y(ω ) − a(θ , ω )x (ω )]

Estimators based on minimizing the criterion (5) are known as Least Absolute Deviation (LAD) estimators [6]. The LAD estimator for the complex wave amplitude is the ML-estimator for the unknown complex signal amplitude in Laplace-like noise. We asssume the noise to be independent across all frequency bins. Then the negative log-likelihood function is the sum over all

(8)

where we introduced the diagonal weighting matrix

W ,ω (θ ) = diag





1

(9)

|yn (ω ) − an (θ , ω )x (ω )|

n=1,...,N

Eq. (8) to zero shows that a stationary point in Lω satisfies

x (ω ) = xˆ,LAD (θ , ω ) =

aH (θ , ω ) W ,ω (θ ) y (ω ) , aH (θ , ω ) W ,ω (θ ) a(θ , ω )

(10)

for all  = 1, . . . , L. Unfortunately, this expression does not allow direct evaluation because the weighting matrix depends on the unknowns. In the following, we assume small amplitude of the incoming plane wave, i.e.,

|an (θ , ω )x (ω )|  min(|yn (ω )| ).

(11)

n

This allows to approximate (9) by

W 0,ω = diag



1

|yn (ω )|



n=1,...,N

.

(12)

A first approximation xˆ(1 ) to the ML estimate (10) is found using (12)

aH (θ , ω ) W 0,ω y (ω ) aH (θ , ω ) W 0,ω a(θ , ω )

(13)

= c (ω ) aH (θ , ω ) W 0,ω y (ω ) = c (ω ) aH (θ , ω ) y˜  (ω ) with

(14)



1 c ( ω ) = H = a (θ , ω ) W 0,ω a(θ , ω ) y˜  (ω ) = W 0,ω y (ω ) =

(5)

(7)

N  ∂ Lω ∂|yn (ω ) − an (θ , ω )x (ω )| = λ ( ω ) ∂ x∗ (ω ) ∂ x∗ (ω ) n=1

(3)

(4)

λ2 (ω ) . 2π

arg min L (θ , x , λ ), x ∈ C1×JL λ  0, θ ∈ 

xˆ(1) (θ , ω ) =

This is the single-frequency likelihood for DOA θ and the vector of complex source amplitudes x(ω ) = [x1 (ω ), . . . , xL (ω )] ∈ C1×L . Thus, the negative single-frequency log-likelihood is

(6)

where  is the feasible set of DOAs. The operator  denotes element-wise > . The numerical solution to (7) is found as follows. First, we optimize all single-frequency log-likelihoods analytically with respect to the complex source amplitudes x(ω ) ∈ CL ∀ω ∈ . The complex conjugate derivative (in the sense of the Wirtinger calculus) of Lω with respect to the complex source amplitude is the direction of largest change [6, Eq. (2.3)].

(1)

where a(θ , ω ) = (a1 (θ , ω ), . . . , aN (θ , ω ))T is the plane wave steering vector for a single wave with unknown complex source amplitude x (ω). We assume a linear sensor array consisting of N sensors

Lω (θ , x(ω ), λ(ω )),

ω ∈

where x ∈ C1×JL is the parameter vector which contains the x(ω) for all ω ∈  and similarly for λ . Here, we perform DOA estimation via minimization of (6), i.e.,

2. Single plane wave with Laplace-like noise A vector-valued time series is observed with a sensor array, lowpass filtered, sampled, analog to digital converted, and segmented into L data pieces of duration T each. Each data piece is transformed via the discrete-time Fourier transform to the frequency domain. We formulate the approximate array data model for multiple measurement vectors for a single plane wave with frequency ω arriving from DOA θ in additive independent identically distributed zero-mean complex-valued Laplace-like noise directly in the frequency domain.



−1

N 

1

n=1

|yn (ω )| T

 y yN 1 ,..., |y1 | |yN |

.

(15) (16)

Note that aH y˜  in (14) resembles a conventional beamformer (CBF), with the steering vector a applied to the phase-only array data vector y˜  . Thus aH y˜  becomes independent of the magnitude of the measurement data which indicates robustness against outliers. The weighting coefficient c in (14) is positive real and independent of DOA θ . The weighting coefficient does, however, depend on the array data magnitudes.

C.F. Mecklenbräuker, P. Gerstoft and E. Zöchmann et al. / Signal Processing 166 (2020) 107262

We average the squared magnitude of (14) over all measurement vectors which gives

B˜L (θ , ω ) =

L 1  (1 ) |xˆ (θ , ω )|2 = aH (θ , ω )Sy˜ (ω )a(θ , ω ) L

(17)

=1

where

Sy˜ (ω ) =

L 1 (c (ω ))2 y˜  (ω )y˜ H (ω ). L

(18)

=1

In the following, we refer to (17) as the weighted phase-only beamformer. We use the first approximation (13) to initialize an iteratively re-weighted least squares (IRWLS) algorithm until convergence [6, Algorithm 3] to the LAD estimate xˆ,LAD (θ , ω ). Substituting x (ω ) = xˆ,LAD (θ , ω ) into (5) gives

Lω (θ , xˆ LAD (θ , ω ), λ ) =

λ

L  =1

(19)

Here xˆ,LAD (θ , ω ) is the th element of xˆ LAD (θ , ω ), the LAD estimate for the complex source amplitude obtained by IRWLS. Analytically minimizing (19) with respect to λ gives the single-frequency ML estimate

λˆ LAD (θ , ω ) = L

=1

2LN

y (ω ) − a(θ , ω )xˆ,LAD (θ , ω )1

(20)

ˆ LAD (θ , ω ) into (19) gives the negative singleSubstituting λ = λ frequency log-likelihood function for the unknown DOA θ ,

ˆ LAD (θ )) Lω (θ , xˆ LAD (θ , ω ), λ = 2LN ln

L   =1

= LN ln

 y (ω ) − a(θ , ω )xˆ,LAD (θ , ω ) + const. 1

L  

 y (ω ) − a(θ , ω )xˆ,LAD (θ , ω )

=1

1

2 + const.

For snapshot  and frequency ω, we factorize the inverse noise covariance as (n (ω ))−1 = W H ,ω W ,ω , where W,ω is chosen to be square and positive definite. For heteroscedastic noise we choose −1 W ,ω = diag(σ1−1 ,,ω , . . . , σN,,ω )

(22)

where σ n,,ω is the standard deviation of the noise at the nth sensor, th measurement vector, at frequency ω. The matrix W,ω is useful for pre-whitening the sensor data. For known diagonal noise covariance nl , the corresponding whitened sensor data is

y (ω ) = W ,ω y (ω ).

(23)

For real-world sensor array data, the noise covariance matrices are unknown and need to be estimated. For the purpose of noise covariance matrix estimation at very low SNR (signals of interest are weak), we may approximate the power of the incoming plane waves as zero. We thus approximate the true noise covariance matrix of the th measurement vector by the sample itself [20], i.e.,

n (ω ) ≈ diag(|y1,l (ω )|2 , . . . , |yN,l (ω )|2 ).

y (ω ) − a(θ , ω )xˆ,LAD (θ , ω )1

λ2 − LN ln . 2π

3

(21)

Finally, all negative single-frequency log-likelihoods are summed (6) and the numerical solution to (7) is found by a global grid search for θ ∈ . The generalization to facilitate multiple DOA estimation is found in [15] for the single-frequency case. 3. Relation to approximate pre-whitening Instead of (9), we may choose alternative weighting matrices W,ω (θ ). This amounts to selecting specific loss functions, other than Laplace’s (6), see [6, Secs. 2.5.2–2.5.3 on pp. 49–53]. The trivial choice W ,ω (θ ) = I N corresponds to the familiar quadratic loss function associated with maximum-likelihood estimation of plane wave amplitudes hidden in additive white Gaussian noise. Popular choices for robust signal processing are, e.g., Huber’s loss function [6, Eq. (2.54] and Tukey’s biweight loss function [6, Eq. (2.55)]. In some situations, many of the important features of a signal in frequency domain are preserved if only the phase is retained. For signals of finite length phase information alone is sufficient to completely reconstruct a signal to within a scale factor [17]. The purpose of this section is to motivate the empirical evidence [18,19] that phase-only beamforming (30) processing with conventional beamforming (CBF) provides improved DOA estimates at low SNR compared to using both amplitude and phase of the array data. Such improvements in DOA estimation are observed when sources are not closely spaced.

(24)

The corresponding approximate pre-whitening matrices are defined in (12) leading to (16). Empirically, it has been found that applying pre-whitening to the data yl only (the dictionary A is non-whitened) is effective in finding the strongest DOA [18,19]. 4. Single-frequency synthetic data simulations Performance is assessed by numerical simulations using synthetic single-frequency data. For the results in Figs. 1 and 2, a single plane wave with DOA −45◦ with additive noise is observed with a uniform linear antenna array with N = 20 elements and spacing λ/2. The feasible set  in (7) is {0◦, 1◦, . . . , 180◦ }. Three types of zero-mean circularly symmetric complex-valued noise n (ω) in (1) are simulated. The first type is i.i.d. Laplace-like noise with probability density (2). Simulations with Laplace-like noise strictly follow the model in Section 2 and allow to assess the gap in peformance between the LAD-based DOA estimator (which minimizes (21)) and its low-cost approximations which maximize (17) or (30). The second type is i.i.d. Gaussian noise with constant parameter σ . Simulations with Gaussian noise serve as a benchmark to compare the gap in performance between the CBF (which arrises as the maximum likelihood DOA estimator for a single plane wave) and the other DOA estimators developed in this work. The third type is heteroscedastic Gaussian noise, cf. [20]. The heteroscedastic noise matrix N = [n1 (ω ), . . . , nL (ω )] is conditionally Gaussian given the parameters σ n for n = 1, . . . , N and  = 1, . . . , L and conditionally independent across sensors and measurement vectors,

p(N |) =

L  N 

1

(2π σk2 ) =1 k=1 ⎡ ⎤ σ11 · · · σ1L .. ⎦. ..  = ⎣ .. . .

σN1

2 2 e−|nk | /σn ,

(25)

(26)

.

···

σNL

The parameter matrix  is an i.i.d. random matrix,

σn = σ sn with sn =  1 NL

tn = 10Un ,

tn N L n=1

,

(27)

2 =1 tn

(28)

where Un is uniformly distributed on [−1, +1]. Simulations with heteroscedastic Gaussian noise allow the comparison with previous results for sparse Bayesian learning in [20]. The deterministic

4

C.F. Mecklenbräuker, P. Gerstoft and E. Zöchmann et al. / Signal Processing 166 (2020) 107262

Fig. 1. RMSE of DOA estimate vs. SNR for a single source at DOA −45◦ (20 element uniform linear antenna array, single measurement vector, L = 1) hidden in additive noise, (a) Laplace-like, (b) Gaussian, (c) heteroscedastic Gaussian.

parameter σ is used to define the SNR in the plots. Fig. 1 shows the performance in terms of the root mean squared error (RMSE) of several DOA estimators for a single plane wave hidden in additive noise. These RMSE results are obtained for the DOA estimators based a single measurement vector, L = 1, as derived in Section 2. When there is just noise in the array data, the investigated DOA estimators become uniformly distributed on [0◦ , 180◦ ] and their stan◦ √ dard deviation becomes 180 ≈ 52◦ . Red curve: RMSE of DOA esti12

mates by maximizing the conventional beamformer (CBF),

BL ( θ , ω ) =

L 2 1   H a ( θ , ω )y ( ω ) L =1





L 1 = aH ( θ , ω ) y  ( ω )y H  ( ω ) a ( θ , ω ). L =1

(29)

Fig. 2. RMSE of DOA estimate vs. SNR for a single source at DOA −45◦ (20 element uniform linear antenna array, multiple measurement vector, L = 4) hidden in additive noise, (a) Laplace-like, (b) Gaussian, (c) heteroscedastic Gaussian. Weighted phase-only beamformer result is for the definition with weighting (17). Phase-only beamformer result is for definition (30).

Purple curve: DOA estimator which maximizes the weighted phase-only beamformer (17), yellow curve with ’◦’-markers: LAD estimator (7), dashed purple curve: sparse Bayesian learning SBL3 [20]. The results for Laplace-like and Gaussian noise in Figs. 1(a) and (b) are all very close to each other they all reach the asymptotic regime above the threshold SNR ≈ 15 dB. For heteroscedastic noise, however, the difference in estimator performance is clearly seen in Fig. 1(c) and the LAD-based DOA estimator has the lowest RMSE over the whole SNR range. The weighted phase-only beamformer approaches the LAD-based DOA estimator for low SNR. Fig. 2 shows RMSE results for several DOA estimators using L = 4 measurement vectors as derived in Section 2. For the CBF, the RMSE of the DOA estimate compared to single measurement vector shows a shift in SNR of 6 dB compared to Fig. 1. The corresponding shifts in SNR for weighted phase-only beamformer (17) and

C.F. Mecklenbräuker, P. Gerstoft and E. Zöchmann et al. / Signal Processing 166 (2020) 107262

5

LAD estimator (7) are either similar or higher than 6 dB. Figs. 1c and 2 c also reveal that CBF and SBL3 suffer from heteroscedastic noise at low SNR in contrast to the weighted phase-only beamformer (17) and LAD-based DOA estimator (21). The phase-only beamformer in Fig. 2 is the weighted phaseonly beamformer made more robust by setting c = 1 for all  = 1, . . . , L in (18). Thus, the phase-only beamformer is defined as

B˜˜L (θ , ω ) =

L (1 ) 1  |xˆ (ω )|2 = aH (θ , ω )S˜ y˜ (ω )a(θ , ω ), L (c (ω ))2

(30)

=1

where we have introduced the sample covariance matrix of the phase-only array data vector y˜  , i.e.,

S˜ y˜ (ω ) =

L 1 H y˜  (ω )y˜  (ω ). L

Fig. 3. Empirical influence function versus sample size L: red: conventional beamformer (29), blue: phase-only beamformer (30). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(31)

=1

The complex sign covariance matrix is closely related to, but different from both, the sign covariance [21] and the spatial sign covariance matrix [6, Eq. (5.6)]. The key difference is in the chosen magnitude normalization of the array data before calculation of the sample covariance. Note that the DOA estimates from the weighted phase-only beamformer (17) are identical to the DOA estimates from the phase-only beamformer when only a single measurement vector is used. In this case L = 1 and the RMSE of the weighted phase-only beamformer (17) in Fig. 1 is identical to the RMSE of the phaseonly beamformer (30). 5. Investigation of qualitative robustness Qualitative robustness can be investigated via the influence function (IF), cf. [6,16,22]. The IF of an estimator T( · ) measures the sensitivity of T(F) to a slight change in the distribution F in the neighborhood of the nominal distribution F0 of the observations. The IF is defined as

IF(z; T, F0 ) = lim

ε →0

T(Fε ) − T(F0 )

ε

,

(32)

where Fε = (1 − ε )F0 + ε G is the ε -contaminated distribution. The contaminating distribution G is defined as the nominal distribution F0 with an additive outlier z ∈ CN in the array data. Thus, (32) measures the effect of an infinitesimal additive contamination z on the estimator T( · ), standardized by the mass of the contamination. A qualitatively robust estimator is characterized by an IF that is continuous and bounded. Continuity implies that small changes in the observed sample cause only small changes in the estimate. The boundedness implies that a small amount of contamination cannot lead to an unbounded error the estimate. Let δ z be the complex point-mass distribution function at z ∈ CN . Hence δ z is chosen as a random point on the N-dimensional real-valued unit-hypersphere scaled by a complex weight rejφ , where φ is uniformly distributed in [0, 2π ]. The magnitude r is denoted as outlier radius in the sequel. All elements of the outlier z have the same phase angle which is similar to a plane wave arriving from broadside at a uniform linear array. Ollila and Koivunen have analytically derived the IF for the MVDR beamformer [23]. Sharif et al. have analytically derived the IF for time-frequency distributions [24]. Furthermore, both compared the IF to the empirical influence functions (EIF), which is a consistent estimator in most cases [16]. We investigate the qualitative robustness of DOA estimators T( · ) via a minor modification of (32) because the evaluation of angular differences should take the circular nature of angles into account. Here, we evaluate the angular difference between two an-

gles θ 1 , θ 2 by

  a θ1 − θ2 = min |θ1 − θ2 | , 360◦ − |θ1 − θ2 |   = 180◦ −  |θ1 − θ2 | −180◦ ,

(33)

which is valid for |θ1 − θ2 | ≤ 360◦ . This definition [25] avoids modulo operations, is non-negative, and bounded by 180◦ . The EIF for DOA estimates is here defined as1

EIF(z; θˆ (· ), Y ) = Eδz

  a   θˆ Y 1:L − θˆ Y 1:L−1 1 L

  a   = L Eδz θˆ Y 1:L − θˆ Y 1:L−1 ,

(34)

where θˆ (· ) is obtained from one of the DOA estimators defined in the previous sections. The multiple measurement vectors are gathered in the matrix Y, where the last column of Y, i.e., yL , contains a contaminated observation. The contaminated observation is sampled from the distribution G which is the nominal distribution F0 with an additive outlier z. The outlier distribution is the complex point-mass δ z distributed uniformly on the hypersphere of radius r. The expectation Eδz {·} is evaluated with respect to the distribution of the outlier z. The EIF according to (34) is the expected angular difference between the DOA estimate for the contaminated sample Y1: L and the DOA estimate for the sample from the nominal distribution Y 1:L−1 . First, we investigate the behavior of the EIF (34) without outlier, i.e. r = 0, versus L and SNR=6 dB. Due to the noise, the EIF is non-zero. The EIF for CBF and the phase-only beamformer are estimated by averaging over 105 realizations and shown in Fig. 3. The EIF is evaluated as an angular difference in degrees per contamination level  = 1L which explains the unusual unit “◦ × Sample Size” in Fig. 3. For small L < 10, the EIF decreases rapidly. This corresponds to the non-asymptotic region of the CBF and phase-only beamformer at low SNR. We see that the EIF reaches a limiting value for large sample size for L > 20, where the CBF approaches the Cramér Rao Lower Bound. Next, we choose to evaluate the EIF at the nominal distribution F0 defined by a single plane wave with additive white Gaussian noise. The SNR of the plane wave signal and the additive white Gaussian noise is 6 dB in a single measurement vector. The DOA estimates are based on L = 100 measurement vectors and the expectation operator in the EIF definition (34) is estimated by averaging over 10 0 0 realizations. The EIF shows L times the expected bias of the contaminated estimator, see Eq. (34). The IF, estimated by the EIF, is independent of the number of snapshots, however. The larger L, the smaller the bias such that the EIF stays constant.

6

C.F. Mecklenbräuker, P. Gerstoft and E. Zöchmann et al. / Signal Processing 166 (2020) 107262

Fig. 4. Empirical influence function for various DOA estimators for L = 100: red: conventional beamformer (29), blue: phase-only beamformer (30), purple: weighted phase-only beamformer (17), yellow: LAD-based DOA estimate (21). Ratio of the outlier sum power to the signal sum power (35) is shown as second x-axis. A zoom into the dashed rectangular box of a) is shown in b). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Figs. 4 (a) and (b) display the EIF versus outlier radius r. For Figs. 4(a) and (b) a second x-axis is plotted which shows the ratio of the outlier sum power to the signal sum power, that is

r2

 = L .  =1 a(θ )x 22

(35)

The outlier shows a strong influence on the conventional beamformer (CBF), see the red curve in Fig. 4(a) according to (29). At outlier radius 80, the EIF of the CBF shows a discontinuity. An outlier shows no influence on the phase-only beamformer (blue curve) according to (30). The qualitative robustness of the LADbased DOA estimator (yellow curve, (21)) is indistinguishable from the weighted phase-only beamformer (purple curve, (17)) since both results are plotted on top of each other. An enlarged view of the dashed region in Fig. 4(a) is shown in Fig. 4(b). It is clearly seen that the phase-only beamformer (30) has a lower EIF compared to all other DOA estimators for outlier radius greater than approx. 10. For very small outlier radius (less than approx. 10) the CBF has the smallest EIF. This is expected as the CBF is optimal for a single DOA in Gaussian noise. 6. Experimental multi-frequency data Here, the weighted phase-only beamformer, phase-only beamformer, and LAD-based DOA estimator are applied to real-world

1

Notation: Y k:l = [yk , . . . , yl ].

measurement data. The simulations have shown that the phaseonly beamformer, weighted phase-only beamformer, and LADbased DOA estimator perform very close to the CBF which is the ML-estimator for a single plane wave in additive Gaussian noise. We therefore expect that the phase-only beamformer, weighted phase-only beamformer, and LAD-based DOA estimator will also perform well for the real-world array data due to asymptotic normality of short-time Fourier transformed array data. However, they should prove to be more robust against additive outliers. The data were collected by Krupp Atlas Elektronik GmbH (Bremen, Germany) in the Bornholm Deep (Baltic Sea) east of Bornholm island with a towed uniform line array (ULA) of 16 hydrophones, dn = (n − 1 ) 2.56 m during October 3–13, 1983. Ocean depth is approximately 50 m and we assume a constant sound speed c = 1475 m/s, [26]. The ULA was towed by the ship “MS Walther von Ledebur.” Shipping noise from the ship “MS Hans Bürkner” served as a maneuvering wideband acoustic source during the experiment. This spacing is larger than the correlation length of the noise, so the noise output from different sensor outputs can be considered uncorrelated. From 16 sensor outputs only the first N = 15 have been at our disposal. We use a sequence of 600 s in which a scenario of multiple broadband sources of shipping noise are present: MS Walther von Ledebur, MS Hans Bürkner, and at least two additional ships. The sequence is recorded at a sample rate of 1024 Hz after low pass filtering with cut off frequency 300 Hz, roll-off of 48 dB/decade, analog-to-digital converted with 12 bit, stored as 16 bit integer format, and partitioned into 150 data pieces of duration T = 4 s each. The data from the 8th hydrophone in the ULA shows 20 dB less power than the others and its statistics deviate significantly from a Gaussian. The phase information of 8th hydrophone, however, was shown to be coherent with the rest of the hydrophones [27]. Further details on data recording, preprocessing, and data conversion are in [28–30]. The processed frequency band up to 250 Hz contains most of the shipping noise [31]. All DOA estimators in this work are based on the a priori assumption that only one plane wave is present in the array data, even though it is known that this simplifying assumption does not match the scenario of four broadband sources. The SNR, as evaluated at a single frequency bin, is too low to allow an acceptable DOA estimate. Therefore, wideband processing is required. We use the set of angular frequencies  = {2π 12, 2π 16, . . . , 2π 252} Hz and the number of processed frequency bins is J = || = 62. The LAD-based DOA estimator (7) sums over all frequency bins. Starting from the singlefrequency CBF (29), the corresponding multi-frequency CBF is here defined as

BL ( θ ) =



BL ( ω , θ ),

(36)

ω ∈

and similarly for the phase-only (30) and weighted phase-only (17) beamformers. Any DOA estimate evaluated for data piece t ∈ {1, . . . , 150} is computed from L = 16 measurement vectors at each frequency ω ∈ . Fig. 5 shows the obtained DOA estimates from multi-frequency processing of the data. Time is shown on the vertical axis as time steps of 4 s durations. The DOA corresponding to the maneuvering ship “MS Hans Bürkner” starts around 105◦ (near broadside of the hydrophone array) and then decreases over time during the first 360 s of the data recording (time steps 0, . . . , 90). During the first 90 time steps, this plane wave arrival is the strongest. Then DOAs appear around 150◦ during time steps 90, . . . , 140 from a distant ship. The DOA estimates are all very similar to each other. This indicates that the investigated DOA estimators perform very similarly to each other.

C.F. Mecklenbräuker, P. Gerstoft and E. Zöchmann et al. / Signal Processing 166 (2020) 107262

7

Fig. 5. DOA estimates for Baltic Sea data (no contamination). All J = 62 frequencies are used. All DOA estimators provide very similar estimates, indicating high SNR gain from wideband processing.

Fig. 6. DOA estimates for Baltic Sea data contaminated by an additive outlier. All J = 62 frequencies are used. DOA estimates by wideband CBF (black dots) and by phase-only beamformer (blue × ). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

8

C.F. Mecklenbräuker, P. Gerstoft and E. Zöchmann et al. / Signal Processing 166 (2020) 107262

Fig. 7. DOA estimates for Baltic Sea data contaminated by an additive outlier. All J = 62 frequency bins are used. DOA estimates by wideband phase-only (red ◦) and by LAD-based estimator (green ∗ ). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Next, we deliberately contaminate the Baltic Sea data set with moderately strong additive outliers at every Lth measurement vector as follows: We simulate a “humm” in the data at angular frequency ω0 = 2π · 100 Hz which is fully coherent over all pairs of sensors. Such humm is commonly caused, e.g., by poorly stabilized AC power supplies and ground loops in the analog parts of the processing chain. The contaminated data are then



yˇ  (ω0 ) =

y  ( ω0 ) y L ( ω0 ) +

√r 1N N

for all  = 1, . . . , (L − 1 ) for  = L,

where y (ω0 ) is the non-contaminated measurement vector, 1N is the all-ones vector, and r is the outlier radius. Fig. 6 shows the DOA estimates from the CBF and phaseonly beamformer from the contaminated data yˇ  .√The chosen contamination is defined by outlier radius r = 1.4 N · yL (ω0 )2 ≈ 5.42yL (ω0 )2 . This outlier radius leads to an outlier to signal sum power ratio (35) of approx. 300%, which is close to the EIF’s discontinuity for the CBF. CBF fails to estimate the DOA because it misinterpretes the contaminating humm as a plane wave arriving from broadside at DOA 90◦ . The phase-only beamformer, however, is unaffected. This is well in agreement with the EIF results from Section 5. Fig. 7 shows the DOA estimates from the weighted phase-only beamformer (red ◦) and the LAD-based estimator (green ∗ ) from the contaminated data. Both estimators suffer severely from the contamination and the DOA estimates are strongly influenced. 7. Conclusion The phase-only beamformer is a low-cost approximation to a DOA estimator based on the Least Absolute Deviation (LAD) criterion. The phase-only beamformer performs well in terms of Root

Mean Squared Error (RMSE) of its DOA estimate for a single DOA at low Signal to Noise Ratio (SNR) for all three investigated noise types. We compare the RMSE of several different DOA estimators versus SNR in simulations: the Conventional Beamformer (CBF), phase-only beamformer, and LAD-based estimator. The simulations indicate that the LAD-based DOA estimator and phase-only beamformer cope well with heteroscedastic noise. The qualitative robustness of the phase-only beamformer, as shown with the Empirial Influence Function (EIF), is also visible in the multi-frequency Baltic Sea data with injected additive outliers. The other three DOA estimates (CBF, weighted phase-only beamformer, LAD-based DOA estimator) suffer from contaminated data. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment The authors thank ATLAS ELEKTRONIK GmbH, Bremen, Germany, for access to the measurements of the Baltic Sea experiment and permission to publish the results. We thank D. Kraus (Hochschule Bremen) for valuable support in this matter. The work at the Inst. of Telecommunications is partially supported by the Christian Doppler Lab “Dependable Wireless Connectivity for the Society in Motion” and the doctoral school “5G Internet of Things”. P. Gerstoft is supported by the Office of Naval Research, Grant No. N0 0 014-18-1-2118.

C.F. Mecklenbräuker, P. Gerstoft and E. Zöchmann et al. / Signal Processing 166 (2020) 107262

References [1] U. Sandkühler, J.F. Böhme, Accuracy of maximum-likelihood estimates for array processing, in: ICASSP ’87. IEEE Int. Conf. on Acoust., Speech, and Sig. Proc., volume 12, 1987, pp. 2015–2018, doi:10.1109/ICASSP.1987.1169416. [2] J. Böhme, Array processing, in: S. Haykin (Ed.), Advances in Spectrum Analysis and Array Processing, Prentice-Hall, Englewood Cliffs, 1991, pp. 1–63. [3] D. Kraus, D. Maiwald, J. Böhme, Maximum likelihood source location estimation via em algorithm, in: Proc. 6th European Signal Processing Conf., Signal Processing VI: Theories and Applications, Brussels, 1992, pp. 649–652. [4] J. Böhme, Statistical array signal processing of measured sonar and seismic data, in: Proc. SPIE 2563 Advanced Signal Processing Algorithms, San Diego, 1995. [5] D. Brillinger, Time series: Data analysis and theory, International Series in Decision Processes, 1st, Holt, Rinehart and Winston, Inc., San Francisco, 1981. [6] A.M. Zoubir, V. Koivunen, E. Ollila, M. Muma, Robust Statistics for Signal Processing, Cambridge University Press, 2018. [7] D. Middleton, A.D. Spaulding, Elements of weak signal detection in non-Gaussian noise, advances in statistical signal processing, in: Signal Detection, volume 2, 1993, pp. 137–215. [8] K.L. Blackard, T.S. Rappaport, C.W. Bostian, Measurements and models of radio frequency noise for indoor wireless communications, IEEE J. Sel. Areas Commun. 11 (7) (1993) 991–1001. [9] S. Jiang, N.C. Beaulieu, Precise ber computation for binary data detection in bandlimited white laplace noise, IEEE Trans. Commun. 59 (6) (2011) 1570–1579. [10] S. Kotz, T.J. Kozubowski, K. Podgorski, The Laplace distribution and generalizations, A revisit with applications to communications, economics, engineering, and finance, Birkhäuser, 2001. [11] T. Eltoft, T. Kim, T. Lee, On the multivariate Laplace distribution, IEEE Sig. Proc. Lett. 13 (5) (2006) 300–303. [12] C.F. Mecklenbräuker, P. Gerstoft, A. Panahi, M. Viberg, Sequential Bayesian sparse signal reconstruction using array data, in: IEEE Trans. Signal Process, volume 61, 2013, pp. 6344–6354. [13] A. Xenaki, P. Gerstoft, K. Mosegaard, Compressive beamforming, J. Acoust. Soc. Am. 136 (1) (2014) 260–271. [14] C.F. Mecklenbräuker, P. Gerstoft, E. Zöchmann, c–LASSO and its dual for sparse signal estimation from array data, Signal Process. 130 (2017) 204–216, doi:10. 1016/j.sigpro.2016.06.029. [15] C.F. Mecklenbräuker, P. Gerstoft, Maximum-likelihood DOA estimation at low SNR in Laplace-like noise, EUSIPCO, A Coruna, Spain, 2019.

9

[16] F.R. Hampel, E.M. Ronchetti, P.J. Rousseeuw, W.A. Stahel, Robust Statistics: The Approach Based on Influence Functions, 196, John Wiley & Sons, 2011. [17] A.V. Oppenheim, J.S. Lim, The importance of phase in signals, Proc. IEEE 69 (5) (1981) 529–541, doi:10.1109/PROC.1981.12022. [18] P. Gerstoft, M.C. Fehler, K.G. Sabra, When katrina hit california, Geophys. Res. Lett. 33 (17) (2006) L17308. [19] Y. Dorfan, S. Gannot, Tree-based recursive expectation-maximization algorithm for localization of acoustic sources, IEEE Trans. Audio, Speech Lang. Process. 23 (10) (2015) 1692–1703. [20] P. Gerstoft, S. Nannuru, C.F. Mecklenbräuker, G. Leus, DOA estimation in heteroscedastic noise, Signal Process. 161 (2019) 63–73, doi:10.1016/j.sigpro.2019. 03.014. [21] E. Ollila, H. Oja, T.P. Hettmansperger, Estimates of regression coefficients based on the sign covariance matrix, J. R. Stat. Soc. Ser. B (Stat. Methodol.) 64 (3) (2002) 447–466, doi:10.1111/1467-9868.00344. [22] E. Ollila, V. Koivunen, Influence functions for array covariance matrix estimators, in: IEEE Workshop on Statistical Signal Processing, volume 20 03, 20 03, pp. 462–465, doi:10.1109/SSP.2003.1289447. [23] E. Ollila, V. Koivunen, Influence function and asymptotic efficiency of scatter matrix based array processors: case MVDR beamformer, in: IEEE Transactions on Signal Processing, volume 57, 2009, pp. 247–259, doi:10.1109/TSP. 20 08.20 07347. [24] W. Sharif, M. Muma, A.M. Zoubir, Robustness analysis of spatial time-frequency distributions based on the influence function, IEEE Trans. Signal Process. 61 (8) (2013) 1958–1971. [25] K.V. Mardia, P.E. Jupp, Directional Statistics, 494, John Wiley & Sons, 20 0 0. [26] T. Seifert, B. Kayser, A high resolution spherical grid topography of the Baltic Sea, Meereswissenschaftlicher Bericht, volume 9, Institut für Ostseeforschung, Warnemünde, 1995. [27] C. Bangert, Analyse und Vorverarbeitung realer Sensorgruppendaten für die Ortsparameterschätzung, Ruhr-Universität Bochum, 1992 Master’s thesis. [28] C. Gierull, Maximum-Likelihood-Schätzung von Wellenparametern mit dem EM-Algorithmus für simulierte und reale Sensorgruppendaten, Ruhr-Universität Bochum, 1990 Master’s thesis. [29] D. Kraus, Approximative Maximum-Likelihood-Schätzung und verwandte Verfahren zur Ortung und Signalschätzung mit Sensorgruppen, Ruhr-Universität Bochum, 1993 Ph.d thesis. [30] D. Maiwald, Breitbandverfahren zur Signalentdeckung und –ortung mit Sensorgruppen in Seismik– und Sonaranwendungen, Ruhr-Universität Bochum, 1995 Ph.d. thesis. [31] R. Urick, Principles of Underwater Sound, 2nd, McGraw-Hill, 1967.