# Joint TOA and DOA Estimation for UWB Systems with Antenna Array via Doubled Frequency Sample Points and Extended Number of Clusters

Authors: Liangliang Gong; Yang Hu; Junyao Zhang; Gaofeng Zhao
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The ultra-wideband (UWB) system, which transmits information using nanosecond or even sub-nanosecond pulses, has been widely applied in precise positioning. In this paper, we investigate the problem of the time of arrival (TOA) estimation and the direction of arrival (DOA) estimation in the UWB systems with antenna array and propose a joint TOA and DOA estimation algorithm with doubled frequency sample points and extended number of clusters. Specifically, the proposed algorithm uses two antennas to receive impinging signals and utilizes the conjugate symmetry characteristic of the delay matrices to extend the sample points as well as the number of clusters. Moreover, in order to obtain TOA estimates with low computational complexity, the proposed algorithm transforms the two-dimensional (2D) spectral search to one-dimensional (1D) searches. The DOA estimates can then be achieved by using the TOA estimation results and the geometric information. Simulation results are given to testify the performance of the proposed algorithm.


---

## Body

## 1. Introduction

The ultra-wideband (UWB) technique is a kind of wireless communication technology which uses nanosecond or even sub-nanosecond pulses as carrier to transmit information. Benefiting from the high transmission rate, low power consumption, and the anti-multipath characteristics, the UWB system has been widely applied in various fields such as radar, imaging, and positioning [1–3]. Due to the extremely narrow pulse, the positioning precision of the UWB system can reach centimeter level or even millimeter level. Even under complex multipath conditions such as indoor environment, the UWB system can still achieve accurate positioning due to its strong anti-multipath and penetration ability [4–6]. At present, IEEE 802.15.4a standard has taken UWB as the preferred technology for positioning application [7].One of the basic problems in the UWB positioning system is the time of arrival (TOA) estimation. TOA estimation algorithms can be divided into two categories: one is the traditional algorithms based on the time domain, and the other is the high-resolution algorithms based on the frequency domain. The former mainly includes the coherent detection method using pulse template matched filter [8] and the incoherent TOA estimation algorithm based on threshold or energy detection [9–11]. The coherent algorithm based on matched filter can obtain TOA estimates with high accuracy but at the same time with high sampling rate, complex receiver structure, and expensive equipment cost. The incoherent TOA estimation algorithm has the advantages of low sampling rate, fast convergence speed, and low hardware resource occupation rate. However, the low sampling rate leads to low time resolution, which reduces the accuracy of the TOA estimation. The traditional time-domain-based algorithms obtain TOA estimates by estimating the arrival time of the direct path (DP) component in the received signal. However, due to the multipath effect and the nonline-of-sight (NLOS) condition, DP may not be the strongest path; thus the algorithm resolution declines. Therefore, super-resolution estimation algorithms based on frequency-domain processing are proposed [12–18]. These algorithms, such as the multiple signal classification (MUSIC) algorithm [12, 13], the propagator method (PM) [14], the estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm [15, 16], and the matrix pencil algorithm [17, 18], model the channel impulse response in frequency domain and realize the TOA estimation using the orthogonality between the signal subspace and the noise subspace, which can achieve high estimation resolution.Researches on DOA estimation in UWB systems were also done. Since the UWB system has high time resolution, some joint TOA and DOA estimation algorithms were proposed. The matrix pencil algorithm was applied to the joint estimation in [19]. In [20], the rough TOA estimates were obtained via energy estimation and the minimum distance criterion, the accurate estimates were then achieved by a low-complexity and high-resolution method based on the signal power delay spectrum, and the DOA estimates were finally obtained by the minimum-variance unbiased estimation using the TOA estimation results. Besides, the electromagnetic vector sensor (EMVS) with multiple-input multiple-output (MIMO) radar [21] also has the potential to be exploited in the UWB systems.In this paper, we investigate the problem of joint TOA and DOA estimation in UWB systems and propose a computationally efficient algorithm. Specifically, the proposed algorithm uses two antennas to receive impinging signals and utilizes the conjugate symmetry characteristic of the delay matrices to double the equivalent frequency sample points and extend the number of clusters. Moreover, to obtain TOA estimates with low computational complexity, the proposed algorithm transforms the two-dimensional (2D) spectral search to one-dimensional (1D) searches. The DOA estimates can then be achieved by using the TOA estimation results and the geometric information. The main contributions of this paper are summarized as follows:(1)
We propose an effective algorithm for UWB systems to jointly obtain the DOA and TOA estimation results.(2)
We transform the time-consuming 2D spectral search to twice 1D searches in order to release the computational burden.(3)
The proposed algorithm can obtain high estimation accuracy due to the doubled sample points and the extended number of clusters.The rest of this paper is organized as follows: Section2 introduces the data model of the UWB system. The proposed joint TOA and DOA estimation algorithm is derived in Section 3. Section 4 then analyses the estimation performance of the proposed algorithm and Section 5 gives the simulation results to testify the algorithm effectiveness. Finally, Section 6 concludes this paper.Notation 1.
In this paper, we use upper (lower) bold characters to represent matrices (vectors).·∗, ·T, and ·H, respectively, denote the operation of conjugate, transpose, and conjugate transpose. ∗ denotes the convolution operation, and ⊗ denotes the Kronecker product. A:,m:n generates a new matrix consisted by the m-th to the n-th columns of the matrix A. E· represents the expectation operation. IM stands for an M×M identity matrix.


## 2. Data Model

### 2.1. Transmit Signals of UWB Systems

In this paper, we utilize the second derivation of the Gaussian pulse as the UWB transmit signal and modulate it with the direct sequence binary phase shift keying (DS-BPSK). The transmit signal of the UWB system can then be expressed as [22](1)st=∑γ=−∞+∞∑n=0Nc−1bγcnpt−γTs−nTc,where bγ∈−1,+1 is the modulated binary data symbol sequence, cn∈−1,+1 is the pseudorandom sequence used to realize the multiple access communication, Tc denotes the repeat period of the pulse, Ts represents the period of the binary data symbol, and Nc is the pulse repetition times of a single binary data symbol. The term pt is the second derivation of the Gaussian pulse and can be further expressed as(2)pt=e−2πt2/Γ21−4πt2/Γ2,where Γ is the pulse forming factor related to the pulse width.


### 2.2. Channel Model of UWB Systems

According to the Saleh–Valenzuela (SV) model [23], we consider that the transmit signal produces multiple multipath components after passing through the channel, and these multipath components arrive at the receiver in the form of clusters. Specifically, we assume there are K clusters and L multipath per cluster in the UWB channel, then the channel impulse response of the kth cluster is given by(3)hkt=∑l=1Lαlkejθlkδt−τlk,where αlk is the channel fading factor of the lth path in the kth cluster, which obeys the Rayleigh distribution. θlk distributes uniformly in the range 0,2π, δ· is the Dirac function, and τlk is the channel delay of the lth path in the kth cluster. Normally, the change rate of the channel is slow compared with the pulse speed of the transmit signal; hence we have τlk=τl. Define βlk=αlkejθlk to denote the random complex fading amplitude, then (3) can be rewritten as(4)hkt=∑l=1Lβlkδt−τl.


### 2.3. Received Signals of UWB Systems

According to the basic theory of the digital signal processing, the received signal of thekth cluster in the time domain can be expressed as [23](5)ykt=st∗hkt+wkt=∑l=1L∑γ=−∞+∞∑n=1Nc−1βlkbγcnpt−γTs−nTc−τl+wkt,where wkt is the additive Gaussian white noise of the kth cluster. Transform (5) into frequency domain, i.e.,(6)Ykω=SωHkω+Wkω=∑l=1LβlkSωe−jωτl+Wkω,where Ykω, Sω, Hkω, and Wkω, respectively, represent the Fourier transformation of ykt, st, hkt, and wkt. Subsequently, by sampling NN>L points with interval Δω=2π/N in the frequency domain, we can obtain the measurement data vector as [23](7)yk=SEτβk+wk,where vector yk=Ykω0,…,YkωN−1T with ωn=nΔωn=0,1,…,N−1. S=diagSω0,…,SωN−1 is an N×N diagonal matrix. Eτ=eτ1,eτ2,…,,eτL is the delay matrix which contains the information of the signal multipath delay, where eτi=1,e−jΔωτi,…,e−jN−1ΔωτiT. Moreover, βk=β1k,β2k,…,βLkT contains the complex channel fading factor of the kth cluster and wk=Wkω0,Wkω1,…,WkωN−1T is the observe vector of the noise.In this paper, we use two array antennas to receive the impinging signals. As shown in Figure1, we assume there are L far-field signals from direction θ1,θ2,…,θL. The TOA of the two antennas is τ=τ1,τ2,…,τL and ς=ς1,ς2,…,ςL, respectively. The distance between the two antennas is d, and c denotes the light speed. According to (7), the received signal of the two antennas in the frequency domain can be expressed as [23](8)Y1=SE1τB+W1,Y2=SE2ςB+W2,where B=β1β2…βK, W1=w11w21…wK1, and W2=w12w22…wK2. E1τ and E2ς are the delay matrices of the two antennas, respectively, which can be expressed as(9)E1τ=e1τ1e1τ2…e1τL=11⋯1e−jΔωτ1e−jΔωτ2⋯e−jΔωτL⋮⋮⋱⋮e−jN−1Δωτ1e−jN−1Δωτ2⋯e−jN−1ΔωτL,(10)E2ς=e2ς1e2ς2…e2ςL=11⋯1e−jΔως1e−jΔως2⋯e−jΔωςL⋮⋮⋱⋮e−jN−1Δως1e−jN−1Δως2⋯e−jN−1ΔωςL.Figure 1
Array model.The estimation of the two channel impulse responses in the frequency domain can be achieved by(11)H^1=Y1S=E1τB+V1,H^2=Y2S=E2ςB+V2,where V1=W1/S and V2=W2/S. In the next section, we propose a joint TOA and DOA estimation algorithm utilizing the channel impulse response H^1 and H^2.


## 2.1. Transmit Signals of UWB Systems

In this paper, we utilize the second derivation of the Gaussian pulse as the UWB transmit signal and modulate it with the direct sequence binary phase shift keying (DS-BPSK). The transmit signal of the UWB system can then be expressed as [22](1)st=∑γ=−∞+∞∑n=0Nc−1bγcnpt−γTs−nTc,where bγ∈−1,+1 is the modulated binary data symbol sequence, cn∈−1,+1 is the pseudorandom sequence used to realize the multiple access communication, Tc denotes the repeat period of the pulse, Ts represents the period of the binary data symbol, and Nc is the pulse repetition times of a single binary data symbol. The term pt is the second derivation of the Gaussian pulse and can be further expressed as(2)pt=e−2πt2/Γ21−4πt2/Γ2,where Γ is the pulse forming factor related to the pulse width.


## 2.2. Channel Model of UWB Systems

According to the Saleh–Valenzuela (SV) model [23], we consider that the transmit signal produces multiple multipath components after passing through the channel, and these multipath components arrive at the receiver in the form of clusters. Specifically, we assume there are K clusters and L multipath per cluster in the UWB channel, then the channel impulse response of the kth cluster is given by(3)hkt=∑l=1Lαlkejθlkδt−τlk,where αlk is the channel fading factor of the lth path in the kth cluster, which obeys the Rayleigh distribution. θlk distributes uniformly in the range 0,2π, δ· is the Dirac function, and τlk is the channel delay of the lth path in the kth cluster. Normally, the change rate of the channel is slow compared with the pulse speed of the transmit signal; hence we have τlk=τl. Define βlk=αlkejθlk to denote the random complex fading amplitude, then (3) can be rewritten as(4)hkt=∑l=1Lβlkδt−τl.


## 2.3. Received Signals of UWB Systems

According to the basic theory of the digital signal processing, the received signal of thekth cluster in the time domain can be expressed as [23](5)ykt=st∗hkt+wkt=∑l=1L∑γ=−∞+∞∑n=1Nc−1βlkbγcnpt−γTs−nTc−τl+wkt,where wkt is the additive Gaussian white noise of the kth cluster. Transform (5) into frequency domain, i.e.,(6)Ykω=SωHkω+Wkω=∑l=1LβlkSωe−jωτl+Wkω,where Ykω, Sω, Hkω, and Wkω, respectively, represent the Fourier transformation of ykt, st, hkt, and wkt. Subsequently, by sampling NN>L points with interval Δω=2π/N in the frequency domain, we can obtain the measurement data vector as [23](7)yk=SEτβk+wk,where vector yk=Ykω0,…,YkωN−1T with ωn=nΔωn=0,1,…,N−1. S=diagSω0,…,SωN−1 is an N×N diagonal matrix. Eτ=eτ1,eτ2,…,,eτL is the delay matrix which contains the information of the signal multipath delay, where eτi=1,e−jΔωτi,…,e−jN−1ΔωτiT. Moreover, βk=β1k,β2k,…,βLkT contains the complex channel fading factor of the kth cluster and wk=Wkω0,Wkω1,…,WkωN−1T is the observe vector of the noise.In this paper, we use two array antennas to receive the impinging signals. As shown in Figure1, we assume there are L far-field signals from direction θ1,θ2,…,θL. The TOA of the two antennas is τ=τ1,τ2,…,τL and ς=ς1,ς2,…,ςL, respectively. The distance between the two antennas is d, and c denotes the light speed. According to (7), the received signal of the two antennas in the frequency domain can be expressed as [23](8)Y1=SE1τB+W1,Y2=SE2ςB+W2,where B=β1β2…βK, W1=w11w21…wK1, and W2=w12w22…wK2. E1τ and E2ς are the delay matrices of the two antennas, respectively, which can be expressed as(9)E1τ=e1τ1e1τ2…e1τL=11⋯1e−jΔωτ1e−jΔωτ2⋯e−jΔωτL⋮⋮⋱⋮e−jN−1Δωτ1e−jN−1Δωτ2⋯e−jN−1ΔωτL,(10)E2ς=e2ς1e2ς2…e2ςL=11⋯1e−jΔως1e−jΔως2⋯e−jΔωςL⋮⋮⋱⋮e−jN−1Δως1e−jN−1Δως2⋯e−jN−1ΔωςL.Figure 1
Array model.The estimation of the two channel impulse responses in the frequency domain can be achieved by(11)H^1=Y1S=E1τB+V1,H^2=Y2S=E2ςB+V2,where V1=W1/S and V2=W2/S. In the next section, we propose a joint TOA and DOA estimation algorithm utilizing the channel impulse response H^1 and H^2.


## 3. Proposed Algorithm

### 3.1. Sample Points and Clusters Extension

The cross-correlation matrix of the channel impulse response can be constructed as(12)RH=EH1H2H≈E1τRBE2Hς,where RB=EBBH∈ℂL×L is a diagonal matrix. We divide the cross-correlation matrix RH∈ℂN×N into two N×N−1 matrices, i.e.,(13)X1=RH:,1:N−1=E1τRBE2aHς,(14)X2=RH:,2:N=E1τRBE2bHς,where X1 and X2 consist of the first and the last N−1 columns of RH and E2aς and E2bς consist of the first and the last N−1 rows of E2ς. According to the form of the delay matrix, we have(15)E2bς=E2aςΨς,where Ψς is an L×L diagonal matrix given by(16)Ψς=e−jΔως1e−jΔως2⋱e−jΔωςL.According to (9) and (10), the two delay matrices of the two antennas are Vandermonde matrices and have the characteristic of conjugate symmetry, i.e.,(17)JNE1∗τ=E1τΦ˜τ,JNE2∗ς=E2ςΦ˜ς,where JN∈ℂN×N is the antidiagonal matrix, and the rotation matrices Φ˜τ and Φ˜ς are, respectively, given by(18)Φ˜τ=e−jΔωτ1e−jΔωτ2e−jΔωτL,Φ˜ς=e−jΔως1e−jΔως2e−jΔωςL.By utilizing the above properties of the delay matrices, we can construct a new matrix usingX1 and X2 as(19)X=X1,JNX2∗X2,JNX1∗=Eeτ,ςBeτ,ς,where(20)Eeτ,ς=E1τE1τΨ∗ς∈ℂ2N×L,Beτ,ς=RBE2aHςΦ˜τRBΨςE2aTς∈ℂL×2N−1.Equation (19) can be seen as the equivalent channel impulse response with doubled sample points and increased number of clusters, which can improve the maximum number of detectable signals as well as the estimation performance.


### 3.2. Reduced-Dimension TOA Estimation

The correlation matrix of the extended observation matrixX can be constructed as RX=XXH. By applying the eigenvalue decomposition, the correlation matrix can be decomposed as(21)RX=UsΛsUsH+UvΛvUvH,where Us and Uv are the signal subspace and the noise subspace, respectively, and Λs=diagλ1λ2…λL and Λv=diagλL+1λL+2…λ2N are diagonal matrices consisting of the largest L eigenvalues and the smallest 2N−L eigenvalues of RX.Similar to the classical MUSIC algorithm, we can construct the 2D-MUSIC spectrum as(22)f2Dτ,ς=1eeHτ,ςUvUvHeeτ,ς,where(23)eeτ,ς=e1τe1τejΔως.Apparently, (22) needs time-consuming 2D spectral search to obtain TOA estimates. Aiming to reduce the computational complexity, we split eeτ,ς as(24)eeτ,ς=1ejΔως⊗e1τ=I2⊗e1τqς,where qς=1,ejΔωςT. Substituting (24) into (23), the spectrum function can be rewritten as(25)f2Dτ,ς=1qHςFτqς,where Fτ=I2⊗e1τHUvUvHI2⊗e1τ and the vector qς satisfies uHqς=1 with u=1,0T. Equation (25) can be regarded as the following optimization problem:(26)minτ,ςqHςFτqςs.t.uHqς=1.According to (26), we can construct the cost function as(27)Lτ,ς=qHςFτqς−ρuHqς−1,where ρ is a constant value. In order to get the extremum, we can construct the partial derivation of Lτ,ς with respect to qς, i.e.,(28)∂Lτ,ς∂qς=2Fτqς+ρu=0.Thus, we haveqς=μF−1τu with μ=−0.5ρ. Considering uHqς=1, the constant μ can be further expressed as(29)μ=1uHF−1τu.Therefore, the vectorqς can be further transformed into(30)qς=F−1τuuHF−1τu.Substituting (30) into (26), then the TOA estimation result of the first antenna is given by(31)τ^=argminτ1uHF−1τu=argmaxτuHF−1τu,which means we can get the TOA estimation of the first antenna by a 1D spectral search with the spectral function(32)fτ=uHI2⊗e1τHUvUvHI2⊗e1τ−1u.Similarly, in order to obtain the TOA estimates of the second antenna, we can exchange the order ofH1 and H2 when constructing the cross-correlation matrix, i.e.,(33)RH=EH2H1H≈E2ςRBE1Hτ.Then, by following the same procedures described above, the 1D spectral function with respect to the TOA of the second antenna is given by(34)fς=uHI2⊗e2ςHUvUvHI2⊗e2ς−1u.After obtaining the TOA estimates of both antennas, we can estimate the DOA of the impinging signals by using the TOA estimates as well as the geometric information,i.e.,(35)θ^l=arcsincτ^l−ς^ld,l=1,2,…,L.The main steps of the proposed algorithm are summarized as follows:(1)
Construct the cross-correlation matrixRH according to (12).(2)
DivideRH into X1 and X2 according to (13) and (14), and construct the extended observation matrix X as (19).(3)
Compute the correlation matrix ofX, and perform the eigenvalue decomposition to obtain the noise subspace Uv.(4)
Construct the spectral function according to (32) and perform the 1D spectral search to obtain the TOA estimates τ^.(5)
ReconstructRH using (33) and repeat Step 2 and Step 3. The TOA estimates ς^ can be obtained by performing 1D spectral search on (34).(6)
Calculate the DOA estimates according to (35).


## 3.1. Sample Points and Clusters Extension

The cross-correlation matrix of the channel impulse response can be constructed as(12)RH=EH1H2H≈E1τRBE2Hς,where RB=EBBH∈ℂL×L is a diagonal matrix. We divide the cross-correlation matrix RH∈ℂN×N into two N×N−1 matrices, i.e.,(13)X1=RH:,1:N−1=E1τRBE2aHς,(14)X2=RH:,2:N=E1τRBE2bHς,where X1 and X2 consist of the first and the last N−1 columns of RH and E2aς and E2bς consist of the first and the last N−1 rows of E2ς. According to the form of the delay matrix, we have(15)E2bς=E2aςΨς,where Ψς is an L×L diagonal matrix given by(16)Ψς=e−jΔως1e−jΔως2⋱e−jΔωςL.According to (9) and (10), the two delay matrices of the two antennas are Vandermonde matrices and have the characteristic of conjugate symmetry, i.e.,(17)JNE1∗τ=E1τΦ˜τ,JNE2∗ς=E2ςΦ˜ς,where JN∈ℂN×N is the antidiagonal matrix, and the rotation matrices Φ˜τ and Φ˜ς are, respectively, given by(18)Φ˜τ=e−jΔωτ1e−jΔωτ2e−jΔωτL,Φ˜ς=e−jΔως1e−jΔως2e−jΔωςL.By utilizing the above properties of the delay matrices, we can construct a new matrix usingX1 and X2 as(19)X=X1,JNX2∗X2,JNX1∗=Eeτ,ςBeτ,ς,where(20)Eeτ,ς=E1τE1τΨ∗ς∈ℂ2N×L,Beτ,ς=RBE2aHςΦ˜τRBΨςE2aTς∈ℂL×2N−1.Equation (19) can be seen as the equivalent channel impulse response with doubled sample points and increased number of clusters, which can improve the maximum number of detectable signals as well as the estimation performance.


## 3.2. Reduced-Dimension TOA Estimation

The correlation matrix of the extended observation matrixX can be constructed as RX=XXH. By applying the eigenvalue decomposition, the correlation matrix can be decomposed as(21)RX=UsΛsUsH+UvΛvUvH,where Us and Uv are the signal subspace and the noise subspace, respectively, and Λs=diagλ1λ2…λL and Λv=diagλL+1λL+2…λ2N are diagonal matrices consisting of the largest L eigenvalues and the smallest 2N−L eigenvalues of RX.Similar to the classical MUSIC algorithm, we can construct the 2D-MUSIC spectrum as(22)f2Dτ,ς=1eeHτ,ςUvUvHeeτ,ς,where(23)eeτ,ς=e1τe1τejΔως.Apparently, (22) needs time-consuming 2D spectral search to obtain TOA estimates. Aiming to reduce the computational complexity, we split eeτ,ς as(24)eeτ,ς=1ejΔως⊗e1τ=I2⊗e1τqς,where qς=1,ejΔωςT. Substituting (24) into (23), the spectrum function can be rewritten as(25)f2Dτ,ς=1qHςFτqς,where Fτ=I2⊗e1τHUvUvHI2⊗e1τ and the vector qς satisfies uHqς=1 with u=1,0T. Equation (25) can be regarded as the following optimization problem:(26)minτ,ςqHςFτqςs.t.uHqς=1.According to (26), we can construct the cost function as(27)Lτ,ς=qHςFτqς−ρuHqς−1,where ρ is a constant value. In order to get the extremum, we can construct the partial derivation of Lτ,ς with respect to qς, i.e.,(28)∂Lτ,ς∂qς=2Fτqς+ρu=0.Thus, we haveqς=μF−1τu with μ=−0.5ρ. Considering uHqς=1, the constant μ can be further expressed as(29)μ=1uHF−1τu.Therefore, the vectorqς can be further transformed into(30)qς=F−1τuuHF−1τu.Substituting (30) into (26), then the TOA estimation result of the first antenna is given by(31)τ^=argminτ1uHF−1τu=argmaxτuHF−1τu,which means we can get the TOA estimation of the first antenna by a 1D spectral search with the spectral function(32)fτ=uHI2⊗e1τHUvUvHI2⊗e1τ−1u.Similarly, in order to obtain the TOA estimates of the second antenna, we can exchange the order ofH1 and H2 when constructing the cross-correlation matrix, i.e.,(33)RH=EH2H1H≈E2ςRBE1Hτ.Then, by following the same procedures described above, the 1D spectral function with respect to the TOA of the second antenna is given by(34)fς=uHI2⊗e2ςHUvUvHI2⊗e2ς−1u.After obtaining the TOA estimates of both antennas, we can estimate the DOA of the impinging signals by using the TOA estimates as well as the geometric information,i.e.,(35)θ^l=arcsincτ^l−ς^ld,l=1,2,…,L.The main steps of the proposed algorithm are summarized as follows:(1)
Construct the cross-correlation matrixRH according to (12).(2)
DivideRH into X1 and X2 according to (13) and (14), and construct the extended observation matrix X as (19).(3)
Compute the correlation matrix ofX, and perform the eigenvalue decomposition to obtain the noise subspace Uv.(4)
Construct the spectral function according to (32) and perform the 1D spectral search to obtain the TOA estimates τ^.(5)
ReconstructRH using (33) and repeat Step 2 and Step 3. The TOA estimates ς^ can be obtained by performing 1D spectral search on (34).(6)
Calculate the DOA estimates according to (35).


## 4. Performance Analysis

### 4.1. Complexity

According to Section3, when estimating the TOA of the first antenna τ, constructing RH needs ON2K, computing RX needs O8N2N−1, and performing eigenvalue decomposition on it requires O8N3. Constructing the spectral function and conducting 1D spectral search needs O4N22N−L+8nsN2+N+1, where ns denotes the search times. Therefore, the total complexity of the proposed algorithm is O2N2K+48N3−8N2L+2+16nsN2+N+1. If we perform 2D spectral search on equation (22), then the complexity is increased to ON2K+24N3−4N2L+2+2ns2N2N+1. The computational complexity of the proposed algorithm and some other algorithms are concluded in Table 1. Besides, the complexity comparison versus the number of frequency sample points of different algorithms is depicted in Figure 2, where K=100, L=3, and ns=100. As shown in Figure 2, the complexity of the proposed algorithm is remarkably reduced compared with the algorithms with 2D spectral search and is close to the complexity of the 2D-PM.Table 1
Complexity comparison of different algorithms.
AlgorithmComplexityProposedO2N2K+48N3−8N2L+2+16nsN2+N+12D proposedON2K+24N3−4N2L+2+2ns2N2N+1Successive PM [24]6K+4L+1N2+3L2N+3L3+4L+2N−1+3L3+22N−L3+4N2N−L2+2ns2N+12N−L2D-PM [14]6K+4L+1N2+2L2N+L3+22N−L3+4N2N−L2+ns22N+12N−LMatrix pencil [19]2K+1N2+6KL2+2KN−2KL3Figure 2
Complexity comparison of different algorithms, whereK=100, L=3, and ns=100.


### 4.2. Algorithm Advantages

We conclude the advantages of the proposed algorithm as follows:(1)
The proposed algorithm reduces the computational complexity by transforming the 2D spectral search to twice 1D spectral search.(2)
The proposed algorithm doubles the equivalent frequency sample points and increases the equivalent number of clusters, which can improve the maximum number of detectable signals as well as the estimation accuracy.(3)
The proposed algorithm can obtain higher estimation accuracy compared with some existed algorithms.


## 4.1. Complexity

According to Section3, when estimating the TOA of the first antenna τ, constructing RH needs ON2K, computing RX needs O8N2N−1, and performing eigenvalue decomposition on it requires O8N3. Constructing the spectral function and conducting 1D spectral search needs O4N22N−L+8nsN2+N+1, where ns denotes the search times. Therefore, the total complexity of the proposed algorithm is O2N2K+48N3−8N2L+2+16nsN2+N+1. If we perform 2D spectral search on equation (22), then the complexity is increased to ON2K+24N3−4N2L+2+2ns2N2N+1. The computational complexity of the proposed algorithm and some other algorithms are concluded in Table 1. Besides, the complexity comparison versus the number of frequency sample points of different algorithms is depicted in Figure 2, where K=100, L=3, and ns=100. As shown in Figure 2, the complexity of the proposed algorithm is remarkably reduced compared with the algorithms with 2D spectral search and is close to the complexity of the 2D-PM.Table 1
Complexity comparison of different algorithms.
AlgorithmComplexityProposedO2N2K+48N3−8N2L+2+16nsN2+N+12D proposedON2K+24N3−4N2L+2+2ns2N2N+1Successive PM [24]6K+4L+1N2+3L2N+3L3+4L+2N−1+3L3+22N−L3+4N2N−L2+2ns2N+12N−L2D-PM [14]6K+4L+1N2+2L2N+L3+22N−L3+4N2N−L2+ns22N+12N−LMatrix pencil [19]2K+1N2+6KL2+2KN−2KL3Figure 2
Complexity comparison of different algorithms, whereK=100, L=3, and ns=100.


## 4.2. Algorithm Advantages

We conclude the advantages of the proposed algorithm as follows:(1)
The proposed algorithm reduces the computational complexity by transforming the 2D spectral search to twice 1D spectral search.(2)
The proposed algorithm doubles the equivalent frequency sample points and increases the equivalent number of clusters, which can improve the maximum number of detectable signals as well as the estimation accuracy.(3)
The proposed algorithm can obtain higher estimation accuracy compared with some existed algorithms.


## 5. Simulation Results

In this section, we perform large number of Monte Carlo trials to examine the effectiveness as well as the superiority of the proposed algorithm. The root mean square error (RMSE) of the trials is used to measure the estimation accuracy, which is given by(36)RMSETOA=1L∑l=1L1Q∑q=1Qτl−τ^l,q2+ςl−ς^l,q2,RMSEDOA=1L∑l=1L1Q∑q=1Qθl−θ^l,q2,where Q is the number of Monte Carlo trials and τ^l,q, ς^l,q, and θ^l,q are the estimation result of τl, ςl, and θl in the qth trial, respectively. In the following simulations, we set Q=500, d=0.09m, and L=3 with τ=0.2ns,0.3ns,0.4ns.Figures3 and 4 show the TOA estimation results of the proposed algorithm, where each point represents a trial, K=100, N=64, and the signal-to-noise ratio (SNR) is −10dB and 10dB, respectively. As shown in the two figures, the proposed algorithm can successfully obtain the paired TOA estimation results and the estimation accuracy is improved with higher SNR.Figure 3
Estimation results of the proposed algorithm, whereSNR=−10dB, K=100, and N=64.Figure 4
Estimation results of the proposed algorithm, whereSNR=10dB, K=100, and N=64.Figures5 and 6, respectively, compare the TOA and DOA estimation accuracy with different number of clusters, where N=64. The simulation results show that the estimation accuracy of both TOA and DOA improves with increased number of clusters.Figure 5
TOA estimation performance of the proposed algorithm in different number of clusters, whereN=64.Figure 6
DOA estimation performance of the proposed algorithm in different number of clusters, whereN=64.Figures7 and 8, respectively, depict the TOA and the DOA estimation performance comparison versus SNR of different algorithms, where K=100 and N=64. It is illustrated clearly that the proposed algorithm can achieve better estimation performance in both TOA and DOA domains compared with the ESPRIT algorithm [16], the successive PM [24], the matrix pencil algorithm [19], and the 2D-PM [14].Figure 7
TOA estimation performance comparison of different algorithms, whereK=100 and N=64.Figure 8
DOA estimation performance comparison of different algorithms, whereK=100 and N=64.


## 6. Conclusion

In this paper, we investigate the problem of joint TOA and DOA estimation in UWB systems and propose a computationally efficient algorithm with doubled frequency sample points and extended number of clusters. The proposed algorithm utilizes the conjugate symmetry characteristic of the delay matrices to extend the sample points as well as the number of clusters and then transforms the 2D spectral search into twice 1D search in order to reduce the computational complexity. Simulations testify that the proposed algorithm can obtain high estimation accuracy and large number of identifiable signals with low complexity.



---

*Source: 7521573-2021-06-07.xml*