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

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. *e 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. *irdly, the proposed scheme achieves enhanced performance by a method of delay compensation, which accounts for observation delay. *rough numerical analysis, we show that DSPCLT achieves performance similar or superior to existing algorithms with only a fraction of computational requirement.


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][5][6], 2-dimensional arrays [7,8], and MIMO radar [9]. However, early subspacebased 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][14][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][19][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. e 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.
irdly, 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. is paper is organized as follows. In Section 2, 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. e 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. e following notations are used throughout this paper. Bold-faced letters are used for matrices and vectors. e complex conjugation, the transposition, and the Hermitian transposition of matrix A are denoted, respectively, by A * , A T , and A H . e M × N zero matrix and N × N identity matrix are denoted, respectively, by O M×N and I N .

System Model and Basic Assumptions
We consider a uniform linear array (ULA) consisting of N sensors with intersensor spacing d and assume that K narrowband far-field signal waves s 1 (t), s 2 (t), · · · , s K (t) of wavelength λ impinge on the array at distinct directions, respectively, making angles θ 1 (t), θ 2 (t), · · · , θ K (t) with the linear array. We denote, by x 1 (t), x 2 (t), · · · , x N (t), the received noisy signals collected at the N sensors. en, if we define the steering vector a(θ k (t)) by a(θ k (t)) � [1, e j] k (t) , · · · , e j(N− 1)] k (t) ] T with ] k (t) � 2πd cos θ k (t)/λ, the system can be modeled with the following equation: where , · · · , s K (t)] T , and n(t) 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) e 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) e vectors s(t) and n(t) are statistically independent complex random vector processes with zero mean satisfying where δ ts is the Kronecker delta function. We further assume that P is a diagonal matrix with positive real numbers p 1 , · · · , p K on its diagonal and that the value of N 0 is assumed to be estimated accurately by the receiver through long-term observation of the noise. (3) e system tracks the DoAs by obtaining angle estimates at time t � 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 − 1)T, nT] for each n � 1, 2, 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. e 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. en, 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. en, in the next two sections, we consider the estimation stage in nonresolvable and resolvable situations, respectively. e 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 by y k [n | n] the state vector at time t � nT and by θ k [n | n] the first component of y k [n | n]. We shall regard θ k [n | n] as the estimate of the DoA θ k (nT) of object k at time nT. Next, we denote by y k [n | n − 1] the predicted state vector obtained in the prediction stage at time t � nT and by θ k [n | n − 1] the first component of y k [n | n − 1].

Initialization Stage.
e initialization stage begins by obtaining the estimates of θ k (− 2T), θ k (− T), and θ k (0) for k � 1, 2, · · · , K, through a certain DoA measurement scheme such as MUSIC [1], ESPRIT [2], or DSPCM [14,15]. e estimates θ k (− 2T), θ k (− T), and θ k (0) shall be denoted, respectively, by en, the initial state vector y k [0 | 0] is computed by We note that the third component of y k [0 | 0] 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. is is because mathematical analysis in later sections is slightly simpler with (3).

Prediction Stage.
After initialization, the state vectors are updated periodically, with period T, in two stages called prediction and estimation. We assume that the state vector y k [n − 1 | n − 1] at previous time (n − 1)T is given, where n ∈ 1, 2, · · · { }. We assume that y k [n − 1 | n − 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 k [n | n − 1], are obtained by the following dynamic matrix equation: where Here, the first, second, and third components of the vector y k [n | n − 1] are interpreted to represent predicted values, based on the previous state vector, for θ k (nT), _ θ k (nT), and € θ k (nT), respectively. After the prediction stage, the distance is computed for l � 1, · · · , k − 1, k + 1, · · · , K. e object k is regarded to be in the nonresolvable situation if d kl 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.

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 vector y k [n | n] is obtained by using quadratic least square regression using the past p estimated DoAs θ k [m | m], m � n − p, n− p + 1, · · · , n − 1. To be more specific, we assume that the true DoA θ [m] satisfies and attempt to obtain the constants a, b, and c by the least square method with the previously estimated DoAs θ k [m | m], m � n − p, n − p + 1, · · · , n − 1. We note that the resultant constants a, b, and c are given by [27] Once the constants a, b, and c are computed by (8), the estimated state vector y k [n | n] is obtained by We note that the second and third components of y k [n | n] can be obtained by the first and the second de-

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 vector y k [n | n − 1] is obtained, in the prediction stage, by simply multiplying the previous estimated state vector y k [n − 1 | n − 1] by a constant matrix F without reflecting the data collected between times (n − 1)T 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 k [n | n] in the resolvable situation.
In particular, the updated estimated state vector y k [n | n] is obtained, using the Luenberger observer, by (10) where g k denotes a 3-dimensional column vector called observer gain [24]. Here, the quantity θ k [n | n] is an estimate of θ k (nT) obtained by a certain algorithm using the data collected at the antenna. In the following, we describe how we obtain θ k [n | n] in the proposed scheme. Before proceeding, we note that the performance of the system can depend on the choice of the observer gain g k . 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 θ k [n | n], which we call the delay compensated direct signal space construction method Newton estimate (DCDSPCMNE) of θ k (nT). e DCDSPCMNE θ k [n | n] is obtained in two steps. In the first step, we obtain a quantity � θ k [n | n] called the direct signal space construction method Newton estimate (DSPCMNE). en, in the second step, we obtain more accurate estimate θ k [n | n] by a method of motion compensation.

Direct Signal Space Construction Method Newton
Estimate. Let us first discuss how we obtain DSPCMNE � θ k [n | n]. In many subspace-based DoA estimation schemes, the DoAs are estimated by minimizing a certain cost function f(θ) defined by where Π is the orthogonal projector into the noise subspace. One major source of computational complexity in subspacebased 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 with f ′ (θ) [11,12,22,25]. We also adopt Newton's iteration method in this paper. Let θ denote the coarse estimate of a zero of f ′ (θ).
en, according to Newton's method, a refined zero θ is given by 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 � θ k [n | n] by where a k [n] and d k [n] are vectors defined, respectively, by and 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.

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 is an estimate of θ k (t), it is not exactly at time t � nT. In this section, we regard � θ k [n | n] as an estimate of θ k (t) at a certain time t � (n − α)T (0 < α < 1) and obtain the DCDSPCMNE θ k [n | n] 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 θ k (t) during this delay.
Before proceeding, we shall assume that T and € θ k (t) are small enough so that _ θ k (t) does not change appreciably over time period [(n − 1)T, 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.
T, and forgetting factor ρ is used to give more weight to more recent data. In the special case in which δ � 1/N s and ρ � 1, we note that the data are collected at times t � (n − 1 + 1/N s )T, (n − 1 + 2/N s )T, · · · , nT, and all the data are equally weighted since ρ � 1. Consequently, we observe that the data are collected symmetrically around time t � (n − (N s − 1)/2N s )T and are equally used. While it is theoretically difficult to justify the validity, we shall assess the value of α as (N s − 1)/2N s in this special case. e assessment of the value of α is more difficult in the general case. In the general case, considering the fact that the data are as the estimated observation time (n − α)T. Next, we note that accounts for the angle change between times (n − 1)T and (n − α)T. e change between times (n − 1)T and nT can be estimated as assuming that angular velocity does not change over [(n − 1)T, nT]. Using this approximation, we obtain the DCDSPCMNE θ k [n | n] by the following equation: Before proceeding, we note that (20) reduces to in the special case of δ � 1/N s .

Selection of Observer Gain
In Section 5, 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 g k . 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 − g k 1 0 0 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.

State Estimation Error.
We first define the state estimation error e k [n] by Mathematical Problems in Engineering 5 where y k (t) denotes the state vector defined as [θ k (t), _ θ k (t), € θ k (t)] T . Also, we define a vector δ k [n], which we call state evolution mismatch at time t � nT, by As discussed in [22], this vector δ k [n] equals to a zero vector if the angular acceleration does not change between times (n − 1)T and nT. However, the angular acceleration may change although the amount may be small, and δ k [n] becomes small but nonzero vector. We define the angle measurement error ε k [n] by In the next proposition, we obtain recurrence relation between e k [n] and e k [n − 1].

Proposition 1.
e state estimation errors satisfy the following recurrent relationship: Proof. We first note, from (10), (4), and (24), that Next, from (22), (23), and (26), it follows that Now, since θ k (t) � c T y k (t) and θ k [n | n − 1] � c T y k [n | n − 1], we next obtain 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 the k th object, which is defined by where N T is the total number of observation intervals. Before proceeding, we note that RMSE k can be expressed, using e k [n], as since First, we diagonalize the matrix F − g k c T F into VΛV − 1 , where Λ is a diagonal matrix whose diagonal elements are eigenvalues of F − g k c T F and V is a matrix whose column vectors are the corresponding eigenvectors. Here, we note that it is possible to represent g k in terms of Λ and V as described in the following lemma.
Here, we note that c T Fc is simply the (1, 1) component of the matrix F, which is 1. Consequently, we obtain 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.
where g ki denotes the i th element of g k .
Proof. We note that the characteristic equation of F − g k c T F is given by 6 Mathematical Problems in Engineering Consequently, for κ 1 , κ 2 e jϕ , and κ 2 e − jϕ to be solutions to (37), the identity equation must hold. Now, comparing the coefficients of the left and right-hand sides of (38) and solving for g ki , we obtain the results. □ □ Although the result in Lemma 1 is simple, it is important in that we can construct the observer gain vector g k by suitably choosing the eigenvectors and eigenvalues that directly affect the system performance. Before proceeding with the performance analysis, let us rewrite (25) as 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 ε k [m] and mismatch δ k [m] with m close to n affect the values of e k [n] since Λ k converges to a zero matrix as k grows. Moreover, we note that the magnitude of δ k [m] is order of magnitudes smaller than that of ε k [n] 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 e k [n] can be approximated as and attempt to compute RMSE k using (31). To obtain reasonably simple and useful result for RMSE k , we shall assume that the angle measurement error ε k [m] has zero mean and that ε k [m] and ε k [n] 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
and that en, where

Mathematical Problems in Engineering
Proof. See Appendix B. □ □ Corollary 1. If we further assume, in eorem 1, that there exists some nonnegative real number σ k such that for all n, then where 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 σ 2 k [n] 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 RMSE k in relation to the choice of Λ. In particular, we note that there is no further information required about the variation of σ 2 k [n] over n, and the RMSE k is represented in the unit of σ k in (47). Secondly, we can regard (47) with as a simple and reasonable approximation of (44) since there is no reason to believe that σ 2 k [n] 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 , κ 2 e jϕ , and κ 2 e − jϕ and compute the right-hand side of (47). Apparently, we need to obtain V for this purpose. However, we can compute F − g k c T F, 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 that e k [n] 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 RMSE k value, which is 0.8646σ k . In the following, we shall express RMSE k in the unit of σ k for simplicity. In Figure 1, we illustrate the RMSE k values as functions of κ 1 and κ 2 for the case of ϕ � 35°. In this figure, we note that the RMSE k value rises as (κ 1 , κ 2 ) moves away from (0.7, 0.5). In particular, we note that the RMSE k value rises relatively mildly if κ 1 or κ 2 is decreased. However, the RMSE k value rises significantly as κ 1 and κ 2 approach 1. We found that the RMSE k 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 RMSE k value changes with the ϕ value.
So far, we discussed how we can choose the optimal observer gain vector with (45) assuming that σ k [n] 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 σ k [n] affects the value of RMSE k . According to our extensive test, the effect of variation of σ k [n] 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 σ k [n] 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 σ k [n], we assume that only object k moves following the trajectory θ k (nT) 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 σ k [n].
For numerical evaluations, we assume that the number N of sensor elements is 12, the total number N T of observation intervals is 60, and the signal-to-noise ratio is 5 dB. With 10 6 simulation run trials, we computed σ 2 k [n], which is illustrated in Figure 5. In this figure, we note that the value of σ 2 k [n] is generally related to the angular acceleration. After obtaining the estimate of σ 2 k [n], we computed RMSE k using (44) in the unit of σ k � To compare the results with the values computed using (47), we evaluated the ratio obtained by dividing the RMSE k value computed using (44) by that using (47)  To understand what happens, let us first consider the case of ϕ � 15°and σ 2 5 [n] 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 σ 2 5 [n]. We observe that the minimum value 0.799 happens near the region κ 1 ≈ 1 and κ 2 ≈ 1 for which the RMSE k values are very large. In contrast, the ratio is very close to 1.0 except for the region in which RMSE k 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,  Mathematical Problems in Engineering the ratio for ϕ � 35°and σ 2 1 [n]. 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 RMSE k value computed with (47) is very close to that with (44) except at κ 1 , κ 2 , and ϕ values for which e k [n] does not converge to zero fast enough, and RMSE k is very large.

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.

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: where θ k (nT) is the estimated DoA of θ k (nT) and N T is the the total number of observation intervals. We note that RMSE k [n] is the RMSE of the DoA of the k th object at time t � nT and that RMSE T represents the total RMSE over all objects and observation intervals. We recall that we regard θ k [n | n] as θ k (nT) in the proposed scheme.
However, we need to be very careful in assessing the performance with RMSE.
is 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. e 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 than 5°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 P F of tracking failure as a performance measure.
For performance evaluations, we first consider a situation in which K � 3 signals impinge upon a ULA with N � 12 sensors separated by half-wavelength from time-varying directions θ 1 (nT), θ 2 (nT), and θ 3 (nT), n � 1, 2, . . ., N T � 60, illustrated in Figure 8. We assume that the signalto-noise ratio of each signal at each sensor is given as 5 dB. During each time interval [nT, (n + 1)T), total N s � 200 data snapshots are uniformly captured at each sensor with δ � 1/N s .
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 ρ N s − 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.75e jπ/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.5e j35π/180 and 0.5e − j35π/180 for DSPCLT. We recall that observer gain is not used in PASTd Kalman.
As the value of M 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 with 10 5 trials. e failure probability P F and RMSE T 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 RMSE T value. In particular, we note that both DSPCLT and PASTd Kalman exhibit not only better RMSE T values than AMEND but also achieve zero failures. e best RMSE T value for each scheme is marked in bold-faced fonts, and we note that DSPCLT achieves the lowest RMSE T value among the three schemes. By the way, the RMSE T value 0.6160 listed on the lowest row for PASTd Kalman corresponds to the forgetting coefficient n/(n + 1) suggested in [21].
To illustrate the performance variation over time, RMSE k [n] 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 RMSE T values. Also, we note DSPCLT achieves slightly better RMSE T value compared with PASTd Kalman. In Figure 9, we observe that DSPCLT exhibits best overall RMSE k [n] 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. is is because the nonresolvable situation is separately treated, in the proposed scheme, using the algorithm described in Section 4.
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 Figure 10, 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.
We next study the performance variation with increased number of objects. For illustration, we consider the situation in which there are K � 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 N s of snapshots is chosen to be 200 with δ � 1/N s , and the value of M is chosen to be the same as K.
e simulations are performed with 10 5 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 RMSE T 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 RMSE T . In particular, we observe that the RMSE T value of DSPCLT is significantly lower than that of AMEND despite the fact that RMSE T is computed excluding the cases in which tracking fails. Next, in Figure 11, the sequence of RMSE k [n] 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.

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.
For more concrete picture of computational complexity, we evaluated the numbers in Table 4 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.
is 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. is 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%

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. e 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. r 1 , · · · , r N− 1 using the relation r 2N− m � r * m , m � 1, 2, · · · , N − 1. Next, we obtain K orthonormal vectors r 1 , · · · , r K by applying Gram-Schmidt procedure on r 1 , · · · , r K , where r k 's are N-dimensional column vectors defined by r k � [r k , r k+1 , · · · , r k+L− 1 ] T . Next, we obtain Π[n] by Π[n] � I L − r 1 r T 1 + · · · + r K r T K . (A.7) the numerical values of θ k (nT) in Figures 4 and 8; (2) 'sigma2_k.xlsx' includes the numerical values of σ 2 k [n] in Figure 5. (Supplementary Materials)