Least-squares based iterative multipath super-resolution technique
aa r X i v : . [ c s . I T ] M a r Least-Squares Based Iterative MultipathSuper-Resolution Technique
Wooseok Nam,
Member, IEEE,
Seung-Hyun Kong*,
Member, IEEE
Dept. of Aerospace Engineering, KAIST,Daejeon, Republic of KoreaE-mail: [email protected], [email protected]
Abstract
In this paper, we study the problem of multipath channel estimation for direct sequence spread spectrum signals.To resolve multipath components arriving within a short interval, we propose a new algorithm called the least-squares based iterative multipath super-resolution (LIMS). Compared to conventional super-resolution techniques,such as the multiple signal classification (MUSIC) and the estimation of signal parameters via rotation invariancetechniques (ESPRIT), our algorithm has several appealing features. In particular, even in critical situations wherethe conventional super-resolution techniques are not very powerful due to limited data or the correlation betweenpath coefficients, the LIMS algorithm can produce successful results. In addition, due to its iterative nature, theLIMS algorithm is suitable for recursive multipath tracking, whereas the conventional super-resolution techniquesmay not be. Through numerical simulations, we show that the LIMS algorithm can resolve the first arrival pathamong closely arriving independently faded multipaths with a much lower mean square error than can conventionalearly-late discriminator based techniques.
Index Terms
Multipath channel, super-resolution, pseudo-noise code, early-late discriminator, least-squares, maximum like-lihood
I. I
NTRODUCTION
A pseudo-noise (PN) code sequence, which is a deterministic but noise-like sequence generated by alinear feedback shift register, is widely used in the areas of radar and sonar signal processing, ranging andpositioning, and digital communication, since it is convenient to estimate the code phase of the receivedsignal using its peaky auto-correlation function (ACF). When the PN code sequence is launched to achannel and received by a remote terminal, the cross-correlation of the received signal and a receiverreplica of the PN code sequence yields a sharp peak at the lag of the propagation delay caused by thechannel. This delay is related to the physical separation between the transmitter and the receiver, andthus gives important information for source localization or target detection. However, in many practicalsituations, the propagation medium suffers from multiple echoes with short delays due to scattering objectsaround the receiver. In such a case, we are frequently interested in estimating the delays of those echoesand, especially, the first arrival path, since it is the most likely to be the line-of-sight path, which containsthe information of the true distance between the transmitter and the receiver.In the multipath channel, the conventional PN code phase estimation method using the ACF, whichfinds the peak of the correlator output, would probably give an incorrect estimate of the first arrivalpath delay unless the first arrival path is the most dominant path. This is because the ACF for the firstarrival path is buried under the sum of other overlapping ACFs of echos, and thus the lag of the largestmagnitude correlator output is often different from the first arrival path delay. We can work around thisproblem to some extent by finding multiple local extrema of the correlator output and selecting the onewith the smallest lag as the estimate of the first arrival path [1]. However, the time resolution of this *Corresponding author. method is limited to the chip duration of the PN code sequence, which is not sufficient in some narrowband applications.In order to improve resolution of the conventional ACF based code phase estimation, super-resolution techniques can be taken into account. The super-resolution techniques, such as the multiple signal classi-fication (MUSIC) [2] and the estimation of signal parameters via rotation invariance techniques (ESPRIT)[3], were originally proposed for direction finding of multiple targets with passive sensor arrays, andfirst applied to the channel multipath resolution in [4]. In [4] and [5], the super-resolution techniqueswere applied directly to the received signal samples without cross-correlation. In particular, in [5], it wasshown that the super-resolution techniques can be applied to received signals sampled at a considerablylower rate than the chip rate of the PN code sequence when the signals have a finite rate of innovation(FRI). On the other hand, in [6], [7], super-resolution techniques for the correlator output signal wereintroduced. These techniques can significantly improve the resolution of multipath delays, compared toconventional techniques, but have some major drawbacks that render them impractical. First, since thesuper-resolution techniques are based on the subspace decomposition of the signal correlation matrix,the signal components contributed by different channel paths should occupy distinct dimensions of thesubspace in order to be distinguishable. For this, each distinct path should be far enough apart fromthe others, and the attenuation of each path should be varying randomly and not be perfectly correlatedwith the others during an observation interval. These conditions are not satisfied when the paths are quitenarrowly spaced, the channel is slowly varying, or the observation interval is not long enough. In addition,it is not so reasonable to assume that only the channel attenuations are randomly varying while the pathdelays are fixed. As a result, additional techniques, such as frequency smoothing [6], [7], are required atthe transmitter to make the assumption valid. As a second drawback, most super-resolution techniques arebatch algorithms that compute a set of estimates from a large number of stationary observations. Althoughthe technique for FRI signals [5] can work with a very small number of observations, it has very low noiseimmunity and also needs a large number of observations to achieve robustness against noise. The radioranging application, however, often requires tracking of the signal as well. That is, a new observation,which can be non-stationary, is periodically given, and a new set of estimates should be immediatelycomputed from the new observation and the previous set of estimates in a recursive manner [8], [9]. Thus,the super-resolution techniques are not adequate for applications where tracking is a concern.In this paper, we consider a least-squares (LS) approach to resolve the short-delay multipath componentswith a correlator output signal. When the noise is Gaussian distributed, our approach includes the maximumlikelihood (ML) approach, which is optimal in performance. The LS and ML approaches for the channelmultipath resolution were presented directly and indirectly in earlier studies [10], [11], [12], but they havenot gathered much interest since the highly nonlinear nature in the parameters to be estimated rendersthe use of these approaches computationally intensive. In [10] and [11], faster algorithms for solving anonlinear LS problem were derived by exploiting the linear predictability of sinusoidal signals, but theyare not applicable to our case. In our work, we propose an iterative algorithm to solve the LS multipathresolution problem, which we will call the least-squares based iterative multipath super-resolution (LIMS)algorithm. Thanks to its iterative and recursive mode of operation, it is appropriate for multipath trackingwith sequential observations of the received signal. Moreover, we show that the tracking of a PN codephase using the early-late (EL) discriminator [8], [9], [13] is a special case of our algorithm. Throughnumerical simulations, we show that our scheme asymptotically achieves the optimal performance givenby the Cramer-Rao bound (CRB) at high carrier-to-noise density ratio. In addition, we show that ouralgorithm outperforms the EL discriminator based technique in the presence of severe multipaths andnoise.This paper is organized as follows. In Section II, we introduce the signal model and the LS formulationfor the multipath resolution. In Section III, we present the LIMS algorithm as a technique to solve the LSmultipath resolution problem, and discuss the application of the LIMS algorithm for multipath tracking.The complexity of the LIMS algorithm and its relationship with the conventional EL discriminator basedtechnique are also discussed in Section III. Results of numerical simulation in various channel models and a comparison with the conventional method are given in Section IV, and Section V concludes thepaper.We will use the following notation throughout this paper. Vectors or matrices are denoted by boldfacesymbols. The p -th element of a vector and the ( p, q ) -th element of a matrix are denoted by [ · ] p and [ · ] p,q ,respectively. The superscripts T and H denote the transpose and conjugate transpose, respectively. TheEuclidean norm of a vector is denoted by k · k , and the infinity norm of a vector, i.e., the largest absolutevalue of elements, is denoted by k · k ∞ . The fields of real and complex numbers are indicated by R and C , respectively, and the N × N identity matrix is denoted by I N . The statistical expectation is denotedby E {·} and the real part of a complex value is denoted by ℜ{·} . Finally, O ( · ) is the big O notation.II. S YSTEM M ODEL
Consider a multipath channel between a transmitter and a receiver with an impulse response h ( t ) = K − X k =0 γ k δ ( t − τ k ) , ≤ t < T i , (1)where K is the number of paths, { γ , . . . , γ K − } are the unknown complex channel coefficients for eachpath, { τ , . . . , τ K − } are the unknown arrival time delays for each path, δ ( t ) is the Dirac delta function,and T i ( ≫ max { τ , . . . , τ K − } ) is a time interval of interest. We assume that the delays { τ , . . . , τ K − } are distinct and τ < τ < · · · < τ K − , without loss of generality. Thus, the delay of the first arrival pathis τ . The transmitter launches a continuous time PN code sequence s ( t ) with a chip duration T c to thechannel (1), and the baseband received signal at the receiver is written as η ( t ) = h ( t ) ∗ s ( t ) + v ( t )= K − X k =0 γ k s ( t − τ k ) + v ( t ) , ≤ t < T i , (2)where v ( t ) is a circularly symmetric, zero-mean, complex Gaussian noise process with auto-correlation N δ ( ζ ) . For the first process at the receiver, the cross-correlation of η ( t ) with s ( t ) is taken over the timeinterval T i , and the correlator output is uniformly sampled at N lags { , T s , . . . , ( N − T s } to yield y n = K − X k =0 γ k R ( nT s − τ k ) + w n , n = 0 , . . . , N − , (3)where R ( ζ ) is the auto-correlation function (ACF) of the PN code sequence given by R ( ζ ) = 1 T i Z T i s ( t ) s ∗ ( t − ζ ) dt , (4)and w n is the filtered noise process w ( ζ ) sampled at ζ = nT s , where w ( ζ ) = 1 T i Z T i v ( t ) s ∗ ( t − ζ ) dt . (5)For a binary phase shift keying (BPSK) modulated PN code sequence with a sufficiently large integrationtime T i ≫ T c , the ACF (4) is ideally given by R ( ζ ) = ζT c + 1 , − T c < ζ ≤ , − ζT c + 1 , < ζ ≤ T c , , otherwise. (6) For notational simplicity in the subsequent expressions, we denote y = [ y , . . . , y N − ] T , c CH = [ γ , . . . , γ K − ] T , t CH = [ τ , . . . , τ K − ] T , and w = [ w , . . . , w N − ] T . Then, from (5), the covariance matrix of the noisevector, E (cid:8) ww H (cid:9) , is given by N T i C , where C has elements [ C ] p,q = T i N E { w p − w ∗ q − } = R (( p − q ) T s ) , ≤ p, q ≤ N . (7)From the above signal model, the LS estimation of the channel coefficients and the path delays isformulated as (cid:8) ˆ c , ˆ t (cid:9) = arg min c , t k Gy − GA ( t ) c k (8a) = arg min c , t k Gy − G ( B ( c , t ) t + b ( c , t )) k , (8b)where c = [ c , . . . , c M − ] T ∈ C M × and t = [ t , . . . , t M − ] T ∈ R M × are the vectors of channelcoefficients and path delays to be estimated, respectively, and M is the assumed number of multipathsfor the estimation. Since the true number of multipaths, K , is generally unknown, we assume that M canbe different from K . Note that there are two equivalent LS formulations (8a) and (8b), and the relatedparameters are defined as [ A ( t )] n +1 ,m +1 = R ( nT s − t m ) , (9a) [ B ( c , t )] n +1 ,m +1 = − c m T c , t m − T c T s < n ≤ t m T s , c m T c , t m T s < n ≤ t m + T c T s , , otherwise, (9b) b ( c , t ) = M − X m =0 b ( c m , t m ) , (9c) [ b ( c m , t m )] n +1 = c m · (cid:16) nT s T c (cid:17) , t m − T c T s < n ≤ t m T s , c m · (cid:16) − nT s T c (cid:17) , t m T s < n ≤ t m + T c T s , , otherwise, (9d) n = 0 , . . . , N − , m = 0 , . . . , M − ,and G ∈ C N × N is a proper nonsingular weighting matrix. For the least-squares estimation (8) to befeasible, we assume that N ≥ M , and A ( t ) and B ( c , t ) have full column rank. The first assumptioncan be reasonable by taking a sufficient number of samples, and so is the second assumption when thedelays { t , . . . , t M − } are distinct and sufficiently spaced by a small fraction of T c [14, 2.4.7]. Remark 1:
The signal model (3) is restricted to using the triangular shaped ACF (6) for simplicity, butit can be extended to general ACFs through a proper linearization around the sampling points. That is,for a general ACF R ( ζ ) , (9b) and (9c) can be replaced by [ B ( c , t )] n +1 ,m +1 = c m · dR ( nT s − t m ) dt m , (10) [ b ( c , t )] n +1 = M − X m =0 c m · (cid:18) R ( nT s − t m ) − dR ( nT s − t m ) dt m t m (cid:19) . (11)In addition, though we assume the uniform sampling with a sampling period T s in (3), all the argumentsin this paper can be immediately extended to non-uniform sampling. Remark 2:
If we let G = C − , the weighting matrix G whitens the noise and thus the LS formulation(8) becomes the maximum likelihood (ML) formulation. In Section IV, we investigate the performancewith and without noise whitening by letting G = C − and G = I N , respectively. III. S
OLVING THE L EAST -S QUARES M ULTIPATH R ESOLUTION P ROBLEM
From the LS formulation (8), we can observe that it is linear in c given t but nonlinear in t due to thedependency of B ( c , t ) and b ( c , t ) on t . Thus, to find the solution to the LS problem, a computationallyexpensive multidimensional search seems to be inevitable. However, there exists a useful workaroundto avoid an exhaustive search: the gradient descent algorithm [15]. The gradient descent algorithm canefficiently find the local minimum of the LS problem in an iterative manner, by stepping from the currentguess of the LS solution to the direction that decreases the LS weight in (8) at every iteration. In thefollowing subsection, we propose a gradient descent algorithm for the LS multipath resolution problem,which will be denoted by the least-squares based iterative multipath super-resolution (LIMS) algorithm. A. Least-squares based iterative multipath super-resolution algorithm
For notational simplicity, let ˜ y , Gy , ˜ A ( t ) = GA ( t ) , ˜ B ( c , t ) = GB ( c , t ) , and ˜ b ( c , t ) = Gb ( c , t ) in(8). The operation flow of the LIMS algorithm is described as follows. First, at the l -th iteration ( l is apositive integer), the intermediate guesses of the channel coefficient and the path delay vectors, c l − and t l − , are given from the previous ( l − -th iteration. Given t l − , a new estimate of the channel coefficientvector, c l , is computed as c l = (cid:16) ˜ A H ( t l − ) ˜ A ( t l − ) (cid:17) − ˜ A H ( t l − )˜ y . (12)Note that the matrix inversion in (12) exists with assumptions that G is nonsingular and A ( t ) has fullcolumn rank. The new estimate c l (12) minimizes the LS weight in (8a) for a given t l − , but it requires amatrix inversion that can sometimes be troublesome in terms of complexity or numerical stability issues.Thus, as an alternative, a refinement by the gradient descent algorithm can be used as c l = c l − + α · ˜ A H ( t l − ) (cid:16) ˜ y − ˜ A ( t l − ) c l − (cid:17) , (13)where α is a properly chosen step size. The refinement c l (13) does not have complexity or stabilityproblems like (12), but may converge more slowly than (12).For the second step, a refinement t l is computed from t l − and c l using the gradient descent algorithm,for which it is required to find the gradient of the LS weight in (8b) with respect to t . However, sincethe LS weight is nonlinear and non-analytic in t , it is not easy to compute the exact gradient. Thus, wecompute an approximate gradient on the basis of the following assumptions: For a small perturbationvector e , B ( c , t + e ) ≃ B ( c , t ) , (14a) b ( c , t + e ) ≃ b ( c , t ) . (14b)From (9b), (9c), and (9d), the assumptions in (14) are quite reasonable for a magnitude of e sufficientlysmall in the sense that k e k ∞ ≪ T s . Under these assumptions, it can be assumed that ˜ B ( c l , t l − ) and ˜ b ( c l , t l − ) are constant in a small region around t l − . As a result, the gradient descent for t reduces to t l = t l − + β · ℜ n ˜ B H ( c l , t l − ) (cid:16) ˜ y − (cid:16) ˜ B ( c l , t l − ) t l − + ˜ b ( c l , t l − ) (cid:17)(cid:17)o = t l − + β · ℜ n ˜ B H ( c l , t l − ) (cid:16) ˜ y − ˜ A ( t l − ) c l (cid:17)o , (15)where β is a proper step size. The step size β should be small enough to satisfy the assumptions in (14),but a systematic criterion for β that guarantees the convergence is not easy to find.Upon getting c l using (12) or (13), and t l using (15), the iteration index l is increased by one and thesame computations are repeated for c l +1 and t l +1 . In this way, the iteration continues until a specific stopcriterion is met. There are a number of stop criteria for this kind of iterative algorithm, e.g., reachinga fixed number of iterations, thresholding the magnitudes of changes of variables, etc., and we omit adetailed discussion of these criteria. The LIMS algorithm is summarized in Table I. B. Complexity issues
The major burden of the LIMS algorithm is the matrix inversion in (12), which is O ( M ) of computationper iteration. The true number of multipaths, K , could be large in practice and, hence, the assumed numberof multipaths, M , should be accordingly large. However, assuming that only a few of the multipaths aredominant, we can keep M small; say, M < . Therefore, when N ≫ M , the matrix inversion maynot be too burdensome compared to the overall complexity. In addition, as noted in Section III-A, thematrix inversion can be avoided by using (13) instead, which is O ( N M ) of computation per iteration. Thetotal number of iterations is another important factor for the complexity. However, as will be discussed inSection III-C, the number of iterations can be considerably smaller when the LIMS algorithm is combinedwith a tracking function such as the delay locked loop. C. Multipath tracking
In many practical realizations of radio ranging systems, a new observation of the received signal isobtained periodically. Successive observations may include different multipath channels when the channelis time varying. Under the assumption that the channels for the successive observations are not statisticallyindependent, the time varying path delays can be efficiently tracked. One such tracking system is the delaylocked loop (DLL) [8], [9]. In the conventional DLL used for PN code phase tracking, the early-late (EL)discriminator is widely used to generate an error signal from the previous estimate of the path delay andthe new observation of the received signal (correlator output). The error signal is then filtered by the loopfilter in the tracking loop and a refined delay estimate is computed. Now, instead of the EL discriminator,we consider using the LIMS algorithm in the tracking loop.Suppose that the ν -th observation of the received signal samples, y ( ν ) , is given. Using the LIMSalgorithm, estimates of the channel coefficient vector ˆ c ( ν ) and the path delay vector ˆ t ( ν ) are computed.Then, when the ( ν + 1) -th observation y ( ν +1) is given, the LIMS algorithm starts a new round of iterationswith ˆ c ( ν ) and ˆ t ( ν ) as initial values. As long as ˆ c ( ν ) and ˆ t ( ν ) are good estimates of the channel coefficientsand the path delays, and the channels for the ν -th and the ( ν + 1) -th observations do not differ much, thenew estimates ˆ c ( ν +1) and ˆ t ( ν +1) might be obtained in a small number of iterations.When the channel is not time-varying but corrupted by different realizations of noise at differentobservations, applying the LIMS algorithm with a small number of iterations for each observation canbe seen as the application of the stochastic gradient method for the channel estimation [15]. Thus, asthe number of observations increases, the accuracy of channel estimation can improve until it reaches thesteady state. TABLE IS
UMMARY OF THE
LIMS
ALGORITHM . Input y = [ y , . . . , y N − ] T Known parameters
Weighting matrix: G ∈ C N × N Step sizes: α and β Initial values c and t Computation : l = 1 , , . . . c l = (cid:16) ˜ A H ( t l − ) ˜ A ( t l − ) (cid:17) − ˜ A H ( t l − )˜ y (cid:16) or c l = c l − + α · ˜ A H ( t l − ) (cid:16) ˜ y − ˜ A ( t l − ) c l − (cid:17)(cid:17) t l = t l − + β · ℜ n ˜ B H ( c l , t l − ) (cid:16) ˜ y − ˜ A ( t l − ) c l (cid:17)o Definitions A ( t ) , B ( c , t ) : given in (9) ˜ y , Gy , ˜ A ( t ) = GA ( t ) , ˜ B ( c , t ) = GB ( c , t ) Fig. 1. Sampling c R ( ζ − t ) at P points. D. Relationship with the conventional early-late discriminator
The EL discriminator computes the difference between the powers (magnitude squares) of two sampleswith different lags at the correlator output. The two samples are called the early and the late samples andare separated by a constant interval. When the peak of an ACF is located in between the two samples,the difference between the two powers serves as a feedback signal that shifts the sampling points so thatthe peak of the ACF is at the midpoint of the sampling points. Intuitively, this operation is establishedon the assumptions that the channel has a single path and that the correlator output maintains the exacttriangular shape of the ACF R ( ζ ) (6). However, these assumptions are not valid when the channel hasmultiple paths and the received signal is corrupted by noise.Basically, the EL discriminator can be thought of as a special case of the LIMS algorithm with thesingle path assumption, i.e., M = 1 , and two samples, i.e., N = 2 . To illustrate this relationship, let usfirst consider P sampling points at ζ = (2 p +1) T s , p = − P, . . . , P − , where P = j T c T s + k . Then, at thesampling points, the sampled values of the correlator output are denoted by y E ( − p ) for p < and y L ( p +1) for p ≥ . Likewise, the sampled values of the scaled ACF, c R ( ζ − t ) , are denoted by r E ( − p ) for p < and r L ( p +1) for p ≥ , as shown in Fig. 1. Now, we temporarily let T s = T c , and thus P = 1 and thenumber of total samples is N = 2 P = 2 . Assuming − T c < t ≤ T c and the number of paths M = 1 , theLS weight is given by (cid:13)(cid:13)(cid:13)(cid:13)(cid:20) y E y L (cid:21) − (cid:20) r E r L (cid:21)(cid:13)(cid:13)(cid:13)(cid:13) (16a) = (cid:13)(cid:13)(cid:13)(cid:13)(cid:20) y E y L (cid:21) − c (cid:20) . − t . t (cid:21)(cid:13)(cid:13)(cid:13)(cid:13) (16b) = (cid:13)(cid:13)(cid:13)(cid:13)(cid:20) y E y L (cid:21) − c (cid:20) − (cid:21) t − c (cid:20) . . (cid:21)(cid:13)(cid:13)(cid:13)(cid:13) . (16c)Note that, since the two samples are one chip apart, the noise levels for the two samples are independentaccording to (7). Hence, by Remark 2, (16) is the ML weight as well. From (9), parameter vectors forthe LS weight (16) are defined as A ( t ) = (cid:20) . − t . t (cid:21) , (17a) B ( c , t ) = c (cid:20) − (cid:21) , (17b) b ( c , t ) = c (cid:20) . . (cid:21) . (17c) −1.5 −1 −0.5 0 0.5 1 1.5−1−0.8−0.6−0.4−0.200.20.40.60.81 Chip offset [chips] D i sc r i m i na t o r ou t pu t Wide EL ( T S = T C )Narrow EL, ( T S = 0.2 T C ) D , T S = T C /2 D , T S = T C /4 D , T S = T C /8 Fig. 2. Discriminator responses.
Applying the LIMS algorithm to the LS weight (16) with the initial value t , = 0 , the first iteration gives c , = (cid:0) A H ( t , ) A ( t , ) (cid:1) − A H ( t , ) (cid:20) y E y L (cid:21) = y E + y L (18)by (12), and t , = t , − βD by (15), where D = −ℜ (cid:26) B H ( c , , t , ) (cid:18)(cid:20) y E y L (cid:21) − c , A ( t , ) (cid:19)(cid:27) = −ℜ { ( y E + y L ) ∗ ( y E − y L ) } = | y E | − | y L | . (19)As expected, D (19), the approximate gradient with respect to t , is the same as the EL discriminator.Similarly, for T s < T c , we can derive D = −ℜ (cid:8) B H ( c , , C − ( y − c , A (0)) (cid:9) , (20)where y = [ y EP , . . . , y E , y L , . . . , y LP ] T , c , = (cid:0) A H (0) C − A (0) (cid:1) − A H (0) C − y , (21)and A (0) and B ( c , , are defined similarly to (17), by letting M = 1 in (9). Note that, in (20) and(21), noise whitening is applied. Then, (20) can be regarded as a generalized P -point discriminatoras an extension of the conventional -point discriminator like the EL discriminator. The quantity D asa function of the code phase offset for different values of T s is evaluated and plotted in Fig. 2. Thediscriminator responses of the EL discriminators with sample spacing T c (wide discriminator) and . T c (narrow discriminator) are also plotted in Fig. 2 for comparison. In the figure, it is observed that, as T s decreases, the discriminator curve gets steeper around the zero offset, which implies stronger resistanceagainst perturbation after the DLL is locked at the peak of the ACF [13].IV. N UMERICAL RESULTS
This section presents simulation results to demonstrate the performance of the LIMS algorithm. Forcomparison, the conventional EL discriminators for PN code phase tracking with wide ( T c ) and narrow( . T c ) spacings are also considered. In particular, performance is evaluated with the estimation accuracyof the delay of the first arrival path and, as the performance measure, the mean square error (MSE) ofthe delay estimate is used. For all simulations, the results are obtained by averaging over trials with independent realizations of noise and channel impulse response. As for the PN code sequence, we usethe global positioning system (GPS) coarse-acquisition (C/A) code sequence [16], which has 1023 chipsper period with a . MHz chip rate; one chip occupies T c = 1 / (1 . × ) seconds. The receivedsignal is generated by superimposing the attenuated, phase rotated, and delayed replicas of the C/A codesequence according to the multipath channel realization, and adding white Gaussian noise with the powerspectral density N . The received carrier power is defined as C = E {| h ( t ) ∗ s ( t ) | } in (2), where theexpectation is over the channel realizations. The carrier power C varies resulting in a carrier to noisedensity ratio ( C/N ) between dBHz and dBHz. At the receiver, the integration time of the correlator, T i , is set to be ms, which corresponds to periods of the GPS C/A code sequence.For all simulations, it is assumed that a window of the correlator output signal of T c width is availablefor the EL discriminator and the LIMS algorithm. The center of the window is defined as the zero lag, i.e., ζ = 0 . It is also assumed that the coarse acquisition of the PN code phase is perfect regardless of C/N so that the true first arrival path is placed at a random lag uniformly distributed over [ − . T c , . T c ) .Note that, since the size of the observation is too small to compute the sample correlation matrix of thereceived signal, conventional super-resolution techniques like MUSIC and ESPRIT cannot be applied tothe considered simulation environment. A. Simulation system
For the simulation of the LIMS algorithm, the sampling period T s = 0 . T c is used. As described inSection III-A, equation (12) is used for the channel coefficient update and, for ease of implementation,the signed gradient descent algorithm is used in place of (15) for the delay update. That is, t l = t l − + β · sgn h ℜ n ˜ B H ( c l , t l − ) (cid:16) ˜ y − ˜ A ( t l − ) c l (cid:17)oi , (22)where sgn[ · ] is the element-wise sign function. For the step size β in (22), we use β = T s / . To startthe iteration, the initial value t = [ t , , . . . , t M − , ] T is given such that t , = − . T c , t M − , = 0 . T c ,and t m, , < m < M − evenly divide the interval ( − . T c , . T c ) . With this initial value, the finalestimates of the channel coefficient vector ˆ c = [ˆ c , . . . , ˆ c M − ] T and the delays vector ˆ t = [ˆ t , . . . , ˆ t M − ] T are obtained by c and t , respectively.Upon completion of the estimation process, a simple path validation process follows to reduce theprobability of false alarm; an estimated path is valid if the estimated power ( | ˆ c m | ) of the path is largerthan % of the total estimated power ( P M − m =0 | ˆ c m | ). The path with the smallest delay among the validpaths is chosen for the first arrival path.Finally, in the simulations of the LIMS algorithm, two different cases, without noise whitening, i.e., G = I M , and with noise whitening, i.e., G = C − / , are included for comparison. B. Two-path non-fading channel
The first test channel considered for simulation is a two-path non-fading channel with c CH = h √ , √ i T and t CH = (cid:2) , T c (cid:3) T . The performance results in terms of the MSE with respect to C/N are shown in Fig.3(a). As a benchmark, the CRB, which is derived in Appendix, is also plotted in the same figure. Withthis channel, the correlator output has a flat top when there is no noise. As a result, the conventional ELdiscriminator, which assumes a single path and the triangular ACF, cannot perform well with the channel.On the other hand, the LIMS algorithm can handle the two paths simultaneously and outperforms theEL discriminator. It is also observed that the LIMS algorithm with noise whitening consistently performsbetter than that without noise whitening. In particular, at high C/N ’s, the LIMS algorithm with noisewhitening asymptotically achieves the CRB.In Fig. 3(b), learning curves of the LIMS algorithm with noise whitening for the first arrival path delayare depicted for independent trials at C/N = 30 dBHz. Observe that, even with a very large initialdeviation, the algorithm successfully converges to the vicinity of the true first arrival path delay within iterations.
20 30 40 50 6010 −6 −5 −4 −3 −2 −1 C / N [dBHz] M ean S qua r e E rr o r [ c h i p s ] Wide EL ( T c )Narrow EL (0.2 T c )LIMS, w/o whitening, M =2LIMS, w/ whitening, M =2CRB (a) Mean square errors F i r s t a rr i v a l pa t h e s t i m a t i on e rr o r ( t − τ ) [ c h i p ] ± CRB
Learning curve (b) Learning curves at
C/N = 30 dBHzFig. 3. Mean square errors and learning curves of the first arrival path delay estimation: c CH = (cid:2) / √ , / √ (cid:3) T and t CH = [0 , T c / T .
20 30 40 50 6010 −4 −3 −2 −1 C / N [dBHz] M ean S qua r e E rr o r [ c h i p s ] Wide EL ( T c )Narrow EL (0.2 T s )LIMS, w/o whitening, M =1LIMS, w/ whitening, M =1LIMS, w/o whitening, M =2LIMS, w/ whitening, M =2LIMS, w/o whitening, M =3LIMS, w/ whitening, M =3LIMS, w/o whitening, M =4LIMS, w/ whitening, M =4 Fig. 4. Mean square errors of the first arrival path delay estimation: Single-path fading channel.
C. Single-path fading channel
For the second simulation, a channel with only one path is considered, where the channel coefficient γ is a zero-mean circular symmetric complex Gaussian random variable with E {| γ | } = 1 . The MSEperformances of the EL discriminators and the LIMS algorithms with different settings for this channel areplotted in Fig. 4. As shown in the figure, the EL discriminators, especially the one with narrow spacing,perform very well for the channel. This is because the assumption that the correlator output has a triangularshape, on which the EL discriminator is based, is valid for this case. Among all performance curves plottedin Fig. 4, the best is achieved by the LIMS algorithm with M = 1 and noise whitening. Compared to theperformance of the EL discriminators, the best performance is an order of magnitude superior. In termsof the generalized discriminator presented in Section III-D, this performance improvement arises fromthe fact that the LIMS algorithm can make use of a larger number of received signal samples and thushas a steeper discriminator response. Finally, in Fig. 4, note that the LIMS algorithm with M = 2 andnoise whitening has a performance comparable to that with M = 1 and noise whitening even though thenumber of assumed paths is not the same as the true value.
20 30 40 50 6010 −2 −1 C / N [dBHz] M ean S qua r e E rr o r [ c h i p s ] Wide EL ( T c )Narrow EL (0.2 T c )LIMS, w/o whitening, M =1LIMS, w/ whitening, M =1LIMS, w/o whitening, M =2LIMS, w/ whitening, M =2LIMS, w/o whitening, M =3LIMS, w/ whitening, M =3LIMS, w/o whitening, M =4LIMS, w/ whitening, M =4 (a) Channel A
20 30 40 50 6010 −2 −1 C / N [dBHz] M ean S qua r e E rr o r [ c h i p s ] Wide EL ( T c )Narrow EL (0.2 T c )LIMS, w/o whitening, M =1LIMS, w/ whitening, M =1LIMS, w/o whitening, M =2LIMS, w/ whitening, M =2LIMS, w/o whitening, M =3LIMS, w/ whitening, M =3LIMS, w/o whitening, M =4LIMS, w/ whitening, M =4 (b) Channel BFig. 5. Mean square errors of the first arrival path delay estimation: Multipath fading channels without band limitation. D. Multipath fading channels
The third test channel models considered for simulation are the multipath fading channels. Two multipathfading channel models, which are denoted by channels A and B, respectively, are considered. Channels Aand B have three and four paths, respectively, and the power-delay profiles are given in Table II. Note that,for both channels, the average power of the first arrival path is dB lower than that of the strongest path.With the given power-delay profiles, the channel coefficients are generated independently as zero-meancircular symmetric complex Gaussian random variables with corresponding variances for each trial. TheMSE performances of the EL discriminators and the LIMS algorithms with different settings are plottedin Figs. 5(a) and 5(b) for channels A and B, respectively. In Figs. 5(a) and 5(b), it is observed that theEL discriminators yield very poor performances. This is because the EL discriminator tends to keep trackof the largest magnitude point at the correlator output, which would be located mostly near the strongestpath and hardly near the first arrival path. Due to the same reason, the performance of the LIMS algorithmwith M = 1 is also bounded. However, the LIMS algorithms with M ≥ show much better performancessince they can track multiple paths simultaneously, and one of the tracked paths is possibly close to thetrue first arrival path.The effect of limited bandwidth on the performance should also be investigated. In practice, the receivedsignal is bandlimited to the pre-correlation bandwidth by a band pass filter before the cross-correlationwith the PN code sequence replica, and thus the ACF may not have as sharp a triangular shape as (6)has [17, 4.2]. This causes a mismatch between the assumed signal model and the true signal, and canentail some performance loss. For the simulation, an ideal band pass filter with a rectangular frequencyresponse centered at the carrier frequency of the signal is assumed. Figs. 6(a) and 6(a) show the resultsobtained for channels A and B, respectively, when an MHz pre-correlation bandwidth is applied. Bycomparing Figs. 5 and 6, it is found that there is no significant loss when using an MHz pre-correlation
TABLE IIP
ARAMETERS OF THE THREE - AND FOUR - PATH FADING CHANNELS .Channel A Channel BPath Relative power (dB) Delay ( T c ) Relative power (dB) Delay ( T c )1 -7.0 0 -7.0 02 0 0.3 -7.0 0.23 -2.2 0.5 0 0.44 - - -2.2 0.6
20 30 40 50 6010 −2 −1 C / N [dBHz] M ean S qua r e E rr o r [ c h i p s ] Wide EL ( T c )Narrow EL (0.2 T c )LIMS, w/o whitening, M =1LIMS, w/ whitening, M =1LIMS, w/o whitening, M =2LIMS, w/ whitening, M =2LIMS, w/o whitening, M =3LIMS, w/ whitening, M =3LIMS, w/o whitening, M =4LIMS, w/ whitening, M =4 (a) Channel A
20 30 40 50 6010 −2 −1 C / N [dBHz] M ean S qua r e E rr o r [ c h i p s ] Wide EL ( T c )Narrow EL (0.2 T c )LIMS, w/o whitening, M =1LIMS, w/ whitening, M =1LIMS, w/o whitening, M =2LIMS, w/ whitening, M =2LIMS, w/o whitening, M =3LIMS, w/ whitening, M =3LIMS, w/o whitening, M =4LIMS, w/ whitening, M =4 (b) Channel BFig. 6. Mean square errors of the first arrival path delay estimation: Multipath fading channels with MHz pre-correlation bandwidth.
20 30 40 50 6010 −2 −1 C / N [dBHz] M ean S qua r e E rr o r [ c h i p s ] Wide EL ( T c )Narrow EL (0.2 T c )LIMS, w/o whitening, M =1LIMS, w/ whitening, M =1LIMS, w/o whitening, M =2LIMS, w/ whitening, M =2LIMS, w/o whitening, M =3LIMS, w/ whitening, M =3LIMS, w/o whitening, M =4LIMS, w/ whitening, M =4 (a) Channel A
20 30 40 50 6010 −2 −1 C / N [dBHz] M ean S qua r e E rr o r [ c h i p s ] Wide EL ( T c )Narrow EL (0.2 T c )LIMS, w/o whitening, M =1LIMS, w/ whitening, M =1LIMS, w/o whitening, M =2LIMS, w/ whitening, M =2LIMS, w/o whitening, M =3LIMS, w/ whitening, M =3LIMS, w/o whitening, M =4LIMS, w/ whitening, M =4 (b) Channel BFig. 7. Mean square errors of the first arrival path delay estimation: Multipath fading channels with MHz pre-correlation bandwidth. bandwidth. The results with a MHz pre-correlation bandwidth are also shown in Figs. 7(a) and 7(b).There is apparent loss in this case, especially for the LIMS algorithms with M = 3 and M = 4 at high C/N ’s. However, the performance loss is still not significant, and, therefore, one can conclude that theLIMS algorithm is quite robust to the band limitation.V. C ONCLUSION
In this paper, the multipath channel estimation problem for direct sequence spread spectrum signalswas investigated. A new technique based on the least-squares criterion, denoted as the least-squaresbased iterative multipath super-resolution (LIMS) algorithm, was proposed for the resolution of shortdelay multipaths. Unlike conventional super-resolution techniques such as MUSIC and ESPRIT, the LIMSalgorithm does not rely on subspace decomposition, and, therefore, has practical advantages. In particular,the LIMS algorithm can work with a short observation of the received signal regardless of the correlationbetween path coefficients. In addition, due to its iterative operation, the LIMS algorithm is suitable forrecursive multipath tracking. We also presented a view of the LIMS algorithm as a generalization of theconventional EL discriminator for PN code phase tracking. Through numerical simulations, it was shownthat the performance of the LIMS algorithm is superior to that of the EL discriminator in estimating the delay of the first arrival path. As a result, the LIMS algorithm can be a useful technique for applicationssuch as radio ranging, global navigation satellite systems (GNSSs), radio propagation analysis, and timesynchronization and multipath tracking in wireless communication systems.A PPENDIX T HE C RAMER -R AO BOUND FOR THE FIRST ARRIVAL PATH DELAY ESTIMATION
From the signal model (3), the inverse Fisher information matrix [18, B.3.3] for the estimation of thechannel coefficients and the path delays is computed as F − = T i N (cid:0) ℜ (cid:8) D H C − D (cid:9)(cid:1) − , (23)where D ∈ C N × K is a concatenated matrix defined as D = [ A ( t CH ) B ( c CH , t CH )] , (24)and A ( · ) and B ( · , · ) are defined in (9). The mean square error (MSE) of each parameter estimation islower bounded by the corresponding diagonal element of the inverse Fisher information matrix. Thus, theMSE of the first arrival path delay estimation is lower bounded by the ( K + 1) -th diagonal element ofthe inverse Fisher information matrix (23). R EFERENCES [1] T. S. Rappaport, “Characterization of UHF multipath radio channels in factory buildings,”
IEEE Trans. Antennas Propagat. , vol. 37,No. 8, pp. 1058–1069, Aug. 1989.[2] R. O. Schmidt, “Multiple emitter location and signal parameter estimation,”
IEEE Trans. Antennas Propagat. , vol. AP-34, pp. 276–280,Mar. 1986.[3] R. Roy and T. Kailath, “ESPIRIT estimation of signal parameters via rotational invariance techniques,”
IEEE Trans. Acoust., Speech,Signal Processing , vol. 37, No. 7, pp. 984–995, July 1989.[4] A. M. Bruckstein, T.-J. Shan, and T. Kailath, “The resolution of overlapping echos,”
IEEE Trans. Acoust., Speech, Signal Processing ,vol. ASSP-33,No. 6, pp. 1357–1367, Dec. 1985.[5] J. Kusuma, I. Maravi´c, and M. Vetterli, “Sampling with finite innovation rate: Channel and timing estimation in UWB and GPS,” in
Proc. ICC , May 2003.[6] T. Manabe and H. Takai, “Superresolution of multipath delay profiles measured by PN correlation method,”
IEEE Trans. AntennasPropagat. , vol. 40, No. 5, pp. 500–509, May. 1992.[7] F. Bouchereau, D. Brady, and C. Lanzl, “Multipath delay estimation using a superresolution PN-correlation method,”
IEEE Trans. onSignal Processing , vol. 49, No. 5, pp. 938–949, May. 2001.[8] J. W. Betz and K. R. Kolodziejski, “Generalized theory of code tracking with an early-late discriminator. Rart I: Lower bound andcoherent processing,”
IEEE Trans. Aerospace and Electronic Systems , vol. 45, No. 4, pp. 1538–1550, Oct. 2009.[9] J. W. Betz and K. R. Kolodziejski, “Generalized theory of code tracking with an early-late discriminator. Part II: Noncoherent processingand numerical results,”
IEEE Trans. Aerospace and Electronic Systems , vol. 45, No. 4, pp. 1551–1564, Oct. 2009.[10] R. Kumaresan and A. K. Shaw, “High resolution bearing estimation without eigendecomposition,” in
Proc. IEEE Int. Conf. Acoust.,Speech, Signal Processing , Tampa, FL, pp. 576–579, Mar. 1985.[11] Y. Bresler and A. Macovski, “Exact maximum likelihood parameter estimation of superimposed exponential signals in noise,”
IEEETrans. Acoust., Speech, Signal Processing , vol. ASSP-34, pp. 1081–1089, Oct. 1986.[12] P. Stoica and A. Nehorai, “MUSIC, maximum likelihood, and Cramer-Rao bound,”
IEEE Trans. Acoust., Speech, Signal Processing ,vol. 37, No. 5, pp. 720–741, May. 1989.[13] M. Irsigler and B. Eissfeller, “Comparison of multipath mitigation techniques with consideration of future signal structures,” in
Proc.ION GPS/GNSS , Portland, OR, Sep. 2003.[14] D. Mu˜noz, F. Bouchereau, C. Vargas, and R. E.-Caldera,
Position Location Techniques and Applications , Academic Press, 2009.[15] S. Haykin,
Adaptive Filter Theory , fourth edition, Prentice Hall, 2002.[16] Pratap Misra and Per Enge,
Global Positioning System: Signals, Measurements, and Performance , second edition, Ganga-Jamuna Press,2001.[17] K. Yu, I. Sharp, and Y. J. Guo,
Ground-Based Wireless Positioning , Wiley, 2009.[18] P. Stoica and R. Moses,