FFT-Based Efficient Algorithms for Time Delay Estimation

Suppose that we have a single sensor receiving a superimposition of attenuated and delayed replicas of a known signal plus noise. From the received data we want to estimate the arrival times of the various replicas and their (complex or real) attenuation coefficients (gains). This is the well-known time delay estimation problem which occurs in many fields including radar, active sonar, wireless communications, Global Navigation Satellite System (GNSS), nondestructive testing, geophysical/seismic exploration, and medical imaging. In this chapter, we will focus on time delay estimation based on one sensor with known probing signal shapes.


Introduction
Suppose that we have a single sensor receiving a superimposition of attenuated and delayed replicas of a known signal plus noise. From the received data we want to estimate the arrival times of the various replicas and their (complex or real) attenuation coefficients (gains). This is the well-known time delay estimation problem which occurs in many fields including radar, active sonar, wireless communications, Global Navigation Satellite System (GNSS), nondestructive testing, geophysical/seismic exploration, and medical imaging. In this chapter, we will focus on time delay estimation based on one sensor with known probing signal shapes.
The most well-known time delay estimator is the matched filter approach. If there is only one signal or the overlapped signals are separated in time by an interval that is much greater than the width of the signal autocorrelation function, then the matched filter is the optimal estimator when the noise is white Gaussian (Ehrenberg et al., 1978) . The resolution capability of the matched filter approach depends on the signal bandwidth and the larger the signal bandwidth, the better the resolution. However, in many situations there exist some practical limitations on increasing the bandwidth of the transmitted signals. How to resolve closely spaced overlapping noisy echoes has attracted the attention of researchers from many fields for several decades.
Several approaches have been suggested for this problem and many of them benefit from the recent development of high resolution sinusoidal frequency estimation and Direction of Arrival (DOA) estimation techniques. Sinusoidal frequency estimation techniques such as MUSIC (Schmidt, 1986), Linear Prediction (Tufts et al., 1982), and MODE (Stoica & Nehorai, 1990) are applied to the time delay estimation problem in (Bian & Last, 1997;Kirsteins, 1987;Kirsteins & Kot, 1990). However, these approaches are only applicable to signals with flat (rectangular) band-limited spectra. Several Maximum Likelihood (ML) approaches have also been suggested for this problem. Multidimensional global optimization algorithms are presented in (Bell & Ewart, 1986;Blackowiak & Rajan, 1995; to analyze a special class of ocean acoustic data that has very oscillatory autocorrelation functions. An efficient approach based on the Expectation Maximization (EM) algorithm (Moon, 1996) is proposed in (Feder & Weinstein, 1988) that decouples the complicated multidimensional optimization problem into a sequence of multiple separate one-dimensional optimization problems. However, its convergence depends highly on the initialization method used and 3 www.intechopen.com no systematic initialization method is given in (Feder & Weinstein, 1988). More recently, time delay estimation lower bounds  and time delay estimation for small-samples (Gedalyahu & Eldar, 2010) are also studied.
In this chapter, a family of relaxation-and FFT-based new algorithms was proposed for different scenarios. The remainder of this chapter is organized as follows. In Section 2, we formulate the problem of interest. Section 3 describes the new algorithms. In Section 4, some applications utilizing the new algorithms are provided. Section 5 concludes the whole chapter.

Data model and problem formulation
Time delay estimation is a well-known traditional problem occurring frequently in radar, active sonar, and many other fields. In this problem, the waveform received at a single sensor consists of delayed replicas of the transmitted signal with different gains. The gains reflect the scattering property of the targets or multipath channel transmission features. The received signal waveform y(t) can be described as ( 1 ) where s(t),0≤ t ≤ T 0 , represents the known transmitted signal (complex or real valued), y(t) denotes the received signal, which is composed of L replicas of s(t) with different (complex or real valued) gains {α l } L l=1 and real valued delays {τ l } L l=1 ,ande(t) is the receiver noise, which is modeled as a zero-mean Gaussian random process.
Usually, the above received analog signal is sampled for digital signal processing. To avoid aliasing, we must sample y(t) according to the bandwidth of s (t).L e tB s denote the double-sided bandwidth of s(t).T h e ny(t) must be sampled with the sampling frequency f s satisfying f s ≥ B s . ( 2 ) After A/D conversion, the sampled received signal has the form y(nT s )= L ∑ l=1 α l s(nT s − τ l )+e(nT s ), n = 0, 1, ··· , N − 1, where T s is the sampling period and is equal to the reciprocal of the sampling frequency f s .
Our problem of interest herein is to estimate {α l , τ l } L l=1 from {y(nT s )} N−1 n=0 with known s(t), 0 ≤ t ≤ T 0 ,or{s(nT s )} N−1 n=0 . Although we could solve the estimation problem in the time-domain (Bell & Ewart, 1986;Blackowiak & Rajan, 1995;Bruckstein et al., 1985;Feder & Weinstein, 1988), we shall consider below solving the problem in the frequency domain and propose a family of relaxation-based algorithms for different scenarios.
Let Y(k), S(k),a n dE(k), k = −N/2, −N/2 + 1, ..., N/2 − 1, denote the discrete Fourier transforms (DFT's) of y(nT s ), s(nT s ),a n de(nT s ), respectively. Provided that aliasing is negligible, then Y(k) canbewrittenas: Note that the time delay estimation problem is similar to the sinusoidal parameter estimation problem except that the exponential signals are weighted by the known signal spectrum. If we divided both sides of (4) by S(k), the problem would become identical to the sinusoidal parameter estimation problem. Yet we should not do so for the following reasons: first, S(k) could be zero for some k; second, the noise E(k)/S(k) will no longer be a white noise even when E(k) is white; third, when E(k) is a white noise, the larger the S(k) at sample k,t h e higher the signal-to-noise ratio (SNR) of the corresponding Y(k) and hence dividing Y(k) by S(k) will de-emphasize those Y(k) ′ s that have high SNRs. Because of this, many well-known sinusoidal parameter estimation algorithms, such as MUSIC (Schmidt, 1986), ESPRIT (Roy et al., 1986), PRONY (Kay, 1988), are not directly applicable to our problem of interest. Using MODE (Stoica & Nehorai, 1990) would require a multidimensional search over a parameter space because we can no longer reparameterize the MODE cost function via the coefficients of a polynomial.

W eighted Fourier transform and RELAX ation (WRELAX) algorithm
We consider below estimating the unknown parameters by minimizing the following nonlinear least squares (NLS) criterion: When e(nT s ) is a zero-mean white Gaussian random process, E(k) is also white since DFT is a unitary transformation. For this white noise case, the NLS approach is the same as the ML method. When E(k) is not white, however, the NLS approach is no longer the ML method. However, it has been shown in (Li & Stoica, 1996) that the NLS approach can still have excellent statistical accuracy.
Minimizing C 1 ({α l , ω l } L l=1 ) with respect to the unknown parameters is a highly nonlinear optimization problem. The cost function has a complicated multimodal shape with a very small attraction domain, which makes it very difficult to find the global minimum. Below, we present a relaxation based optimization algorithm to obtain the NLS parameter estimates. Before we present our approach, let us consider the following preparations. Let and a(ω l )= e jω l (−N/2) e jω l (−N/2+1) ··· e jω l (N/2−1) T , where (·) T denotes the transpose. Denote Then (6) becomes where · denotes the Euclidean norm. Minimizing C 2 (α l , ω l ) with respect to α l yields the estimateα l of α lα where (·) H represents the conjugate transpose and · F denotes the Frobenius norm (Stewart, 1973). Then the estimateω l of ω l is obtained as follows: w h e r ew eh a v eu s e dt h ef a c tt h a tb H (ω l )b(ω l )= S 2 F and hence is independent of ω l .H e n c eω l is obtained as the location of the dominant peak of the magnitude squared of the Fourier transform, |a H (ω l )(S * Y l )| 2 , which can be efficiently computed by using the fast Fourier transform (FFT) with the weighted data vector S * Y l padded with zeros. An alternative scheme to zero-padding FFT is to find an approximate peak location first by using FFT without much zero-padding and then perform a fine search nearby the approximate peak location by, for example, the fminfunction in MATLAB, which uses the Golden section search algorithm. With the estimate of ω l at hand,α l can be easily computed from the corresponding complex height:α With the above simple preparations, we now present the WRELAX algorithm.

www.intechopen.com
Iterate the previous two substeps until "practical convergence" is achieved (to be discussed later on).
Iterate the previous three substeps until "practical convergence".
Remaining Steps: Continue similarly until L is equal to the desired or estimated number of signals. (Whenever L is unknown, it can be estimated from the available data, for instance, by using the generalized Akaike information criterion (AIC) rules which are particularly tailored to the WRELAX method of parameter estimation. See, for example, (Li & Stoica, 1996).) The "practical convergence" in the iterations of the above WRELAX method may be determined by checking the relative change of the cost function C 1 ({ω l ,α l } L l=1 ) in (6) between two consecutive iterations. The algorithm is bound to converge to at least some local minimum point (Karmanov, 1977). The convergence speed depends on the time delay spacing of the signals. If the spacing between any two signals is larger than the reciprocal of the signal bandwidth, the algorithm converges in a few steps. As the spacing of the signals becomes closer, the convergence speed becomes slower.
Once we have obtained the estimates {ω l } L l=1 , the estimates {τ l } L l=1 of {τ l } L l=1 can be determined by using (5).
At this point, we would like to point out the relationship between WRELAX and the conventional matched filter approach. The matched filter approach can also be formulated in the frequency domain. Let The matched filter method searches for the L largest peak positions of F(ω) as the estimates of {ω l } L l=1 , and then the gains are determined as followŝ Hence when there is only one signal, this one-dimensional matched filter approach is equivalent to the WRELAX algorithm. However, when there are multiple signals that are not well separated, this conventional matched filter approach will perform poorly. In this case, a multidimensional matched filter method (Bell & Ewart, 1986) could be used and the method is equivalent to the NLS fitting approach (Bell & Ewart, 1986). The WRELAX algorithm decouples the multidimensional matched filters into a sequence of one-dimensional matched filters. Thus the excellent parameter estimation performance of the NLS fitting approach can be achieved at a much lower implementation cost.
Similar to the WRELAX algorithm, the EM algorithm proposed in (Feder & Weinstein, 1988) also transforms the multidimensional optimization problem into a series of one-dimensional optimization problems. The detailed implementations of the algorithms, however, are quite different. The EM algorithm consists of two steps, the E (Estimate) step and the M (Maximize) step. The idea is to decompose the observed data into their signal components (the E step) and then to estimate the parameters of each signal component separately (the M step). The algorithm is iterative, using the current parameter estimates to decompose the observed data.
At each E step, the residue error corresponding to the current estimates is also decomposed among different signal components. Although initial conditions are needed by EM, no systematic initialization method is given in (Feder & Weinstein, 1988). We have also found that the performance of EM is very sensitive to the initial conditions used. Even with the same initial conditions, our numerical examples show that the convergence speed of EM can be much slower than the last step of WRELAX. Further, WRELAX does not require any initial conditions before its iterations and the first L − 1 steps of WRELAX can provide an excellent initial condition for Step L.
Consider next the case where {α l } L l=1 are real-valued. Minimizing C 2 (α l , ω l ) with respect to α l and ω l yieldsα where Re(X) denotes the real part of X,and The WRELAX algorithm could also be implemented in the time domain, which is based on the correlations. However, we prefer to use the frequency domain version of WRELAX. For the time domain version, we could be restricted to use the discrete values of {τ l } L l=1 if we only know the sampled version of s(t). For this case, if a more accurate delay estimate is required, then one has to resort to interpolation (Bell & Ewart, 1986). This inconvinence can be avoided by transforming the problem to the frequency domain, where {τ l } L l=1 can take on a continuum of values. Even without considering the additional interpolation cost, the computational load of the time domain correlation-based WRELAX is heavier than that of the frequency domain WRELAX.
We present several numerical examples to demonstrate the performance of the proposed algorithms.
In the following examples, we use a windowed chirp signal, s(t)= w(t)e jβ(t− T 0 2 ) 2 ,0 ≤ t ≤ T 0 ,w h e r eβ is the chirp rate and w(t) is a bell-shaped window function. We use N = 64, β = π × 10 12 , the signal bandwidth B s = βT 0 /π, and the sampling frequency f s = 2B s . T 0 is chosen in such a way that T 0 =(N/2 − 1)T s . The resolution limit of the conventional matched filter method is around τ e = 1/B s . L = 2 signals are assumed to be superimposed together with α 1 = e jπ/8 , α 2 = e jπ/4 . In the first example, the additive noise is zero-mean white Gaussian and SNR=SNR1=SNR2=10dB. The time delay spacing between the two signals is τ 2 − τ 1 = τ e . Even in this case, the conventional matched filter method fails to resolve the two signals, as can be seen from Figure 1(a), where the horizontal axis denotes the normalized time delay τ/T and the two vertical lines indicate the true time delays of the two signals. However, using WRELAX we can resolve them very well. As pointed out before, the WRELAX algorithm can be viewed as transforming the multi-dimensional matched filters into a sequence of one-dimensional matched filters. The outputs of the two matched filters for all iterations are plotted in Figure 1 is smaller than the corresponding true gains). After several steps, they converge to the true time delays and gains.
In this example, the time delay spacing between the two signals is τ 2 − τ 1 = 0.5τ e .Th eMS Es for the first signal using WRELAX are compared with the corresponding CRBs in Figure 2 and the MSE and CRB curves for the other signal are similar. From Figure 2, it can be noted that the MSEs obtained by using WRELAX approach the corresponding CRBs as the SNR increases.

T oeplitz property based W eighted Fourier transform and RELAX ation (TWRELAX) algorithm
Here, we extend the above WRELAX algorithm to the case of multiple looks. Two scenarios will be considered, which include 1) fixed delays but arbitrary gains and 2) fixed delays and fixed gains. In radar applications, the two cases correspond to two different target fluctuation models, e.g., Swerling I and Swerling II (Barton, 1988).

Fixed delays but arbitrary gains
Consider the case where multiple pulses are transmitted and the ranges of target scatterers remain the same but their gains change randomly during the observation interval.
Let Y (m) be the DFT of the received vector due to the mth pulse. Then where α (m) l denotes the gain of the lth scatterers due to the mth pulse and the noise vectors .
We now extend the WRELAX algorithm to this multiple look case. The extended WRELAX algorithm minimizes the following NLS criterion: where b(ω l ) is defined in (12).
Before we present the extended WRELAX algorithm, let us consider the following preparations. Let where α ..,L,i =l,m=1,...,M are assumed given. Then the cost function Minimizing Minimizing and ω l yieldŝ

Fixed delays and fixed gains
When both the delays and gains of the target scatterers remain the same during the multiple look interval, we can derive an ML estimator when the noise is assumed to be a zero-mean colored Gaussian noise with an unknown covariance matrix Q. Note that although we could continue to use the NLS approach for the current problem, we prefer to take the noise statistics into account since we will show below that doing so in this case introduces little difficulties for sufficiently large M. For the former problems, modeling the noise with an unknown covariance matrix Q makes the ML approach ill-defined due to too many unknowns (Li et al. b, 2002).
Let Y (m) be the DFT of the received data vector due to the mth pulse which can be written as where the noise vectors {E (m) } M m=1 are assumed to be zero-mean colored Gaussian random vectors with an unknown covariance matrix Q that are independent of each other. Let and Then The log-likelihood function of Y (m) is proportional to (within an additive constant): where det(·) denotes the determinant of a matrix and tr(·) denotes the trace of a matrix. Consider first the estimate of Q and the unstructured estimate of C = bα.I ti se a s yt os h o w that the estimateQ of Q isQ whereĈ may be obtained by minimizing the following cost function: Then To minimize det(G),wehaveĈ Then using theĈ in (39) to replace the C in (34) yieldŝ With these notations, the above C 5 can be rewritten as where we have used the fact that det(I + ab)=det(I + ba) if the dimensions of a and b permit. Hence minimizing C 5 is equivalent to minimizing which is again a highly nonlinear optimization problem.
We consider below using the relaxation based approach to minimize C 6 ({α l , ω l } L l=1 ).Let ,i =l are assumed given. Then minimizing C 6 becomes minimizing Consider first the case of complex-valued {α l } L l=1 . Minimizing C 7 (α l , ω l ) with respect to α l and ω l yields:α Consider next the case of real-valued {α l } L l=1 . Minimizing C 7 (α l , ω l ) with respect to α l and ω l yieldsα In the above derivations,Q −1 plays the role of whitening the noise. A good estimate of Q requires a large number of independent data vectors (i.e., M should be large enough as compared with N). When M is small, the noise covariance matrix estimated from (40) is singular or near singular. At least N data vectors are needed to guarantee that the matrixQ is non-singular with probability one.
Usually, the receiver noise e(t) in (1) can be modeled as a zero-mean stationary and ergodic Gaussian stochastic process. Let the covariance matrix corresponding to the sampled noise vector [e(0), e(T s ) ··· e((N − 1)T s )] T be Q t ,wherethesubscript"t" represents the covariance matrix of the noise in the time domain. Then Q t is a Hermitian and Toeplitz matrix. The frequency domain noise covariance matrix Q is related to Q t as follows where γ is the DFT matrix, It can be shown that, in general, Q is no longer a Toeplitz matrix. However, we can use the Toeplitz property of Q t to improve the estimation performance. First, we can obtainQ by using (40). Then the estimateQ t of Q t can be obtained by using (49), which iŝ Duetoafinitenumberofdatavectors,Q t is no longer a Toeplitz matrix. Although there are many ways to to modifyQ t to obtain a Toeplitz matrixQ (T) t , in this paper we use the following simple approach. Letq t (i, j) be the (i, j)th element ofQ t .D e fi n ê ThenQ UsingQ (T) instead ofQ in (45)-(48), wherê we obtain a new algorithm referred to as TWRELAX. The TWRELAX algorithm can greatly improve the estimation performance of WRELAX, especially when M is small as compared with N, as can be seen from the numerical examples below.
For notational convenience, we denote the new algorithms with and without exploiting the Toeplitz property of the noise covariance matrix by TWRELAX and WRELAX, respectively. Here, the time delay spacing between the two signals is τ 2 − τ 1 = 0.5τ e .We have used ǫ = 0.001 to test the convergence of TWRELAX and WRELAX. The colored noise is modeled as a first-order autoregressive (AR) process with coefficient a 1 =-0.85.
The MSEs of TWRELAX ("•") and WRELAX ("×") are compared with the corresponding CRBs (solid line) in Figures 3 and 4 as a function of the SNR and the normalized look number log 2 (M/N), respectively. We fix M = 4N in Figure 3 and SNR=−5d Bi nF i g u r e4 . F r o m Figure 3, it can be noted that, when M is sufficiently large, both TWRELAX and WRELAX can approach the CRBs for a wide range of SNRs and the former performs slightly better than the latter. However, when M is small as compared to N, TWRELAX outperforms WRELAX significantly, as can be seen from approach the CRBs. Note also that in Figure 4, due to the inversion of poorly estimated noise covariance matrices, some points of the MSE curves of WRELAX are beyond the scope of the axis limits.

Hybrid-WRELAX algorithm and EX tended Invariance Principle Based WRELAX (EXIP-WRELAX) algorithm
When bandpass real-valued probe signals are used, the correlation function between the received and the known transmitted signals oscillates near the carrier frequency of the transmitted signal. In this case, many existing time delay estimation algorithms perform poorly due to converging to local optimum points. Here, two efficient algorithms are proposed to deal with the above problem. First we assume that the signal amplitudes are complex-valued and use WRELAX to obtain the initial estimates of the delays and the amplitudes of the superimposed signals by minimizing a much smoother NLS cost function. Then the initial estimates are refined with two approaches. One approach (referred to as Hybrid-WRELAX) uses the last step of the WRELAX algorithm to minimize the true NLS cost function corresponding to the real-valued signal amplitudes. The other approach (referred to as EXIP-WRELAX) uses the extended invariance principle (EXIP). For Hybrid-WRELAX, the refinement step is iterative, while it is not for EXIP-WRELAX.
Real-valued signals are often bandpass signals that occur, for example, in underwater sonar and ultra wideband ground penetrating radar applications. Bandpass signals have highly oscillatory correlation functions, which makes the super resolution time delay estimation problem more difficult. The larger the center frequency of the pass band, the sharper the oscillation of the correlation function.
The same data model as that in the Section 2 is adopted. However, the transmitted signal s(t) and the received signal y(t) are real-valued, and the gains {α l } L l=1 and the noise e(t) are also real-valued.

Hybrid-WRELAX
Since both the transmitted signal s(t) and the received signal y(t) are real-valued, their Fourier transforms are conjugate symmetric, i.e., Y(−k)=Y * (k) and S(−k)=S * (k), k = 1, 2, ···, N/2 − 1, where (·) * denotes the complex conjugate, and Y(−N/2), Y(0), S(−N/2),andS(0) are real-valued. It can be readily shown that the cost function (6) is equivalent to where {W(k)=1} −1 k=−N/2+1 and W(−N/2)=W(0)=1/ √ 2. We assume that e(nT s ) is a real-valued zero-mean white Gaussian random process with variance σ 2 .Y e tE(k) will not be a circularly symmetric complex-valued zero-mean white Gaussian random process since E(−k)=E * (k), k = 1, 2, ··· , N/2 − 1. (The circularly symmetric assumption on the noise is widely used in the literature (Stoica & Moses, 1997).) The cost function C 8 ({α l , ω l } L l=1 ) in (56) with {α l } L l=1 being real-valued is referred to as the true cost function. Minimizing C 8 ({α l , ω l } L l=1 ) with respect to the unknown parameters is a highly nonlinear optimization problem. For narrowband transmitted signals, the cost function is highly oscillatory and have numerous closely spaced local minima, which makes it very difficult to find the global minimum. By assuming the real-valued amplitudes {α l } L l=1 to be complex-valued, a much smoother cost function can be obtained. This is equivalent to formulate the original time delay estimation problem in its complex analytic signal form. Since the analytic signal of the transmitted signal is lowpass, its autocorrelation function is no longer oscillatory. This is the conventional complex demodulation process and is widely used in practice. Although it is much easier to find the global minimum of the cost function corresponding to complex-valued amplitudes, the so-obtained estimates can be much less accurate than those obtained by minimizing the true cost function. The two cost functions share the same global minimum only when there is no noise. However, as suggested in Vaccaro et al., 1992), we can minimize the cost function associated with complex-valued amplitudes to obtain the initial conditions needed to minimize the true cost function. Below, we present a relaxation based global minimizer of the NLS criterion based on this idea. The algorithm is referred to as the Hybrid-WRELAX algorithm. It simply requires a sequence of weighted Fourier transforms.
Before we present our approach, let us consider the following preparations. Let Then (56) becomes where · denotes the Euclidean norm. Minimizing C 9 (α l , ω l ) with respect to the real-valued α l yields the estimateα l of α lα where (·) H denotes the conjugate transpose, Re(Z) represents the real part of Z,a n d · F denotes the Frobenius norm (Stewart, 1973). (More specifically, S F = ∑ 0 n=−N/2 |W(n)S(n)| 2 .) Then the estimateω l of ω l is obtained as follows: w h e r ew eh a v eu s e dt h ef a c tt h a tb H (ω l )b(ω l )=∑ 0 n=−N/2 |W(n)S(n)| 2 and hence is independent of ω l .H e n c ê ω l is obtained as the location of the dominant peak of Re 2 ā H (ω l )(S * Ȳ l ) . With the estimate of ω l at hand,α l is easily computed from the corresponding complex height by usingω l to replace ω l in (64).
Similarly, minimizing C 9 (α l , ω l ) with respect to ω l and the complex-valued α l , respectively, yields the estimatesω l of ω l andα l of α l , whereω l can also be found via FFT with the weighted data vector.
With the above simple preparations, we now present the Hybrid-WRELAX algorithm.
Step 2: Refine the estimates obtained in Step 1 with the last step of the WRELAX algorithm (i.e., the last substep of Step 1 above) by using Equations (64) and (65) derived for the real-valued {α l } L l=1 and using {ω l } L l=1 and the real parts of {α l } L l=1 obtained in Step 1 as initial conditions. Iteratively update {ω l ,α l }, l = 1, 2, ··· , L, until "practical convergence".
Note that WRELAX can be used directly for signals with real-valued amplitudes. For this case, the approach would consist of the substeps of Step 1 above except that (64) and (65) will be used instead of (66) and (67), respectively. We will use a numerical example in Section 5 to show the problem encountered by the direct use of WRELAX when the cost function is highly oscillatory.

EXIP-WRELAX
The Invariance Principle (IP) of ML estimators is well known in the estimation theory (Zehna, 1966). The invariance principle gives a simple answer to the relationship between the minimizers of a given cost function parameterized in two different ways in some special cases. By appropriately reparameterizing the original cost function and enlarging the supporting domain of the parameter space, less accurate estimates can be obtained from this simple data model. These estimates may be refined to asymptotically achieve the performance available using the original data model. This is the basic idea behind the Extended Invariance Principle (EXIP) proposed in  and (Söderström & Stoica b, 1989) for the purpose of achieving some computational advantages. In this section, we present an EXIP based algorithm, referred to as the EXIP-WRELAX algorithm, that avoids dealing with the highly oscillatory true cost function entirely.
By using (58) and (62), the cost function (56) with {α l } L l=1 being real-valued can be written in the following vector form where By replacing the real-valued amplitudes {α l } L l=1 with the complex-valued amplitudes {α l } L l=1 (notations introduced for the sake of clarity) in (68), we obtain the following cost function: with Im(Z) denotes the imaginary part of Z,and Let where with I and 0 denote the L × L identity matrix and the L × L matrix with zero elements, respectively. Using the EXIP principle Söderström & Stoica b, 1989), we can obtain a new estimateη fromη by solving the following weighted least squares where It has been shown in Söderström & Stoica b, 1989) thatη is asymptotically (for large N or high SNR) statistically equivalent toη. The weighting matrix W EXIP is simply the Fisher In f ormation Matrix (possibly scaled by a constant) for the complex-valued {α l } L l=1 withη replaced by its estimateη . It can be easily shown that The EXIP-WRELAX algorithm is composed of two steps. The first step is the same as Step 1 of the Hybrid-WRELAX algorithm and the second step is to refine the initial conditions obtained in Step 1 by using (81). Compared to the Hybrid-WRELAX algorithm, the second step of the EXIP-WRELAX algorithm is non-iterative and avoids dealing with the highly oscillatory true NLS cost function entirely.
Our numerical examples show that at low SNR, the former tends to outperform the latter. The MSEs of the WRELAX ("+") for assuming {α l } L l=1 being complex-valued, Hybrid-WRELAX ("•"), and EXIP-WRELAX ("×") are compared with the CRBs obtained by assuming {α l } L l=1 being complex-valued (dashed line) and real-valued (solid line) in Figure  5. Note that both Hybrid-WRELAX and EXIP-WRELAX achieve the corresponding CRB. Note also that the threshold effect is obvious in Figure 5, where the MSEs deviate away from the CRBs at low SNR. Although the WRELAX for assuming {α l } L l=1 being complex-valued also attains its corresponding CRB (dashed line) at high SNR, this wrong CRB can be larger than the true CRB by approximately 30 dB. (Note that the former CRB is expected to be worse than the latter CRB due to the parsimony principle In this example, Hybrid-WRELAX outperforms EXIP-WRELAX at low SNR.

Method Of Direction Estimation based WRELAX (MODE-WRELAX) algorithm
WRELAX was extended in to deal with the real-valued signals with highly oscillatory correlation functions (Hybrid-WRELAX and EXIP-WRELAX). The resolution of WRELAX are much higher than that of the conventional matched filter approach. However, when the signals are very closely spaced in arrival times, the convergence speed of WRELAX decreases rapidly. Here, we study how MODE can be used with our efficient WRELAX algorithm for super resolution time delay estimation. The new algorithm is referred to as MODE-WRELAX. Although MODE can provide very poor amplitude estimates and WRELAX has the slow convergence problem, MODE-WRELAX outperforms both MODE and WRELAX. MODE-WRELAX can be used for both complex-and real-valued signals (including those with highly oscillatory correlation functions).
The same data model as that in the Section 2 is adopted. Here, the transmitted signal s(t) represents an arbitrary known transmitted signal. We assume that s(t) , y(t), e(t) and {α l } L l=1 are either all complex-valued or all real valued, which will be dealt with in the following respectively.

MODE-WRELAX for complex-valued signals
Assume that (·) T denote the transpose and let where a(ω l ) is same as that in (10).
Then the data model (4) can be written in the following vector form: where Y, S and E are same as that in (7)-(9), α is same as that in (70).
When S is an identity matrix, then the above time delay estimation issue becomes a sinusoidal parameter estimation problem and MODE is an asymptotically statistically efficient estimator of {ω l } L l=1 for complex-valued signals Stoica & Sharman b, 1990). The MODE algorithm Stoica & Sharman b, 1990) can be easily extended to the data model in (83) where S is an arbitrary diagonal matrix as follows. The MODE estimates {ω l } L l=1 of {ω l } L l=1 can be obtained by minimizing the following cost function and To avoid the search over the parameter space, C 10 ({ω l } L l=1 ) can also be reparametrized in terms of another parameter vectorb = b 0b1 ···b L T ,w h er e{b l } L l=0 are the coefficients of the following polynomial: Since the polynomialb(z) in (87) has all of its zeros on the unit circle, its coefficients {b l } satisfy the conjugate symmetry constraint : where (·) * denotes the complex conjugate. Let Assume that the diagonal elements of S are nonzero (see Remark 1 for more discussions). Let It can be readily verified that B H A = 0 and henceB (84) is equivalent to minimizing Note thatB HB in (91) can be replaced by a consistent estimate without affecting the asymptotically statistical efficiency of the minimizer of (91). Henceb can be obtained computationally efficiently as follows: whereB 0 is the initial estimate of B obtained by replacingb withb (0) in (89). The initial valuễ b (0) is obtained by settingB HB in (91) to I: are obtained, the amplitudes α are estimated by applying the linear least-squares approach to whereÂ is formed by replacing {ω l } L l=1 with {ω l } L l=1 in (82). Remark 1: MODE cannot be implemented efficiently to avoid the search over the parameter space when S(k)=0f o rs o m ek. The most commonly used complex analytic signal s(t) is low-pass. For this case, we can select a contiguous segment of Y satisfying |S(k)| > 0, K 1 ≤ k ≤ K 2 , and preferrably with |S(k)| above a certain threshold to avoid numerical problems. We can then apply MODE to the segment {Y(k)} K 2 k=K 1 to estimate {ω l } L l=1 . Remark 2: The amplitude estimates given above can be very poor when the SNR is not sufficiently high. This is because some of the MODE estimates {ω l } L l=1 can be so closely spaced thatÂ in (94) is seriously ill-conditioned. We use a simple spacing adjustment scheme to avoid this problem. After obtaining the MODE estimates {ω l } L l=1 of {ω l } L l=1 ,wefirstsort them in the ascending order and then check the spacing between two adjacent estimates. If the distance between any two estimates, sayω 1 andω 2 (ω 1 ≤ω 2 ), is smaller than a predefined threshold, say ω t , we adjust the estimates by replacingω 1 withω 1 − 0.5 ω t andω 2 witĥ ω 2 + 0.5 ω t . The amplitudes are then estimated using the adjusted estimates of {ω l } L l=1 .Th is spacing adjustment step is ad hoc but can be used to provide good initial delay and amplitude estimates to replace the first L − 1 steps of WRELAX.
The MODE estimates {ω l } L l=1 of {ω l } L l=1 and {α l } L l=1 of {α l } L l=1 , which may not be optimal, especially for real-valued signals, can be refined by using the last step of the WRELAX algorithm.
WRELAX is a relaxation-based minimizer of the following nonlinear least-squares (NLS) criterion: When the signals are not spaced very closely, WRELAX usually converges in a few steps. However, when the signals are very closely spaced, the convergence speed of WRELAX is very slow. Yet by using the above MODE algorithm to obtain the initial conditions and then using the last step of WRELAX to refine them, super resolution time delay estimation can be achieved with a fast convergence speed.
Before we present the MODE-WRELAX algorithm, let us consider the following preparations. Let where {α i ,ω i } L i=1,i =l are assumed to be given . Then (95) becomes Minimizing C 15 (α l , ω l ) with respect to ω l and the complex-valued α l yieldŝ With the above preparations, we now present the steps of the MODE-WRELAX algorithm for complex-valued signals.
Step (1): Select a contiguous segment of data vector Y (for MODE use only) so that |S(k)| > 0, K 1 ≤ k ≤ K 2 . Apply MODE to the segment to obtain {ω l } L l=1 .A d j u s t{ω l } L l=1 so that the minimum spacing of {ω l } L l=1 is at least ω t . Obtain the estimates {α l } L l=1 of {α l } L l=1 by using (94).
Similarly, we can use MODE as an initialization method for the EM time delay estimation algorithm (Moon, 1996), which is referred to as MODE-EM. However, we have found through numerical simulations that the convergence speed of MODE-EM is slower than that of MODE-WRELAX.

MODE-WRELAX for real-valued signals
For bandpass real-valued signals, the cost function in (95) is a highly oscillatory cost function and is very difficult to find its global minimum. Although MODE is derived for complex-valued signals, we can apply it by assuming the real-valued amplitudes {α l } L l=1 to be complex-valued. These initial estimates are then refined by the WRELAX algorithm. Since the attraction domain of the cost function (95) is extremely small, a very good initial condition is required to achieve the global convergence of any minimizer of (95). The MODE estimates are first refined by WRELAX by assuming {α l } L l=1 to be complex-valued since the attraction domain of (95) becomes much larger when assuming the real-valued {α l } L l=1 to be complex-valued . The so-obtained estimates are refined again by WRELAX by using the fact that {α l } L l=1 are real-valued. With the above preparations, we now present the steps of the MODE-WRELAX algorithm for real-valued signals.
Step (1): Select a contiguous segment of data vector Y so that |S(k)| > 0, K 1 ≤ k ≤ K 2 . By assuming the real-valued {α l } L l=1 to be complex-valued, obtain the estimates {ω l } L l=1 and {α l } L l=1 in the same way as Step (1) of the MODE-WRELAX algorithm for complex-valued signals.
Step (2): Refine the estimates obtained in Step (1) above by using the last step of WRELAX by assuming complex-valued signals. Take the real parts of the so-obtained amplitude estimates as the amplitude estimates {α l } L l=1 of {α l } L l=1 .
Step (3): Refine the estimates obtained in Step (2) above by using the last step of WRELAX and the fact that the signals are real-valued. Now we consider an example where τ 2 − τ 1 = 0.2τ e . The MSEs of MODE ("•"), WRELAX ("×"), and MODE-WRELAX (" * ") are compared with the corresponding CRBs (solid line) in Figure 6. Note that due to the highly oscillatory cost functions and very closely spaced signals, WRELAX converges to some local minimum instead of the global one, which yields very poor estimates. Since the MODE amplitude estimates are obtained without spacing adjustment, they are so poor at low SNR that some of their MSEs are above the axis limit due to the inversion of ill-conditioned matrices corresponding to very closely spaced delay estimates. Although the MSEs of the MODE estimates are close to the CRBs corresponding to the complex-valued amplitudes when the SNR is high, the wrong CRBs (not shown to avoid too many lines in the figure) can be larger than the true CRBs, which correspond to the real-valued amplitudes, by approximately 30 dB.

Pavement profiling via Ground Penetrating Radar (GPR)
The detection and classification of roadway subsurface anomalies are very important for the design and quality evaluation of highways. Ultra wideband ground penetrating radar emits nonsinusoidal impulses with extremely large bandwidth (several GHz) and is very suitable for this application because of its high range resolution (on the order of several centimeters). The returned echoes of the ultra wideband ground penetrating radar are superimposed real-valued signals reflected from the boundaries of different media (layers, voids, etc.), which can be described by (1). Both the delays and gains are very useful for the detection and classification of roadway subsurface anomalies. The delays can be used to determine the layer thickness or anomaly location and the gains can be used to classify the type of media because the gains are related to the reflection coefficient at the boundary between two media with different dielectric constants. Once we get the estimates of the media dielectric constants, we can judge the type of the media.
WRELAX has been applied to layer stripping, based on a single-scan basis using experimental data and has been proven to be very useful. However, the straightforward extension of WRELAX from one scan to multi-scans is not feasible since, under moving conditions, the antenna can change its height relative to the road surface. As the height above the road decreases (or increases), the amplitude and the time delay of the received signal changes as well, which can lead to incorrect estimates of layer thickness and permittivity. A motion compensation and parameter estimation algorithm for GPR is presented for simultaneous motion compensation and time delay estimation. More details can be found in our papers Su & Wu, 2000;.

Acoustic proximity ranging system
Proximity ranging is very important in a wide range of remote sensing applications, including level detection, robot manipulation, process control, nondestructive testing, and cavity thickness monitoring. Various types of sensors based on different physical principles such as capacitive or inductive proximity sensors, laser displacement sensors, and ultrasonic sensors, can be used to perform proximity ranging. Among these sensors, ultrasonic sensors have many important advantages over the others. However, due to the presence of secondary echoes, using ultrasonic sensors for very close proximity ranging is very difficult. When ultrasonic sensors face a sound-hard reflective surface, the reflected sound wave can bounce back and forth several times between the sensors and the reflection surface before decaying to zero, which results in unwanted strong and overlapping secondary echoes in the received signal. The time delays for the secondary echoes are approximately integer multiples of the time delay of the first echo. The matched filter based methods cannot resolve two echoes with a time spacing less than the reciprocal of the signal bandwidth. Hence, for most of the very short distance measurement scenarios, the matched filter based algorithms tend to fail or suffer from severe performance degradations due to their poor resolutions. WRELAX and its extended version as a super-resolution time delay estimation approaches are developed for general purposes. They do not exploit the a priori information of the integer multiple time delays and the nonnegative amplitude due to the acoustically hard reflections. More details can be found in our papers Li et al. b, 2003).

Conclusion
A family of relaxation-and FFT-based efficient time-delay estimation algorithms were presented for different scenarios in this chapter. By avoiding the computationally demanding multidimensional search over the parameter space, the proposed algorithms minimize the NLS criterion at a much lower implementation cost. They are more efficient and systematic than existing algorithms. Some practical applications utilizing the proposed algorithms are also presented. Theoretical analysis and simulations demonstrate the efficiency of the proposed new algorithms. The field of signal processing has seen explosive growth during the past decades; almost all textbooks on signal processing have a section devoted to the Fourier transform theory. For this reason, this book focuses on the Fourier transform applications in signal processing techniques. The book chapters are related to DFT, FFT, OFDM, estimation techniques and the image processing techqniques. It is hoped that this book will provide the background, references and the incentive to encourage further research and results in this area as well as provide tools for practical applications. It provides an applications-oriented to signal processing written primarily for electrical engineers, communication engineers, signal processing engineers, mathematicians and graduate students will also find it useful as a reference for their research activities.