Next Article in Journal
Additive Manufacturing of a Miniaturized X-Band Single-Ridge Waveguide Magic-T for Monopulse Radar Applications
Next Article in Special Issue
Hybrid T-Shaped Sensor Array Composed of Acoustic Vector Sensors and Scalar Sensors
Previous Article in Journal
Space-Borne System-in-Package Based on High Reliability Microwave Interconnections
Previous Article in Special Issue
A Target Detection Method of Distributed Passive Radar without Direct-Path Signal
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Wideband Direction-of-Arrival Estimation Based on Hierarchical Sparse Bayesian Learning for Signals with the Same or Different Frequency Bands

1
School of Marine Science and Technology, Northwestern Polytechnical University, Xi’an 710072, China
2
Shaanxi Key Laboratory of Underwater Information Technology, Xi’an 710072, China
*
Author to whom correspondence should be addressed.
Electronics 2023, 12(5), 1123; https://doi.org/10.3390/electronics12051123
Submission received: 7 January 2023 / Revised: 14 February 2023 / Accepted: 21 February 2023 / Published: 25 February 2023
(This article belongs to the Special Issue Advances in Array Signal Processing)

Abstract

:
Wideband sparse Bayesian learning (WSBL) based on joint sparsity achieves high direction-of-arrival (DOA) estimation precision when the signals share the same frequency band. However, when the signal frequency bands are non-overlapped or partially overlapped, i.e., the frequency bands are different, the performance of the method degrades due to the improper prior on signal. This paper aims at extending the WSBL to a more general version, which is also suitable for the cases where the signal frequency bands are non-overlapped or partially overlapped. Given that the signals are sparsely distributed in the space, the signal matrix whose column is composed of the signal in each frequency bin is row-sparse. Moreover, the signal vectors in some frequency bins have different sparse supports when the signals occupy the different frequency bands. Therefore, a hierarchical sparse prior is assigned to the signal matrix, where a set of hyperparameters are used to ensure the row-sparsity and the other set are used to adjust the signal sparsity in each frequency bin. The DOAs are finally estimated in the Bayesian framework. The simulation results verify that the proposed method achieves good performance on estimation precision in both the same and different frequency band scenarios.

1. Introduction

Direction-of-arrival (DOA) estimation is an important research subject in sonar array signal processing [1,2]. Recently, sparse reconstruction methods [3,4,5,6,7,8,9] have attracted considerable attention. In comparison with traditional high-resolution methods [10,11,12], they can be applied in the lower signal-to-noise ratio (SNR) cases [3]. This kind of method divides the space into the discretized grid points and the steering vectors of the grid points form the dictionary matrix. Then, a linear combination among these bases is selected to fit the received data. These methods can be roughly classified into two categories: lp-norm-based method [3,4,5] and sparse Bayesian learning (SBL) method [13,14]. The lp-norm-based method needs to well tune a regularization parameter to control the tradeoff between the sparsity levels and the data fitting error. Different value of this parameter leads to different estimates [3]. However, the optimal parameter is hard obtained [15,16]. In contrast, the SBL method can automatically tune the hyperparameters from the received data. Therefore, it can be easily applied in practical signal processing.
Nowadays, wideband signal processing is widely used in many applications [17]. A direct way to achieve wideband DOA estimation is to transform the data into frequency domain through Fourier transformation. As such, the signal in each frequency bin is narrowband signal. The narrowband sparse reconstruction methods can be directly applied in each frequency bin, and the estimation results are then summed to obtain the wideband spatial spectrum, which is an independent processing (IP) among frequency bins [5]. The estimation precision of incoherent processing is seriously affected by a large estimation bias in a frequency bin. Hence, it needs high SNR to maintain good performance [18]. To coherently process the data in each frequency bin, the focusing matrix [18,19,20] is used to transform the signal subspace in each frequency bin into the focusing frequency bin. In this manner, the wideband DOA estimation problem is transformed into the narrowband one at the focusing frequency. However, the DOA information is needed as a prior in this kind of method to design the focusing matrix, which is known as the focusing angles. A large bias between the focusing angles and true DOAs leads to inaccurate DOA estimates [21].
To avoid above problems, He et al. [22] proposed a wideband sparse spectrum fitting (W-SpSF) method. It expresses the received data in all frequency bins into a matrix, and the signal vectors in all frequency bins are also written into a matrix, which forms a signal matrix. Clearly, this matrix is a row-sparse matrix since the signal vectors are sparse. As such, this method uses l2,1-norm to the signal matrix as the sparse penalty term to control the row-sparsity, thus jointly processing the data in all frequency bins. However, the W-SpSF is a kind of lp-norm-based method. As mentioned above, the optimal regularization parameter is hard selected [15,16]. Moreover, the method suffers from low computational efficiency since the dimension of the dictionary matrix increases by expressing all frequency bins into a matrix. Das et al. [23] proposed the wideband SBL (WSBL) method. This method assigns the same prior to the signal vectors in all frequency bins, where the common hyperparameters control the signal sparsity in each frequency bin. In this manner, the sparsity of the estimate is determined by all frequency bins jointly, and the estimated signal matrix must be a row-sparse matrix. This thought is also used in other works [24,25,26]. However, this kind of method only considers the condition where the signals occupy the same frequency band and the signal vectors share the same sparse support in all frequency bins. It does not consider the condition where the signals from different angles have different frequency bands and the sparse supports in different frequency bins are not totally the same. Under such conditions, the performances of these methods are obviously degraded.
This paper extends the WSBL to a more general version, which is also suitable for the case where the signal frequency bands are different. As mentioned above, the signal matrix is a row-sparse matrix. The proposed method firstly assigns a set of hyperparameters to control the row-sparsity of the signal matrix. Then, considering that the sparse supports of the signal vectors in different frequency bins are not totally the same if the signal frequency bands are different, other hyperparameters are used to adjust the signal sparsity in each frequency bin. In this manner, a hierarchical sparse prior is assigned to the signal matrix and a Bayesian probabilistic model is established. The variational Bayesian inference (VBI) [27,28] is used to estimate the hyperparameters in this model, thus achieving wideband DOA estimation. The proposed method can be regarded as a combination of the IP and joint processing (JP). Numerical simulation results verify that the proposed method achieves high DOA estimation precision under the conditions where the signal frequency bands are non-overlapped, partially overlapped and common.
The following notations are used throughout this paper: ( · ) H , ( · ) T and ( · ) * represent the conjugate transpose, transpose and conjugate, respectively; t r { · } and denote the trace operator and Hadamard product, respectively; · 2 and | · | denote the l2-norm and modulus, respectively; d i a g ( x ) is a diagonal matrix with entries in vector x being its diagonal elements; d i a g ( X ) represents a vector composed of the diagonal entries of matrix X. X 1 means the inverse of the matrix X r o u n d ( x ) rounds x to the nearest integer; I M is an M × M identity matrix; 1 Q × 1 and 0 Q are Q × 1 vectors whose entries are all 1 and 0, respectively.

2. Signal Model Establishment

Assume that K signals impinge on an M-element array from θ = [ θ 1 , , θ K ] T . The sampling frequency is f s . The received data are divided into N segments and are transformed into frequency domain in each segment through L ˜ -point discrete Fourier transformation (DFT). In this manner, the frequency resolution is Δ f = f s / L ˜ . Record the analysis frequency range as [ f L ,   f H ] . The array output model x l ( n ) M × 1 in the lth sub-band is established as
x l ( n ) = A l ( θ ) s l ( n ) + e l ( n ) ,             n = 1 , , N ,   l = 1 , , L
where the center frequency of the lth sub-band is f l = r o u n d ( f L / Δ f + l 1 ) Δ f , L = r o u n d [ ( f H f L ) / Δ f ] + 1 , x l ( n ) M × 1 , s l ( n ) K × 1 , and e l ( n ) M × 1 denote the DFT coefficients of the received data, the signals, and the additive noise at the nth segment, respectively. A l ( θ ) = [ a l ( θ 1 ) , , a l ( θ K ) ] , where a l ( θ ) is the steering vector from the direction θ . Clearly, s l k ( n ) = 0 if the frequency f l is not contained in the signal at θ k , where s l k ( n ) is the kth entry of s l ( n ) .
With the whole space divided into the discretized grid θ = [ θ 1 , , θ Q ] T , Equation (1) can be rewritten as
x l ( n ) = A l ( θ ) s l ( n ) + e l ( n ) ,             n = 1 , , N ,   l = 1 , , L ,
where A l ( θ ) = [ a l ( θ 1 ) , , a l ( θ Q ) ] , and s l ( n ) Q × 1 is a zero-padded version of s l ( n ) . For any k = 1 , , K , s l q ( n ) = s l k ( n ) holds if θ q = θ k , where s l q ( n ) is the qth entry of s l ( n ) . Otherwise, s l q ( n ) = 0 .
Figure 1 shows three different cases of frequency bands. In Figure 1a, the signals share the same band. As such, signal vectors s l ( n ) ,   l = 1 , , L share the same sparsity, and the same prior is assigned to them to make full use of common sparsity profile [23]. Figure 1b,c show the non-overlapped bands and partially overlapped bands, respectively. In these cases, the sparse supports of signals in some frequency bins are different. Hence, it is unsuitable to assign the same prior to them. A direct way to solve this problem is to estimate the DOAs in each frequency bin independently. However, if a large estimation bias occurs in a frequency bin, it will seriously affect the final estimation precision.
As shown in Figure 1, the row of [ s 1 ( n ) , , s L ( n ) ] is sparse since the number of grid points is much larger than that of signals, and the sparse supports of s l ( n ) ,   l = 1 , , L are not totally the same if signals occupy different frequency bands. Hence, two sets of hyperparameters are used in this paper: one set control the row-sparsity of the signal matrix, and the other set adjust the signal sparsity in each frequency bin. In this manner, one can ensure that the signal matrix estimate is row-sparse and the sparse supports of estimates of s l ( n ) ,   l = 1 , , L are different when the signals occupy the different frequency bands.

3. Proposed Method

3.1. Bayesian Criteria

(1) Likelihood: Under an assumption of white complex Gaussian noise, the likelihood of x l ( n ) ,             n = 1 , , N ,   l = 1 , , L is written as
l = 1 L n = 1 N p ( x l ( n ) | s l ( n ) ; σ l ) = l = 1 L n = 1 N C N ( A l ( θ ) s l ( n ) , σ l I M ) ,
where C N ( μ , Σ ) represents the complex Gaussian distribution with mean of μ and variance matrix of Σ , σ l represents the noise power.
(2) Prior: As analyzed above, the signal matrix [ s 1 ( n ) , , s L ( n ) ] is a row-sparse matrix, and the sparse supports of s l ( n ) ,   l = 1 , , L are different as shown in Figure 1b,c. As such, a hierarchical sparse prior is assigned to [ s 1 ( n ) , , s L ( n ) ] to control the row-sparsity of this matrix and the sparsity of s l ( n ) ,   l = 1 , , L , as given by
n = 1 N l = 1 L p ( s l ( n ) ; α l , γ ) = n = 1 N l = 1 L C N ( 0 Q , Γ l ) ,
where Γ l = d i a g ( α l γ ) . The hyperparameter γ = [ γ 1 , , γ Q ] T controls the row-sparsity of the signal matrix [ s 1 ( n ) , , s L ( n ) ] . If γ q 0 , the qth row [ s 1 q ( n ) , , s L q ( n ) ] 0 , i.e., there exists no signal at θ q . The other hyperparameter α l = [ α l 1 , , α l Q ] T adjusts the sparsity of s l ( n ) adaptively. If α l q 0 , s l q ( n ) 0 , i.e., frequency f l is not contained in the signal at θ q , which ensures that the sparse supports of signal vector estimates are different when the signals occupy different frequency bands.
The inverse Gamma distributions are assigned to α l and γ , that is
p ( α l | ρ l ) = q = 1 Q p ( α l q | ρ l q ) = q = 1 Q I G ( 1 , ρ l q ) = q = 1 Q ρ l q α l q 2 exp ( ρ l q α l q ) ,   l = 1 , , L ,
p ( γ | ρ ) = q = 1 Q p ( γ q | ρ q ) = q = 1 Q I G ( 1 , ρ q ) = q = 1 Q ρ q γ q 2 exp ( ρ q γ q ) ,
where I G ( · ) represents the inverse Gamma distribution, ρ l = [ ρ l 1 , , ρ l Q ] T and ρ = [ ρ 1 , , ρ Q ] T are the scale parameters of the inverse Gamma distribution. Hence, the prior over s l ( n ) for a given γ is expressed as
n = 1 N l = 1 L p ( s l ( n ) | γ , ρ l ) = n = 1 N l = 1 L p ( s l ( n ) | γ , α l ) p ( α l | ρ l ) d α l = n = 1 N l = 1 L q = 1 Q p ( s l q ( n ) | γ q , α l q ) p ( α l q | ρ l q ) d α l q , = n = 1 N l = 1 L q = 1 Q 2 ρ l q π γ q 1 ( | s l q ( n ) | 2 / γ q + ρ l q ) 2
which is a generalized t-distribution. Similarly, the prior over s l ( n ) for a given α l is also a generalized t-distribution:
n = 1 N l = 1 L p ( s l ( n ) | α l , ρ ) = n = 1 N l = 1 L p ( s l ( n ) | γ , α l ) p ( γ | ρ ) d γ = n = 1 N l = 1 L q = 1 Q p ( s l q ( n ) | γ q , α l q ) p ( γ q | ρ q ) d γ q , = n = 1 N l = 1 L q = 1 Q 2 ρ q π α l q 1 ( | s l q ( n ) | 2 / α l q + ρ q ) 2
As shown in [29], in comparison with complex Gaussian distribution, the generalized t-distribution has a heavier tail, thus promoting the sparsity. Moreover, the tail becomes heavier as the decreasing of ρ l q and ρ q [30].
To adaptively adjust the result sparsity from the received data, the Gamma distributions are assigned to the parameters ρ l ,   l = 1 , , L and ρ , as given by
{ p ( ρ l ; a , b ) = q = 1 Q p ( ρ l q ; a , b ) = q = 1 Q G ( a , b ) = q = 1 Q b a Γ ( a ) ρ l q a 1 exp ( b ρ l q ) p ( ρ ; a , b ) = q = 1 Q p ( ρ q ; a , b ) = q = 1 Q G ( a , b ) = q = 1 Q b a Γ ( a ) ρ q a 1 exp ( b ρ q ) ,
where G ( · ) represents the Gamma distribution, Γ ( a ) = 0 t a 1 exp ( t ) d t is the Gamma function, a and b are known as the shape and rate parameters of the Gamma distribution, respectively. a and b are often set to small values, thus inducing non-informative prior on ρ l ,   l = 1 , , L and ρ [27]. They are set to 10 6 in this paper.
Figure 2 illustrates the directed acyclic graph of hierarchical Bayesian model. It can be observed that s l ( n ) , α l , γ , ρ l and ρ are the latent variables, while σ l is treated as a parameter. Based on Equations (3)~(6) and Equation (9), the posterior distribution is written as follows:
l = 1 L n = 1 N p ( s l ( n ) , α l , γ , ρ l , ρ | x l ( n ) ; σ l , a , b ) = l = 1 L n = 1 N p ( x l ( n ) , s l ( n ) , α l , γ , ρ l , ρ ; σ l , a , b ) p ( x l ( n ) ) = l = 1 L n = 1 N p ( x l ( n ) | s l ( n ) ; σ l ) p ( s l ( n ) | α l , γ ) p ( α l | ρ l ) p ( γ | ρ ) p ( ρ l ; a , b ) p ( ρ ; a , b ) p ( x l ( n ) ) ,
where p ( x l ( n ) , s l ( n ) , α l , γ , ρ l , ρ ; σ l , a , b ) is the joint distribution of all unknown and observed quantities.

3.2. Proposed Method

It is difficult to obtain the posterior p ( s l ( n ) , α l , γ , ρ l , ρ | x l ( n ) ; σ l , a , b ) directly. Hence, the VBI [27,28] is used to approximate Equation (10). Under the assumption that the approximated posterior is factorizable, the approximated posterior distribution with respect to latent variables is formulated as
q ( Θ ) = l = 1 L n = 1 N q ( s l ( n ) ) q ( α l ) q ( γ ) q ( ρ l ) q ( ρ ) .
In this equation, q ( · ) is the approximated posterior, and
Θ = [ s 1 T ( 1 ) , , s 1 T ( N ) , s 2 T ( 1 ) , , s L T ( N ) , α 1 T , , α L T , γ T , ρ 1 T , , ρ L T , ρ T ] T
is a vector containing all latent variables.
The VBI iteratively updates the latent variables in Θ , thus minimizing the KL divergence between the true and the approximated posteriors, i.e., [27,28]
q ^ ( Θ ) = arg min q ( Θ ) K L ( q ( Θ ) l = 1 L n = 1 N p ( s l ( n ) , α l , γ , ρ l , ρ | x l ( n ) ; σ l , a , b ) ) ,
where q ^ ( Θ ) is the estimation result for q ( Θ ) . From Equation (12), one can compute the approximated posterior with respect to the individual variable, as given by [27,28]
ln q ^ ( Θ i ) = ln l = 1 L n = 1 N p ( x l ( n ) , s l ( n ) , α l , γ , ρ l , ρ ; σ l , a , b ) q ^ ( Θ \ Θ i ) .
In this equation, Θ i is the ith entry of Θ , · q ( · ) denotes the expectation operator with respect to q ( · ) , and q ^ ( Θ \ Θ i ) represents the approximated posterior of Θ except Θ i .
In the VBI framework, the posteriors of the entries in Θ are updated alternately while the rest are fixed. Once the approximated posterior q ^ ( Θ i ) is obtained from Equation (13), the mean of Θ i with respect to q ^ ( Θ i ) is regarded as the estimate of Θ i , which is recorded as Θ i . The parameter σ l is estimated via maximizing the expected log-likelihood function since no informative prior is assigned to it. The updating for each unknown parameter in an iteration is provided as follows:
(1) Update for s l ( n ) : As shown in Equation (13), the approximated posterior with respect to s l ( n ) is calculated by
ln q ^ ( s l ( n ) ) ln p ( x l ( n ) | s l ( n ) ; σ l ) p ( s l ( n ) | α l , γ ) q ^ ( α l ) q ^ ( γ ) .
With the substitution of Equations (3) and (4) into Equation (14), the approximated posterior distribution of s l ( n ) follows a complex Gaussian distribution, as given by
q ^ ( s l ( n ) ) = C N ( s l ( n ) , Σ s l ) ,
where
{ Σ s l = ( Γ l 1 + σ ^ l 1 A l H ( θ ) A l ( θ ) ) 1 s l ( n ) = σ ^ l 1 Σ s l A l H ( θ ) x l ( n ) .
In the above equation, σ ^ l is the estimate of σ l which is provided in the following;
(2) Update for α l : The approximated posterior of α l is given by
ln q ^ ( α l ) n = 1 N ln p ( s l ( n ) | α l , γ ) + ln p ( α l | ρ l ) q ^ ( s l ( n ) ) q ^ ( γ ) q ^ ( ρ l ) .
With the substitution of Equations (4) and (5) into Equation (17), the following relation is obtained:
ln q ^ ( α l q ) n = 1 N [ ln α l q + | s l q ( n ) | 2 γ q 1 α l q ] 2 ln α l q ρ l q α l q = ( N + 2 ) ln α l q γ q 1 n = 1 N | s l q ( n ) | 2 + ρ l q α l q ,
which follows an inverse Gamma distribution:
α l q ~ I G ( N + 1 , γ q 1 n = 1 N | s l q ( n ) | 2 + ρ l q ) .
Thus, the mean of α l q is given by
α l q = γ q 1 n = 1 N | s l q ( n ) | 2 + ρ l q N ;
(3) Update for γ : The approximated posterior distribution with respect to γ is given as follows:
ln q ^ ( γ ) l = 1 L n = 1 N ln p ( s l ( n ) | α l , γ ) + ln p ( γ | ρ ) q ( s l ( n ) ) q ( α l ) q ( ρ ) q = 1 Q l = 1 L n = 1 N [ ln γ q + | s l q ( n ) | 2 α l q 1 γ q ] q = 1 Q [ 2 ln γ q + ρ q γ q ] . = q = 1 Q [ ( L N + 2 ) ln γ q + l = 1 L α l q 1 n = 1 N | s l q ( n ) | 2 + ρ q γ q ]
It can be observed that γ q ,   q = 1 , , Q follows an inverse Gamma distribution:
γ q ~ I G ( L N + 1 , l = 1 L α l q 1 n = 1 N | s l q ( n ) | 2 + ρ q ) .
Thus, the mean of γ q is as follows:
γ q = l = 1 L α l q 1 n = 1 N | s l q ( n ) | 2 + ρ q L N ;
(4) Update for ρ l : The approximated posterior of ρ l is obtained by
ln q ^ ( ρ l ) ln p ( α l | ρ l ) + ln p ( ρ l | a , b ) q ( α l ) q = 1 Q [ ln ρ l q + α l q 1 ρ l q ] q = 1 Q [ ( a 1 ) ln ρ l q + b ρ l q ] = q = 1 Q [ a ln ρ l q + ( α l q 1 + b ) ρ l q ]
i.e., ρ l q ~ G ( a + 1 , < α l q 1 > + b ) . The mean of ρ l q is computed by
ρ l q = a + 1 α l q 1 + b ;
(5) Update for ρ : Similarly, the approximated posterior of ρ is obtained by
ln q ^ ( ρ ) ln p ( γ | ρ ) + ln p ( ρ | a , b ) q ( α l ) q = 1 Q [ ln ρ q + γ q 1 ρ q ] q = 1 Q [ ( a 1 ) ln ρ q + b ρ q ] , = q = 1 Q [ a ln ρ q + ( γ q 1 + b ) ρ q ]
which also obeys a Gamma distribution. The mean of ρ q is expressed as
ρ q = a + 1 γ q 1 + b ;
(6) Update for σ l : σ l is updated via maximizing the expected log-likelihood function, that is,
arg max σ l n = 1 N ln p ( x l ( n ) | s l ( n ) ; σ l ) q ( s l ( n ) ) = arg min σ l n = 1 N [ M ln σ l + 1 σ l x l ( n ) A l ( θ ) s l ( n ) 2 2 q ( s l ( n ) ) ] . = arg min σ l   M N   ln σ l + 1 σ l n = 1 N x l ( n ) A l ( θ ) s l ( n ) 2 2 q ( s l ( n ) )
After differentiating Equation (28) with respect to σ l and setting it to zero, the estimate of σ l is obtained by
σ ^ l = 1 M N n = 1 N x l ( n ) A l ( θ ) s l ( n ) 2 2 q ( s l ( n ) ) .
In the above equation, x l ( n ) A l ( θ ) s l ( n ) 2 2 q ( s l ( n ) ) is obtained by
x l ( n ) A l ( θ ) s l ( n ) 2 2 q ( s l ( n ) ) = x l ( n ) A l ( θ ) s l ( n ) 2 2 + t r ( A l H ( θ ) A l ( θ ) Σ s l ) .
The initial values of the iteration are set to α l ( 0 ) = d i a g ( A l H ( θ ) R ^ l A l ( θ ) / M 2 ) , where R ^ l = n = 1 N x l ( n ) x l H ( n ) / N is the sample covariance matrix and superscript (i) represents the ith iteration, γ ( 0 ) = [ l = 1 L a l H ( θ 1 ) R ^ l a l ( θ 1 ) / M 2 , , l = 1 L a l H ( θ Q ) R ^ l a l ( θ Q ) / M 2 ] T , ρ l ( 0 ) =   ρ ( 0 ) = 1 Q × 1 , and σ ^ l ( 0 ) = 1 . The signal power estimate p ^ q at θ q is computed by
p ^ q = 1 N n = 1 N | s l q ( n ) | 2 ,   q = 1 , , Q .
The stop criterion of the iteration is that | | p ^ ( i + 1 ) p ^ ( i ) | | 2 / | | p ^ ( i ) | | 2 τ for a small tolerance τ or the iteration reaches the maximum number I t e r max , where p ^ = [ p ^ 1 , , p ^ Q ] T . The proposed method is summarized in Algorithm 1.
Algorithm 1. Summary of the proposed algorithm
Input: α l ( 0 ) , γ ( 0 ) , ρ l ( 0 ) , ρ ( 0 ) , σ ^ l ( 0 ) ,   l = 1 , , L , and i = 0
Do
Update s l ( n ) ( i + 1 ) and Σ s l ( i + 1 ) using Equation (16).
Update α l ( i + 1 ) and γ ( i + 1 ) via Equations (20) and (23), respectively.
Update ρ l ( i + 1 ) and ρ ( i + 1 ) via Equations (25) and (27), respectively.
Compute σ ^ l ( i + 1 ) and p ^ ( i + 1 ) using Equations (29) and (31), respectively.
i = i + 1
while | | p ^ ( i ) p ^ ( i 1 ) | | 2 / | | p ^ ( i ) | | 2 > τ and i < I t e r max
Output: p ^ ( i )
As shown in Equations (20) and (23), it can be seen that the update of α l ( i + 1 ) only depends on the data in the lth frequency bin while the update of γ ( i + 1 ) joints all frequency bins. Therefore, the proposed method can be regarded as a combination of IP and JP.

4. Simulation Results

Consider a 32-element uniform linear array (ULA) placed under the sea surface in the simulations. In general, the sound speed in the seawater is considered as 1500 m/s [31]. The array element spacing is 4 m (i.e., the half wavelength of 187.5 Hz). Two far-field sound signals (recorded as Signal1 and Signal2) impinge on the array from −10° and −8° with the assumption of plane wave propagation, where ±90° are defined as the endfire direction. Three different frequency cases are considered as follows:
(1)
Case I: Signal1: [90, 134] Hz, Signal2: [136, 180] Hz. The bandwidth of signals is 44 Hz;
(2)
Case II: Signal1: [90, 160] Hz, Signal2: [110, 180] Hz. The bandwidth of signals is 70 Hz;
(3)
Case III: Signal1: [90, 180] Hz, Signal2: [90, 180] Hz. The bandwidth of signals is 90 Hz.
In the first case, the frequency bands of signals are non-overlapped. In the second case, the signals have partially overlapped bands. In the last case, the signals share the same frequency band. The data are sampled by the acquisition system and the sample frequency is 5 kHz. The received data are divided into 30 segments, i.e., N = 30, with the length of 2500. Then, a 2500-point DFT is applied in each segment, i.e., L ˜ = 2500 . The frequency resolution is 2 Hz. The analysis frequency range is from 90 Hz to 180 Hz. The total number of the sub-bands is L = 46.
The performance of the proposed method is compared with the methods in [23,27], and W-SpSF [22]. The method in [27] is a narrowband SBL method. It is used to estimate DOAs in each frequency bin and the results are summed incoherently. Hence, it is recorded as IP-WSBL in the simulation. The method in [23] assigns the same prior in all frequency bins and jointly processes the data in all frequency bins. It is recorded as JP-WSBL. The results of traditional method, i.e., conventional beamforming (CBF) [17], are also shown for comparison purpose. The observed space [−90, 90]° is uniformly discretized into 181 grid points, i.e., Q = 181. The tolerance τ and I t e r max are set to 10 3 and 3000, respectively. All results are performed in MATLAB on a PC with an Intel Core i7-9700F CPU and an 8 GB RAM.
Figure 3 illustrates the space–frequency map estimation results of the above methods under three different cases by fixing the SNR to −5 dB. The color bar represents the normalized signal power. The row of the bright area gives the DOA information, and the corresponding column provides the signal frequency information in this direction. The results of the proposed method and the W-SpSF are similar to the true values, but the proposed method provides cleaner space–frequency maps. The JP-WSBL achieves a comparable result to the proposed method under the Case III. However, it cannot work well under the first two cases since the sparse supports of the signal vectors in some frequency bins are different and the assumption of a common sparse profile is unsuitable under such conditions. The performance of the IP-WSBL degrades when two signals are in the overlapping frequency bands. Furthermore, the CBF cannot resolve two signals in the overlapping frequency band due to the poor resolution.
Figure 4 illustrates the spectra of the above methods under three different cases. It can be seen that the proposed method and the W-SpSF can estimate two signals under all considered conditions. The IP-WSBL can hardly estimate two signals for Case III since the DOA estimation precision at low frequency bins is poor, affecting the final DOA estimation result. The performance of JP-WSBL degrades under Case I because of the unsuitable assumption on signal sparsity in each frequency bin. Furthermore, the CBF cannot resolve two closely-spaced signals under all considered conditions due to poor resolution.
Then, the performances of above methods on estimation precision are compared from the perspective of the root-mean-square error (RMSE):
RMSE = 1 K W k = 1 K w = 1 W ( θ ^ k w θ k ) 2 ,
where θ ^ k w is the estimated DOA of the kth signal in the wth trial, and W is the total of the Monte Carlo runs. The value of W is set to 200. Figure 5 illustrates the RMSE versus the SNR for Case I. The performances of the proposed method, the IP-WSBL and the W-SpSF are similar to each other and can provide reasonable estimates when SNR ≥ −15 dB. However, the W-SpSF needs to tune a regularization parameter well, which is a difficult task. The JP-WSBL cannot work well at low SNR cases since the assumption in this method is unsuitable under such conditions. The RMSE of the CBF is larger than those of other methods as the SNR increases because it cannot resolve two signals.
Figure 6 shows the RMSE versus the SNR under the condition where the signals have partially overlapped bands. As the overlap of the signal bands increases, two signals are more difficult to be resolved. It is observed that the performances of the proposed method and W-SpSF are nearly unchanged in comparison with the results in Figure 5. The IP-WSBL estimates the DOAs in each frequency bin independently, and the estimation results are not satisfying at the low frequency bins. Hence, it needs higher SNR than the proposed method to provide stable estimates. The RMSE of JP-WSBL is still larger than the proposed method since the sparse supports in different frequency bins are not totally the same. The CBF suffers from low resolution, resulting in large RMSEs in all considered situations.
The RMSE versus the SNR under Case III is shown in Figure 7. Once again, the performances of the proposed method and the W-SpSF are similar to each other. Unlike the W-SpSF, the proposed method gets rid of tuning the hyperparameters manually. Given that the signal vectors share the same sparsity in this simulation, the JP-WSBL assigns the same prior to signal vectors in all frequency bins to make full use of common sparse profile, thus processing the data in all frequency bins jointly. In this manner, the JP-WSBL achieves a comparable performance to the proposed method. The estimation precision of IP-WSBL is affected by large estimation errors in some frequency bins. Hence, it needs higher SNR to achieve high estimation precision, and its RMSE approaches that of the proposed method when the SNR is larger than −5 dB. Furthermore, the RMSE of the CBF is larger than those of other methods due to the low resolution.
Finally, the computational efficiency of each method is compared from the perspective of the mean running time over 200 simulations by fixing SNR to −5 dB. Table 1 shows the mean running time of each method under three different frequency cases. The proposed method, JP-WSBL and IP-WSBL iteratively estimate the parameters. The workload in each iteration depends on the updating of Σ s l , which is O ( Q 3 ) . However, the proposed method achieves higher computational efficiency than the other two methods mainly because it requires a smaller number of the iterations to converge. The W-SpSF expresses the data in all frequency bins into a matrix, which increases the dimension of the matrix, and the workload is O ( L 3 Q 3 ) [22]. Hence, it suffers from low computational efficiency. The CBF needs no iteration and avoids matrix inversion. Hence, it achieves the highest computational efficiency, and the workload is O ( Q M 2 ) . However, it suffers from low resolution.

5. Conclusions

This paper modifies the WSBL to improve its performance under the condition where the signals from different angles have different frequency bands. Given that the signal matrix [ s 1 ( n ) , , s L ( n ) ] is a row-sparse matrix, and some columns share the different sparse supports when the signals occupy the different frequency bands, different hyperparameters are applied to control the row-sparsity of the signal matrix and the sparsity of the signal vector in each frequency bin. As such, a hierarchical sparse prior is assigned to the signal matrix. The hyperparameters are estimated through VBI to achieve wideband DOA estimation in the Bayesian framework, avoiding the choice of regularization parameter. The proposed method can be regarded as a combination of the JP-WSBL and IP-WSBL. The numerical simulations are performed under the conditions where the signal bands are the same, partially overlapped and non-overlapped. The simulation results verify that the proposed method achieves better performance than the IP-WSBL when the bands of the signals are overlapped, and achieves better performance than the JP-WSBL under the condition that the signals share the different frequency bands.

Author Contributions

Methodology, Y.Y. and Y.Z.; software, L.Y.; validation, Y.Y. and Y.Z.; formal analysis, Y.Y.; writing–original draft preparation, Y.Y. and Y.Z.; writing–review and editing, L.Y. and Y.W.; supervision, Y.Y. and Y.W. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported in part by the National Natural Science Foundation of China under Grant 11974286 and 61971353, in part by the National Key R&D Program of China under Grant 2021YFB3203001.

Data Availability Statement

Data sharing not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Johnson, D.H.; Degraaf, S.R. Improving the resolution of bearing in passive sonar arrays by eigenvalue analysis. IEEE Trans. Acoust. Speech Signal Process. 1982, 30, 638–647. [Google Scholar] [CrossRef] [Green Version]
  2. Ren, S.; Ge, F.; Guo, X.; Guo, L. Eigenanalysis-based adaptive interference suppression and its application in acoustic source range estimation. IEEE J. Ocean. Eng. 2015, 40, 903–916. [Google Scholar] [CrossRef]
  3. Gerstoft, P.; Xenaki, A.; Mecklenbräuker, C.F. Multiple and single snapshot compressive beamforming. J. Acoust. Soc. Am. 2015, 138, 2003–2014. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Zheng, J.; Kaveh, M. Sparse spatial spectral estimation: A covariance fitting algorithm, performance and regularization. IEEE Trans. Signal Process. 2013, 61, 2767–2777. [Google Scholar] [CrossRef]
  5. Malioutov, D.; Cetin, M.; Willsky, A.S. A sparse signal reconstruction perspective for source localization with sensor arrays. IEEE Trans. Signal Process. 2005, 53, 3010–3022. [Google Scholar] [CrossRef] [Green Version]
  6. Wen, C.; Xie, X.; Shi, G. Off-grid DOA estimation under nonuniform noise via variational sparse Bayesian learning. Signal Process. 2017, 137, 69–79. [Google Scholar] [CrossRef]
  7. Zhang, Y.; Yang, Y.; Yang, L.; Guo, X. Root sparse asymptotic minimum variance for off-grid direction-of-arrival estimation. Signal Process. 2019, 163, 225–231. [Google Scholar] [CrossRef]
  8. Chen, P.; Chen, Z.; Zhang, X.; Liu, L. SBL-based direction finding method with imperfect array. Electronics 2018, 7, 426. [Google Scholar] [CrossRef] [Green Version]
  9. Ling, Y.; Gao, H.; Ru, G.; Chen, H.; Li, B.; Cao, T. Grid reconfiguration method for off-grid DOA estimation. Electronics 2019, 8, 1209. [Google Scholar] [CrossRef] [Green Version]
  10. Schmidt, R.O. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propag. 1986, 34, 276–280. [Google Scholar] [CrossRef] [Green Version]
  11. Roy, R.; Kailath, T. ESPRIT-estimation of signal parameters via rotational invariance techniques. IEEE Trans. Acoust. Speech Signal Process. 1989, 37, 984–995. [Google Scholar] [CrossRef] [Green Version]
  12. Rao, B.D.; Hari, K.V.S. Performance analysis of root-MUSIC. IEEE Trans. Acoust. Speech Signal Process. 1989, 37, 1939–1949. [Google Scholar] [CrossRef]
  13. Tipping, M.E. Sparse Bayesian learning and the relevance vector machine. J. Mach. Learn. Res. 2001, 1, 211–244. [Google Scholar]
  14. Yang, J.; Liao, G.; Li, J. An efficient off-grid DOA estimation approach for nested array signal processing by using sparse Bayesian learning strategies. Signal Process. 2016, 128, 110–122. [Google Scholar] [CrossRef]
  15. Huang, M.; Huang, L. Sparse recovery assisted DOA estimation utilizing sparse Bayesian learning. In Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Calgary, AB, Canada, 15–20 April 2018. [Google Scholar]
  16. Li, S.; Xie, D. Compressed symmetric nested arrays and their application for direction-of-arrival estimation of near-field sources. Sensors 2016, 16, 1939. [Google Scholar] [CrossRef] [Green Version]
  17. Kim, H.; Viberg, M. Two decades of array signal parameter estimation. IEEE Signal Process. Mag. 1996, 13, 67–94. [Google Scholar]
  18. Wang, H.; Kaveh, M. Coherent signal-subspace processing for the detection and estimation of angles of arrival of multiple wideband sources. IEEE Trans. Acoust. Speech Signal Process. 1985, 33, 823–831. [Google Scholar] [CrossRef] [Green Version]
  19. Hung, H.; Kaveh, M. Focussing matrices for coherent signal-subspace processing. IEEE Trans. Acoust. Speech Signal Process. 1988, 36, 1272–1281. [Google Scholar] [CrossRef]
  20. Valaee, S.; Kabal, P. Wideband array processing using a two-sided correlation transformation. IEEE Trans. Signal Process. 1995, 43, 160–172. [Google Scholar] [CrossRef] [Green Version]
  21. Sellone, F. Robust auto-focusing wideband DOA estimation. Signal Process. 2006, 86, 17–37. [Google Scholar] [CrossRef]
  22. He, Z.; Shi, Z.; Huang, L.; So, H.C. Underdetermined DOA estimation for wideband signals using robust sparse covariance fitting. IEEE Signal Process. Lett. 2015, 22, 435–439. [Google Scholar] [CrossRef]
  23. Das, A.; Sejnowski, T.J. Narrowband and wideband off-grid direction-of-arrival estimation via sparse Bayesian learning. IEEE J. Ocean. Eng. 2018, 43, 108–118. [Google Scholar] [CrossRef]
  24. Das, A. Real-valued sparse Bayesian learning for off-grid direction-of-arrival (DOA) estimation in ocean acoustics. IEEE J. Ocean. Eng. 2021, 46, 172–182. [Google Scholar] [CrossRef]
  25. Jiang, Y.; He, M.; Liu, W.; Feng, M. Underdetermined wideband DOA estimation for off-grid targets: A computationally efficient sparse Bayesian learning approach. IET Radar Sonar Navig. 2020, 14, 1583–1591. [Google Scholar] [CrossRef]
  26. Hu, N.; Sun, B.; Zhang, Y.; Dai, J.; Wang, J.; Chang, C. Underdetermined DOA estimation method for wideband signals using joint nonnegative sparse Bayesian leaning. IEEE Signal Process. Lett. 2017, 24, 535–539. [Google Scholar] [CrossRef]
  27. Tzikas, D.G.; Likas, A.C.; Galatsanos, N.P. The variational approximation for Bayesian inference. IEEE Signal Process. Mag. 2008, 25, 131–146. [Google Scholar] [CrossRef]
  28. Themelis, K.E.; Rontogiannis, A.A.; Koutroumbas, K.D. A variational Bayes framework for sparse adaptive estimation. IEEE Trans. Signal Process. 2014, 62, 4723–4736. [Google Scholar] [CrossRef] [Green Version]
  29. Yang, J.; Yang, Y.; Lu, J. A variational Bayesian strategy for solving the DOA estimation problem in sparse array. Digit. Signal Process. 2019, 90, 28–35. [Google Scholar] [CrossRef]
  30. Yang, L.; Fang, J.; Cheng, H.; Li, H. Sparse Bayesian dictionary learning with a Gaussian hierarchical model. Signal Process. 2017, 130, 93–104. [Google Scholar] [CrossRef] [Green Version]
  31. Yang, L.; Hou, X.; Yang, Y. Self-calibration for sparse uniform linear arrays with unknown direction-dependent sensor phase by deploying an individual standard sensor. Electronics 2023, 12, 60. [Google Scholar] [CrossRef]
Figure 1. The signals with (a) the same, (b) non-overlapped and (c) partially overlapped frequency bands. The white and green squares correspond to the absence and presence of a signal, respectively.
Figure 1. The signals with (a) the same, (b) non-overlapped and (c) partially overlapped frequency bands. The white and green squares correspond to the absence and presence of a signal, respectively.
Electronics 12 01123 g001
Figure 2. Directed acyclic graph of the hierarchical Bayesian model.
Figure 2. Directed acyclic graph of the hierarchical Bayesian model.
Electronics 12 01123 g002
Figure 3. (a) The true space–frequency map, as well as the estimation results of (b) the proposed method, (c) JP-WSBL, (d) IP-WSBL, (e) W-SpSF and (f) CBF.
Figure 3. (a) The true space–frequency map, as well as the estimation results of (b) the proposed method, (c) JP-WSBL, (d) IP-WSBL, (e) W-SpSF and (f) CBF.
Electronics 12 01123 g003aElectronics 12 01123 g003b
Figure 4. The comparison between CBF and (a) the proposed method, (b) JP-WSBL, (c) IP-WSBL and (d) W-SpSF under three different cases.
Figure 4. The comparison between CBF and (a) the proposed method, (b) JP-WSBL, (c) IP-WSBL and (d) W-SpSF under three different cases.
Electronics 12 01123 g004aElectronics 12 01123 g004b
Figure 5. The RMSE versus the SNR under Case I.
Figure 5. The RMSE versus the SNR under Case I.
Electronics 12 01123 g005
Figure 6. The RMSE versus the SNR under Case II.
Figure 6. The RMSE versus the SNR under Case II.
Electronics 12 01123 g006
Figure 7. The RMSE versus the SNR under Case III.
Figure 7. The RMSE versus the SNR under Case III.
Electronics 12 01123 g007
Table 1. Mean running time of each method under three different frequency cases.
Table 1. Mean running time of each method under three different frequency cases.
Case ICase IICase III
the proposed method6.07 s 9.73 s11.15 s
JP-WSBL8.30 s 11.76 s12.36 s
IP-WSBL7.56 s10.88 s12.65 s
W-SpSF26.49 s26.28 s27.14 s
CBF0.02 s0.02 s0.02 s
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Yang, Y.; Zhang, Y.; Yang, L.; Wang, Y. Wideband Direction-of-Arrival Estimation Based on Hierarchical Sparse Bayesian Learning for Signals with the Same or Different Frequency Bands. Electronics 2023, 12, 1123. https://doi.org/10.3390/electronics12051123

AMA Style

Yang Y, Zhang Y, Yang L, Wang Y. Wideband Direction-of-Arrival Estimation Based on Hierarchical Sparse Bayesian Learning for Signals with the Same or Different Frequency Bands. Electronics. 2023; 12(5):1123. https://doi.org/10.3390/electronics12051123

Chicago/Turabian Style

Yang, Yixin, Yahao Zhang, Long Yang, and Yong Wang. 2023. "Wideband Direction-of-Arrival Estimation Based on Hierarchical Sparse Bayesian Learning for Signals with the Same or Different Frequency Bands" Electronics 12, no. 5: 1123. https://doi.org/10.3390/electronics12051123

APA Style

Yang, Y., Zhang, Y., Yang, L., & Wang, Y. (2023). Wideband Direction-of-Arrival Estimation Based on Hierarchical Sparse Bayesian Learning for Signals with the Same or Different Frequency Bands. Electronics, 12(5), 1123. https://doi.org/10.3390/electronics12051123

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop