# Direct Signal Space Construction Luenberger Tracker with Quadratic Least Square Regression for Computationally Efficient DoA Tracking

Authors: Bukeun Byeon; Ho-Kyoung Lee; Do-Sik Yoo
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, we propose a computationally efficient direction-of-arrival (DoA) tracking scheme called the direct signal space construction Luenberger tracker (DSPCLT) with quadratic least square (QLS) regression. Also, we study analytically how the choice of observer gain affects the algorithm performance and illustrate how we can use the theoretical results in determining optimal observer gain value. The proposed scheme (DSPCLT) has several distinct features compared with existing algorithms. First, it requires only a fraction of computational complexity compared with other schemes. Secondly, it maintains robustness by treating separately the special case of object overlap in which subspace-based algorithms often suffer from lack of resolvability. Thirdly, the proposed scheme achieves enhanced performance by a method of delay compensation, which accounts for observation delay. Through numerical analysis, we show that DSPCLT achieves performance similar or superior to existing algorithms with only a fraction of computational requirement.


---

## Body

## 1. Introduction

In this paper, we propose a computationally efficient direction-of-arrival (DoA) tracking scheme called the direct signal space construction Luenberger tracker (DSPCLT) with quadratic least square (QLS) regression. During the last several decades, various subspace-based schemes have been proposed for DoA estimation such as multiple signal classification algorithms [1], estimation of signal parameter via rotational invariance technique (ESPRIT) [2], and root-MUSIC method [3]. Recently, because of their high resolution, these algorithms have been studied in various applications such as nonlinear arrays [4–6], 2-dimensional arrays [7, 8], and MIMO radar [9]. However, early subspace-based algorithms generally suffer from the issue of large computational complexity, typically due to the need for eigenvalue decomposition (ED) or singular-value decomposition (SVD). For this reason, various schemes have been proposed to reduce the system complexity, such as the propagator method (PM) [10], orthogonal propagator method (OPM) [11], subspace method without eigendecomposition (SWEDE) [12], direct signal space construction method (DSPCM) [13–15], and projection approximation subspace tracking-based deflation (PASTd) algorithm [16].However, these algorithms exhibit difficulty in estimating the DoAs when one or more signals suddenly disappear or when DoAs overlap. [17] As means to overcome such issues, a number of DoA tracking schemes have been proposed in [18–20]. Recently, authors in [21, 22] proposed to combine existing tracking techniques with computational efficient subspace techniques to avoid the issues of overlapping or suddenly disappearing targets while keeping the complexity minimal. For example, the authors in [21] proposed an algorithm by combining PASTd and Kalman filtering [23]. We shall refer to this algorithm as PASTd Kalman in this paper. PASTd Kalman reduces the computational complexity of constructing the subspace by avoiding ED and SVD with recursive least square (RLS) technique. However, the algorithm still requires relatively high complexity due to the adaptation of Kalman filtering. More recently, the authors in [22] proposed an algorithm called the adaptive method of estimating DoA (AMEND) which replaces the Kalman filter with Luenberger state observer [24]. However, the complexity of AMEND is not smaller than PASTd Kalman because it employs, for subspace construction, a relatively computationally expensive algorithm called the subspace-based method without eigendecomposition (SUMWE) [25].In this paper, we propose DSPCLT that achieves similar or better performance with computational complexity significantly lower than existing algorithms. The proposed scheme has several distinct features compared with existing schemes. First, it achieves very low complexity by employing DSPCM [14, 15] for subspace construction and Luenberger observer in place of the Kalman filter for estimation and filtering. Secondly, it achieves robustness by treating separately the special case of object overlap in which subspace-based algorithms suffer from the lack of resolvability. Thirdly, the proposed scheme achieves enhanced performance by a method of delay compensation, which takes into account the fact that the DoA estimation at present time is performed with observed data collected previously. We show, through numerical analysis, that the proposed scheme achieves performance similar or superior to existing algorithms with only a fraction of computational requirement.Another important topic studied in this paper is the choice of observer gains used to implement Luenberger observer. In most existing schemes employing Luenberger observer, optimal observer gain values are determined through simulations. For example, a certain observer gain value was suggested in [22] after simulations. However, we note that there is infinite variety of situations depending on the number of antenna sensor elements, the numbers of objects, the speeds of movements, the signal-to-noise ratios of the target, the shape of object trajectory, and so on. Hence, for a selection of a certain observer gain value to be accepted as a valid choice, there must be some degree of guarantee that a particular observer gain value will provide optimal or near-optimal performance in other situations, which is generally very difficult with pure simulations. For this reason, we study analytically how the choice of observer gain affects the algorithm performance and how we can choose optimal observer gain value.This paper is organized as follows. In Section2, the system model and basic assumptions are described together with the notational definitions used throughout this paper. In Sections 3∼5, the proposed scheme (DSPCLT) is delineated, which consists of three stages called initialization, prediction, and estimation. The stages of initialization and prediction are introduced in Section 3 and that of estimation is presented in Sections 4 and 5 depending on whether objects are resolvable by the subspace-based algorithm employed. In Section 6, we study theoretically how the selection of observer gain value affects the algorithm performance and describe how we can make optimal observer gain value selection using the analytical results. Next, we show, in Section 7, that the proposed scheme achieves similar or better performance with only a fraction of computational requirements compared with other algorithms. Finally, we draw conclusion in Section 8.The following notations are used throughout this paper. Bold-faced letters are used for matrices and vectors. The complex conjugation, the transposition, and the Hermitian transposition of matrixA are denoted, respectively, by A∗, AT, and AH. The M×N zero matrix and N×N identity matrix are denoted, respectively, by OM×N and IN.


## 2. System Model and Basic Assumptions

We consider a uniform linear array (ULA) consisting ofN sensors with intersensor spacing d and assume that K narrowband far-field signal waves s1t,s2t,⋯,sKt of wavelength λ impinge on the array at distinct directions, respectively, making angles θ1t,θ2t,⋯,θKt with the linear array. We denote, by x1t,x2t,⋯,xNt, the received noisy signals collected at the N sensors. Then, if we define the steering vector aθkt by aθkt=1,ejνkt,⋯,ejN−1νktT with νkt=2πdcosθkt/λ, the system can be modeled with the following equation:(1)xt=Aθ1t,⋯,θKtst+nt,where xt=x1t,x2t,⋯,xNtT, Aθ1t,⋯,θKt=aθ1t,⋯,aθKt, st=s1t,⋯,sKtT, and nt is the vector consisting of the noises added at the sensors. In this paper, we make, on the system model, the following basic assumptions:(1)
The signal wavelengthλ is known a priori and that d is chosen to be λ/2. And the number K of impinging signals is predetermined by appropriate methods such as the one in [26].(2)
The vectorsst and nt are statistically independent complex random vector processes with zero mean satisfying(2)EstssH=Pδts,EstssT=OK×K,EntnsH=N0δtsIN,EntnsT=ON×N,whereδts is the Kronecker delta function. We further assume that P is a diagonal matrix with positive real numbers p1,⋯,pK on its diagonal and that the value of N0 is assumed to be estimated accurately by the receiver through long-term observation of the noise.(3)
The system tracks the DoAs by obtaining angle estimates at timet=0,T, 2T,3T,…, where T is a predefined positive real number. We assume further that T is chosen to be small enough that the DoAs are essentially constant during the time duration n−1T,nT for each n=1,2,3,….


## 3. Initialization and Prediction

In this paper, we propose a computationally efficient DoA tracking scheme called the direct signal space construction Luenberger tracker (DSPCLT) with quadratic least square regression. As in AMEND, the DoA is tracked by sequential update of state that is defined as a three-dimensional vector consisting of direction-of-arrival and its first and second derivatives. The process of update is carried out in two stages, which we call prediction and estimation. In the prediction stage, the state is tentatively updated purely based on the previous state. Then, more accurate updated state is obtained in the estimation stage, which is the major part of the proposed scheme.As in AMEND, we propose to use, in the estimation stage, the Luenberger observer using a subspace-based algorithm. In particular, we use Newton iteration method to obtain the improved estimate of the DoA for use in the Luenberger observer. While such subspace-based tracking schemes provide excellent performance in general, we need to take special care when the DoAs of objects overlap, in which case the Newton iteration may not provide robust estimate of the DoA. Consequently, it is wise to consider such a special case separately.Since it is rather complicated to describe the estimation stage, we describe only the initialization stage and the prediction stage in this section. Then, in the next two sections, we consider the estimation stage in nonresolvable and resolvable situations, respectively. The proposed scheme employs, for each object, a 3-dimensional state vector whose components represent the direction-of-arrival and its first and second derivatives, respectively. At the initialization stage, the state vector is assumed to be obtained by a certain existing method such as MUSIC [1]. After the initialization stage, the state vector is updated at a regular interval T in two stages called prediction and estimation described in the following.Before proceeding, let us clarify the notations used in the following. First, we denote byy^knn the state vector at time t=nT and by θ^knn the first component of y^knn. We shall regard θ^knn as the estimate of the DoA θknT of object k at time nT. Next, we denote by y^knn−1 the predicted state vector obtained in the prediction stage at time t=nT and by θ^knn−1 the first component of y^knn−1.


### 3.1. Initialization Stage

The initialization stage begins by obtaining the estimates ofθk−2T, θk−T, and θk0 for k=1,2,⋯,K, through a certain DoA measurement scheme such as MUSIC [1], ESPRIT [2], or DSPCM [14, 15]. The estimates θk−2T, θk−T, and θk0 shall be denoted, respectively, by θ˜k−2−2,θ˜k−1−1, and θ˜k00. Then, the initial state vector y^k00 is computed by(3)y^k00=θ˜k00,θ˜k00−θ^k−1−1/T,θ˜k00−2θ˜k−1−1+θ˜k−2−2/T2T.We note that the third component ofy^k00 is set to be zero in AMEND [22] and PASTd Kalman [21]. In fact, this initialization provides better performance in the proposed scheme as in AMEND and PASTd Kalman according to our simulations. Consequently, we shall use this method of initialization previously used in AMEND and PASTd Kalman when we evaluate the performance in a later section. However, we shall use initialization (3) for the proposed scheme despite slightly worse performance. This is because mathematical analysis in later sections is slightly simpler with (3).


### 3.2. Prediction Stage

After initialization, the state vectors are updated periodically, with periodT, in two stages called prediction and estimation. We assume that the state vector y^kn−1n−1 at previous time n−1T is given, where n∈1,2,⋯. We assume that y^kn−1n−1 has been obtained in the initialization stage if n=1 and in the previous estimation stage if n>1. In this stage, namely, in the prediction stage, rough estimates of the next state, which we denote by y^knn−1, are obtained by the following dynamic matrix equation:(4)y^knn−1=F⋅y^kn−1n−1,where(5)F=1T0.5T201T001.Here, the first, second, and third components of the vectory^knn−1 are interpreted to represent predicted values, based on the previous state vector, for θknT, θ˙knT, and θ¨knT, respectively. After the prediction stage, the distance(6)dkl=θ^knn−1−θ^lnn−1,is computed for l=1,⋯,k−1,k+1,⋯,K. The object k is regarded to be in the nonresolvable situation if dkl is less than or equal to a certain predetermined threshold value for some i,j with l≠k. Naturally, the object k is regarded to be resolvable if it is not in the nonresolvable situation.


## 3.1. Initialization Stage

The initialization stage begins by obtaining the estimates ofθk−2T, θk−T, and θk0 for k=1,2,⋯,K, through a certain DoA measurement scheme such as MUSIC [1], ESPRIT [2], or DSPCM [14, 15]. The estimates θk−2T, θk−T, and θk0 shall be denoted, respectively, by θ˜k−2−2,θ˜k−1−1, and θ˜k00. Then, the initial state vector y^k00 is computed by(3)y^k00=θ˜k00,θ˜k00−θ^k−1−1/T,θ˜k00−2θ˜k−1−1+θ˜k−2−2/T2T.We note that the third component ofy^k00 is set to be zero in AMEND [22] and PASTd Kalman [21]. In fact, this initialization provides better performance in the proposed scheme as in AMEND and PASTd Kalman according to our simulations. Consequently, we shall use this method of initialization previously used in AMEND and PASTd Kalman when we evaluate the performance in a later section. However, we shall use initialization (3) for the proposed scheme despite slightly worse performance. This is because mathematical analysis in later sections is slightly simpler with (3).


## 3.2. Prediction Stage

After initialization, the state vectors are updated periodically, with periodT, in two stages called prediction and estimation. We assume that the state vector y^kn−1n−1 at previous time n−1T is given, where n∈1,2,⋯. We assume that y^kn−1n−1 has been obtained in the initialization stage if n=1 and in the previous estimation stage if n>1. In this stage, namely, in the prediction stage, rough estimates of the next state, which we denote by y^knn−1, are obtained by the following dynamic matrix equation:(4)y^knn−1=F⋅y^kn−1n−1,where(5)F=1T0.5T201T001.Here, the first, second, and third components of the vectory^knn−1 are interpreted to represent predicted values, based on the previous state vector, for θknT, θ˙knT, and θ¨knT, respectively. After the prediction stage, the distance(6)dkl=θ^knn−1−θ^lnn−1,is computed for l=1,⋯,k−1,k+1,⋯,K. The object k is regarded to be in the nonresolvable situation if dkl is less than or equal to a certain predetermined threshold value for some i,j with l≠k. Naturally, the object k is regarded to be resolvable if it is not in the nonresolvable situation.


## 4. Estimation Stage in the Nonresolvable Situation

We note that the nonresolvable situation becomes less likely if the number of antenna elements is increased. In other words, the nonresolvable situation does not happen very often and does not last very long even when it happens if the number of antenna elements is large enough. Although the nonresolvable situation may not happen very often and the algorithm used in this situation is not the major contribution of this paper, we consider this situation first because the proposed algorithm for the resolvable situation is much more complicated to describe.In the nonresolvable situation, the estimated state vectory^knn is obtained by using quadratic least square regression using the past p estimated DoAs θ^kmm, m=n−p,n−p+1,⋯,n−1. To be more specific, we assume that the true DoA θm satisfies(7)θkm=am−n2+bm−n+c,and attempt to obtain the constants a,b, and c by the least square method with the previously estimated DoAs θ^kmm, m=n−p,n−p+1,⋯,n−1. We note that the resultant constants a,b, and c are given by [27](8)abc=p∑m=−p−1m∑m=−p−1m2∑m=−p−1m∑m=−p−1m2∑m=−p−1m3∑m=−p−1m2∑m=−p−1m3∑m=−p−1m4−1∑m=n−pn−1θ^kmm∑m=n−pn−1m−nθkmm∑m=n−pn−1m−n2θ^kmm.Once the constantsa,b, and c are computed by (8), the estimated state vector y^knn is obtained by(9)y^knn=cb2a.We note that the second and third components ofy^knn can be obtained by the first and the second derivative of θkm=am−n2+bm−n+c at m=n.


## 5. Estimation Stage in the Resolvable Situation

In this section, we describe how we obtain the estimated state vector in the resolvable situation. First, we recall that the predicted state vectory^knn−1 is obtained, in the prediction stage, by simply multiplying the previous estimated state vector y^kn−1n−1 by a constant matrix F without reflecting the data collected between times n−1T and nT. While the data collected, captured at the antenna, are still not used to obtain the updated estimated state vector in the nonresolvable situation, the data are used to obtain the updated estimated state vector in the resolvable situation. In the following, we describe how the data are used to obtain y^knn in the resolvable situation.In particular, the updated estimated state vectory^knn is obtained, using the Luenberger observer, by(10)y^knn=y^knn−1+θ˜knn−θ^knn−1gk,where gk denotes a 3-dimensional column vector called observer gain [24]. Here, the quantity θ˜knn is an estimate of θknT obtained by a certain algorithm using the data collected at the antenna. In the following, we describe how we obtain θ˜knn in the proposed scheme. Before proceeding, we note that the performance of the system can depend on the choice of the observer gain gk. In many cases as in AMEND, no serious discussion is provided regarding how the choice of the observer gain affects the performance or how we can achieve optimal or near-optimal performance. We discuss this issue in the next section. Hence, we focus, in this section, solely on how we obtain the quantity θ˜knn, which we call the delay compensated direct signal space construction method Newton estimate (DCDSPCMNE) of θknT. The DCDSPCMNE θ˜knn is obtained in two steps. In the first step, we obtain a quantity θˇknn called the direct signal space construction method Newton estimate (DSPCMNE). Then, in the second step, we obtain more accurate estimate θ˜knn by a method of motion compensation.


### 5.1. Direct Signal Space Construction Method Newton Estimate

Let us first discuss how we obtain DSPCMNEθˇknn. In many subspace-based DoA estimation schemes, the DoAs are estimated by minimizing a certain cost function fθ defined by(11)fθ=aHθΠaθ,where Π is the orthogonal projector into the noise subspace. One major source of computational complexity in subspace-based schemes is in the construction of the orthogonal projector Π. In this paper, we propose to use the direct signal space construction method [14, 15] to reduce the computational complexity significantly. In Appendix A.1, it was described how we can construct the orthogonal projector Π in the proposed scheme.Since the minimum search process also takes some nonnegligible computational complexity, various methods have been proposed to reduce the burden. One such method, which is particularly useful when some coarse estimate of the DoA is available, is to employ Newton’s iteration method of searching zeroes withf′θ [11, 12, 22, 25]. We also adopt Newton’s iteration method in this paper. Let θ¯ denote the coarse estimate of a zero of f′θ. Then, according to Newton’s method, a refined zero θ˜ is given by(12)θ˜=θ¯−f′θ¯f″θ¯=θ¯−Rea′Hθ¯Πaθ¯Rea″Hθ¯Πaθ¯+a′Hθ¯Πa′θ¯,where Re⋅ denotes the real part of ⋅. Here, we note that the first term of the denominator is much smaller than the second term of the denominator since Πaθ¯≈0. For this reason, the first term of the denominator is open dropped in implementations.Using (12) and neglecting the first term of the denominator, we obtain DSPCMNE θˇknn by(13)θˇknn=θ^knn−1−RedHnΠ^nandHnΠ^ndn,where akn and dkn are vectors defined, respectively, by(14)akn=aθ^knn−1,and(15)dkn=daθdθθ=θ^knn−1.Also, the matrixΠ^n in (13) denotes an estimate of the propagator π. We discuss, in Appendix A.2, how we can construct the estimate Π^n in the proposed scheme.


### 5.2. Delay Compensated Direct Signal Space Construction Method Newton Estimate

We note that the orthogonal projectorΠ^n is obtained by using the data collected during the time interval n−1T,nT. Hence, while DSPCMNE θˇknn is an estimate of θkt, it is not exactly at time t=nT. In this section, we regard θˇknn as an estimate of θkt at a certain time t=n−αT (0<α<1) and obtain the DCDSPCMNE θ˜knn by compensating the movement during the delay αT. For this purpose, we first need to assess the value of α and then estimate the amount of change in θkt during this delay.Before proceeding, we shall assume thatT and θ¨kt are small enough so that θ˙kt does not change appreciably over time period n−1T,nT. However, we note, even with this simplifying assumption, that it is very difficult to determine theoretically justifiable value for αT since the process of obtaining Π^n does not depend linearly on the observation times. As described in Appendix A.2, the data are collected at times t=n−Ns−1δT,n−Ns−2δT,⋯,n−Ns−NsT, and forgetting factor ρ is used to give more weight to more recent data. In the special case in which δ=1/Ns and ρ=1, we note that the data are collected at times t=n−1+1/NsT,n−1+2/NsT,⋯,nT, and all the data are equally weighted since ρ=1. Consequently, we observe that the data are collected symmetrically around time t=n−Ns−1/2NsT and are equally used. While it is theoretically difficult to justify the validity, we shall assess the value of α as Ns−1/2Ns in this special case.The assessment of the value ofα is more difficult in the general case. In the general case, considering the fact that the data are collected at t=n−Ns−1δT,n−Ns−2δT,⋯,n−Ns−NsT and that data collected at time t=n−Ns−lδT are weighted by 1−ρρNs−l/1−ρNs in the computation of r˜m in (A.3) and (A.4), we shall regard the weighted time average(16)n−Ns−1δT⋅1−ρρNs−11−ρNs+n−Ns−2δT⋅1−ρρNs−21−ρNs+⋯+n−Ns−NsδT⋅1−ρρNs−Ns1−ρNs=n−Ns−1ρNs+1−NsρNs+ρ1−ρ1−ρNs⋅δT,as the estimated observation time n−αT.Next, we note that(17)θˇknn−θ^kn−1n−1,accounts for the angle change between times n−1T and n−αT. The change between times n−1T and nT can be estimated as(18)11−αθˇknn−θ^kn−1n−1,assuming that angular velocity does not change over n−1T,nT. Using this approximation, we obtain the DCDSPCMNE θ˜knn by the following equation:(19)θ˜knn=θ^kn−1n−1+11−αθˇknn−θ^kn−1n−1,(20)=θ^kn−1n−1+1−ρ1−ρNsθˇknn−θ^kn−1n−11−ρ1−ρNs−Ns−1ρNs+1−NsρNs+ρδ.Before proceeding, we note that (20) reduces to(21)θ˜knn=θ^kn−1n−1+NS1−ρ1−ρNsNs−Ns+1ρ+ρNs+1θˇknn−θ^kn−1n−1,in the special case of δ=1/Ns.


## 5.1. Direct Signal Space Construction Method Newton Estimate

Let us first discuss how we obtain DSPCMNEθˇknn. In many subspace-based DoA estimation schemes, the DoAs are estimated by minimizing a certain cost function fθ defined by(11)fθ=aHθΠaθ,where Π is the orthogonal projector into the noise subspace. One major source of computational complexity in subspace-based schemes is in the construction of the orthogonal projector Π. In this paper, we propose to use the direct signal space construction method [14, 15] to reduce the computational complexity significantly. In Appendix A.1, it was described how we can construct the orthogonal projector Π in the proposed scheme.Since the minimum search process also takes some nonnegligible computational complexity, various methods have been proposed to reduce the burden. One such method, which is particularly useful when some coarse estimate of the DoA is available, is to employ Newton’s iteration method of searching zeroes withf′θ [11, 12, 22, 25]. We also adopt Newton’s iteration method in this paper. Let θ¯ denote the coarse estimate of a zero of f′θ. Then, according to Newton’s method, a refined zero θ˜ is given by(12)θ˜=θ¯−f′θ¯f″θ¯=θ¯−Rea′Hθ¯Πaθ¯Rea″Hθ¯Πaθ¯+a′Hθ¯Πa′θ¯,where Re⋅ denotes the real part of ⋅. Here, we note that the first term of the denominator is much smaller than the second term of the denominator since Πaθ¯≈0. For this reason, the first term of the denominator is open dropped in implementations.Using (12) and neglecting the first term of the denominator, we obtain DSPCMNE θˇknn by(13)θˇknn=θ^knn−1−RedHnΠ^nandHnΠ^ndn,where akn and dkn are vectors defined, respectively, by(14)akn=aθ^knn−1,and(15)dkn=daθdθθ=θ^knn−1.Also, the matrixΠ^n in (13) denotes an estimate of the propagator π. We discuss, in Appendix A.2, how we can construct the estimate Π^n in the proposed scheme.


## 5.2. Delay Compensated Direct Signal Space Construction Method Newton Estimate

We note that the orthogonal projectorΠ^n is obtained by using the data collected during the time interval n−1T,nT. Hence, while DSPCMNE θˇknn is an estimate of θkt, it is not exactly at time t=nT. In this section, we regard θˇknn as an estimate of θkt at a certain time t=n−αT (0<α<1) and obtain the DCDSPCMNE θ˜knn by compensating the movement during the delay αT. For this purpose, we first need to assess the value of α and then estimate the amount of change in θkt during this delay.Before proceeding, we shall assume thatT and θ¨kt are small enough so that θ˙kt does not change appreciably over time period n−1T,nT. However, we note, even with this simplifying assumption, that it is very difficult to determine theoretically justifiable value for αT since the process of obtaining Π^n does not depend linearly on the observation times. As described in Appendix A.2, the data are collected at times t=n−Ns−1δT,n−Ns−2δT,⋯,n−Ns−NsT, and forgetting factor ρ is used to give more weight to more recent data. In the special case in which δ=1/Ns and ρ=1, we note that the data are collected at times t=n−1+1/NsT,n−1+2/NsT,⋯,nT, and all the data are equally weighted since ρ=1. Consequently, we observe that the data are collected symmetrically around time t=n−Ns−1/2NsT and are equally used. While it is theoretically difficult to justify the validity, we shall assess the value of α as Ns−1/2Ns in this special case.The assessment of the value ofα is more difficult in the general case. In the general case, considering the fact that the data are collected at t=n−Ns−1δT,n−Ns−2δT,⋯,n−Ns−NsT and that data collected at time t=n−Ns−lδT are weighted by 1−ρρNs−l/1−ρNs in the computation of r˜m in (A.3) and (A.4), we shall regard the weighted time average(16)n−Ns−1δT⋅1−ρρNs−11−ρNs+n−Ns−2δT⋅1−ρρNs−21−ρNs+⋯+n−Ns−NsδT⋅1−ρρNs−Ns1−ρNs=n−Ns−1ρNs+1−NsρNs+ρ1−ρ1−ρNs⋅δT,as the estimated observation time n−αT.Next, we note that(17)θˇknn−θ^kn−1n−1,accounts for the angle change between times n−1T and n−αT. The change between times n−1T and nT can be estimated as(18)11−αθˇknn−θ^kn−1n−1,assuming that angular velocity does not change over n−1T,nT. Using this approximation, we obtain the DCDSPCMNE θ˜knn by the following equation:(19)θ˜knn=θ^kn−1n−1+11−αθˇknn−θ^kn−1n−1,(20)=θ^kn−1n−1+1−ρ1−ρNsθˇknn−θ^kn−1n−11−ρ1−ρNs−Ns−1ρNs+1−NsρNs+ρδ.Before proceeding, we note that (20) reduces to(21)θ˜knn=θ^kn−1n−1+NS1−ρ1−ρNsNs−Ns+1ρ+ρNs+1θˇknn−θ^kn−1n−1,in the special case of δ=1/Ns.


## 6. Selection of Observer Gain

In Section5, we discussed how the updated estimated state vector is obtained in the normal situation, namely, in the resolvable situation. In particular, the estimated state vector is updated, in the proposed scheme, using the Luenberger observer as described in (10). One important issue that was not discussed in Section 5 was the selection of the observer gain gk. If certain conditions are met, the estimation errors approach zero as time progresses. For example, the state error vectors converge to zero if and only if all the eigenvalues of the matrix F−gk100 have magnitudes strictly less than one in the case of the algorithm in [22]. However, not so much information is available about the effect of observer gain selection beyond such theoretical basics. Hence, observer gains are generally selected through simulations. However, the question arises whether an observer gain selected to be optimal through simulations in certain situations will lead to reasonably good performance in other situations. In this section, we study the effect of observer gain selection through theoretical and numerical analysis and describe how we can obtain reasonably optimal observer gains.


### 6.1. State Estimation Error

We first define the state estimation errorekn by(22)ekn=y^knn−yknT,where ykt denotes the state vector defined as θkt,θ˙kt,θ¨ktT. Also, we define a vector δkn, which we call state evolution mismatch at time t=nT, by(23)δkn=yknT−F⋅ykn−1T.As discussed in [22], this vector δkn equals to a zero vector if the angular acceleration does not change between times n−1T and nT. However, the angular acceleration may change although the amount may be small, and δkn becomes small but nonzero vector. We define the angle measurement error εkn by(24)εkn=θ˜knn−θknT.In the next proposition, we obtain recurrence relation betweenekn and ekn−1.Proposition 1.
The state estimation errors satisfy the following recurrent relationship:(25)ekn=F−gkcTF⋅ekn−1+gkεkn+cTδkngk−δkn.Proof.
We first note, from (10), (4), and (24), that(26)y^knn=Fy^n−1n−1+θknT+εkn−θ^knn−1gk.
Next, from (22), (23), and (26), it follows that(27)ekn=Fy^n−1n−1+θknT+εkn−θ^knn−1gk−F⋅ykn−1T+δkn=Fy^n−1n−1−ykn−1T+θknT+εkn−θ^knn−1gk−δkn.
Now, sinceθkt=cTykt and θ^knn−1=cTy^knn−1, we next obtain(28)ekn=F⋅ekn−1+cTyknT−y^knn−1gk+εkngk−δkn.(29)=F⋅ekn−1+cTF⋅ykn−1T+δkn−F⋅ykn−1n−1gk+εkngk−δkn.
Finally, we obtain (25) by rearranging (29). □


### 6.2. Root Mean Squared Error and Observer Gain

In this section, we attempt to study theoretically how the observer gain affects the performance of angle estimation. In particular, we study how the selection of observer gain is related to the root mean squared error (RMSE) of the DoA of thekth object, which is defined by(30)RMSEk=1NT∑n=1NTEθ^knn−θknT2,where NT is the total number of observation intervals. Before proceeding, we note that RMSEk can be expressed, using ekn, as(31)RMSEk=1NT∑n=1NTEcTekn2,since(32)θ^knn−θknT=cTy^knn−cTynT=cTekn.First, we diagonalize the matrixF−gkcTF into VΛV−1, where Λ is a diagonal matrix whose diagonal elements are eigenvalues of F−gkcTF and V is a matrix whose column vectors are the corresponding eigenvectors. Here, we note that it is possible to represent gk in terms of Λ and V as described in the following lemma.Lemma 1.
(33)gk=VI−ΛV−1c.Proof.
SinceF−gkcTF=VΛV−1, it follows that(34)gkcTFc=Fc−VΛV−1c.
Here, we note thatcTFc is simply the 1,1 component of the matrix F, which is 1. Consequently, we obtain(35)gk=Fc−VΛV−1c=c−VΛV−1c,from which (33) readily follows. □Moreover, the observer gain vector is determined solely by the eigenvalues, namely, the diagonal elements ofΛ, as described in the following lemma.Lemma 2.
Letκ1, κ2ejϕ, and κ2e−jϕ denote the diagonal elements of Λ, where κ1 and κ2 are nonnegative real numbers and ϕ is a real number. (Since the characteristic equation of F−gkcTF is a third-order equation with real coefficients, at least one of the three eigenvalues must be real-valued and the other two must form a complex conjugate pair). Then,(36)gk1=1−κ1κ22,gk2=12T3κ1κ22−κ22−2κ1κ2cosϕ−κ1−2κ2cosϕ+3,gk3=1T2−κ1κ22+κ22+2κ1κ2cosϕ−κ1−2κ2cosϕ+1,where gki denotes the ith element of gk.Proof.
We note that the characteristic equation ofF−gkcTF is given by(37)λ3+gk1+gk2T+0.5gk3T2−3λ2+−2gk1−gk2T+0.5gk3T2+3λ+gk1−1=0.
Consequently, forκ1, κ2ejϕ, and κ2e−jϕ to be solutions to (37), the identity equation(38)λ3+gk1+gk2T+0.5gk3T2−3λ2+−2gk1−gk2T+0.5gk3T2+3λ+gk1−1=λ3−κ1+2κ2cosϕκ2+κ2κ2+2κ1cosϕλ−κ1κ22=0,must hold. Now, comparing the coefficients of the left and right-hand sides of (38) and solving for gki, we obtain the results. □Although the result in Lemma1 is simple, it is important in that we can construct the observer gain vector gk by suitably choosing the eigenvectors and eigenvalues that directly affect the system performance. Before proceeding with the performance analysis, let us rewrite (25) as(39)ekn=VΛV−1ekn−1+VI−ΛV−1cεkn+VI−ΛV−1ccT−Iδkn,by using (33). By repeatedly using (39), we can represent the state estimation error ekn in terms of its initial value ek0, the angle measurement error εkm, and the state evolution mismatch δkm as described in the following:(40)ekn=VΛnV−1ek0+∑m=1nVΛn−mI−ΛV−1cεkm+∑m=1nVΛn−mI−ΛV−1ccT−V−1δkm.We note that the second and third terms of (38) involve summation from m=1 to m=n. Consequently, the second and third terms may grow larger as n grows. However, if we choose the magnitudes of all the diagonal elements of Λ to be sufficiently less than one, only the error εkm and mismatch δkm with m close to n affect the values of ekn since Λk converges to a zero matrix as k grows. Moreover, we note that the magnitude of δkm is order of magnitudes smaller than that of εkn unless the observation time interval T is chosen to be very large, in which the tracking scheme in this paper does not work properly. Consequently, if we choose the magnitudes of all the diagonal elements of Λ to be small enough compared with 1, the third term of (40) is negligibly small in usual situations. Consequently, in the following, we shall assume that ekn can be approximated as(41)ekn=VΛnV−1ek0+∑m=1nVΛn−mI−ΛV−1cεkm,and attempt to compute RMSEk using (31). To obtain reasonably simple and useful result for RMSEk, we shall assume that the angle measurement error εkm has zero mean and that εkm and εkn are uncorrelated if n≠m, which are reasonably accurate assumptions in real situations. In the following theorem, we summarize the result.Theorem 1.
Assume that(42)Eεkn=0,and that(43)Eεkmεkn=δmnσk2n.Then,(44)RMSEk=1NT∑n=1NTcTVΛnV−1TVΛnV−1Tc+∑m=1ncTVΛn−mI−ΛV−1c2σk2m,where(45)T=σk20σk20Tσk20σk20Tσk20+σk2−1T2σk20+2σk2−1T3σk20T2σk20+2σk2−1T3σk20+4σk2−1+σk2−2T4.Proof.
See Appendix B. □Corollary 1.
If we further assume, in Theorem1, that there exists some nonnegative real number σ˜k such that(46)σkn=σ˜k,for all n, then(47)RMSEk=σ˜k1NT∑n=1NTcTVΛnV−1TVΛnV−1Tc+∑m=1ncTVΛn−mI−ΛV−1c2,where(48)T=11T1T21T2T23T31T23T36T4.According to our test with various simulation environments, assumptions (42) and (43) hold very accurately. However, assumption (46) does not hold, and the value of σk2n usually fluctuates between 1/2 and 2 times of the average value. Despite the invalidity of the assumption, the result in Corollary 1 is very useful. First, we note the simplicity of (47) compared with (44). Consequently, it is much easier, with (47), to grasp intuition about the behavior of RMSEk in relation to the choice of Λ. In particular, we note that there is no further information required about the variation of σk2n over n, and the RMSEk is represented in the unit of σ˜k in (47). Secondly, we can regard (47) with σ˜k=1/NT∑nNTσk2n as a simple and reasonable approximation of (44) since there is no reason to believe that σk2n will fluctuate extremely over n.Now, we shall show how we can select eigenvalues for optimal or near-optimal performance with (47), which is not environment dependent. Later, in this section, we shall also discuss how much the results with (47) are different from those with (44) using a number of simulation environment setups. We first choose eigenvalues κ1,κ2ejϕ, and κ2e−jϕ and compute the right-hand side of (47). Apparently, we need to obtain V for this purpose. However, we can compute F−gkcTF, which equals to VΛV−1, using Lemma 2. We further note that the right-hand side of (44) can be evaluated once we obtain VΛV−1.We recall thatekn converges to zero as n grows if and only if both κ1 and κ2 are smaller than 1. We consider 0.05,0.10,⋯,0.95 as candidate values for κ1 and κ2 and 0°,5°,⋯,175° as the values for ϕ. Consequently, we consider a total of 19×19×36 eigenvalue set combinations. Out of all considered values, the case of κ1=0.7, κ2=0.5, and ϕ=35° provides the lowest RMSEk value, which is 0.8646σ˜k. In the following, we shall express RMSEk in the unit of σ˜k for simplicity. In Figure 1, we illustrate the RMSEk values as functions of κ1 and κ2 for the case of ϕ=35°. In this figure, we note that the RMSEk value rises as κ1,κ2 moves away from 0.7,0.5. In particular, we note that the RMSEk value rises relatively mildly if κ1 or κ2 is decreased. However, the RMSEk value rises significantly as κ1 and κ2 approach 1. We found that the RMSEk values generally rise slowly as the ϕ value moves away from 35°. Although it is not possible to include all results with different ϕ values, we include the results with ϕ=15° and ϕ=55° in Figures 2 and 3 to illustrate how the RMSEk value changes with the ϕ value.Figure 1
RMSEk for ϕ=35° as a function of κ1 and κ2.Figure 2
RMSEk for ϕ=15° as a function of κ1 and κ2.Figure 3
RMSEk for ϕ=55° as a function of κ1 and κ2.So far, we discussed how we can choose the optimal observer gain vector with (45) assuming that σkn is a constant independent of n. Hence, the validity of the proposed method of obtaining the observer gain will depend on how the effect of variations in σkn affects the value of RMSEk. According to our extensive test, the effect of variation of σkn is very negligible except for extreme situations that do not affect the optimal selection of the observer gain. To illustrate in a space economical way, we consider a situation in which there are five objects with DoA trajectories illustrated in Figure 4. In realistic environments, the DoA trajectories may not change as fast as in Figure 4. However, if we choose a slowly moving environment, the σkn value would not change very much, and hence the effect would also be negligible. For this reason, we have chosen a very fast changing environment shown in Figure 4. To obtain σkn, we assume that only object k moves following the trajectory θknT and perform simulations with the proposed scheme (DSPCLT). We recall that we divided the situations into resolvable and nonresolvable situations to avoid overlapping situations. For this reason, we assume that only object k exists in computing σkn.Figure 4
DoA trajectories used for simulations.For numerical evaluations, we assume that the numberN of sensor elements is 12, the total number NT of observation intervals is 60, and the signal-to-noise ratio is 5 dB. With 106 simulation run trials, we computed σk2n, which is illustrated in Figure 5. In this figure, we note that the value of σk2n is generally related to the angular acceleration. After obtaining the estimate of σk2n, we computed RMSEk using (44) in the unit of σ˜k=1/NT∑n=1NTσk2n. To compare the results with the values computed using (47), we evaluated the ratio obtained by dividing the RMSEk value computed using (44) by that using (47) for each set of κ1,κ2, and ϕ values. Some results of the evaluations are summarized in Table 1. In this table, the number pair a,b means the minimum and maximum values of the ratios over all κ1 and κ2 values for given ϕ value and sequence σk2n. For example, for ϕ=15° and σ12n, the ratio falls between 0.965 and 1.141. With this set of numbers, we can conclude how different the two RMSEk values computed by (42) and (45) are to the extreme. However, the extreme values do not tell the whole story since the extreme values, namely, the maximum and minimum ratio values, occur at κ1 and κ2 values close to 1 so that the RMSEk value is very large. In other words, the extreme values happen with eigenvalue selections for which ekn do not converge to zero very fast.Figure 5
σk2n in the unit of σ˜k.Table 1
(Minimum, maximum) pair of theRMSEk value ratio.
σ12nσ22nσ32nσ42nσ52nϕ=15°0.965,1.1410.973,1.0220.842,0.9840.908,0.9800.799,0.987ϕ=35°0.965,1.1240.973,1.0170.859,0.9840.916,0.9800.822,0.987ϕ=55°0.967,1.0500.973,0.9950.922,0.9840.949,0.9800.908,0.987To understand what happens, let us first consider the case ofϕ=15° and σ52n in which the minimum value is particularly small. In Figure 6, we plotted the ratios as a function of κ1 and κ2 for this case, namely, for ϕ=15° and σ52n. We observe that the minimum value 0.799 happens near the region κ1≈1 and κ2≈1 for which the RMSEk values are very large. In contrast, the ratio is very close to 1.0 except for the region in which RMSEk remains to be acceptably small. Moreover, we recall that the results in Figure 6 correspond to some extreme case in which the ratio fluctuates a lot. As a more typical case, we depict, in Figure 7, the ratio for ϕ=35° and σ12n. We note that the ratio is again very close to 1.0 except for the region in which κ1≈1 and κ2≈1. Consequently, we observe that the RMSEk value computed with (47) is very close to that with (44) except at κ1, κ2, and ϕ values for which ekn does not converge to zero fast enough, and RMSEk is very large.Figure 6
RMSE5 computed with time-varying σ5n for ϕ=15°.Figure 7
RMSE1 computed with time-varying σ1n for ϕ=35°.


## 6.1. State Estimation Error

We first define the state estimation errorekn by(22)ekn=y^knn−yknT,where ykt denotes the state vector defined as θkt,θ˙kt,θ¨ktT. Also, we define a vector δkn, which we call state evolution mismatch at time t=nT, by(23)δkn=yknT−F⋅ykn−1T.As discussed in [22], this vector δkn equals to a zero vector if the angular acceleration does not change between times n−1T and nT. However, the angular acceleration may change although the amount may be small, and δkn becomes small but nonzero vector. We define the angle measurement error εkn by(24)εkn=θ˜knn−θknT.In the next proposition, we obtain recurrence relation betweenekn and ekn−1.Proposition 1.
The state estimation errors satisfy the following recurrent relationship:(25)ekn=F−gkcTF⋅ekn−1+gkεkn+cTδkngk−δkn.Proof.
We first note, from (10), (4), and (24), that(26)y^knn=Fy^n−1n−1+θknT+εkn−θ^knn−1gk.
Next, from (22), (23), and (26), it follows that(27)ekn=Fy^n−1n−1+θknT+εkn−θ^knn−1gk−F⋅ykn−1T+δkn=Fy^n−1n−1−ykn−1T+θknT+εkn−θ^knn−1gk−δkn.
Now, sinceθkt=cTykt and θ^knn−1=cTy^knn−1, we next obtain(28)ekn=F⋅ekn−1+cTyknT−y^knn−1gk+εkngk−δkn.(29)=F⋅ekn−1+cTF⋅ykn−1T+δkn−F⋅ykn−1n−1gk+εkngk−δkn.
Finally, we obtain (25) by rearranging (29). □


## 6.2. Root Mean Squared Error and Observer Gain

In this section, we attempt to study theoretically how the observer gain affects the performance of angle estimation. In particular, we study how the selection of observer gain is related to the root mean squared error (RMSE) of the DoA of thekth object, which is defined by(30)RMSEk=1NT∑n=1NTEθ^knn−θknT2,where NT is the total number of observation intervals. Before proceeding, we note that RMSEk can be expressed, using ekn, as(31)RMSEk=1NT∑n=1NTEcTekn2,since(32)θ^knn−θknT=cTy^knn−cTynT=cTekn.First, we diagonalize the matrixF−gkcTF into VΛV−1, where Λ is a diagonal matrix whose diagonal elements are eigenvalues of F−gkcTF and V is a matrix whose column vectors are the corresponding eigenvectors. Here, we note that it is possible to represent gk in terms of Λ and V as described in the following lemma.Lemma 1.
(33)gk=VI−ΛV−1c.Proof.
SinceF−gkcTF=VΛV−1, it follows that(34)gkcTFc=Fc−VΛV−1c.
Here, we note thatcTFc is simply the 1,1 component of the matrix F, which is 1. Consequently, we obtain(35)gk=Fc−VΛV−1c=c−VΛV−1c,from which (33) readily follows. □Moreover, the observer gain vector is determined solely by the eigenvalues, namely, the diagonal elements ofΛ, as described in the following lemma.Lemma 2.
Letκ1, κ2ejϕ, and κ2e−jϕ denote the diagonal elements of Λ, where κ1 and κ2 are nonnegative real numbers and ϕ is a real number. (Since the characteristic equation of F−gkcTF is a third-order equation with real coefficients, at least one of the three eigenvalues must be real-valued and the other two must form a complex conjugate pair). Then,(36)gk1=1−κ1κ22,gk2=12T3κ1κ22−κ22−2κ1κ2cosϕ−κ1−2κ2cosϕ+3,gk3=1T2−κ1κ22+κ22+2κ1κ2cosϕ−κ1−2κ2cosϕ+1,where gki denotes the ith element of gk.Proof.
We note that the characteristic equation ofF−gkcTF is given by(37)λ3+gk1+gk2T+0.5gk3T2−3λ2+−2gk1−gk2T+0.5gk3T2+3λ+gk1−1=0.
Consequently, forκ1, κ2ejϕ, and κ2e−jϕ to be solutions to (37), the identity equation(38)λ3+gk1+gk2T+0.5gk3T2−3λ2+−2gk1−gk2T+0.5gk3T2+3λ+gk1−1=λ3−κ1+2κ2cosϕκ2+κ2κ2+2κ1cosϕλ−κ1κ22=0,must hold. Now, comparing the coefficients of the left and right-hand sides of (38) and solving for gki, we obtain the results. □Although the result in Lemma1 is simple, it is important in that we can construct the observer gain vector gk by suitably choosing the eigenvectors and eigenvalues that directly affect the system performance. Before proceeding with the performance analysis, let us rewrite (25) as(39)ekn=VΛV−1ekn−1+VI−ΛV−1cεkn+VI−ΛV−1ccT−Iδkn,by using (33). By repeatedly using (39), we can represent the state estimation error ekn in terms of its initial value ek0, the angle measurement error εkm, and the state evolution mismatch δkm as described in the following:(40)ekn=VΛnV−1ek0+∑m=1nVΛn−mI−ΛV−1cεkm+∑m=1nVΛn−mI−ΛV−1ccT−V−1δkm.We note that the second and third terms of (38) involve summation from m=1 to m=n. Consequently, the second and third terms may grow larger as n grows. However, if we choose the magnitudes of all the diagonal elements of Λ to be sufficiently less than one, only the error εkm and mismatch δkm with m close to n affect the values of ekn since Λk converges to a zero matrix as k grows. Moreover, we note that the magnitude of δkm is order of magnitudes smaller than that of εkn unless the observation time interval T is chosen to be very large, in which the tracking scheme in this paper does not work properly. Consequently, if we choose the magnitudes of all the diagonal elements of Λ to be small enough compared with 1, the third term of (40) is negligibly small in usual situations. Consequently, in the following, we shall assume that ekn can be approximated as(41)ekn=VΛnV−1ek0+∑m=1nVΛn−mI−ΛV−1cεkm,and attempt to compute RMSEk using (31). To obtain reasonably simple and useful result for RMSEk, we shall assume that the angle measurement error εkm has zero mean and that εkm and εkn are uncorrelated if n≠m, which are reasonably accurate assumptions in real situations. In the following theorem, we summarize the result.Theorem 1.
Assume that(42)Eεkn=0,and that(43)Eεkmεkn=δmnσk2n.Then,(44)RMSEk=1NT∑n=1NTcTVΛnV−1TVΛnV−1Tc+∑m=1ncTVΛn−mI−ΛV−1c2σk2m,where(45)T=σk20σk20Tσk20σk20Tσk20+σk2−1T2σk20+2σk2−1T3σk20T2σk20+2σk2−1T3σk20+4σk2−1+σk2−2T4.Proof.
See Appendix B. □Corollary 1.
If we further assume, in Theorem1, that there exists some nonnegative real number σ˜k such that(46)σkn=σ˜k,for all n, then(47)RMSEk=σ˜k1NT∑n=1NTcTVΛnV−1TVΛnV−1Tc+∑m=1ncTVΛn−mI−ΛV−1c2,where(48)T=11T1T21T2T23T31T23T36T4.According to our test with various simulation environments, assumptions (42) and (43) hold very accurately. However, assumption (46) does not hold, and the value of σk2n usually fluctuates between 1/2 and 2 times of the average value. Despite the invalidity of the assumption, the result in Corollary 1 is very useful. First, we note the simplicity of (47) compared with (44). Consequently, it is much easier, with (47), to grasp intuition about the behavior of RMSEk in relation to the choice of Λ. In particular, we note that there is no further information required about the variation of σk2n over n, and the RMSEk is represented in the unit of σ˜k in (47). Secondly, we can regard (47) with σ˜k=1/NT∑nNTσk2n as a simple and reasonable approximation of (44) since there is no reason to believe that σk2n will fluctuate extremely over n.Now, we shall show how we can select eigenvalues for optimal or near-optimal performance with (47), which is not environment dependent. Later, in this section, we shall also discuss how much the results with (47) are different from those with (44) using a number of simulation environment setups. We first choose eigenvalues κ1,κ2ejϕ, and κ2e−jϕ and compute the right-hand side of (47). Apparently, we need to obtain V for this purpose. However, we can compute F−gkcTF, which equals to VΛV−1, using Lemma 2. We further note that the right-hand side of (44) can be evaluated once we obtain VΛV−1.We recall thatekn converges to zero as n grows if and only if both κ1 and κ2 are smaller than 1. We consider 0.05,0.10,⋯,0.95 as candidate values for κ1 and κ2 and 0°,5°,⋯,175° as the values for ϕ. Consequently, we consider a total of 19×19×36 eigenvalue set combinations. Out of all considered values, the case of κ1=0.7, κ2=0.5, and ϕ=35° provides the lowest RMSEk value, which is 0.8646σ˜k. In the following, we shall express RMSEk in the unit of σ˜k for simplicity. In Figure 1, we illustrate the RMSEk values as functions of κ1 and κ2 for the case of ϕ=35°. In this figure, we note that the RMSEk value rises as κ1,κ2 moves away from 0.7,0.5. In particular, we note that the RMSEk value rises relatively mildly if κ1 or κ2 is decreased. However, the RMSEk value rises significantly as κ1 and κ2 approach 1. We found that the RMSEk values generally rise slowly as the ϕ value moves away from 35°. Although it is not possible to include all results with different ϕ values, we include the results with ϕ=15° and ϕ=55° in Figures 2 and 3 to illustrate how the RMSEk value changes with the ϕ value.Figure 1
RMSEk for ϕ=35° as a function of κ1 and κ2.Figure 2
RMSEk for ϕ=15° as a function of κ1 and κ2.Figure 3
RMSEk for ϕ=55° as a function of κ1 and κ2.So far, we discussed how we can choose the optimal observer gain vector with (45) assuming that σkn is a constant independent of n. Hence, the validity of the proposed method of obtaining the observer gain will depend on how the effect of variations in σkn affects the value of RMSEk. According to our extensive test, the effect of variation of σkn is very negligible except for extreme situations that do not affect the optimal selection of the observer gain. To illustrate in a space economical way, we consider a situation in which there are five objects with DoA trajectories illustrated in Figure 4. In realistic environments, the DoA trajectories may not change as fast as in Figure 4. However, if we choose a slowly moving environment, the σkn value would not change very much, and hence the effect would also be negligible. For this reason, we have chosen a very fast changing environment shown in Figure 4. To obtain σkn, we assume that only object k moves following the trajectory θknT and perform simulations with the proposed scheme (DSPCLT). We recall that we divided the situations into resolvable and nonresolvable situations to avoid overlapping situations. For this reason, we assume that only object k exists in computing σkn.Figure 4
DoA trajectories used for simulations.For numerical evaluations, we assume that the numberN of sensor elements is 12, the total number NT of observation intervals is 60, and the signal-to-noise ratio is 5 dB. With 106 simulation run trials, we computed σk2n, which is illustrated in Figure 5. In this figure, we note that the value of σk2n is generally related to the angular acceleration. After obtaining the estimate of σk2n, we computed RMSEk using (44) in the unit of σ˜k=1/NT∑n=1NTσk2n. To compare the results with the values computed using (47), we evaluated the ratio obtained by dividing the RMSEk value computed using (44) by that using (47) for each set of κ1,κ2, and ϕ values. Some results of the evaluations are summarized in Table 1. In this table, the number pair a,b means the minimum and maximum values of the ratios over all κ1 and κ2 values for given ϕ value and sequence σk2n. For example, for ϕ=15° and σ12n, the ratio falls between 0.965 and 1.141. With this set of numbers, we can conclude how different the two RMSEk values computed by (42) and (45) are to the extreme. However, the extreme values do not tell the whole story since the extreme values, namely, the maximum and minimum ratio values, occur at κ1 and κ2 values close to 1 so that the RMSEk value is very large. In other words, the extreme values happen with eigenvalue selections for which ekn do not converge to zero very fast.Figure 5
σk2n in the unit of σ˜k.Table 1
(Minimum, maximum) pair of theRMSEk value ratio.
σ12nσ22nσ32nσ42nσ52nϕ=15°0.965,1.1410.973,1.0220.842,0.9840.908,0.9800.799,0.987ϕ=35°0.965,1.1240.973,1.0170.859,0.9840.916,0.9800.822,0.987ϕ=55°0.967,1.0500.973,0.9950.922,0.9840.949,0.9800.908,0.987To understand what happens, let us first consider the case ofϕ=15° and σ52n in which the minimum value is particularly small. In Figure 6, we plotted the ratios as a function of κ1 and κ2 for this case, namely, for ϕ=15° and σ52n. We observe that the minimum value 0.799 happens near the region κ1≈1 and κ2≈1 for which the RMSEk values are very large. In contrast, the ratio is very close to 1.0 except for the region in which RMSEk remains to be acceptably small. Moreover, we recall that the results in Figure 6 correspond to some extreme case in which the ratio fluctuates a lot. As a more typical case, we depict, in Figure 7, the ratio for ϕ=35° and σ12n. We note that the ratio is again very close to 1.0 except for the region in which κ1≈1 and κ2≈1. Consequently, we observe that the RMSEk value computed with (47) is very close to that with (44) except at κ1, κ2, and ϕ values for which ekn does not converge to zero fast enough, and RMSEk is very large.Figure 6
RMSE5 computed with time-varying σ5n for ϕ=15°.Figure 7
RMSE1 computed with time-varying σ1n for ϕ=35°.


## 7. System Performance and Computational Complexity

In the section, we compare the performance and the computational complexity of the proposed scheme with those of two existing algorithms, namely, AMEND and PASTd Kalman. Most of all, we observe that the proposed scheme (DSPCLT) exhibits excellent performance with significantly low computational complexity in comparison with AMEND and PASTd Kalman.


### 7.1. Performance Evaluation of Algorithms

First, let us compare the performance of DSPCLT with that of AMEND and PASTd Kalman. As measures of performance, we use the following quantities:(49)RMSEkn=Eθ^knT−θknT2,(50)RMSET=1KNT∑k=1K∑n=1NTEθ^knT−θknT2,where θ^knT is the estimated DoA of θknT and NT is the the total number of observation intervals. We note that RMSEkn is the RMSE of the DoA of the kth object at time t=nT and that RMSET represents the total RMSE over all objects and observation intervals. We recall that we regard θ^knn as θ^knT in the proposed scheme.However, we need to be very careful in assessing the performance with RMSE. This is because there are some possibilities of tracking wrong objects after passing through a crossover point or of object tracking divergence due to lack of sensitivity to object movements, which can result in huge estimation errors. The value of RMSE can be meaningless when such catastrophic tracking failure happens. To make the value of RMSE more meaningful, we shall regard the situation as tracking failure if and only if the estimated DoA of any of the signals is more than5° away from the actual DoA and compute the RMSE value disregarding such cases of tracking failure. (In this paper, we shall measure DoAs and RMSE in degrees rather than in radians.) Since the RMSE values do not account for the cases of tracking failure, we shall also consider, in addition to RMSE, the probability PF of tracking failure as a performance measure.For performance evaluations, we first consider a situation in whichK=3 signals impinge upon a ULA with N=12 sensors separated by half-wavelength from time-varying directions θ1nT, θ2nT, and θ3nT, n=1,2,…, NT=60, illustrated in Figure 8. We assume that the signal-to-noise ratio of each signal at each sensor is given as 5 dB. During each time interval nT,n+1T, total Ns=200 data snapshots are uniformly captured at each sensor with δ=1/Ns.Figure 8
DoA trajectories used for simulations.We note that the choice of forgetting factorρ affects the tracking performance. Hence, we compare the performance with various ρ values from 0.91 to 0.99. In the case of PASTd Kalman, we also consider, instead of ρNs−n, the forgetting coefficient n/n+1, which was used in [21]. For AMEND, we use the observer gain vector obtained with the eigenvalues 0.75,0.75ejπ/18, and 0.75e−jπ/18 as suggested in [22]. However, we select the observer gain vector obtained, using Lemma 2, with eigenvalues 0.7,0.5ej35π/180 and 0.5e−j35π/180 for DSPCLT. We recall that observer gain is not used in PASTd Kalman.As the value ofM is used to compute r˜m in (A.3) and (A.4), we use M=K since the RMSE performance enhances rapidly as M increases from 1 to K and then becomes virtually invariant as it further grows. Also, we separate the resolvable and nonresolvable situations with threshold values 5°,4.5°, and 4° for N=12,16, and 20, respectively.Simulations are performed with105 trials. The failure probability PF and RMSET obtained through the simulations are summarized in Table 2. We may say that AMEND exhibits the most reliable performance in that it does not experience failure for any choice of ρ value. However, AMEND is worst in terms of RMSET value. In particular, we note that both DSPCLT and PASTd Kalman exhibit not only better RMSET values than AMEND but also achieve zero failures. The best RMSET value for each scheme is marked in bold-faced fonts, and we note that DSPCLT achieves the lowest RMSET value among the three schemes. By the way, the RMSET value 0.6160 listed on the lowest row for PASTd Kalman corresponds to the forgetting coefficient n/n+1 suggested in [21].Table 2
PF and RMSET (K=3,N=12,M=3,Ns=200).
ρDSPCLTAMENDPASTd KalmanPF%RMSET°PF%RMSET°PF%RMSET°0.910.400.435400.47190.830.34500.920.240.415500.46340.570.32760.930.120.374200.45610.300.30840.940.090.354100.45100.100.29280.950.050.318100.45030.030.27920.960.030.301400.46020.020.28310.9700.271500.490500.32470.9800.292300.571300.43080.9900.385100.749600.6328—————00.6160To illustrate the performance variation over time,RMSEkn is plotted in Figure 9. In this figure, we used ρ=0.95 for AMEND and 0.97 for DSPCLT and PASTd Kalman. In Table 2, we note that AMEND exhibits the worst RMSET values. Also, we note DSPCLT achieves slightly better RMSET value compared with PASTd Kalman. In Figure 9, we observe that DSPCLT exhibits best overall RMSEkn performance except for the region where the DoAs of objects overlap. In particular, we note that the performance of DSPCLT is significantly better than others for the first object during n∈1,2,⋯,10. From Figure 8, we note that the DoA of the first object changes very fast over n∈1,2,⋯,10. Consequently, we can conclude that the delay compensation algorithm introduced in Section 5.2 indeed provides improved performance in fast changing environments. We also observe that DSPCLT is less susceptible to object crossover compared with AMEND. This is because the nonresolvable situation is separately treated, in the proposed scheme, using the algorithm described in Section 4.Figure 9
Performance variation over time.N=12,K=3.One of the main reasons for the overall performance superiority of the proposed scheme is the comparatively larger effective aperture size with the same number of sensor elements. For illustration of the results of effective aperture size, we compare, in Figure10, realizations of the function fθ, defined in (11), at time index n=19. In this figure, we observe that fθ is much sharper for DSPCLT than for AMEND and PASTd Kalman. For this reason, DSPCLT provides better overall RMSE performance than AMEND rivaling PASTd Kalman that employs much more sophisticated Kalman filter algorithm.Figure 10
Comparison offθ functions at time index n=19.We next study the performance variation with increased number of objects. For illustration, we consider the situation in which there areK=5 objects with trajectories illustrated in Figure 4. Here, we note that AMEND can track only K objects satisfying K<N/2, while DSPCLT and PASTd Kalman can track upto N−1 objects. Consequently, at least N=11 sensor elements are needed for AMEND to be able to track K=5 objects. To see the effect of the number of sensor elements, we consider three values 12,16, and 20 as the number N of antenna elements of the ULA. Again, the received SNR of each of the 5 objects is set to be 5 dB, the number Ns of snapshots is chosen to be 200 with δ=1/Ns, and the value of M is chosen to be the same as K.The simulations are performed with105 runs again, and the results are summarized in Table 3 and Figure 11. As in the case of K=3, we obtained simulation results with various forgetting factor values and included the best results in the table and the figure. In particular, the value ρ indicates the forgetting factor value that leads to the best RMSET value. In Table 3 and Figure 11, we observe that AMEND exhibits the worst overall performance in almost every aspect. We observe that both AMEND and DSPCLT experienced nonzero failure probabilities, which decrease as the number of sensor elements increases. Consequently, we can conclude that PASTd Kalman achieves the most reliable performance among the three considered schemes. However, we note that DSPCLT exhibits performance better than PASTd Kalman with a sufficiently large number of sensor elements. We also note that DSPCLT, compared with AMEND, not only achieves lower failure probability but also exhibits lower RMSET. In particular, we observe that the RMSET value of DSPCLT is significantly lower than that of AMEND despite the fact that RMSET is computed excluding the cases in which tracking fails. Next, in Figure 11, the sequence of RMSEkn values for each of the five objects is depicted for the case of N=20. As in the previous case of K=3 objects, AMEND exhibits the worst performance, and DSPCLT achieves similar or better performance, compared with PASTd Kalman, except around crossover points.Table 3
PF and RMSET (K=5,M=5,Ns=200,N=12,16,20).
NDSPCLTAMENDPASTd KalmanρPF%RMSET°ρPF%RMSET°ρPF%RMSET°120.974.410.50800.9722.161.15890.9700.3787160.970.040.27910.950.310.43600.9700.3135200.9500.18840.9300.29370.9600.2433Figure 11
Performance variation over time.N=20,K=5.


### 7.2. Computational Complexity

In this section, we compare the computational complexities of DSPCLT, AMEND, and PASTd Kalman. Before proceeding, we note that the complexity of DSPCLT varies depending on whether the object is in the resolvable or nonresolvable situation. Since the complexity of DSPCLT in the nonresolvable situation is significantly smaller than that in the nonresolvable situation and nonresolvable situation happens only very rarely when the number of sensor elements is large enough, we shall assume that only resolvable situation happens when computing the complexity of DSPCLT. Before proceeding, we note that all the three tracking schemes DSPCLT, AMEND, and PASTd Kalman can be divided into three processes, namely, “subspace construction,” “projector construction,” and “filtering.” Here, the subspace construction is in essence the construction of the vectors that generate the signal subspace.In PASTd Kalman [21] and AMEND [22], subspace construction is carried out using the recursive least square (RLS) method, while the vectors in DSPCLT are constructed directly from the estimated correlation matrix using (A.1). In the projector construction stage, an orthogonal projector matrix such as Πn is constructed. And state vectors are updated using Kalman filter or Luenberger observer in the final stage of filtering. As described in Section 1, Kalman filter is used in PASTd Kalman, and Luenberger observer is used in DSPCLT and AMEND. In Table 4, the numbers of real floating point operations needed to carry out each of the three processes for each of the three schemes as well as the total numbers are presented.Table 4
Required floating point operations.
DSPCLTAMENDPASTd KalmanSubspace construction5M2N−K+1Ns64K2N−2KNsK28N+7NsProjector construction8K2N−K216N−KN−2K28KN2Filtering162N−K216N−K28N2Ns+8K+2N2Total5M2N−K+1Ns+8K+162N−K264K2N−2KNs+16N−KN−2K2+N−K28KN+8N2+7KNs+161+KN2For more concrete picture of computational complexity, we evaluated the numbers in Table4 with the choice of parameters used in obtaining Figures 9 and 11, the results of which are presented, respectively, in Tables 5 and 6. We note that a significant portion of computational complexity is needed in subspace construction. We note that AMEND requires particularly large numbers of operations in the subspace construction compared with DSPCLT. This is because AMEND uses matrix inversion and QR decomposition even though it does not use ED nor SVD. We note that DSPCM does not require matrix inversion nor QR decomposition. We note that PASTd Kalman, compared with DSPCLT and AMEND, requires larger amount of operations in the filtering stage since it uses computationally expensive Kalman filtering. Finally, we note that DSPCLT requires 2∼6 times more operations to carry out projector construction and filtering compared with AMEND. This is because DSPCLT has larger effective aperture size despite that the same number of sensors is used. However, the overall computational complexity of DSPCLT is significantly smaller than that of AMEND and PASTd Kalman because the computational complexity of the second stage is much smaller than that of either the first or the third stage. Overall, DSPCLT requires only about 11.99% (18.78%) and 7.70% (19.96%) real floating point operations compared with AMEND and PASTd Kalman, respectively, for the case of N=12,M=K=3,Ns=200 and N=20,M=K=5,Ns=200.Table 5
Required floating point operations:N=12,M=K=3,Ns=200.
DSPCLTAMENDPASTd KalmanSubspace construction66,000691,200205,800Projector construction10,5845,1843,456Filtering7,0561,296236,160Total83,640697,680445,416Table 6
Required floating point operations:N=20,M=K=5,Ns=200.
DSPCLTAMENDPASTd KalmanSubspace construction180,0003,200,000567,000Projector construction49,00024,00016,000Filtering19,6003,600662,400Total248,6003,227,6001,245,400


## 7.1. Performance Evaluation of Algorithms

First, let us compare the performance of DSPCLT with that of AMEND and PASTd Kalman. As measures of performance, we use the following quantities:(49)RMSEkn=Eθ^knT−θknT2,(50)RMSET=1KNT∑k=1K∑n=1NTEθ^knT−θknT2,where θ^knT is the estimated DoA of θknT and NT is the the total number of observation intervals. We note that RMSEkn is the RMSE of the DoA of the kth object at time t=nT and that RMSET represents the total RMSE over all objects and observation intervals. We recall that we regard θ^knn as θ^knT in the proposed scheme.However, we need to be very careful in assessing the performance with RMSE. This is because there are some possibilities of tracking wrong objects after passing through a crossover point or of object tracking divergence due to lack of sensitivity to object movements, which can result in huge estimation errors. The value of RMSE can be meaningless when such catastrophic tracking failure happens. To make the value of RMSE more meaningful, we shall regard the situation as tracking failure if and only if the estimated DoA of any of the signals is more than5° away from the actual DoA and compute the RMSE value disregarding such cases of tracking failure. (In this paper, we shall measure DoAs and RMSE in degrees rather than in radians.) Since the RMSE values do not account for the cases of tracking failure, we shall also consider, in addition to RMSE, the probability PF of tracking failure as a performance measure.For performance evaluations, we first consider a situation in whichK=3 signals impinge upon a ULA with N=12 sensors separated by half-wavelength from time-varying directions θ1nT, θ2nT, and θ3nT, n=1,2,…, NT=60, illustrated in Figure 8. We assume that the signal-to-noise ratio of each signal at each sensor is given as 5 dB. During each time interval nT,n+1T, total Ns=200 data snapshots are uniformly captured at each sensor with δ=1/Ns.Figure 8
DoA trajectories used for simulations.We note that the choice of forgetting factorρ affects the tracking performance. Hence, we compare the performance with various ρ values from 0.91 to 0.99. In the case of PASTd Kalman, we also consider, instead of ρNs−n, the forgetting coefficient n/n+1, which was used in [21]. For AMEND, we use the observer gain vector obtained with the eigenvalues 0.75,0.75ejπ/18, and 0.75e−jπ/18 as suggested in [22]. However, we select the observer gain vector obtained, using Lemma 2, with eigenvalues 0.7,0.5ej35π/180 and 0.5e−j35π/180 for DSPCLT. We recall that observer gain is not used in PASTd Kalman.As the value ofM is used to compute r˜m in (A.3) and (A.4), we use M=K since the RMSE performance enhances rapidly as M increases from 1 to K and then becomes virtually invariant as it further grows. Also, we separate the resolvable and nonresolvable situations with threshold values 5°,4.5°, and 4° for N=12,16, and 20, respectively.Simulations are performed with105 trials. The failure probability PF and RMSET obtained through the simulations are summarized in Table 2. We may say that AMEND exhibits the most reliable performance in that it does not experience failure for any choice of ρ value. However, AMEND is worst in terms of RMSET value. In particular, we note that both DSPCLT and PASTd Kalman exhibit not only better RMSET values than AMEND but also achieve zero failures. The best RMSET value for each scheme is marked in bold-faced fonts, and we note that DSPCLT achieves the lowest RMSET value among the three schemes. By the way, the RMSET value 0.6160 listed on the lowest row for PASTd Kalman corresponds to the forgetting coefficient n/n+1 suggested in [21].Table 2
PF and RMSET (K=3,N=12,M=3,Ns=200).
ρDSPCLTAMENDPASTd KalmanPF%RMSET°PF%RMSET°PF%RMSET°0.910.400.435400.47190.830.34500.920.240.415500.46340.570.32760.930.120.374200.45610.300.30840.940.090.354100.45100.100.29280.950.050.318100.45030.030.27920.960.030.301400.46020.020.28310.9700.271500.490500.32470.9800.292300.571300.43080.9900.385100.749600.6328—————00.6160To illustrate the performance variation over time,RMSEkn is plotted in Figure 9. In this figure, we used ρ=0.95 for AMEND and 0.97 for DSPCLT and PASTd Kalman. In Table 2, we note that AMEND exhibits the worst RMSET values. Also, we note DSPCLT achieves slightly better RMSET value compared with PASTd Kalman. In Figure 9, we observe that DSPCLT exhibits best overall RMSEkn performance except for the region where the DoAs of objects overlap. In particular, we note that the performance of DSPCLT is significantly better than others for the first object during n∈1,2,⋯,10. From Figure 8, we note that the DoA of the first object changes very fast over n∈1,2,⋯,10. Consequently, we can conclude that the delay compensation algorithm introduced in Section 5.2 indeed provides improved performance in fast changing environments. We also observe that DSPCLT is less susceptible to object crossover compared with AMEND. This is because the nonresolvable situation is separately treated, in the proposed scheme, using the algorithm described in Section 4.Figure 9
Performance variation over time.N=12,K=3.One of the main reasons for the overall performance superiority of the proposed scheme is the comparatively larger effective aperture size with the same number of sensor elements. For illustration of the results of effective aperture size, we compare, in Figure10, realizations of the function fθ, defined in (11), at time index n=19. In this figure, we observe that fθ is much sharper for DSPCLT than for AMEND and PASTd Kalman. For this reason, DSPCLT provides better overall RMSE performance than AMEND rivaling PASTd Kalman that employs much more sophisticated Kalman filter algorithm.Figure 10
Comparison offθ functions at time index n=19.We next study the performance variation with increased number of objects. For illustration, we consider the situation in which there areK=5 objects with trajectories illustrated in Figure 4. Here, we note that AMEND can track only K objects satisfying K<N/2, while DSPCLT and PASTd Kalman can track upto N−1 objects. Consequently, at least N=11 sensor elements are needed for AMEND to be able to track K=5 objects. To see the effect of the number of sensor elements, we consider three values 12,16, and 20 as the number N of antenna elements of the ULA. Again, the received SNR of each of the 5 objects is set to be 5 dB, the number Ns of snapshots is chosen to be 200 with δ=1/Ns, and the value of M is chosen to be the same as K.The simulations are performed with105 runs again, and the results are summarized in Table 3 and Figure 11. As in the case of K=3, we obtained simulation results with various forgetting factor values and included the best results in the table and the figure. In particular, the value ρ indicates the forgetting factor value that leads to the best RMSET value. In Table 3 and Figure 11, we observe that AMEND exhibits the worst overall performance in almost every aspect. We observe that both AMEND and DSPCLT experienced nonzero failure probabilities, which decrease as the number of sensor elements increases. Consequently, we can conclude that PASTd Kalman achieves the most reliable performance among the three considered schemes. However, we note that DSPCLT exhibits performance better than PASTd Kalman with a sufficiently large number of sensor elements. We also note that DSPCLT, compared with AMEND, not only achieves lower failure probability but also exhibits lower RMSET. In particular, we observe that the RMSET value of DSPCLT is significantly lower than that of AMEND despite the fact that RMSET is computed excluding the cases in which tracking fails. Next, in Figure 11, the sequence of RMSEkn values for each of the five objects is depicted for the case of N=20. As in the previous case of K=3 objects, AMEND exhibits the worst performance, and DSPCLT achieves similar or better performance, compared with PASTd Kalman, except around crossover points.Table 3
PF and RMSET (K=5,M=5,Ns=200,N=12,16,20).
NDSPCLTAMENDPASTd KalmanρPF%RMSET°ρPF%RMSET°ρPF%RMSET°120.974.410.50800.9722.161.15890.9700.3787160.970.040.27910.950.310.43600.9700.3135200.9500.18840.9300.29370.9600.2433Figure 11
Performance variation over time.N=20,K=5.


## 7.2. Computational Complexity

In this section, we compare the computational complexities of DSPCLT, AMEND, and PASTd Kalman. Before proceeding, we note that the complexity of DSPCLT varies depending on whether the object is in the resolvable or nonresolvable situation. Since the complexity of DSPCLT in the nonresolvable situation is significantly smaller than that in the nonresolvable situation and nonresolvable situation happens only very rarely when the number of sensor elements is large enough, we shall assume that only resolvable situation happens when computing the complexity of DSPCLT. Before proceeding, we note that all the three tracking schemes DSPCLT, AMEND, and PASTd Kalman can be divided into three processes, namely, “subspace construction,” “projector construction,” and “filtering.” Here, the subspace construction is in essence the construction of the vectors that generate the signal subspace.In PASTd Kalman [21] and AMEND [22], subspace construction is carried out using the recursive least square (RLS) method, while the vectors in DSPCLT are constructed directly from the estimated correlation matrix using (A.1). In the projector construction stage, an orthogonal projector matrix such as Πn is constructed. And state vectors are updated using Kalman filter or Luenberger observer in the final stage of filtering. As described in Section 1, Kalman filter is used in PASTd Kalman, and Luenberger observer is used in DSPCLT and AMEND. In Table 4, the numbers of real floating point operations needed to carry out each of the three processes for each of the three schemes as well as the total numbers are presented.Table 4
Required floating point operations.
DSPCLTAMENDPASTd KalmanSubspace construction5M2N−K+1Ns64K2N−2KNsK28N+7NsProjector construction8K2N−K216N−KN−2K28KN2Filtering162N−K216N−K28N2Ns+8K+2N2Total5M2N−K+1Ns+8K+162N−K264K2N−2KNs+16N−KN−2K2+N−K28KN+8N2+7KNs+161+KN2For more concrete picture of computational complexity, we evaluated the numbers in Table4 with the choice of parameters used in obtaining Figures 9 and 11, the results of which are presented, respectively, in Tables 5 and 6. We note that a significant portion of computational complexity is needed in subspace construction. We note that AMEND requires particularly large numbers of operations in the subspace construction compared with DSPCLT. This is because AMEND uses matrix inversion and QR decomposition even though it does not use ED nor SVD. We note that DSPCM does not require matrix inversion nor QR decomposition. We note that PASTd Kalman, compared with DSPCLT and AMEND, requires larger amount of operations in the filtering stage since it uses computationally expensive Kalman filtering. Finally, we note that DSPCLT requires 2∼6 times more operations to carry out projector construction and filtering compared with AMEND. This is because DSPCLT has larger effective aperture size despite that the same number of sensors is used. However, the overall computational complexity of DSPCLT is significantly smaller than that of AMEND and PASTd Kalman because the computational complexity of the second stage is much smaller than that of either the first or the third stage. Overall, DSPCLT requires only about 11.99% (18.78%) and 7.70% (19.96%) real floating point operations compared with AMEND and PASTd Kalman, respectively, for the case of N=12,M=K=3,Ns=200 and N=20,M=K=5,Ns=200.Table 5
Required floating point operations:N=12,M=K=3,Ns=200.
DSPCLTAMENDPASTd KalmanSubspace construction66,000691,200205,800Projector construction10,5845,1843,456Filtering7,0561,296236,160Total83,640697,680445,416Table 6
Required floating point operations:N=20,M=K=5,Ns=200.
DSPCLTAMENDPASTd KalmanSubspace construction180,0003,200,000567,000Projector construction49,00024,00016,000Filtering19,6003,600662,400Total248,6003,227,6001,245,400


## 8. Conclusion

In this paper, we proposed the direct signal space construction method Luenberger tracker with quadratic least square regression for computationally efficient DoA tracking. Also, we studied analytically how the selection of observer gain value affects the system performance and described how we can use the theoretical results to optimize the observer gain value. We note that the proposed algorithm combines the direct signal space construction method and Luenberger observer to achieve computational complexity significantly lower than other schemes. Also, the algorithm achieves enhanced performance by employing a method of delay compensation which accounts for the difference in the times of estimation and data collection. The proposed scheme avoids the possible issue of nonresolvability by treating separately the situation of object overlap. As a result, the proposed scheme achieves performance similar or superior to existing algorithms with significantly lower computational complexity, which was verified with a variety of simulation results.



---

*Source: 4697381-2020-06-17.xml*