# DOA Estimation for Highly Correlated and Coherent Multipath Signals with Ultralow SNRs

Authors: Li Cheng; Yang Li; Lianying Zou; Yong Qin
Journal: International Journal of Antennas and Propagation
 (2019)

---

## Abstract

In a typical multipath propagation environment, there exists a strong direct path signal accompanying with several weak multipath signals. Due to the strong direct path interference and other masking effects, the Direction-of-Arrival (DOA) of a weak multipath signal is hard to be estimated. In this paper, a novel method is proposed to estimate the DOA of multipath signals with ultralow signal-to-noise ratio (SNR). The main idea is to increase the SNR and signal-to-interference ratio (SIR) of the desired multipath signal in time-delay domain before DOA estimation processing. Firstly, the cross-correlation functions of the direct path signal and the received array signal are calculated. Then, they are combined and constructed to an enhanced array signal. Under certain conditions, the SNR and SIR of the desired signal can be significantly increased. Finally, the DOAs of multipath signals can be estimated by conventional technologies, and the associated time delays can be measured on the DOA-time-shift map. The SNR and SIR gains of the desired signal are analyzed theoretically, and theoretical analysis also indicates that the Cramer–Rao bound can be reduced. Simulation examples are presented to verify the advantages of the proposed method.


---

## Body

## 1. Introduction

Direction-of-Arrival (DOA) plays an important role in wireless communication, radar, sonar, and other areas. In the past few decades, many methods, such as subspace-based approaches, maximum likelihood methods, and sparse-representation-based methods [1–5] have been studied extensively on this topic.In some multipath environments, such as ionospheric sounding [6], global navigation satellite system [7], and airborne radar [8], some of the incident signals are highly correlated with each other, and some of them can be coherent. In this scenario, the performance of conventional DOA estimation methods is affected dramatically. For example, the popular subspace-based methods MUltiple SIgnal Classification (MUSIC) fails to span the signal subspace in the case of correlated waves. MUSIC can be extended to be applicable by adopting decorrelation preprocessing, such as the spatial smoothing [9] and Toeplitz approximation techniques [10]. Also, the maximum likelihood methods [11] and sparse-representation-based methods [12] can be applicable to the coherent scenario. However, they do not consider the case that the signal-to-noise ratio (SNR) of the desired multipath signal whose DOA needs to be estimated can be very low, and there exists strong direct path which can mask the adjacent weak path. For example, in passive bistatic radar system, the power difference between the echo signal and direct path signal can be larger than 100 dB, which indicates the signal-to-interference ratio (SIR) of the desired echo signal is very low [13], thus estimating the DOA of echo signal is a challenging problem [14].In low SNR case, the number of sources is hard to estimate correctly, and the signal subspace may swap with the noise subspace [15]. Some research studies have been conducted: a MUSIC-like DOA estimation method without estimating the number of sources is proposed in [16], its performance holds stable when SNR is low; a sparse representation of array covariance vectors in an overcomplete basis is proposed in [17], which is statistically robust in low SNR cases; a two-stage DOA estimation method for low SNR signals is proposed in [18], which enhances the angular resolution and estimation accuracy, particularly for the case when the array antenna elements received the low SNR signals; a linear prediction orthogonal propagator method is proposed in [19], which is more robust to the noise, especially in low SNR scenarios. Although these methods improved the performance of DOA estimation in low SNR scenarios, the performance is constrained by the Cramer–Rao bound (CRB). The CRB relates to many factors, such as the number of snapshots, the array aperture, angular separations, and correlations of the signals, no doubt that SNR is one of the most important factors among them.This work mainly focuses on the DOA estimation for partly correlated or coherent signals with ultralow SNR in multipath environment, and the precondition that a strong direct path signal is included in the multipath is required, or the reference signal is available in advance. To overcome the SNR limitation, we provide a novel method that can enhance the desired multipath signal before the DOA estimation procedure. We realize this goal through four steps: firstly, the direct path signal is obtained by the conventional DOA estimation and beamforming technologies with the received original array signal (OAS); secondly, the cross-correlation of the direct path signal and OAS are calculated, and then they are combined to a single snapshot array signal; thirdly, the single snapshot is extended to multisnapshots-enhanced array signal (EAS), under certain conditions, the EAS with particular time shift can be enhanced. At last, by using the time-shift scanning and conventional DOA estimation approaches, the DOAs of multipath signals can be estimated and the associated time delays can be measured on a two-dimensional DOA-time-shift map.The main contributions of the work are as follows: (a) based on the OAS, a new array signal EAS is constructed, in which the desired multipath signal can be enhanced; (b) the SNR gain and SIR gain in the EAS are derived in theory; (c) the DOA and the associated time delay of the multipath signal with ultralow SNR can be estimated or measured with the EAS; (d) the CRB for DOA estimation can be reduced with the proposed method starting with the OAS.The rest of the paper is organized as follows: The basic signal model and array signal processing are provided in Section2. The proposed method and its implementation are developed in Section 3. Some simulation examples are presented in Section 4. A brief conclusion appears in Section 5. In the paper, E⋅, ⋅T, ⋅H, ⋅−1, ⋅, ⊙, and superscript ∗ denote the expectation, transpose, Hermitian transpose, inverse, absolute value, Hadamard product, and complex conjugate, respectively. Boldfaced variable denotes matrix or vector.


## 2. Problem Formulation

### 2.1. Array Signal Model

Consider a multipath propagation environment, assuming that there is one direct path signal with its SNR significantly larger than other multipath signals, and there is no strong interference with the same frequency. Multipath signals incident on a uniform linear array (ULA) withM omnidirectional antennas, and the source is in the far field of the array. The received array signal can be expressed as [9](1)xk=∑i=0Lηisk+τiai+nk,where k denotes the snapshot point, the number of multipath signals is L+1, the ith multipath signal is with steering vector ai, ηi denotes amplitude attenuation, and τi denotes time delay. Let i=0 represents the direct path; we can make η0=1 and τ0=0 for reference. Array signal xk, steering vector ai, and array noise nk are M×1 dimensional complex vector. The signals are stationary Gaussian random process, and the additive noise is a spatially white Gaussian process. The power of direct path signal and noise is σs2 and σn2, respectively, and the power of ith multipath signal is ηi2σs2. xk is called OAS in the following.The data covariance matrix can be expressed asExkxHk. Since the actual covariance matrix is unknown in practice, it is often replaced by K snapshots sample data covariance matrix R=1/K∑k=1KxkxHk. The data covariance matrix can be further expressed in different forms which depend on whether the signals are coherent or not [20].The eigen decomposition ofR is(2)R=∑i=1MγiuiuiH=USΓSUSH+UNΓNUNH,where γi and ui are the eigenvalues and the corresponding eigenvectors of R. The eigenvalues are sorted in descending order, γ1>⋯>γd≫γd+1=⋯=γM=σn2, where d is the number of the distinguishable signals. The eigenvectors corresponding to the largest d eigenvalues span the signal subspace US=u1,…,ud, and the other eigenvectors corresponding to the smaller eigenvalues span the noise subspace UN=ud+1,…,uM.


### 2.2. DOA Estimation

There are numerous methods that can be used to estimate the DOA of directional signal with OAS. One of the classic methods is Capon spatial spectrum, which is calculated by(3)PCaponθ=1aHθR−1aθ,where θ is the DOA of the incident signal. The peaks of the spatial spectrum correspond to the DOAs of distinguishable signals. Although the Capon spatial spectrum is implemented easily, its resolution is limited, and the performance is decreased in low SNR scenarios.The MUSIC is one of the very popular methods for super-resolution DOA estimation [21], and it is given by(4)PMUSICθ=1aHθUNUNHaθ.Resolution and accuracy of an eigen structure method like MUSIC depend on the number of snapshots, array aperture, angular separations, SNRs, and correlations of the signals, etc. The CRB provides an algorithm-independent benchmark against which various algorithms can be compared. The stochastic CRB of DOA estimation with OAS is computed as [22](5)CRBOAS=12KReDHI−AAHA−1AHD⊙BxT−1,where D=∂aθ0/∂θ0,∂aθ1/∂θ1,…,∂aθL/∂θL, A=a0,a1,…,aL, Bx=diagSNROASi, and SNROASi=ηi2σs2/σn2, i=0,1,…,L. When other parameters are fixed, the snapshots number K and the SNR (relate to R and σn2) decide the CRB of OAS. In the following, we will develop a novel method to reduce the CRB with the same OAS.


## 2.1. Array Signal Model

Consider a multipath propagation environment, assuming that there is one direct path signal with its SNR significantly larger than other multipath signals, and there is no strong interference with the same frequency. Multipath signals incident on a uniform linear array (ULA) withM omnidirectional antennas, and the source is in the far field of the array. The received array signal can be expressed as [9](1)xk=∑i=0Lηisk+τiai+nk,where k denotes the snapshot point, the number of multipath signals is L+1, the ith multipath signal is with steering vector ai, ηi denotes amplitude attenuation, and τi denotes time delay. Let i=0 represents the direct path; we can make η0=1 and τ0=0 for reference. Array signal xk, steering vector ai, and array noise nk are M×1 dimensional complex vector. The signals are stationary Gaussian random process, and the additive noise is a spatially white Gaussian process. The power of direct path signal and noise is σs2 and σn2, respectively, and the power of ith multipath signal is ηi2σs2. xk is called OAS in the following.The data covariance matrix can be expressed asExkxHk. Since the actual covariance matrix is unknown in practice, it is often replaced by K snapshots sample data covariance matrix R=1/K∑k=1KxkxHk. The data covariance matrix can be further expressed in different forms which depend on whether the signals are coherent or not [20].The eigen decomposition ofR is(2)R=∑i=1MγiuiuiH=USΓSUSH+UNΓNUNH,where γi and ui are the eigenvalues and the corresponding eigenvectors of R. The eigenvalues are sorted in descending order, γ1>⋯>γd≫γd+1=⋯=γM=σn2, where d is the number of the distinguishable signals. The eigenvectors corresponding to the largest d eigenvalues span the signal subspace US=u1,…,ud, and the other eigenvectors corresponding to the smaller eigenvalues span the noise subspace UN=ud+1,…,uM.


## 2.2. DOA Estimation

There are numerous methods that can be used to estimate the DOA of directional signal with OAS. One of the classic methods is Capon spatial spectrum, which is calculated by(3)PCaponθ=1aHθR−1aθ,where θ is the DOA of the incident signal. The peaks of the spatial spectrum correspond to the DOAs of distinguishable signals. Although the Capon spatial spectrum is implemented easily, its resolution is limited, and the performance is decreased in low SNR scenarios.The MUSIC is one of the very popular methods for super-resolution DOA estimation [21], and it is given by(4)PMUSICθ=1aHθUNUNHaθ.Resolution and accuracy of an eigen structure method like MUSIC depend on the number of snapshots, array aperture, angular separations, SNRs, and correlations of the signals, etc. The CRB provides an algorithm-independent benchmark against which various algorithms can be compared. The stochastic CRB of DOA estimation with OAS is computed as [22](5)CRBOAS=12KReDHI−AAHA−1AHD⊙BxT−1,where D=∂aθ0/∂θ0,∂aθ1/∂θ1,…,∂aθL/∂θL, A=a0,a1,…,aL, Bx=diagSNROASi, and SNROASi=ηi2σs2/σn2, i=0,1,…,L. When other parameters are fixed, the snapshots number K and the SNR (relate to R and σn2) decide the CRB of OAS. In the following, we will develop a novel method to reduce the CRB with the same OAS.


## 3. The Proposed Method

### 3.1. Direct Path Signal Beamforming

Firstly, the direct path signal is required. In some applications, the reference signal such as the pilot signal in communication and the transmitted signal in active radar are priori known, and the direct path signal can be obtained from the reference signal. However, the reference signal is unavailable in many applications. The beamforming technology which is used to filter the directional signal in space can be used to obtain the direct path signal:yDk=wHxk, where w is the beamformer’s weight vector with its main beam pointing to the DOA of direct path signal. In order to suppress the interference adaptively as well as prevent desired signal self-nulling, the robust adaptive beamforming technologies are suggested to be used. The robust Capon beamformer [23], the robust adaptive beamformer [24], or their improved technologies [25, 26] can be used. After the adaptive beamforming, we assume the obtained direct path signal is pure enough, which means yDk=wHxk≈sk.


### 3.2. Array Cross-Correlation

The self-correlation function ofsk is calculated by(6)cKτ=∑k=1Ksks∗k+τ.When the time shiftτ=0, cK0≈Kσs2. The self-correlation function has a property that the response is maximized when the time shift τ equals to 0:(7)cK0≥cKτ.The cross-correlation is a measure of similarity of two series as a function of the displacement of one relative to the other. The discrete form of cross-correlation function of themth antenna’s received signal xmk and direct path signal yDk≈sk is defined as(8)rKmτ=∑k=1KxmkyD∗k+τ≈∑k=1Kxmks∗k+τ,m=1,2,…,M.The peaks of cross-correlation function correspond to the multipath signals that are matched inτ=τi.Taking a combination ofrKmτ as rKτ=rK1τ,rK2τ,…,rKMτT, it can be decomposed as(9)rKτ=∑k=1Kxks∗k+τ=∑i=0L∑k=1Kηisk+τis∗k+τai+∑k=1Ks∗k+τnk=∑i=0LηicKτ−τiai+∑k=1Ks∗k+τnk.The matrixrKτ is M×1 dimensional, and its form is similar with OAS: the signal component is ∑i=0LηicKτ−τiai, and the noise component is ∑k=1Ks∗k+τnk.Assuming that the path-1 signal witha1,τ1 is the desired signal whose DOA needs to be estimated. If τ=τ1, (9) can be decomposed as(10)rKτ1≈η1cK0a1+∑i=0,2LηicKτ1−τiai+∑k=1Ks∗k+τ1nk.Through the Cauchy–Schwarz inequality, we have(11)∑k=1Ks∗k+τ1nmk2≤∑k=1Ksk+τ12∑k=1Knmk2=K2σs2σn2.The SIR and the SNR of the path-1 signal inxk and rKτ1 are listed in Table 1. Results show that the SIR and SNR of the path-1 is enhanced in rKτ1.Table 1
The comparison of SIR and SNR.
x
k
=
∑
i
=
0
L
η
i
s
k
+
τ
i
a
i
+
n
k
r
K
τ
1
≈
η
1
c
K
0
a
1
+
∑
i
=
0,2
L
η
i
c
K
τ
1
−
τ
i
a
i
+
∑
k
=
1
K
s
∗
k
+
τ
1
n
k
SIR
=
η
1
2
/
η
i
2
=
η
1
2
c
2
0
/
η
i
2
c
2
τ
1
−
τ
i
≥
η
1
2
/
η
i
2
SNR
=
η
1
2
σ
s
2
/
σ
n
2
≥
η
1
2
c
2
0
/
K
2
σ
s
2
σ
n
2
=
η
1
2
σ
s
2
/
σ
n
2


### 3.3. Enhanced Array Signal

Unfortunately,rKτ1 in (10) has only one snapshot; thus, it is hard to use conventional DOA estimation approaches directly. Although there are many methods that can be used to estimate the DOA with a single snapshot, their performance is decreased compared with multisnapshots. In the following, the single snapshot array signal rKτ1 is extended to multisnapshots array signal.As shown in Table2, the K-length snapshots in xk, k = 1, 2, …, K are divided into K2 uniform segments, and each segment has K1=K/K2 subsnapshots.Table 2
Multisnapshots array signal.
Segments
1
2
…
K
2
Subsnapshots
1∼K1
K
1+1∼2K1
…
1 (K2–1)K1+1∼K2 K1
OAS
x (1∼K1)
x (K1+1∼2K1)
…
x[(K2–1)K1+1∼K2 K1]
Direct path
s (1∼K1)
s (K1+1∼2K1)
…
s[(K2–1)K1+1∼K2 K1]
Correlation
r
K1 (1, τ)
r
K
1
2
,
τ
…
r
K1 (K2, τ)The cross-correlation of the direct path signal and the array signal for subsnapshotsk2 is(12)rK1k2,τ=∑k=k2−1K1+1k2K1xks∗k+τ=∑i=0L∑k=k2−1K1+1k2K1ηisk+τis∗k+τai+∑k=k2−1K1+1k2K1s∗k+τnk=∑i=0LηicK1τ−τiai+∑k=k2−1K1+1k2K1s∗k+τnk,k2=1,2,…,K2.Comparing withrKτ in (9), the form of rK1k2,τ in (12) is also similar with OAS. We call the combination rK11,τ,rK12,τ,…,rK1K2,τ as the EAS, in which the path with time delay τ is enhanced. The difference and equivalence of OAS in (1) and EAS in (12) are described in Table 3.Table 3
Comparison of OAS and EAS.
x
k
r
K
1
k
2
,
τ
Path-i array signal
η
i
s
k
+
τ
i
a
i
∑
k
=
k
2
−
1
K
1
+
1
k
2
K
1
η
i
s
k
+
τ
i
s
∗
k
+
τ
a
i
Array noise
n
k
∑
k
=
k
2
−
1
K
1
+
1
k
2
K
1
s
∗
k
+
τ
n
k
Snapshot sequence
k
=
1
,
…
,
K
k
2
=
1,2
,
…
,
K
2Assuming that the snapshots are infinite, according to the discrete time matched filters theory, the noise element atmth antenna is(13)E∑k=1Ks∗k+τ1nmk2=E∑ka=1Ks∗ka+τ1nmka∑kb=1Kskb+τ1nm∗kb=∑ka=1K∑kb=1Ks∗ka+τ1skb+τ1Enmkanm∗kb=∑ka=1K∑kb=1Ks∗ka+τ1skb+τ1σn2δka,kb=σn2∑k=1Ksk+τ12=Kσs2σn2,m=1,2,…,M,where δ is an impulse function. Hence, for finite snapshots, we can make(14)1K2∑k2=1K2∑k=k2−1K1+1k2K1s∗k+τnmk2≈K1σs2σn2.The power of signal in path-i is(15)1K2∑k2=1K2∑k=k2−1K1+1k2K1ηisk+τis∗k+τ2=ηi2K2∑k2=1K2cK12τ−τi=ηicK12τ−τi.Ifτ=τ1, the relationship of SIR and SNR in the OAS and EAS are(16)SIREAS=η12cK120ηi2cK12τ1−τi=cK10cK1τ1−τi2SIROAS,(17)SNREAS≈η12cK120K1σs2σn2=η12K12σs4K1σs2σn2=K1η12σs2σn2=K1SNROAS.Hence, compared with the OAS, the SIR in EAS is enhanced withcK10/cK1τ1−τi2 times, and the SNR in EAS is enhanced with K1 times. Notice that the subsnapshots can be set partly overlapped, and then the snapshots K2 or the SNR gain K1 can be further increased.It is worth noting that the noise∑k=k2−1K1+1k2K1s∗k+τnk may not be white. However, in the detection and estimation theory, if the SNR is high, the noise covariance matrix can be assumed to be inappreciable [27]. In EAS, the SNR of the desired signal has been increased about K1 times, which can be extremely large; hence, we can treat the noise in EAS as white noise for simplification.


### 3.4. CRB Analysis

There are many existing methods which can be used to estimate the DOAs of multipath signals with EAS. Because the array signal structures in OAS and EAS are similar, the CRB for DOA estimation with EAS is similar with (5)(18)CRBEAS=12K2ReDHI−AAHA−1AHD⊙BrT−1,where Br=diagSNREASi, i=0,1,…,L. To compute the CRB of a single path-i, we can set SNREASi=K1SNROASi. Comparing the CRBEAS and the CRBOAS, although the snapshots number K2<K, the SNR satisfies SNREASi>SNROASi. Therefore, the CRBEAS can be lower than CRBOAS under certain conditions.


### 3.5. Spatial Spectrum Estimation

The MUSIC approach is used to estimate the spatial spectrum of multipath signals. The data covariance matrix of EAS is(19)Rτ=1K2∑k2=1K2rk2,τrHk2,τ.Taking the eigenvalue decomposition ofRt=ESΛSESH+ENΛNENH, where ES and EN are the signal subspace and noise subspace, respectively. When τ=τi, since the SIR is enlarged in EAS, which indicates other paths are attenuated, it is foreseeable that there are two signals having the same lager power: the path-i signal and the direct path signal. Hence, it is reasonable to set the dimension of EN to M−2. When τ≠τi, the M−2 dimensional noise subspace may cause another false alarm peak; fortunately, it is so low that can be omitted. Therefore, the dimension of EN can be set equal to M−2 directly.Finally, the MUSIC spatial spectrum is calculated by(20)PMUSICθ,τ=1aHθENENHaθ.Assuming the time delays of multipath signals lie betweenτmin and τmax, the range τmin,τmax can be divided to discrete time-shift sequences τ1, τ2,…, and then the spatial spectra with different time-shift sequences are calculated. By stacking all the spatial spectrums along the time shifts, we can obtain a two dimensional DOA-time-shift map. Only if the time shift equals to the time delay of a multipath signal, the weak multipath signal can be enhanced, and a peak corresponding to its DOA will appear on the spatial spectrum. Therefore, the DOA of the ith multipath signal with time delay τ=τi corresponds to a peak on the DOA-time-shift map. That is, to say, we can estimate the DOAs and measure their associated time delays on the DOA-time-shift map for all the distinguishable multipath signals.The MUSIC approach can be replaced, and other DOA estimation methods can be used with the EAS, so do the technologies can deal with coherent signals, such as the spatial smoothing technology [9].


### 3.6. Implementation

The implementation of the proposed method is summarized as follows:Step 1: estimating the DOA of direct path signal withK-snapshots OAS; any DOA estimation technology is suitable.Step 2: obtaining the direct path signal with beamforming technology.Step 3: dividing theK-snapshots OAS and direct path signal into k2 = 1, 2, …, K2 segments, and each segment has the same length K1. For each segment, calculating the cross-correlation function of the direct path signal and the OAS by (12).Step 4: scanning in the time shift to estimate the DOAs with MUSIC approach:(21)Forτ=τ1,τ2,…,Rτ=1K2∑k2=1K2rk2,τrHk2,τ,Rt=ESΛSESH+ENΛNENH,Pθ,τ=1aHθENENHaθ,End.The spectral peaks of two-dimensional DOA-time-shift mapPθ,τ are corresponding to the DOAs and associated time delays of multipath signals. Step 4 can be replaced by the MUSIC approach with forward spatial smoothing technology, if coherent multipath signal exists [9]. In case of two coherent signals, the dimension of the noise subspace should be set to Ms−3, where Ms is the dimension of subarrays in the forward spatial smoothing technology.The computational complexity of the main operations are as follows: (a) beamforming in step 2:OM3, if the method in [24] is used; (b) DOA estimation in steps 1 and 4: OM3, if the MUSIC is used; and (c) cross-correlation in step 3: OK1logK1, if the fast Fourier transform is used.


## 3.1. Direct Path Signal Beamforming

Firstly, the direct path signal is required. In some applications, the reference signal such as the pilot signal in communication and the transmitted signal in active radar are priori known, and the direct path signal can be obtained from the reference signal. However, the reference signal is unavailable in many applications. The beamforming technology which is used to filter the directional signal in space can be used to obtain the direct path signal:yDk=wHxk, where w is the beamformer’s weight vector with its main beam pointing to the DOA of direct path signal. In order to suppress the interference adaptively as well as prevent desired signal self-nulling, the robust adaptive beamforming technologies are suggested to be used. The robust Capon beamformer [23], the robust adaptive beamformer [24], or their improved technologies [25, 26] can be used. After the adaptive beamforming, we assume the obtained direct path signal is pure enough, which means yDk=wHxk≈sk.


## 3.2. Array Cross-Correlation

The self-correlation function ofsk is calculated by(6)cKτ=∑k=1Ksks∗k+τ.When the time shiftτ=0, cK0≈Kσs2. The self-correlation function has a property that the response is maximized when the time shift τ equals to 0:(7)cK0≥cKτ.The cross-correlation is a measure of similarity of two series as a function of the displacement of one relative to the other. The discrete form of cross-correlation function of themth antenna’s received signal xmk and direct path signal yDk≈sk is defined as(8)rKmτ=∑k=1KxmkyD∗k+τ≈∑k=1Kxmks∗k+τ,m=1,2,…,M.The peaks of cross-correlation function correspond to the multipath signals that are matched inτ=τi.Taking a combination ofrKmτ as rKτ=rK1τ,rK2τ,…,rKMτT, it can be decomposed as(9)rKτ=∑k=1Kxks∗k+τ=∑i=0L∑k=1Kηisk+τis∗k+τai+∑k=1Ks∗k+τnk=∑i=0LηicKτ−τiai+∑k=1Ks∗k+τnk.The matrixrKτ is M×1 dimensional, and its form is similar with OAS: the signal component is ∑i=0LηicKτ−τiai, and the noise component is ∑k=1Ks∗k+τnk.Assuming that the path-1 signal witha1,τ1 is the desired signal whose DOA needs to be estimated. If τ=τ1, (9) can be decomposed as(10)rKτ1≈η1cK0a1+∑i=0,2LηicKτ1−τiai+∑k=1Ks∗k+τ1nk.Through the Cauchy–Schwarz inequality, we have(11)∑k=1Ks∗k+τ1nmk2≤∑k=1Ksk+τ12∑k=1Knmk2=K2σs2σn2.The SIR and the SNR of the path-1 signal inxk and rKτ1 are listed in Table 1. Results show that the SIR and SNR of the path-1 is enhanced in rKτ1.Table 1
The comparison of SIR and SNR.
x
k
=
∑
i
=
0
L
η
i
s
k
+
τ
i
a
i
+
n
k
r
K
τ
1
≈
η
1
c
K
0
a
1
+
∑
i
=
0,2
L
η
i
c
K
τ
1
−
τ
i
a
i
+
∑
k
=
1
K
s
∗
k
+
τ
1
n
k
SIR
=
η
1
2
/
η
i
2
=
η
1
2
c
2
0
/
η
i
2
c
2
τ
1
−
τ
i
≥
η
1
2
/
η
i
2
SNR
=
η
1
2
σ
s
2
/
σ
n
2
≥
η
1
2
c
2
0
/
K
2
σ
s
2
σ
n
2
=
η
1
2
σ
s
2
/
σ
n
2


## 3.3. Enhanced Array Signal

Unfortunately,rKτ1 in (10) has only one snapshot; thus, it is hard to use conventional DOA estimation approaches directly. Although there are many methods that can be used to estimate the DOA with a single snapshot, their performance is decreased compared with multisnapshots. In the following, the single snapshot array signal rKτ1 is extended to multisnapshots array signal.As shown in Table2, the K-length snapshots in xk, k = 1, 2, …, K are divided into K2 uniform segments, and each segment has K1=K/K2 subsnapshots.Table 2
Multisnapshots array signal.
Segments
1
2
…
K
2
Subsnapshots
1∼K1
K
1+1∼2K1
…
1 (K2–1)K1+1∼K2 K1
OAS
x (1∼K1)
x (K1+1∼2K1)
…
x[(K2–1)K1+1∼K2 K1]
Direct path
s (1∼K1)
s (K1+1∼2K1)
…
s[(K2–1)K1+1∼K2 K1]
Correlation
r
K1 (1, τ)
r
K
1
2
,
τ
…
r
K1 (K2, τ)The cross-correlation of the direct path signal and the array signal for subsnapshotsk2 is(12)rK1k2,τ=∑k=k2−1K1+1k2K1xks∗k+τ=∑i=0L∑k=k2−1K1+1k2K1ηisk+τis∗k+τai+∑k=k2−1K1+1k2K1s∗k+τnk=∑i=0LηicK1τ−τiai+∑k=k2−1K1+1k2K1s∗k+τnk,k2=1,2,…,K2.Comparing withrKτ in (9), the form of rK1k2,τ in (12) is also similar with OAS. We call the combination rK11,τ,rK12,τ,…,rK1K2,τ as the EAS, in which the path with time delay τ is enhanced. The difference and equivalence of OAS in (1) and EAS in (12) are described in Table 3.Table 3
Comparison of OAS and EAS.
x
k
r
K
1
k
2
,
τ
Path-i array signal
η
i
s
k
+
τ
i
a
i
∑
k
=
k
2
−
1
K
1
+
1
k
2
K
1
η
i
s
k
+
τ
i
s
∗
k
+
τ
a
i
Array noise
n
k
∑
k
=
k
2
−
1
K
1
+
1
k
2
K
1
s
∗
k
+
τ
n
k
Snapshot sequence
k
=
1
,
…
,
K
k
2
=
1,2
,
…
,
K
2Assuming that the snapshots are infinite, according to the discrete time matched filters theory, the noise element atmth antenna is(13)E∑k=1Ks∗k+τ1nmk2=E∑ka=1Ks∗ka+τ1nmka∑kb=1Kskb+τ1nm∗kb=∑ka=1K∑kb=1Ks∗ka+τ1skb+τ1Enmkanm∗kb=∑ka=1K∑kb=1Ks∗ka+τ1skb+τ1σn2δka,kb=σn2∑k=1Ksk+τ12=Kσs2σn2,m=1,2,…,M,where δ is an impulse function. Hence, for finite snapshots, we can make(14)1K2∑k2=1K2∑k=k2−1K1+1k2K1s∗k+τnmk2≈K1σs2σn2.The power of signal in path-i is(15)1K2∑k2=1K2∑k=k2−1K1+1k2K1ηisk+τis∗k+τ2=ηi2K2∑k2=1K2cK12τ−τi=ηicK12τ−τi.Ifτ=τ1, the relationship of SIR and SNR in the OAS and EAS are(16)SIREAS=η12cK120ηi2cK12τ1−τi=cK10cK1τ1−τi2SIROAS,(17)SNREAS≈η12cK120K1σs2σn2=η12K12σs4K1σs2σn2=K1η12σs2σn2=K1SNROAS.Hence, compared with the OAS, the SIR in EAS is enhanced withcK10/cK1τ1−τi2 times, and the SNR in EAS is enhanced with K1 times. Notice that the subsnapshots can be set partly overlapped, and then the snapshots K2 or the SNR gain K1 can be further increased.It is worth noting that the noise∑k=k2−1K1+1k2K1s∗k+τnk may not be white. However, in the detection and estimation theory, if the SNR is high, the noise covariance matrix can be assumed to be inappreciable [27]. In EAS, the SNR of the desired signal has been increased about K1 times, which can be extremely large; hence, we can treat the noise in EAS as white noise for simplification.


## 3.4. CRB Analysis

There are many existing methods which can be used to estimate the DOAs of multipath signals with EAS. Because the array signal structures in OAS and EAS are similar, the CRB for DOA estimation with EAS is similar with (5)(18)CRBEAS=12K2ReDHI−AAHA−1AHD⊙BrT−1,where Br=diagSNREASi, i=0,1,…,L. To compute the CRB of a single path-i, we can set SNREASi=K1SNROASi. Comparing the CRBEAS and the CRBOAS, although the snapshots number K2<K, the SNR satisfies SNREASi>SNROASi. Therefore, the CRBEAS can be lower than CRBOAS under certain conditions.


## 3.5. Spatial Spectrum Estimation

The MUSIC approach is used to estimate the spatial spectrum of multipath signals. The data covariance matrix of EAS is(19)Rτ=1K2∑k2=1K2rk2,τrHk2,τ.Taking the eigenvalue decomposition ofRt=ESΛSESH+ENΛNENH, where ES and EN are the signal subspace and noise subspace, respectively. When τ=τi, since the SIR is enlarged in EAS, which indicates other paths are attenuated, it is foreseeable that there are two signals having the same lager power: the path-i signal and the direct path signal. Hence, it is reasonable to set the dimension of EN to M−2. When τ≠τi, the M−2 dimensional noise subspace may cause another false alarm peak; fortunately, it is so low that can be omitted. Therefore, the dimension of EN can be set equal to M−2 directly.Finally, the MUSIC spatial spectrum is calculated by(20)PMUSICθ,τ=1aHθENENHaθ.Assuming the time delays of multipath signals lie betweenτmin and τmax, the range τmin,τmax can be divided to discrete time-shift sequences τ1, τ2,…, and then the spatial spectra with different time-shift sequences are calculated. By stacking all the spatial spectrums along the time shifts, we can obtain a two dimensional DOA-time-shift map. Only if the time shift equals to the time delay of a multipath signal, the weak multipath signal can be enhanced, and a peak corresponding to its DOA will appear on the spatial spectrum. Therefore, the DOA of the ith multipath signal with time delay τ=τi corresponds to a peak on the DOA-time-shift map. That is, to say, we can estimate the DOAs and measure their associated time delays on the DOA-time-shift map for all the distinguishable multipath signals.The MUSIC approach can be replaced, and other DOA estimation methods can be used with the EAS, so do the technologies can deal with coherent signals, such as the spatial smoothing technology [9].


## 3.6. Implementation

The implementation of the proposed method is summarized as follows:Step 1: estimating the DOA of direct path signal withK-snapshots OAS; any DOA estimation technology is suitable.Step 2: obtaining the direct path signal with beamforming technology.Step 3: dividing theK-snapshots OAS and direct path signal into k2 = 1, 2, …, K2 segments, and each segment has the same length K1. For each segment, calculating the cross-correlation function of the direct path signal and the OAS by (12).Step 4: scanning in the time shift to estimate the DOAs with MUSIC approach:(21)Forτ=τ1,τ2,…,Rτ=1K2∑k2=1K2rk2,τrHk2,τ,Rt=ESΛSESH+ENΛNENH,Pθ,τ=1aHθENENHaθ,End.The spectral peaks of two-dimensional DOA-time-shift mapPθ,τ are corresponding to the DOAs and associated time delays of multipath signals. Step 4 can be replaced by the MUSIC approach with forward spatial smoothing technology, if coherent multipath signal exists [9]. In case of two coherent signals, the dimension of the noise subspace should be set to Ms−3, where Ms is the dimension of subarrays in the forward spatial smoothing technology.The computational complexity of the main operations are as follows: (a) beamforming in step 2:OM3, if the method in [24] is used; (b) DOA estimation in steps 1 and 4: OM3, if the MUSIC is used; and (c) cross-correlation in step 3: OK1logK1, if the fast Fourier transform is used.


## 4. Simulation Examples

In the following simulation examples, a ULA withM=16 antennas and half-wavelength spacing is considered. Assuming each antenna is omnidirectional, the array has been calibrated. The signal is stationary Gaussian random process, and the additive noise is a spatially white Gaussian process. The baseband signal is generated by a complex white noise sequence filtered with a low-pass filter, and the bandwidth is 80 kHz (also 40 kHz in Example 1). The sampling rate is 100 kHz, and the segments number K2 is 50.The proposed method is compared with the following three methods: (a) Capon: the Capon spatial spectrum estimator in (3), which is an easy implantation and widely used method; (b) MUSIC: the MUSIC spatial spectrum estimator in (4), which is a well-known super resolution method. (c) OMP-LP: the linear prediction orthogonal propagator method in [19], which is robust to the noise, especially in low SNR scenarios.


### 4.1. SIR and SNR Gain

In the first example, we evaluate the influence of signal bandwidth and snapshot number on the gain of SIR and SNR. The SNR and time delay parameters are set the same as Example 2, while the bandwidth is set to 40 kHz and 80 kHz, and the snapshot numberK varies from 5000 to 50,000 which means K1 varies from 100 to 1000. Assume that path-2 is the desired signal and path-0 and path-1 are interference.The SIR gain is calculated by(22)GainSIRi=SIREASSIROAS=∑k2=1K2∑k=k2−1K1+1k2K1sk+τ2s∗k+τ22∑k2=1K2∑k=k2−1K1+1k2K1sk+τis∗k+τ22,i=0,1.The SNR gain is calculated by(23)GainSNR=SNREASSNROAS=∑k2=1K2∑k=k2−1K1+1k2K1sk+τ2s∗k+τ22∑k2=1K2∑k=k2−1K1+1k2K1s∗k+τ2n1k2×σn2.Example 1.
The SIR and SNR gain simulation.
Figure1 shows the SIR and SNR gain with different bandwidths and snapshot numbers. Results show that (a) the gain SIR0 and gain SIR1 are almost equal for the same bandwidth and K. It is because the path-0 and path-1 interference are attenuated at almost the same level for the simulated signal. Actually, it can be concluded from (16) that the attenuation level of interference relates to the shape of cross-correlation function of the signal. (b) For the same K, when the bandwidth increases from 40 kHz to 80 kHz, the SIR gain doubled. It is because the value of cross-correlation processing depends on the time bandwidth product; if time K is fixed and bandwidth doubled, the gain doubled correspondingly. (c) The SIR gain and SNR gain increase almost linearly as the snapshot number K increases. (d) The SNR gain is independent of bandwidth, and it is very close to the theoretical result in (17), which indicates the SNR gain equals to K1.Figure 1
The SIR and SNR gains with different bandwidths and snapshots.


### 4.2. Performance Comparison

In Examples 2 and 3, the multipath signals’ parameters are listed in Table4. Other parameters are: K = 50000, K1 = 1000, and K2 = 50.Table 4
Simulated parameters for examples 2 and 3.
Multipath
Path-0
Path-1
Path-2
DOA (degree)
60
90
120
SNR (dB)
10
−20
−30
Time delay (ms)
0
0.5
0.8Example 2.
Spatial spectrum comparison.
Figure2 displays the DOA-time-shift spectrum map. Results show that (a) there is a line at θ = 30° throughout all the time delays, and it is the residual direct path extended to all the time delays. (b) The peak of the line is with DOA = 30°, and time delay = 0 is corresponding to the direct path signal path-0. (c) The DOAs and their associated time delays of two weak signals path-1 and path-2 can be distinguished correctly and clearly.
Further study shows that the signal with larger bandwidth has higher resolution in time delay. Also, it is shown in Figure1 that the larger the bandwidth, the larger the SIR gain. However, the MUSIC is a narrowband approach. To apply the proposed method in both the DOA and time delay estimation, the baseband signal should be narrowband compared with the carrier signal, and then the larger the bandwidth, the better the performance, if other parameters are fixed.
To compare the performance with other methods, we extract the maximum value along the time delay 0∼1 ms at each direction to construct a one-dimensional spatial spectrum for the proposed method. Results in Figure3 show that (a) all the methods can estimate the DOA of the direct path signal with SNR = 10 dB, the OMP-LP generates the sharpest peak at the direct path. (b) Both the MUSIC and the OMP-LP and the proposed method can estimate the DOA of path-1 with SNR = –20 dB, and the proposed method generates the sharpest peak. (c) Only the proposed method generates a peak at the DOA of path-2 with SNR = –30 dB.Figure 2
DOA-time-shift spectrum.Figure 3
Results comparison for the four methods.Example 3.
The RMSE comparison.
The SNR of path-2 varies from −40 dB to 0 dB. The root mean square error (RMSE) is calculated for path-2,RMSE=1/1000∑i=11000θ^2i−θ22, where θ2 and θ^2i are the actual and the ith time estimated DOA of path-2, respectively, and θ^2i is obtained by searching the maximum value in 110° : 0.01° : 130°. 1000 independent runs are performed. Results in Figure 4 show that (a) benefiting from the desired signal is enhanced; the CRB of EAS is lower than the CRB of OAS, and as the SNR becomes lower, the CRB-EAS gets much lower than the CRB-OAS. (b) The RMSE of the proposed method is lower than other methods, especially at much lower SNR scenario. (c) The RMSE of the proposed method is lower than the CRB of OAS.
The SNR of path-2 is fixed at −20 dB, while the number of snapshotsK varies from 5000 to 50,000, which indicates K2 = 50, and K1 varies from 100 to 1000. Results in Figure 5 show that (a) the CRB of EAS is lower than the CRB of OAS. (b) As the number of snapshots increase, the RMSE of the proposed method gets lower, and when the number of subsnapshots K1 surplus 200, the performance of the proposed method is stable. (c) The RMSE of the proposed method is lower than other methods and is lower than the CRB of OAS when K1 surplus 150.Figure 4
RMSE versus SNR.Figure 5
RMSE versus snapshot number.


### 4.3. The Same DOA and Coherent Signals with Ultralow SNRs

In Examples 4 and 5, the multipath signals’ parameters are listed in Table5. Other parameters are: K = 500000, K1 = 10000, and K2 = 50. The signals in path-4 and path-5 have the same DOA, and the signals in path-3 and path-4 are coherent.Table 5
Simulated parameters for Examples 4 and 5.
Multipath
Path-0
Path-1
Path-2
Path-3
Path-4
Path-5
Path-6
Path-7
DOA (degree)
30
70
130
120
100
100
140
80
SNR (dB)
10
−36
−40
−37
−35
−38
−34
−39
Time delay (ms)
0
0.2
0.1
0.4
0.4
0.6
0.8
0.9Example 4.
The proposed method applied in ultralow SNR
To our knowledge, the existing DOA estimation methods cannot work at ultralow SNR scenario, such as the SNR is below −40 dB, and there exists a strong direct path at the same time. Figure6 shows the DOA-time-shift spectrum calculated by the proposed method. Results show that (a) the DOAs and the associated time delays of path-1, -2, -5, -6, and -7 can be distinguished clearly. (b) Path-3 and path-4 cannot be distinguished because they are coherent.Figure 6
DOA-time-shift spectrum.Example 5.
The proposed method with spatial smoothing.
Then, we apply the MUSIC approach with forward spatial smoothing technology [9] in the proposed method. The number of subarrays is set to 2, which means the dimension of subarray is M−1=15. The dimension of noise subspace is set to 15 − 3 = 12. Results in Figure 7 show that the DOAs as well as their associated time delays of two coherent signals path-3 and path-4 can be distinguished correctly and clearly.Figure 7
DOA-time-shift spectrum with spatial smoothing.


## 4.1. SIR and SNR Gain

In the first example, we evaluate the influence of signal bandwidth and snapshot number on the gain of SIR and SNR. The SNR and time delay parameters are set the same as Example 2, while the bandwidth is set to 40 kHz and 80 kHz, and the snapshot numberK varies from 5000 to 50,000 which means K1 varies from 100 to 1000. Assume that path-2 is the desired signal and path-0 and path-1 are interference.The SIR gain is calculated by(22)GainSIRi=SIREASSIROAS=∑k2=1K2∑k=k2−1K1+1k2K1sk+τ2s∗k+τ22∑k2=1K2∑k=k2−1K1+1k2K1sk+τis∗k+τ22,i=0,1.The SNR gain is calculated by(23)GainSNR=SNREASSNROAS=∑k2=1K2∑k=k2−1K1+1k2K1sk+τ2s∗k+τ22∑k2=1K2∑k=k2−1K1+1k2K1s∗k+τ2n1k2×σn2.Example 1.
The SIR and SNR gain simulation.
Figure1 shows the SIR and SNR gain with different bandwidths and snapshot numbers. Results show that (a) the gain SIR0 and gain SIR1 are almost equal for the same bandwidth and K. It is because the path-0 and path-1 interference are attenuated at almost the same level for the simulated signal. Actually, it can be concluded from (16) that the attenuation level of interference relates to the shape of cross-correlation function of the signal. (b) For the same K, when the bandwidth increases from 40 kHz to 80 kHz, the SIR gain doubled. It is because the value of cross-correlation processing depends on the time bandwidth product; if time K is fixed and bandwidth doubled, the gain doubled correspondingly. (c) The SIR gain and SNR gain increase almost linearly as the snapshot number K increases. (d) The SNR gain is independent of bandwidth, and it is very close to the theoretical result in (17), which indicates the SNR gain equals to K1.Figure 1
The SIR and SNR gains with different bandwidths and snapshots.


## 4.2. Performance Comparison

In Examples 2 and 3, the multipath signals’ parameters are listed in Table4. Other parameters are: K = 50000, K1 = 1000, and K2 = 50.Table 4
Simulated parameters for examples 2 and 3.
Multipath
Path-0
Path-1
Path-2
DOA (degree)
60
90
120
SNR (dB)
10
−20
−30
Time delay (ms)
0
0.5
0.8Example 2.
Spatial spectrum comparison.
Figure2 displays the DOA-time-shift spectrum map. Results show that (a) there is a line at θ = 30° throughout all the time delays, and it is the residual direct path extended to all the time delays. (b) The peak of the line is with DOA = 30°, and time delay = 0 is corresponding to the direct path signal path-0. (c) The DOAs and their associated time delays of two weak signals path-1 and path-2 can be distinguished correctly and clearly.
Further study shows that the signal with larger bandwidth has higher resolution in time delay. Also, it is shown in Figure1 that the larger the bandwidth, the larger the SIR gain. However, the MUSIC is a narrowband approach. To apply the proposed method in both the DOA and time delay estimation, the baseband signal should be narrowband compared with the carrier signal, and then the larger the bandwidth, the better the performance, if other parameters are fixed.
To compare the performance with other methods, we extract the maximum value along the time delay 0∼1 ms at each direction to construct a one-dimensional spatial spectrum for the proposed method. Results in Figure3 show that (a) all the methods can estimate the DOA of the direct path signal with SNR = 10 dB, the OMP-LP generates the sharpest peak at the direct path. (b) Both the MUSIC and the OMP-LP and the proposed method can estimate the DOA of path-1 with SNR = –20 dB, and the proposed method generates the sharpest peak. (c) Only the proposed method generates a peak at the DOA of path-2 with SNR = –30 dB.Figure 2
DOA-time-shift spectrum.Figure 3
Results comparison for the four methods.Example 3.
The RMSE comparison.
The SNR of path-2 varies from −40 dB to 0 dB. The root mean square error (RMSE) is calculated for path-2,RMSE=1/1000∑i=11000θ^2i−θ22, where θ2 and θ^2i are the actual and the ith time estimated DOA of path-2, respectively, and θ^2i is obtained by searching the maximum value in 110° : 0.01° : 130°. 1000 independent runs are performed. Results in Figure 4 show that (a) benefiting from the desired signal is enhanced; the CRB of EAS is lower than the CRB of OAS, and as the SNR becomes lower, the CRB-EAS gets much lower than the CRB-OAS. (b) The RMSE of the proposed method is lower than other methods, especially at much lower SNR scenario. (c) The RMSE of the proposed method is lower than the CRB of OAS.
The SNR of path-2 is fixed at −20 dB, while the number of snapshotsK varies from 5000 to 50,000, which indicates K2 = 50, and K1 varies from 100 to 1000. Results in Figure 5 show that (a) the CRB of EAS is lower than the CRB of OAS. (b) As the number of snapshots increase, the RMSE of the proposed method gets lower, and when the number of subsnapshots K1 surplus 200, the performance of the proposed method is stable. (c) The RMSE of the proposed method is lower than other methods and is lower than the CRB of OAS when K1 surplus 150.Figure 4
RMSE versus SNR.Figure 5
RMSE versus snapshot number.


## 4.3. The Same DOA and Coherent Signals with Ultralow SNRs

In Examples 4 and 5, the multipath signals’ parameters are listed in Table5. Other parameters are: K = 500000, K1 = 10000, and K2 = 50. The signals in path-4 and path-5 have the same DOA, and the signals in path-3 and path-4 are coherent.Table 5
Simulated parameters for Examples 4 and 5.
Multipath
Path-0
Path-1
Path-2
Path-3
Path-4
Path-5
Path-6
Path-7
DOA (degree)
30
70
130
120
100
100
140
80
SNR (dB)
10
−36
−40
−37
−35
−38
−34
−39
Time delay (ms)
0
0.2
0.1
0.4
0.4
0.6
0.8
0.9Example 4.
The proposed method applied in ultralow SNR
To our knowledge, the existing DOA estimation methods cannot work at ultralow SNR scenario, such as the SNR is below −40 dB, and there exists a strong direct path at the same time. Figure6 shows the DOA-time-shift spectrum calculated by the proposed method. Results show that (a) the DOAs and the associated time delays of path-1, -2, -5, -6, and -7 can be distinguished clearly. (b) Path-3 and path-4 cannot be distinguished because they are coherent.Figure 6
DOA-time-shift spectrum.Example 5.
The proposed method with spatial smoothing.
Then, we apply the MUSIC approach with forward spatial smoothing technology [9] in the proposed method. The number of subarrays is set to 2, which means the dimension of subarray is M−1=15. The dimension of noise subspace is set to 15 − 3 = 12. Results in Figure 7 show that the DOAs as well as their associated time delays of two coherent signals path-3 and path-4 can be distinguished correctly and clearly.Figure 7
DOA-time-shift spectrum with spatial smoothing.


## 5. Conclusion

A novel method was developed to estimate the DOAs of multipath signals with ultralow SNRs. By calculating the cross-correlation function of the direct path signal and the received OAS and then constructing an EAS, the desired multipath signal can be enhanced. The DOA of desired multipath signal can be estimated by the MUSIC approach with the EAS. Theoretical analysis and simulation results indicate that the SNR gain of the desired multipath signal equals to the snapshots number, so the proposed method is able to estimate the DOAs of multipath signals with ultralow SNRs, if the number of snapshots is large enough. Coherent signals with ultralow SNRs can also be measured with spatial smoothing technology by the proposed method. Furthermore, the CRB can be reduced in the EAS.



---

*Source: 2837315-2019-10-13.xml*