Adaptive realization of a maximum likelihood time delay estimator

ABSTRACT This paper presents an adaptive maximum likelihood method for estimating the time difference of arrival of a source signal at two spatially separate sensors. It is well-known that the maximum likelihood technique achieves the Cramer-Rao lower bound for time delay estimation error for certain signal conditions. The a-β tracker is a heuristic mechanism that is heavily used in target tracking applications. In this work, we combine an adaptive realization of the maximum likelihood time delay estimator with the a-β tracker to obtain significant improvement in the performance of the tracker. Experimental results showing 2 to 8 dB improvement in the mean-square estimation error over the conventional a-β tracker for various signal-to-noise ratios are also included in the paper.

ADAPTIVE REALIZATION OF A MAXIMUM LIKELIHOOD TIME DELAY ESTIMATOR Peter J. Hahnl V. John Mathewsl Thao D. Tran2 lDepartment of Electrical Engineering, University of Utah, Salt Lake City, UT 84112, USA 2Digital Systems Resources, 12450 Fair Lakes Circle, Suite 500, Fairfax, VA 22030, USA ABSTRACT This paper presents an adaptive maximum likelihood method for estimating the time difference of arrival of a source signal at two spatially separate sensors. It is well-known that the maximum likelihood technique achieves the Cramer-Rao lower bound for time delay estimation error for certain signal conditions. The a-(3 tracker is a heuristic mechanism that is heavily used in target tracking appli­cations. In this work, we combine an adaptive realization of the maximum likelihood time delay estimator with the a-(3 tracker to obtain significant improvement in the per­formance of the tracker. Experimental results showing 2 to 8 dB improvement in the mean-square estimation error over the conventional a-(3 tracker for various signal-to-noise ratios are also included in the paper. 1. INTRODUCTION This paper presents an adaptive maximum likelihood method for estimating the time difference of arrival of a source signal at two spatially separate sensors. The signals arriving at the two sensors from a target are modeled in this work as (1) and X2(n) = s(n - D(n)) + 172(n), (2) where s(n) is the source signal and 17l(n) and 172(n) are noise components that have zero-mean values and are un­correlated with each other and the source signal. The com­ponent of the source signal arriving at the second sensor is delayed by an amount D(n). We assume in this work that the amplitudes of the signals arriving at the two sensors are approximately the same. Note that the delay can vary with time, especially when the receiver array or the target is moving. The basic approach for estimating the time delay is to estimate the cross-correlation function of the two signals and then to find the time lag at which this estimate peaks [1]. Using the model for the received signals in (1) and (2) and assuming stationary source signals with constant delay D, we see that the cross-correlation function Rl2 (m) 'This work was supported in part by an IBM Departmental Grant and a University of Utah Research Fellowship. 0-7803-3192-3/96 $5.00©1996 IEEE is given by Rl2 (m) = t'{aol(n)x2(n - m)} = t'{s(n)s(n + D - m)} (3) = Rs,,(m - D) , where t' {-} denotes the statistical expectation of {.} and Rss (') is the autocorrelation function of s(n). Since the autocorrelation function peaks at zero lag, we can see that the time delay can be estimated by finding the location of the peak of R12(m). The generalized cross-correlation (GCC) method [1, 6] of estimating the time delays involves smoothing the cross­correlation function using certain weighting functions in the frequency domain. Let G]2 (J) denote the cross-spectral density function of xl(n) and x2(n). Then, the generalized cross-correlation function of Xl (n) and X2 (n) is defined as (4) where W(J) is a weighting function and F-l {.} denotes the inverse-Fourier transform of {-}. The weighting function is selected to reduce the eJfect of noise on the estimates. The maximum likelihood method for finding the estimate of the time-delay between two spatially separate sensors was developed in [2, 3, 6]. This method was found to achieve the Cramer-Rao lower bound for Gaussian input statistics. The weighting function W(M)(J) for this estimator is (5) where Afl2 (J) is the complex coherence function of Xl (n) and x2(n) defined as (6) The generalized cross-correlation method with this maxi­mum likelihood weighting function is refered to as the ML­GCC method in this paper, The a-(3 tracker [7, 8] is: a heuristic mechanism that is heavily used in target tracking applications. In this paper, we combine an adaptive realization of the maximum likeli­hood time delay estimator with the 0l-f3 tracker to obtain significant improvement in the performance of the tracker. The adaptive realization of the ML-GCC function is similar 3121 to the method proposed by Youn and Mathews [9], and is described in Section 3. One novel feature of our realization is that the step size of the adaptive filters is updated online based on the input signal-to-noise measurement to provide as close to the best possible performance without manual su pervision. The rest of this paper is organized as follows. The next section discusses the 0I-{3 tracker. The third section dis­cusses the adaptive realization of the ML-GCC function. The adaptive maximum likelihood 0I-{3 tracker algorithm is developed in Section 4. Experimental results are given in Section 5. The concluding remarks are made in Section 6. 2. THE 0I-j3 TRACKER The 0I-{3 tracker [7, 8] is a block adaptive algorithm that attempts to estimate the time delay from the peak of an estimate of the cross-correlation function and tracks the changes in the time delay over time. The 0I-{3 tracker has three major components: (1) the cross-correlation function estimation, (2) the asymmetry error estimation, and (3) the parameter update. The cross-correlation function estima­tion block generates an estimate of the normalized cross­correlation function between the two input signals over a contiguous block of data. The cross-correlation function for the kth block is estimated as (7) where N is the number of samples in each block and L repre­sents the offset between the first samples of successive input blocks. Let f( k) denote the estimate of the time delay for the kth block. Assuming that this estimate is reasonably ac­curate, we can expect that the normalized cross-correlation function estimate for this block given by R' ( k) _ R12 (m, k) 12 m, - , J Rll (0, k)R22(O' k) (8) peaks in the vicinity of f( k). The asymmetry error c(k) = RIZ(f(k) + 1, k) - R12(f(k) - 1, k) (9) provides a measure of the error between the actual time delay r(k) and the estimate of the time delay f(k) [3, 4]. The 0/-{3 tracker iteratively estimates the time delay f(k) and its derivative +( k) as i-(k + 1) = r(k) + (31~) e:(k) (10) and f(k + 1) = f(k) + t::.tr(k + 1) + OI(k)c(k) , (11) after every block of N samples. The parameter t::.t in the above equations is the time between updates. The param­eters OI(k) and (3(k) control the rate of convergence and tracking of the system, and are updated heuristically [7, 8]. The 01 and (3 update equations are adapted on the basis of an estimate of the input signal-to-noise ratio (SNR) [7, 8]. The estimate of the SNR is found from the peak of the normalized cross-correlation estimate [7, 8] as (12) where f( k) is the location of the peak of the normal­ized cross-correlation function and we have assumed that the SNRs of the two received signals are very close to each other. The cross-correlation function estimate R12 (f( k), k) is found by interpolating from the original es­timates R12 (m, k) with a polynomial interpolation scheme [7, 8]. The heuristic argument employed for updating OI( k) and (3( k) is that when the SNR is large, the algorithm should have a small effective integration time which allows track­ing of higher bearing rate targets. The parameter (3( k) is updated as (3(k) = ;3ogK(K(k»g.(c:(k),K(k» , (13) where ;30 is a weighting constant, K(k) is a smoothed estimate of the value of the peak of the normalized cross-correlation function R12(f(k), k), c(k) is a smoothed estimate of the asymmetry error, and g}((K(k» and g.(c:(k), K(k» are gain parameters. Both K(k) and c:(k) are smoothed with a single pole filter. The gain g}(K(k» is given by I K(k) gK(K(k» = V Ko ' (14) where Ko is the value of the peak of the normalized cross­correlation function for the lowest useful SNR expected. The gain g.(c:(k), K(k» is given by g.(c(k), K(k» = 1 + 31c:(k)1 (15) J100K~ + K2(k) The sequence OI(k) is a function of j3(k) and K(k) and is given by OI(k) = (16) The details of the derivations as well as the heuristics of the above equations are described in [7, 8]. Our experience is that this method, in spite of its heuristic nature, outper­forms almost all adaptive time delay estimation techniques available in the literature for the types of signals we have tested the procedures on. 3. THE ADAPTIVE ML-GCC FUNCTION ESTIMATION The adaptive realization of the ML-GCC function [9] in­volves the use of two normalized LMS adaptive filters. The first filter estimates Xl (n) using X2 (n) and the second es­timates x2(n) using xl(n). Let k12 (m,n) and h21 (m,n) represent the coefficients of the two adaptive filters at time n. These coefficients are updated as [4] 3122 Adaptive ML-GCC Estimator Asymmetry Error E(k) u-/3 Update ~k) Figure 1. Block diagram of the adaptive ML algo­rithm. where f.L( n) is the step size sequence that controls the speed of convergence and steady-state characteristics of the adaptive filter, hij (n) is the vector of the adaptive fil­ter coefficients given by hij(n) = [kij(-M,n),kij(-M + 1, n), ... , hij(M, n)V, and Xj( n) is the input vector defined as Xj(n) = [xj(n + M), xj(n + M - 1), ... , xj(n - M)f. The error eij(n) between the desired signal xi(n) and the input block Xj(n) is given by (18) It was shown in [9) that represents an estimate of the magnitude-squared coher­ence (MSC) function of the input signals at time n. In (19), Hij(f, n) represents the discrete-time Fourier trans­form (dtFt) of the NLMS filter coefficients hij (m, n) at time n. The ML-GCC function may now be estimated as R(M\ )-F-1{ 1-'Y12(f,nW H12(f,n)} (20) 12 m,n - (1-1-'Y12(f,nW)IH12(f,n)I' 4. THE ADAPTIVE ML a-{3 TRACKER The adaptive ML a-{3 tracker combines the ML-GCC func­tion estimation algorithm with the a-{3 tracker algorithm. The system first finds an adaptive estimate of the ML-GCC function as described in the previous section. The estimate of the ML-GCC function is obtained using (19) and (20) once every L samples. The peak location of the ML-GCC function estimate f(k) for the kth iteration of the a-{3 tracker is found using an iterative bisection al­gorithm [5) that finds the zero of Instead of using the estimate of the time delay error given by (9), our method computes the estimate of the error as e(k) = r(k) - r(k) . (22) This estimate of the error is used in place of the asymmetry error in the a and (3 update equations (10) and (11). The sequences a(k) and (3(k) are still updated using (13) and (16). A block diagram of the system is shown in Figure 1. One major difficulty in implementing this algorithm is in the selection of the step size for the adaptive filters. The problem of choosing an optimal step size value or an appro­priate step size sequence is not easy since there are no ana­lytical results that relate the step size to the performance of Figure 2. A three-dimensional plot of the optimal step size f.L as a function of SNR and r. the time delay estimators. In stationary operating environ­ments, a smaller step size results in a better steady-state estimate, but slower adaptation. A larger step size con­versely results in a poorer steady-state estimate, but faster adaptation. In stationary environments, one desires an ini­tially large step size and then a smaller step size after the transient stage. A non-stationary environment requires that the value of step size be large enough to track the changes in the operating environment, yet small enough that the estimate is accurate. Our system employs an empirically derived step-size up­date algorithm for the adaptive filters. A large number of experiments were conducted to find the value of the step size that minimizes the mean-square error (MSE) between the estimate of the time delay and the actual time delay for a large range of SNRs aIlld rates of variation of the time delay. The different rates of change in the actual time delay were simulated using linearly varying time delay functions with different slopes. For each set of SNR and i-, simu­lations were performed for a large set of step sizes in the range [4.0 xl0-5 , 5.0 x 10-3 )., The step size resulting in the smallest average MSE in the time delay estimate was tabu­lated for each of the representative environments. A three­dimensional plot of the "optimal" f.L value for the different SNR and rate of change in the time delay ("taudot") are shown in Figure 2. The ML adaptive algorithm uses the tabulated experimental information by choosing the step size that provided the best performance for the estimate of the SNR and the rate of change of the time delay function. The SNR was estimated using the technique from (12) on the same block of L samples using the estimate of the cross­correlation function given in (7) and (8) separately from the estimation of the ML-GCC function. The estimate of the rate of change in the time delay is given in the difference equation (10). The step size parameter is kept constant between updates of the a-{3 tracker. 5. SIMULATION RESULTS The conventional and adapHve ML a-{3 trackers were sim­ulated using 30 Monte-Carlo runs with a. ba.ndpass filtered Gaussian signal and additive white Gaussian noise at 6 dif- 3123 Oti9 Adapt ML "0-._ o SNR (dB) ,0 Figure 3. Comparison of the original a-f3 tracker and the maximum likelihood technique for a con­stant delay. ferent SNR levels. The source signal belonged to a zero­mean Gaussian process with a fiat spectrum in the normal­ized frequency range [0.06,0.25] cycles/sample. The a-f3 tracker used a block size of 1000 samples with no overlap of the samples. The time-varying parameters of the adaptive ML 0'-13 tracker and the estimate of the time delay were up­dated every 1000 samples. The two normalized LMS filters each had 51 coefficients. The value for 130 = 7.22 X 10-3 was choosen to result in a maximum possible f3( k) of 0.3. The value for Ko = 0.0245 was choosen to match the peak of the cross-correlation function when in the environment with the smallest expected SNR of -16 dB. The experiments were conducted for SNR values ranging from 10 dB to -15 dB in steps of -5 dB. Both a linear and constant delay function were used to evaluate the "steady-state" response of the systems. Each experiment used 390,000 samples of data. The constant delay function used in the experiments was at two sampling periods. The linear delay function started at two sampling periods, and decreased by one sampling pe­riod every 100,000 samples. In each experiment the system was initialized to the proper time delay and rate of change of the time delay. The average mean-square error between the estimate of the time delay and the actual time delay was calculated over the 30 runs for the last 300,000 sam­ples. The MSE of the time delay estimates are ploted as a function of the input SNR in Figure 3 for the constant delay and in Figure 4 for the linear delay. The curves labeled as "Orig" refer to the original a-f3 tracker. The curves labeled as "Adapt ML" refer to the adaptive maximum likelihood realization. Our technique improved the MSE by at least 2 dB and by as much as 5 dB in the non-stationary environ­ment, and it improved the MSE by at least 5 dB and by as much as 8 dB in the stationary environment. 6. CONCLUDING REMARKS This paper discussed the estimation of the time-difference of arrival of a signal between two spatially separate sensors. A new time-delay estimation algorithm based on the a-f3 tracking mechanism and the maximum likelihood estimator was presented. Experimental comparisons showed that the Orig AdaptML ·5 0 SNR (dB) 'a. Figure 4. Comparison of the original a-f3 tracker and the maximum likelihood technique for a linear delay function. new tracker provides significant performance improvement over the conventional a-f3 tracker. A novel component of the system is an empirically derived algorithm for selecting the step size of the adaptive filters online and without manual supervision. 3124 REFERENCES [1] G. C. Carter, Ed., "Special Issue on Time Delay Esti­mation", IEEE Trans. on Acoustics, Speech, and Sig­nal Proc., Vol. ASSP-29, June 1981. [2] E. J. Hannan and P. J. Thomson, "The Estimation of Coherence and Group Delay", Biometrica, Vol. 58, No.3, 1971. [3] E. J. Hannan and P. J. Thomson, "Estimating Group Delay", Biometrica, VoL 60, No.2, 1973. [4] S. Haykin, Adaptive Filter Theory (Second Edition), Prentice Hall, Englewood Cliffs, N J, 1991. [5] D. Kincaid and W. Cheney, Numerical Analysis, Brooks/Cole, Pacific Grove, CA, 1991. [6] C. H. Knapp and G. C. Carter, "The generalized cor­relation method for estimation of time delay," IEEE Trans. on Acoustics, Speech, and Signal Proc., VoL ASSP-24, pp. 320-327, Aug. 1976. [7] C. J. Wenk, "A Design Methodology for Continuous Gain Adaptive PBB Time Delay Tracking", Analysis and Technology, Inc. Report No. P-1717-2-84, Mar. 1984. [8] C. J. Wenk, "Continuous Gain Adaptive Time Delay Tracking", Analysis and Technology, Inc. Report No. P-2078-3-86, May 1986. [9] D. H. Youn and V. J. Mathews, "Adaptive Realization of the Maximum Likelihood Processor for Time Delay Estimation", IEEE Trans. on Acoustics, Speech, and Signal Proc., Vol. ASSP-32, Aug. 1984.