**Efficient FFT-based Algorithms for Multipath Interference Mitigation in GNSS**

Renbiao Wu, Qiongqiong Jia, Wenyi Wang and Jie Li

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/59940

## **1. Introduction**

[11] Lee, MH., *Jacket Matrices-Construction and Its Application for Fast Cooperative Wireless Signal Processin*g. LAP LAMBERT Academic publishing, Germany, November, 2012.

[12] Wang, CL. and Chen, CY., "High-throughput VLSI architectures for the 1-D and 2-D discrete cosine transform," *IEEE Trans. Circuits Syst. Video Technol*., vol. 5, pp. 31–40,

[13] Wang, Z., "Fast Algorithm for the Discrete W Transform and for the Discrete Fourier Transform," *IEEE Trans. on Acoustics, Speech and Signal Process*., vol. 32, No. 4, pp. 803

[14] Lee, MH., "High speed multidimensional systolic arrays for discrete Fourier trans‐

[15] Kim, KJ., Fan, Y.,Iltis, R. A., Poor, H. V., and Lee, M. H., "A reduced feedback pre‐ coder for MIMO-OFDM cooperative diversity system," *IEEE Trans. Veh. Technol.*, vol.

[16] Jang, U.,Cho, K.,Ryu, W., and Lee, HJ., "Interference management with block diago‐ nalization for macro/femto coexisting networks," *ETRI Journal*, vol. 34, pp. 297–307,

[17] Hou, HS., "A fast recursive algorithm for computing the discrete cosine transform," *IEEE Trans. Acoust., Speech, Signal Process*., vol. 35, no. 10, pp. 1455–1461, Oct. 1987. [18] Chen, WH., Smith, CH., and Fralick, SC., "A fast computational algorithm for the discrete cosine transform," *IEEE Trans. Commun*., vol. 25, no. 9, pp. 1004–1009, Sep.

[19] Spencer, Q. H., Lee, A. Swindlehurst, M. Haardt, "Zero-forcing methods for down‐ link spatial multiplexing in multiuser MIMO channels", *IEEE Trans. Signal Process*.,

[20] Andrews, HC., Caspari, KL., "A generalized technique for spectral analysis," *IEEE*

[21] Reju, VG.,Koh, SN.,Soon, IY., "Convolution using discrete sine and cosine trans‐

[22] Lee, MH., Khan, MHA., Sarker, MA. L., Guo, Y. and Kim, KJ., "A MIMO LTE precod‐ ing based on fast diagonal weighted Jacket matrices", *Fiber and Integrated Optics*, *Tay‐*

[23] Khan, MHA., Li, J., Lee, MH., "A block diagonal Jacket matrices for MIMO broadcast channel" *IEEE International Symposium on Broadband Multimedia Systems and Broadcast‐*

forms", IEEE Signal Processing Letters, vol. 14, no. 7, July 2007.

*lor and Francis*, Invited paper, vol. 31, no. 2, pp. 111-132, March 2012.

form," *IEEE Trans. Circuits Syst. II,* vol. 39, no. 12, pp. 876–879, Dec. 1992.

Feb. 1995.

June 2012.

1977.

– 816, Aug. 1984.

61, pp. 584–596, Feb. 2012.

24 Fourier Transform - Signal Processing and Physical Sciences

vol. 52,no. 2, Feb. 2004.

*Trans. Computers*, vol. 19, no. 1, pp.16-17, 1970.

*ing*, Brunel University, June 4-7th, 2013, UK.

GNSS (Global Navigation Satellite System) has been found application in many areas. In some cases, the performance requirements for GNSS are very high. There are many error sources that would degrade the positioning performance of GNSS, e.g., clock errors, ephemeris errors, tropospheric propagation delay, and multipath. Many positioning errors mentioned above are constant for all GNSS receivers in a given small area and can be removed or reduced by using the popular differential technique. However, due to the geographical position difference between the reference station and the receiver, the multipath environment of receivers, such as amplitudes and number of multipath signals, is totally different with that of the reference station. Thus differential technique can not eliminate the multipath error. Many studies have shown that the multipath interference will lead to a position error around several meters which endangers the reliability and accuracy of GNSS. Therefore, multipath interference mitigation has been a hot topic in the field of satellite navigation receiver design.

Multipath interferences are the signals reflected by the objects around the GNSS receiver. Then the multipath interference and the LOS (line of sight) signal are simultaneously received by antenna which brings a phase distortion in the tracking loops of receivers. Finally, the phase distortion results in the tracking and positioning error.

There have been many studies about the effect of multipath interference. Kalyanaraman et al. in [1] analyzed the multipath effects on code tracking loop. Kos et al. [2] provided a detailed analysis of the multipath effects on the positioning error. Main multipath interference mitigation techniques are based on antenna technique and signal processing algorithms.

© 2015 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited.

By mounting the antenna in a well-designed place based on the multipath environment, Maqsood et al. in [3] compared the performances of multipath interference mitigation for different antennas. Ray et al. in [4] proposed a multipath mitigation algorithm based on an antenna array. A design principle of antenna was proposed by Alfred et al. in [5] for satellite navigation landing system. The multipath mitigation technique based on the special antenna can only suppress the multipath signal coming from the ground below the antenna. However, it is useless for the multipath interference signal reflected by the objects above the antenna.

Popular signal processing techniques are the narrow correlator and MEDLL (Multipath Estimated Delay Locked Loop). Narrow correlator is adopted in many GNSS receivers which suppresses multipath interference by reducing the early-late correlator spacing. The signal model for the narrow correlator is provided by Michael et al. in [6]. Cannon et al. in [7] analyzed the performance of narrow correlator in the satellite navigation system. The performance of narrow correlator can be improved by decreasing the early-late correlator spacing. However, the narrow correlator assumes that the bandwidth of the received signal is infinite which is invalid in the practical applications. Thus, when the early-late correlation spacing is less than the reciprocal of the channel bandwidth, the tracking error cannot be further decreased by reducing the early-late correlator spacing. MEDLL is a multipath interference mitigation algorithm based on the statistical theory as in [8], which estimates the time delays via the maximum likelihood criterion. Therefore, the complexity of MEDLL is much higher than narrow correlator.

This chapter firstly presents a new acquisition algorithm for GNSS. Then two multipath interference mitigation algorithms based on the DPE (Decoupled Parameter Estimation) parameter estimation algorithms are presented. The FFT algorithm is utilized to reduce the computational complexity in all proposed algorithms. Numerical results are provided to demonstrate the performances of the proposed algorithms. The remainder of this chapter is arranged as follows. Section 2 presents the new acquisition algorithm. Two multipath inter‐ ference mitigation algorithms are separately presented in Section 3. Section 4 concludes the chapter.

## **2. A novel GNSS acquisition algorithm**

## **2.1. Data model and problem formulation**

GNSS signal is composed of three parts, the carrier, the ranging codes and the data codes. In order to facilitate the presentation, GPS is taken as an example. However, the proposed algorithm is also suitable for other GNSS system. The carrier of GPS (Global Navigation System) consists of the three different frequency bands, i.e., L1 (1575.42 MHz), L2 (1227.6 MHz) and L5 (1176.45MHz), see [9-10]. The ranging codes are pseudo random noise (PRN) codes including C/A code, P code and M code which are known in advance. Data codes or the navigation data are binary codes containing the satellite ephemeris, satellite status, clock system, the orbit perturbation correction of the satellite clock motion state and the atmospheric refraction correction, etc..

Suppose that one receiver obtained the signals from P satellites and the received data can be expressed as

$$\mathbf{y}(t) = \sum\_{p=1}^{p} \alpha\_{p} d\_{p}(t - \tau\_{p}) c\_{p}(t - \tau\_{p}) e^{i\alpha\_{\Phi}(t - \tau\_{p})} + \mathbf{e}(t) \tag{1}$$

where *dp*(*t*) represents the navigation data of the *<sup>p</sup> th* satellite, *cp*(*t*) is the C/A code of the *<sup>p</sup> th* satellite, *e*(*t*) is the thermal noise, *αp*, *τp* and *ωdp* denote the amplitude, time delay and Doppler frequency of the *p th* satellite, respectively. It should be noted that only the L1 signal is considered in this chapter.

After A/D conversion, the data can be written as

By mounting the antenna in a well-designed place based on the multipath environment, Maqsood et al. in [3] compared the performances of multipath interference mitigation for different antennas. Ray et al. in [4] proposed a multipath mitigation algorithm based on an antenna array. A design principle of antenna was proposed by Alfred et al. in [5] for satellite navigation landing system. The multipath mitigation technique based on the special antenna can only suppress the multipath signal coming from the ground below the antenna. However, it is useless for the multipath interference signal reflected by the objects above the antenna. Popular signal processing techniques are the narrow correlator and MEDLL (Multipath Estimated Delay Locked Loop). Narrow correlator is adopted in many GNSS receivers which suppresses multipath interference by reducing the early-late correlator spacing. The signal model for the narrow correlator is provided by Michael et al. in [6]. Cannon et al. in [7] analyzed the performance of narrow correlator in the satellite navigation system. The performance of narrow correlator can be improved by decreasing the early-late correlator spacing. However, the narrow correlator assumes that the bandwidth of the received signal is infinite which is invalid in the practical applications. Thus, when the early-late correlation spacing is less than the reciprocal of the channel bandwidth, the tracking error cannot be further decreased by reducing the early-late correlator spacing. MEDLL is a multipath interference mitigation algorithm based on the statistical theory as in [8], which estimates the time delays via the maximum likelihood criterion. Therefore, the complexity of MEDLL is much higher than

This chapter firstly presents a new acquisition algorithm for GNSS. Then two multipath interference mitigation algorithms based on the DPE (Decoupled Parameter Estimation) parameter estimation algorithms are presented. The FFT algorithm is utilized to reduce the computational complexity in all proposed algorithms. Numerical results are provided to demonstrate the performances of the proposed algorithms. The remainder of this chapter is arranged as follows. Section 2 presents the new acquisition algorithm. Two multipath inter‐ ference mitigation algorithms are separately presented in Section 3. Section 4 concludes the

GNSS signal is composed of three parts, the carrier, the ranging codes and the data codes. In order to facilitate the presentation, GPS is taken as an example. However, the proposed algorithm is also suitable for other GNSS system. The carrier of GPS (Global Navigation System) consists of the three different frequency bands, i.e., L1 (1575.42 MHz), L2 (1227.6 MHz) and L5 (1176.45MHz), see [9-10]. The ranging codes are pseudo random noise (PRN) codes including C/A code, P code and M code which are known in advance. Data codes or the navigation data are binary codes containing the satellite ephemeris, satellite status, clock system, the orbit perturbation correction of the satellite clock motion state and the atmospheric

narrow correlator.

chapter.

**2. A novel GNSS acquisition algorithm**

**2.1. Data model and problem formulation**

26 Fourier Transform - Signal Processing and Physical Sciences

refraction correction, etc..

$$\text{cy}(n) = \sum\_{p=1}^{p} \alpha\_p d\_p (n - \tau\_p) c\_p (n - \tau\_p) e^{\langle \alpha\_{\phi} \cdot (n - \tau\_{\rho}) \rangle} + \text{e}(n) \tag{2}$$

The conventional acquisition procedure is a two-dimensional searching algorithm, which is time consuming even when using parallel searching algorithm. Moreover, the frequency resolution of the conventional acquisition algorithm cannot satisfy the requirements of the tracking loop. Hence, a more accurate frequency estimation is required before tracking in the conventional receiver.

Therefore, a new GNSS acquisition algorithm is proposed in this section. The Doppler frequency is firstly estimated. After that, the initial code phase is then obtained via NLS (Nonlinear Least Square) fitting. Compared with the conventional acquisition algorithm, the proposed algorithm can obtain a comparative performance with a lower computational complexity.

#### **2.2. Principle of the novel acquisition algorithm**

It can be noted from the data model in section 2.1 that there are three unknown parameters in the acquisition process, which are the PRN index of the *p th* satellite, the corresponding Doppler frequency *ωdp* and its initial code phase *τp*.

#### *2.2.1. Doppler frequency estimation*

Due to the navigation data and C/A code are ±1 in GPS system, it can be seen from equation (2) that the Fourier spectrum of each satellite contains multiple frequency components. Therefore, the Doppler frequency cannot be extracted from the spectrum directly. To eliminate the influence on frequency estimation brought by the navigation data and C/A code, we square equation (2) as in [11]

$$\begin{split} \mathbf{y}^{2}(n) &= \sum\_{\rho=1}^{p} \alpha\_{\rho}^{2} d\_{\rho}^{2} (n - \tau\_{\rho}) c\_{\rho}^{2} (n - \tau\_{\rho}) e^{\langle 2a\_{\phi} (n - \tau\_{\rho}) \rangle} + 2e(n) \sum\_{\rho=1}^{p} \alpha\_{\rho} d\_{\rho} (n - \tau\_{\rho}) c\_{\rho} (n - \tau\_{\rho}) e^{\langle a\_{\phi} (n - \tau\_{\rho}) \rangle} \\ &+ 2 \sum\_{\rho=1}^{p-1} \sum\_{\rho=1}^{p} \alpha\_{\rho} a\_{\rho} d\_{\rho} (n - \tau\_{\rho}) d\_{\rho} (n - \tau\_{\rho}) c\_{\rho} (n - \tau\_{\rho}) c\_{\rho} (n - \tau\_{\rho}) e^{\langle a\_{\phi} (n - \tau\_{\rho}) + a\_{\phi} (n - \tau\_{\iota}) \rangle} + e^{2}(n) \end{split} \tag{3}$$

Since the value of navigation data and C/A code is ±1 in GPS system, then the above equation can be simplified as

$$\text{y}^2(n) = \sum\_{\rho=1}^p \alpha\_\rho^2 e^{\sqrt{2}a a\_\theta \left(n - \varepsilon\_\rho\right)} + e\_1(n) \tag{4}$$

where,

$$\begin{split} e\_i(n) &= 2\sum\_{p=1}^{P-1} \sum\_{r=p+1}^{P} \alpha\_p \alpha\_r d\_p(n - \tau\_p) d\_r(n - \tau\_r) c\_p(n - \tau\_p) c\_r(n - \tau\_r) e^{\int \left[ a\_{\theta^\*}(n - \tau\_\rho) + a\_{\theta^\*}(n - \tau\_r) \right]} \\ &+ 2e(n) \sum\_{p=1}^{P} \alpha\_p d\_p(n - \tau\_p) c\_p(n - \tau\_p) e^{\int a\_{\theta^\*}(n - \tau\_\rho)} + e^2(n) \end{split} \tag{5}$$

It can be seen from equation (4) that the spectrum of each satellite only contains a single dominant frequency component corresponding to the Doppler frequency. For the reason that not only the correlation of C/A code of different satellites is very small but also the signal and noise are uncorrelated, it is obviously that in equation (5) the product terms of different C/A codes as well as the product of noise and signal are close to zero. Fourier analysis method can be directly used to estimate the signal frequency.

The Fourier transform of equation (4) can be expressed as

$$F(\alpha \mathbf{e}) = \left| \sum\_{n=-N/2}^{N/2-1} \mathbf{y}^2(n) \mathbf{e}^{-j\alpha n} \right| \tag{6}$$

where *ω* = 2*ωdp* can be obtained based on the locations of the first *p* dominant peaks of *F*(*ω*). A more accurate estimation result can be achieved by using FFT (Fast Fourier Transform) padding with zeros.

The estimated value of *ω* ^ is in the range of −*π*, *π* , thus from equation (6) we know that the estimated value of *ω* ^ *dp* −*π* / 2, *π* / 2 . As *ωdp* −*π*, *π* , the estimated frequency could be *ω* ^ *dp* or *ω* ^ *dp* + *π*. Therefore, the frequency ambiguity will lead to a mistake. Due to the IF (Intermediate Frequency) of the satellite signals is definitely known after down conversion and the range of Doppler shift is between ±10kHz [12], the difference between true frequency and the ambiguity (6)

. A

Chapter 1

*n N F y ne*

where 2 

zeros.

/2 1 <sup>2</sup> / 2 ( ) ()

*<sup>N</sup> <sup>j</sup> <sup>n</sup>*

*dp* can be obtained based on the locations of the first *p* dominant peaks of *F*( )

Figure 1 Frequency spectrum before and after square operation **Figure 1.** Frequency spectrum before and after square operation

2 2 2 2 2( ) ( )

*P P j n j n pp p p p pp p p p*

*d n dn c n cn e e n*

Since the value of navigation data and C/A code is ±1 in GPS system, then the above equation

<sup>1</sup> [ ( ) ( )]

 t


*P P jn n prp P r r p P r r*

( ) 2

It can be seen from equation (4) that the spectrum of each satellite only contains a single dominant frequency component corresponding to the Doppler frequency. For the reason that not only the correlation of C/A code of different satellites is very small but also the signal and noise are uncorrelated, it is obviously that in equation (5) the product terms of different C/A codes as well as the product of noise and signal are close to zero. Fourier analysis method can

> /2 1 <sup>2</sup> / 2 ( ) ()

=- <sup>=</sup> å

w

*n N F y ne*


w

*<sup>N</sup> j n*

where *ω* = 2*ωdp* can be obtained based on the locations of the first *p* dominant peaks of *F*(*ω*). A more accurate estimation result can be achieved by using FFT (Fast Fourier Transform)

*dp* + *π*. Therefore, the frequency ambiguity will lead to a mistake. Due to the IF (Intermediate Frequency) of the satellite signals is definitely known after down conversion and the range of Doppler shift is between ±10kHz [12], the difference between true frequency and the ambiguity


*dp p*

 t

w t  t


*dp p dr r*

 t

(4)

w tw

(6)

^ is in the range of −*π*, *π* , thus from equation (6) we know that the

*dp* −*π* / 2, *π* / 2 . As *ωdp* −*π*, *π* , the estimated frequency could be *ω*

 t *dp p dr r*

 t

 t w t

(3)

(5)

^ *dp* or

*dp p dp p*

at

w tw

1

<sup>1</sup> [ ( ) ( )] <sup>2</sup>


*P P jn n prp P r r p P r r*

+ - --- +

w t

2 ( ) ( ) ( )( ) ( )

2 2 2( )

a

( ) ( ) -

<sup>=</sup> å <sup>+</sup> *dp p <sup>P</sup> j n p p*

*y n e en* w t

1

=

() 2 ( ) ( ) ( ) ( )

*e n d n dn c n cn e*

 t

> t

= - ---

2( ) ( ) ( ) ( )

*en d n c n e e n*

+ -- +

aa

be directly used to estimate the signal frequency.

The Fourier transform of equation (4) can be expressed as

*<sup>P</sup> j n pp p p p*

 t

( ) ( ) ( ) 2( ) ( ) ( )

*y n d n c n e en d n c n e*

= -- + --

1 1

 t

= =

*p p*

 t

å å

1 1

can be simplified as

1

padding with zeros.

estimated value of *ω*

*ω* ^

The estimated value of *ω*

^

where,

aa

at

28 Fourier Transform - Signal Processing and Physical Sciences

1 1

at

= = +

å å

*p rp*

1

=

*p*

å

 t

= = +

å å

*p rp*

one is generally greater than the Doppler shift. Therefore, ambiguity frequency can be excluded by determining whether the frequency is in the range of the Doppler shift. estimated value of ˆ [ / 2, / 2] *dp* . As [ , ] *dp* , the estimated frequency could be ˆ*dp* or ˆ*dp* . Therefore, the frequency ambiguity will lead to a mistake. Due to the IF (Intermediate

, thus from equation (6) we know that the

ˆ is in the range of [ , ]

The frequency spectrum of the received data before and after the square operation is shown in Figure 1. It can be noted that before square operation, there are multiple dominant frequency components. Thus the Doppler frequency cannot be extracted directly. However, after being squared, there is only one dominant Doppler frequency present in the spectrum, so an accurate Doppler frequency can be estimated. Frequency) of the satellite signals is definitely known after down conversion and the range of Doppler shift is between 10kHz [12], the difference between true frequency and the ambiguity one is generally greater than the Doppler shift. Therefore, ambiguity frequency can be excluded by determining whether the frequency is in the range of the Doppler shift.

#### *2.2.2. Code phase estimation and matching*

The estimated value of

After obtaining the Doppler frequency, the initial phases and the PRN index of the C/A code from all satellites can be estimated by using its autocorrelation characteristic. Let *sp*(*n*)=*dp*(*n*)*cp*(*n*)*<sup>e</sup> <sup>j</sup>ωdpn* , then equation (2) can be expressed as

$$\mathbf{y}(n) = \sum\_{p=1}^{P} \alpha\_p \mathbf{s}\_p(n - \tau\_p) + \mathbf{e}(n) \tag{7}$$

As the Doppler frequency *ω* ^ *dq* has been estimated, *ω* ^ *dq* can be used to take the place of *ωdq*. Then the signal sent by the *q th* satellite can be reconstructed as

$$
\hat{s}\_q(n) = d\_q(n)c\_q(n)e^{\hat{\imath}\hat{\alpha}\_{\hat{\omega}}n} \tag{8}
$$

Suppose the time delay from the *<sup>q</sup> th* satellite to the receiver is *τq*, then the received signal can be written as

$$\mathbf{y}\_q(n) = \alpha\_q \hat{\mathbf{s}}\_q(n - \tau\_q) + \mathbf{e}\_2(n) \tag{9}$$

where *e*2(*n*)=∑ *p*=1 *p*≠*q P αpsp*(*n* −*τp*) + *αqs*(*n* −*τq*)−*αqs* ^(*<sup>n</sup>* <sup>−</sup>*τq*) <sup>+</sup> *<sup>e</sup>*(*n*). The DFT (Discrete Fourier Transform) of

equation (9) can be written as

$$\mathbf{y}\_{\;/\;}(k) = \alpha\_q \hat{\mathbf{s}}\_{\;/\;\theta}(k) e^{\beta \mathbf{e}\_q k} + e\_{\;/\;2}(k) \tag{10}$$

where *yf* (*k*), *s* ^ *fq*(*k*) and *e <sup>f</sup>* 2(*k*) are the DFT of *y*(*n*), *s* ^ *<sup>q</sup>*(*n*) and *e*2(*n*) respectively, *ω<sup>q</sup>* = −2*πτ<sup>q</sup> f <sup>s</sup>* / *N* with *f <sup>s</sup>* is the sampling frequency. Note, here *xf* is used to represent the DFT of *x*. Hence, the estimation of *τq* is transformed to the estimation of *ωq*.

To estimate *ωq*, we define the following NLS cost function [12-15]

$$\left| z\_f \left( \hat{\alpha}\_q \right) = \sum\_{k=-N/2}^{N/2-1} \left| y\_f(k) - \hat{s}\_{\hat{\rho}}(k) \alpha\_q e^{\beta \alpha\_q k} \right|^2 \tag{11}$$

Let

$$\hat{\mathbf{S}}\_{\hat{\rho}} = \text{diag}\left\{ \hat{\mathbf{s}}\_{\hat{\rho}}(-N/2), \hat{\mathbf{s}}\_{\hat{\rho}}(-N/2 + \mathbf{l}), \dots, \hat{\mathbf{s}}\_{\hat{\rho}}(N/2 - \mathbf{l}) \right\} \tag{12}$$

$$\mathbf{a}(\alpha o\_q) = \left[ e^{j\alpha\_\circ(-N/2)}, e^{j\alpha\_\circ(-N/2+1)}, \dots, e^{j\alpha\_\circ(N/2-1)} \right]^T \tag{13}$$

$$\mathbf{y}\_f = \left[ \mathbf{y}\_f(-N/2), \mathbf{y}\_f(-N/2+1), \dots, \mathbf{y}\_f(N/2-1) \right]^r \tag{14}$$

where ( )<sup>T</sup> denotes the transpose operation. Then the cost function in equation (11) can be transformed as

$$Z\left(\hat{\boldsymbol{\alpha}}\_{q}\right) = \left\|\mathbf{y}\_{\boldsymbol{f}} - \boldsymbol{\alpha}\_{q}\hat{\mathbf{S}}\_{\boldsymbol{\beta}q}\mathbf{a}(\boldsymbol{\alpha}\_{q})\right\|^{2} \tag{15}$$

where denotes the Euclidean norm. Minimizing *Z*(*ω* ^ *<sup>q</sup>*) with respect to *ωq* yields the estimated *ω* ^ *<sup>q</sup>* as

$$\hat{\boldsymbol{\alpha}}\_q = \arg\max\_{\boldsymbol{\alpha}\_q} \left| \mathbf{a}^H(\boldsymbol{\alpha}\_q) \hat{\mathbf{S}}\_{\boldsymbol{\beta}q}^H \mathbf{y}\_f \right|^2 \tag{16}$$

The time delay estimation can be further obtained by

1

*dq* has been estimated, *ω*

=

*p*

^

30 Fourier Transform - Signal Processing and Physical Sciences

the signal sent by the *q th* satellite can be reconstructed as

*αpsp*(*n* −*τp*) + *αqs*(*n* −*τq*)−*αqs*

estimation of *τq* is transformed to the estimation of *ωq*.

w

As the Doppler frequency *ω*

be written as

where *e*2(*n*)=∑

where *yf* (*k*), *s*

Let

where ( )<sup>T</sup>

transformed as

*p*=1 *p*≠*q*

^

equation (9) can be written as

*P*

*P*

() ( ) ()

*pp p*

 t

^

w

(7)

(9)

*<sup>q</sup>*(*n*) and *e*2(*n*) respectively, *ω<sup>q</sup>* = −2*πτ<sup>q</sup> f <sup>s</sup>* / *N*

(11)

^(*<sup>n</sup>* <sup>−</sup>*τq*) <sup>+</sup> *<sup>e</sup>*(*n*). The DFT (Discrete Fourier Transform) of

(10)

*dq* can be used to take the place of *ωdq*. Then

(8)

= -+ å

<sup>ˆ</sup> ˆ () () () = *dq j n q qq s n d nc ne*

<sup>2</sup> () ( ) () = -+ ˆ *q qq q yn sn en* a

<sup>2</sup> () () () = ˆ + *<sup>q</sup> j k <sup>f</sup> q fq <sup>f</sup> y k s ke e k* w

( ) <sup>2</sup> /2 1

( / 2) ( / 2 1) ( / 2 1) ( ) , ,..., - -+ - = é ù ë û *qq q <sup>T</sup> j N j N jN*

=- - + - é ù ( / 2), ( / 2 1),..., ( / 2 1) ë û

*<sup>q</sup> ee e* ww

**a**

<sup>=</sup> å - *<sup>q</sup> <sup>N</sup> j k*

/ 2 <sup>ˆ</sup> () () <sup>ˆ</sup> -

*f q f fq q*

*z yk s k e*

=-

*k N*

with *f <sup>s</sup>* is the sampling frequency. Note, here *xf* is used to represent the DFT of *x*. Hence, the

^

w

{ } <sup>ˆ</sup> = - -+ - ˆˆ ˆ ( / 2), ( / 2 1),..., ( / 2 1) *fq fq fq fq* **<sup>S</sup>** *diag s N s N s N* (12)

 w

denotes the transpose operation. Then the cost function in equation (11) can be

*ff f <sup>f</sup>* **y** *y N y N yN* (14)

(13)

*T*

 a

a

*fq*(*k*) and *e <sup>f</sup>* 2(*k*) are the DFT of *y*(*n*), *s*

To estimate *ωq*, we define the following NLS cost function [12-15]

w

Suppose the time delay from the *<sup>q</sup> th* satellite to the receiver is *τq*, then the received signal can

 t

*yn s n en* a

$$
\hat{\tau}\_q = -\hat{\alpha}\_q N / (2\pi f\_s) \tag{17}
$$

Equation (16) can be solved by using FFT algorithm with a low computation burden. Since the PRN index of the received signal is unknown, consequently, it is not feasible to reconstruct *s* ^ *<sup>q</sup>*(*n*) directly by equation (8). Therefore, we should reconstruct signal of each satellite with the estimated Doppler frequency *ω* ^ *<sup>q</sup>*. Then the cross correlation of the reconstructed signal and the received one can be calculated. Furthermore, the corresponding time delay can be estimated from equation (16).

It is well-known that the cross correlation coefficients of different C/A codes are very small. In addition, the signal and noise are uncorrelated. Thus Fourier analysis results of the above mentioned component are close to zero. The maximum correlation value can only be obtained in the case that the C/A code of the reconstructed signal is the same as the received one. Therefore, the PRN index can be obtained. Furthermore, the code delay can be estimated by the location of the maximum correlation value. Then, we obtained all the unknown parameters.

#### **2.3. Comparison of acquisition algorithms**

Suppose the searching time for Doppler frequency estimation is *Q* in the conventional acquisition algorithm, *P* is the number of satellite received by the receiver. In GPS, the typical values of *Q* and *P* are 21 and 4-8, respectively, as in [9, 11]. The total number of satellites *M* is 32. Since Doppler frequency and code delay are unknown in the conventional acquisition algorithm, *Q* frequency points should be searched over certain frequency range. What's more, *M* satellites are also to be searched for each of the *Q* frequency points. While in the proposed algorithm, the Doppler frequencies of *P* satellites have been estimated, only *P* accurate frequencies need to be searched.

Figure 2 compares the searching process of the proposed algorithm and the conventional one. It can be seen that the computation complexity of the proposed algorithm is significantly lower than that of the conventional algorithm in that the typical value of *Q* is obviously greater than *P*.

Comparison of computation burden between new acquisition algorithm and conventional acquisition algorithm are shown Table 1. It can be clearly seen from Table 1 that, since *Q* is

the conventional algorithm in that the typical value of *Q* is obviously greater than *P* .

**Figure 2.** Comparison of searching process between new algorithm and conventional algorithm

always greater than *P*, the computation complexity of the proposed algorithm is lower than the conventional acquisition procedure. In practices, only 4 satellites are enough to determine the user position, which is less than 8, hence the computation burden should be further decreased. As the number of satellite is equal to 4, only half of the original computation complexity is required. Comparison of computation burden between new acquisition algorithm and conventional acquisition algorithm are shown Table 1. It can be clearly seen from Table 1 that, since *Q* is always greater than *P* , the computation complexity of the proposed algorithm is lower than the conventional acquisition procedure. In practices, only 4 satellites are enough to determine the user position, which is less than

Figure 2. Comparison of searching process between new algorithm and conventional algorithm


**Table 1.** Comparison of computation burden between new algorithm and conventional one

#### **2.4. Numerical results**

To verify the performance of the proposed algorithm, GPS signals from 8 satellites with PRN indexes 1, 2, 13, 20, 22, 24, 25 and 27 are simulated. Considering the IF is 1.25MHz, the sampling rate is 7.5MHz, the pre integration time is 1ms and the signal to noise ratio (SNR) is -20dB. The frequency searching step of the conventional acquisition algorithm is 1kHz. The acquisition results are shown in Figure 3. Chapter 1

(b) proposed acquisition algorithm

The term multipath is derived from the fact that a signal transmitted from a GNSS satellite can follow a 'multiple' number of propagation 'paths' to the receiving antenna. This is possible because the signal can be reflected back to the antenna from surrounding objects, including the earth's surface.

() () () *<sup>c</sup> <sup>j</sup> <sup>t</sup> x t dtcte*

where *d t* is the navigation data, *c t* is the C/A code in GPS. Here, we take the two path signal model as an example, which can be generalized to multiple reflection path signals directly. The

(18)

Figure 3 Compassion of the acquisition results **Figure 3.** Compassion of the acquisition results

**3.1 Multipath Data Model** 

**3 Two Novel Methods for Multipath Mitigation** 

Suppose the signal transmitted by GNSS satellites can be written as

received signal including the LOS signal with a single reflection can be written as

always greater than *P*, the computation complexity of the proposed algorithm is lower than the conventional acquisition procedure. In practices, only 4 satellites are enough to determine the user position, which is less than 8, hence the computation burden should be further decreased. As the number of satellite is equal to 4, only half of the original computation

Figure 2. Comparison of searching process between new algorithm and conventional algorithm Comparison of computation burden between new acquisition algorithm and conventional acquisition algorithm are shown Table 1. It can be clearly seen from Table 1 that, since *Q* is always greater than *P* , the computation complexity of the proposed algorithm is lower than the conventional acquisition procedure. In practices, only 4 satellites are enough to determine the user position, which is less than

**Figure 2.** Comparison of searching process between new algorithm and conventional algorithm

(b) searching process of conventional algorithm

Chapter 1

1 2 3 4 5 6 7 8

 

*d d d d d d d d*

(a) searching process of new algorithm 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

*dd d d d d d d d d d d d d d d d d d d d*

 

 

the conventional algorithm in that the typical value of *Q* is obviously greater than *P* .

*c c c c c c c c c c c*

0 31 32 *c c*

*c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c*  

 

*c c c c c c c c c c c c c c c c c c c c c*

32 Fourier Transform - Signal Processing and Physical Sciences

complexity is required.

In Figure 3, the horizontal axis stands for the PRN index of the satellites and the vertical axis denotes the acquisition metric. Since the proposed algorithm only calculates the correlation of the received satellite signal and the local reference signal, correlation results of the other satellite are not shown in Figure 3(b). It can be noted that acquisition results are identical, and the peak value of the proposed algorithm is higher than the conventional one. Therefore, the proposed algorithm has a comparative performance with the conventional one.

## **3. Two novel methods for multipath mitigation**

#### **3.1. Multipath data model**

The term multipath is derived from the fact that a signal transmitted from a GNSS satellite can follow a 'multiple' number of propagation 'paths' to the receiving antenna. This is possible because the signal can be reflected back to the antenna from surrounding objects, including the earth's surface.

Suppose the signal transmitted by GNSS satellites can be written as

$$\mathbf{x}(t) = d(t)\mathbf{c}(t)e^{\langle \alpha\_i t \rangle} \tag{18}$$

where *d*(*t*) is the navigation data, *c*(*t*) is the C/A code in GPS. Here, we take the two path signal model as an example, which can be generalized to multiple reflection path signals directly. The received signal including the LOS signal with a single reflection can be written as

$$\overline{\hat{y}}(t) = \sum\_{p=1}^{2} \alpha\_{p}^{\prime} d(t - \tau\_{p}) c(t - \tau\_{p}) e^{i a\_{i}(t - t\_{p})} + e(t) \tag{19}$$

where *α*<sup>1</sup> ′, *τ*<sup>1</sup> is amplitude and code delay of the LOS signal, *α*<sup>2</sup> ′, *τ*<sup>2</sup> is amplitude and code delay of the reflect signal, *e*(*t*) is the thermal noise. Equation (19) can also be deployed as

$$\begin{split} \overline{\mathbf{y}}(t) &= \sum\_{p=1}^{2} \alpha\_{p}^{\prime} d(t - \tau\_{p}) c(t - \tau\_{p}) e^{j a\_{t} (t - \tau\_{p})} + e(t) \\ &= \left( \alpha\_{1}^{\prime} d\left(t - \tau\_{1}\right) c\left(t - \tau\_{1}\right) e^{-j a\_{t} \tau\_{1}} + \alpha\_{2}^{\prime} d\left(t - \tau\_{2}\right) c\left(t - \tau\_{2}\right) e^{-j a\_{t} \tau\_{2}} \right) e^{j a\_{t} t} + e(t) \end{split} \tag{20}$$

To obtain the Doppler frequency of the LOS signal, we use the following relation

$$
\phi = \alpha \iota\_c \tau\_1 = \alpha \iota\_c \frac{R\_1(t)}{c} \tag{21}
$$

Assume the radial velocity of the receiver relative to the satellite is *v*1, then the range between them can be given by

$$R\_t(t) = R\_0 + \nu\_t t \tag{22}$$

where **R**<sup>0</sup> is the initial range between receiver and satellite, then equation (21) can be further given by

$$\boldsymbol{\phi} = \boldsymbol{\alpha}\_c \boldsymbol{\tau}\_1 = \boldsymbol{\alpha}\_c \frac{R\_0 + \mathbf{v}\_i t}{c} = 2\pi \left(\frac{R\_0}{\lambda} + \frac{\mathbf{v}\_i t}{\lambda}\right) \tag{23}$$

As the relative radial velocity is constant, the Doppler frequency can be obtained by

$$
\rho \alpha\_{d1} = \frac{d\phi}{dt} = \frac{\mathbf{v}\_1}{\mathcal{A}} \tag{24}
$$

The propagation range of the multipath signal can be written as

$$R\_2(t) = R\_0 + \nu\_\text{l} t + \Delta R(t) \tag{25}$$

where *Δ***R**(*t*) is the propagation range difference between the LOS signal and the reflected one. Since satellites are far away from the receiver, the range difference *Δ***R**(*t*) can be considered as constant in the short integration time. Then equation (21) can be represented as

$$
\rho \alpha\_c \tau\_p = \alpha\_c \frac{R\_0 + \mathbf{v}\_i t}{c} = 2\pi \left(\frac{R\_0}{\lambda} + \frac{\mathbf{v}\_i t}{\lambda} + \frac{\Delta R}{\lambda}\right) \tag{26}
$$

Further simplify of equation (20) we can get

In Figure 3, the horizontal axis stands for the PRN index of the satellites and the vertical axis denotes the acquisition metric. Since the proposed algorithm only calculates the correlation of the received satellite signal and the local reference signal, correlation results of the other satellite are not shown in Figure 3(b). It can be noted that acquisition results are identical, and the peak value of the proposed algorithm is higher than the conventional one. Therefore, the

The term multipath is derived from the fact that a signal transmitted from a GNSS satellite can follow a 'multiple' number of propagation 'paths' to the receiving antenna. This is possible because the signal can be reflected back to the antenna from surrounding objects, including

() () () = *<sup>c</sup> j t xt dtcte*

The received signal including the LOS signal with a single reflection can be written as

<sup>2</sup> ( )

() ( ) ( ) () -

 t


*R t c*

*d t ct e d t ct e e et*

*c c c*

 t wt

 w

(21)

*j j jt*

= -- + å ¢ *c p j t pp p*

*yt dt ct e et*

of the reflect signal, *e*(*t*) is the thermal noise. Equation (19) can also be deployed as

( ( ) ( ) ( ) ( ) ) 1 2

To obtain the Doppler frequency of the LOS signal, we use the following relation

 at

= -- +- - + ¢ ¢

<sup>1</sup> ( ) = = *c c* <sup>1</sup>


*j t*

w t

11 1 2 2 2

w t w

w t

where *d*(*t*) is the navigation data, *c*(*t*) is the C/A code in GPS. Here, we take the two path signal model as an example, which can be generalized to multiple reflection path signals directly.

(18)

(19)

′, *τ*<sup>2</sup> is amplitude and code delay

( )

(20)

proposed algorithm has a comparative performance with the conventional one.

**3. Two novel methods for multipath mitigation**

34 Fourier Transform - Signal Processing and Physical Sciences

Suppose the signal transmitted by GNSS satellites can be written as

1

at

′, *τ*<sup>1</sup> is amplitude and code delay of the LOS signal, *α*<sup>2</sup>

=

*p*

<sup>2</sup> ( )

() ( ) ( ) ()

 t  t

> f wt w

= -- + ¢

å *c p*

*yt dt ct e et*

*pp p*

at

at

1

=

*p*

**3.1. Multipath data model**

the earth's surface.

where *α*<sup>1</sup>

$$\overline{\mathbf{y}}(t) = \left(\alpha\_1' d\left(t - \tau\_1\right) c\left(t - \tau\_1\right) e^{-j\phi\_1} + \alpha\_2' d\left(t - \tau\_2\right) c\left(t - \tau\_2\right) e^{-j\left(\theta\_1 \ast \Lambda \phi\right)}\right) e^{j2\pi d/t} + e(t) \tag{27}$$

where *ω* ′ *<sup>d</sup>* =*ω<sup>d</sup>* <sup>1</sup> + *ωc* is the Doppler frequency of the received signal, *φ*1=**R**<sup>0</sup> / *λ* and *φ*2=*φ*<sup>1</sup> + *Δφ* are the phase of the LOS signal and the reflect one respectively, with *Δφ* an extra phase caused by the extra propagation range *Δ***R**. Accordingly, the LOS signal and reflect signal share the same Doppler shift.

After A/D conversion, the transformed digital signal can be written as

$$\overline{\mathbf{y}}(n) = \sum\_{p=1}^{2} \alpha'\_{p} d(n - \tau\_{p}) \mathbf{c}(n - \tau\_{p}) e^{-j\phi\_{p}} e^{j\alpha'\_{p} n} + e(n) \tag{28}$$

Assume the Doppler frequency has been estimated accurately, then the complex phase *<sup>e</sup> <sup>j</sup><sup>ω</sup>* ′ *d n* can be compensated after down-conversion. Hence, equation (28) can be written as

$$\mathbf{y}(n) = \sum\_{\rho=1}^{2} \alpha\_{\rho}^{\prime} d(n - \tau\_{\rho}) \mathbf{c}(n - \tau\_{\rho}) e^{-j\phi\_{\rho}} + \mathbf{z}^{\prime}(n) \tag{29}$$

where *e*˜(*n*) is equal to *e*(*n*)*<sup>e</sup>* <sup>−</sup> *<sup>j</sup><sup>ω</sup>* ′ *d n* that share the same statistics characteristic with *e*(*n*). So we still use *e*(*n*) to represent the noise data.

The in-phase component of equation (29) can be given by

$$\text{Cov}(n) = \sum\_{\rho=1}^{2} \alpha\_{\rho}^{\prime} d(n - \tau\_{\rho}) \mathbf{c}(n - \tau\_{\rho}) \cos \left(\varphi\_{\rho} - \varphi\_{c}\right) + e(n) \tag{30}$$

where *φc* is the synthetic phase of the LOS signal and the reflect one. Assume *d* is the early-late correlator spacing in the classical DLL, the estimated code delay of the direct signal is *τ* ^ <sup>1</sup>. Then the early and late code can be written as

$$\mathbf{s}\_E(t) = \mathbf{y}\_I(t - \hat{\mathbf{r}}\_1 - d / \mathcal{D}) \tag{31}$$

$$s\_L(t) = y\_I(t - \hat{\tau}\_1 + d / \, 2) \tag{32}$$

In classical DLL, the in-phase early and late correlation value is given as

$$R\_{\varepsilon}(\varepsilon) = \sum\_{p=1}^{2} \alpha\_{\rho} R(\varepsilon + \Delta \tau\_{\rho} - d \,/\ 2) \cos(\varphi\_{\rho} - \varphi\_{\varepsilon}) \tag{33}$$

$$R\_L(\boldsymbol{\varepsilon}) = \sum\_{p=1}^{2} \alpha\_p R(\boldsymbol{\varepsilon} + \Delta \boldsymbol{\tau}\_p + d \wedge \mathcal{D}) \cos(\boldsymbol{\varphi}\_p - \boldsymbol{\varphi}\_c) \tag{34}$$

where *ε* =*τ*1− *τ* ^ <sup>1</sup> is the estimation error of the LOS signal, namely tracking error. *Δτ* =*τ*2−*τ*1 is the relative delay of the multipath signal to the LOS. The discrimination function of the classical DLL (Delay Locked Loop) can be written as

**Figure 4.** Correlation peak when multipath is present

2 1

36 Fourier Transform - Signal Processing and Physical Sciences

=

2 1

=

*d n*

The in-phase component of equation (29) can be given by

1

*y n dn cn* a

 t

=

*p*

*p*

where *e*˜(*n*) is equal to *e*(*n*)*<sup>e</sup>* <sup>−</sup> *<sup>j</sup><sup>ω</sup>* ′

still use *e*(*n*) to represent the noise data.

I

the early and late code can be written as

at

*p*

( ) ( )( ) ( ) - ¢

 t

Assume the Doppler frequency has been estimated accurately, then the complex phase *<sup>e</sup> <sup>j</sup><sup>ω</sup>* ′

j w

> j

(28)

(29)

*e n* (30)

(31)

(32)

(33)

(34)

^ <sup>1</sup>. Then

that share the same statistics characteristic with *e*(*n*). So we

*d n*

= -- + å ¢ *<sup>p</sup> <sup>d</sup> <sup>j</sup> j n pp p*

*yn d n cn e e en*

can be compensated after down-conversion. Hence, equation (28) can be written as

( ) ( )( ) ( ) -

( ) <sup>2</sup>

( ) ( ) ( )cos ( )

 t

where *φc* is the synthetic phase of the LOS signal and the reflect one. Assume *d* is the early-late

 jj

= - - -+ å ¢ *p p p pc*

correlator spacing in the classical DLL, the estimated code delay of the direct signal is *τ*

<sup>1</sup> ( ) ( / 2) = -- ˆ *E I st yt d* t

<sup>1</sup> ( ) ( / 2) = -+ ˆ *L I st yt d* t

( ) ( / 2)cos( )

( ) ( / 2)cos( )

relative delay of the multipath signal to the LOS. The discrimination function of the classical

*L pp* = +D + å *p c* -

 t

 jj

 jj

<sup>1</sup> is the estimation error of the LOS signal, namely tracking error. *Δτ* =*τ*2−*τ*1 is the

*E pp* = +D - å *p c* -

 t

In classical DLL, the in-phase early and late correlation value is given as

*R Rd*

*R Rd*

2 1

=

2 1

=

 ae

*p*

 ae

*p*

e

e

DLL (Delay Locked Loop) can be written as

where *ε* =*τ*1− *τ*

^

 t

= -- + å ¢ *<sup>p</sup> <sup>j</sup> pp p*

*yn dn cn e e n*

at

$$\begin{split} D(\boldsymbol{\varepsilon}) &= R\_{\boldsymbol{\varepsilon}}(\boldsymbol{\varepsilon}) - R\_{\boldsymbol{\varepsilon}}(\boldsymbol{\varepsilon}) \\ &= \sum\_{\rho=1}^{2} \alpha\_{\rho} \Big[ R(\boldsymbol{\varepsilon} + \boldsymbol{\Lambda}\boldsymbol{\tau}\_{\rho} - \boldsymbol{d} \;/\, 2) - R(\boldsymbol{\varepsilon} + \boldsymbol{\Lambda}\boldsymbol{\tau}\_{\rho} + \boldsymbol{d} \;/\, 2) \Big] \cos(\rho\_{\rho} - \rho\_{\boldsymbol{c}}) \end{split} \tag{35}$$

In GNSS navigation, the receiver is concerned to maximize the correlation function between the received and locally generated signals. This can be accomplished by determining the locations of the zero output of the discriminator which corresponds to the maximum of the correlation function. Here, the early-late correlator is used to determine the position of this zero. However, the presence of multipath introduces some bias in the position of the first arrival peak and has an impact in the user's position, which can be clearly seen in Figure 4. Then the classical DLL failed to cope with multipath propagation, see [7-8].

#### **3.2. Code delay estimation**

#### *3.2.1. Code delay estimation in the correlation domain via NLS*

In this subsection, a code delay estimation algorithm based on DPE in correlation domain is proposed. To deploy the proposed algorithm, we combine the complex constant phase *e* <sup>−</sup> *<sup>j</sup>φ<sup>p</sup>* and the real amplitude *α* ′ *<sup>p</sup>* in equation (29) together as a new complex variable *α* ″ *<sup>p</sup>*. Then equation (29) can be written as

$$\mathbf{y}(n) = \sum\_{\rho=1}^{2} \alpha\_{\rho}'' d(n - \tau\_{\rho}) \mathbf{c}(n - \tau\_{\rho}) + \mathbf{e}(n) \tag{36}$$

The navigation data cycle is much longer than that of the C/A code, but only one cycle of the C/A code is used in the code delay estimation. Consequently, the navigation data jump can be neglected in the signal reconstruction. Then the navigation data +1 or −1 is written into *α* ″ *p* and denoted as *αp*. Therefore equation (36) can be written as

$$\mathbf{y}(n) = \sum\_{\rho=1}^{2} \alpha\_{\rho} \mathbf{c}(n - \tau\_{\rho}) + \mathbf{e}(n) \tag{37}$$

To simplify the following expression, *s*() is introduced to represent *c*(). Then equation (37) can be represented as

$$\mathbf{y}(n) = \sum\_{p=1}^{2} \alpha\_p \mathbf{s}(n - \tau\_p) + \mathbf{e}(n) \tag{38}$$

Multipath interference mitigation based on code delay estimation focus on the estimation of the LOS delay *τ*<sup>1</sup> and the reflect delay *τ*2, where the estimated result can be used to migrate the multipath. After acquisition, the correlation between the received satellite signal and the local reference signal is used to estimate the parameter of both the LOS signal and the reflected one. The proposed code delay estimation algorithm can take the place of the classical DLL, which can be described in details as follows.

To obtain the initial code delay, the correlation function between the received satellite signal and the reference signal is represented as

$$r(n) = \text{corr}\left(s(n), \mathbf{y}(n)\right) \tag{39}$$

where corr() represents the cross correlation operation. Since the noise is uncorrelated with the signal, equation (39) can be further expressed as

$$r(n) = \sum\_{p=1}^{2} \alpha\_p r^s (n - \tau\_p) + \mathcal{w}(n) \tag{40}$$

where *r <sup>s</sup>* (*n*) is the autocorrelation of *s*(*n*), *w*(*n*)=corr(*s*(*n*), *e*(*n*)) can still be thought as noise that sharing the same statistic characteristics with *e*(*n*). Apply DFT to equation (40)

Efficient FFT-based Algorithms for Multipath Interference Mitigation in GNSS http://dx.doi.org/10.5772/59940 39

$$p\_f(k) = r\_f^s(k) \sum\_{\rho=1}^2 \alpha\_\rho e^{j\alpha\_\rho k} + \omega\_f(k) \tag{41}$$

where *pf* (*k*), *rf s* (*k*), *wf* (*k*) are the DFT of *p*(*n*), *<sup>r</sup> <sup>s</sup>* (*n*) and *w*(*n*) respectively, *ω<sup>p</sup>* = −2*πτ<sup>p</sup>* / *NTs* with *Ts* is the sampling interval.

Define a NLS cost function as

2 1

38 Fourier Transform - Signal Processing and Physical Sciences

=

2 1

=

2 1

=

*p*

*p*

*p*

and denoted as *αp*. Therefore equation (36) can be written as

be represented as

where *r <sup>s</sup>*

can be described in details as follows.

and the reference signal is represented as

signal, equation (39) can be further expressed as

( ) ( )( ) ( )

() ( ) ()

() ( ) ()

= -+ å *p p*

*yn sn en* a

 t

To simplify the following expression, *s*() is introduced to represent *c*(). Then equation (37) can

 t

Multipath interference mitigation based on code delay estimation focus on the estimation of the LOS delay *τ*<sup>1</sup> and the reflect delay *τ*2, where the estimated result can be used to migrate the multipath. After acquisition, the correlation between the received satellite signal and the local reference signal is used to estimate the parameter of both the LOS signal and the reflected one. The proposed code delay estimation algorithm can take the place of the classical DLL, which

To obtain the initial code delay, the correlation function between the received satellite signal

where corr() represents the cross correlation operation. Since the noise is uncorrelated with the

 t

(*n*) is the autocorrelation of *s*(*n*), *w*(*n*)=corr(*s*(*n*), *e*(*n*)) can still be thought as noise that

() ( ) ()

= -+ å *<sup>s</sup> p p*

*rn r n wn* a

2 1

=

sharing the same statistic characteristics with *e*(*n*). Apply DFT to equation (40)

*p*

= -+ å *p p*

*yn cn en* a

The navigation data cycle is much longer than that of the C/A code, but only one cycle of the C/A code is used in the code delay estimation. Consequently, the navigation data jump can be neglected in the signal reconstruction. Then the navigation data +1 or −1 is written into *α* ″

 t

(36)

(37)

(38)

*rn sn yn* ( ) corr ( ), ( ) = ( ) (39)

(40)

*p*

= - -+ å ¢¢ *pp p*

*yn d n cn en* at

$$C\_1\left(\left\{\alpha\_p, \alpha\_p\right\}\_{p=1}^2\right) = \sum\_{k=-N/2}^{N/2} \left\| p\_f(k) - r\_f^s(k) \sum\_{\rho=1}^2 \alpha\_\rho e^{\prime \alpha\_\rho k} \right\|^2 \tag{42}$$

The unknown parameters {*ωp*, *αp*} *<sup>p</sup>*=1 <sup>2</sup> can be obtained by minimizing the cost function *C*1({*ωp*, *αp*} *<sup>p</sup>*=1 <sup>2</sup> ). Whereas, searching over a multidimensional space to solve the NLS problem requires higher computational complexity. To reduce the complexity, the Weighted Fourier Transform and RELAXation (WRELAX) algorithm (see [12]) is utilized here to solve the NLS problem. Let

$$\mathbf{a}(\alpha o\_{\rho}) = \left[ e^{j o\_{\rho}(-N/2)}, e^{j o\_{\rho}(-N/2 \times 1)}, \dots, e^{j o\_{\rho}(N/2 - 1)} \right]^{\mathrm{T}} \tag{43}$$

$$\mathbf{p}\_f = \left[ p\_f(-N/2), p\_f(-N/2+1), \dots, p\_f(N/2-1) \right]^\top \tag{44}$$

$$\mathbf{R}\_f = \text{diag}\left\{ r\_f^s(-N \;/\ \mathcal{D}), r\_f^s(-N \;/\ \mathcal{D}+1), \dots, r\_f^s(N \;/\ \mathcal{D}-1) \right\} \tag{45}$$

Consequently the cost function in equation (42) can be rewritten as

$$C\_1(\{o o\_p, \alpha\_p\}\_{p=1}^2) = \left\| \mathbf{p}\_f - \sum\_{p=1}^2 \alpha\_p \mathbf{R}\_f \mathbf{a}(o\_p) \right\|^2 \tag{46}$$

Further, we assume {*ω* ^ *<sup>q</sup>*, *α* ^ *<sup>q</sup>*}*q*=1,*q*<sup>≠</sup> *<sup>p</sup>* <sup>2</sup> have been estimated, then

$$\mathbf{p}\_{\boldsymbol{\beta}^{\boldsymbol{\rho}}} = \mathbf{p}\_{\boldsymbol{\beta}^{\boldsymbol{\epsilon}}} - \sum\_{q=1 \atop q \neq p}^{2} \hat{\alpha}\_{q} [\mathbf{R}\_{\boldsymbol{\beta}} \mathbf{a}(\hat{\boldsymbol{\alpha}}\_{p})] \tag{47}$$

Substitute equation (47) into equation (46) we have

$$C\_2(\left\{\alpha o\_\rho, \alpha\_\rho\right\}\_{p=1}^2) = \left\|\mathbf{p}\_{\not\!\! \rho} - \alpha\_\rho \mathbf{R}\_f \mathbf{a}(o\_\rho)\right\|^2 \tag{48}$$

By minimizing the cost function *C*2, estimation of *ω* ^ *<sup>p</sup>* and *<sup>α</sup>* ^ *<sup>p</sup>* can be obtained from

$$\hat{\boldsymbol{\phi}}\_{\boldsymbol{\rho}} = \arg\max\_{\boldsymbol{\phi}\_{\boldsymbol{\rho}}} \left| \mathbf{a}^{\mathrm{H}} (\boldsymbol{\phi}\_{\boldsymbol{\rho}}) \left( \mathbf{R}\_{f}^{\ast} \mathbf{p}\_{\boldsymbol{\rho}\boldsymbol{\rho}} \right) \right|^{2} \tag{49}$$

and

$$\hat{\alpha}\_{\rho} = \frac{\mathbf{a}^{\mathrm{H}}(o\nu\_{\rho}) \left(\mathbf{R}\_{\prime}^{\ast}\mathbf{p}\_{\rho\mathrm{e}}\right)}{\left\|\mathbf{R}\_{\prime}\right\|\_{F}^{2}}\Bigg|\_{o\nu\_{\rho} = \hat{\alpha}\_{\rho}}\tag{50}$$

where *F* represents the Frobenius norm. Hence, *ω* ^ *<sup>p</sup>* can be obtained as the location of the periodogram <sup>|</sup> **<sup>a</sup>***Η*(*ωp*)(**<sup>R</sup>** *<sup>f</sup>* **<sup>p</sup>** *fp*)| <sup>2</sup> , which can be efficiently computed by using FFT to the data sequence **R** *<sup>f</sup>* **p** *fp* padded with zeros. (Note the padding with zeros is necessary to determine *ω* ^ *<sup>p</sup>* with more accuracy.) Then *<sup>α</sup>* ^ *<sup>p</sup>* is easily computed from the complex height of the peak of **a***Η*(*ω* ^ *<sup>p</sup>*)(**<sup>R</sup>** *<sup>f</sup>* **<sup>p</sup>** *fp*) **R** *<sup>f</sup> <sup>F</sup>* <sup>2</sup> .

With the above simple preparations, we now proceed to present the correlation domain algorithm for the minimization of the NLS cost function. The proposed algorithm comprises the following steps.

Step 1. Assume *p* =1. Obtain {*ω* ^ 1, *α* ^ 1} from equation (49) and equation (50).

Step 2. Assume *p* =2 (the LOS signal and one reflected path). Compute **p** *<sup>f</sup>* <sup>2</sup> with (47) by using {*ω* ^ 1, *α* ^ 1} obtained in Step (1). Obtain {*ω* ^ 2, *α* ^ 2} from **p** *<sup>f</sup>* <sup>2</sup> described above.

Next, compute **p** *<sup>f</sup>* <sup>1</sup> with (47) by using {*ω* ^ 2, *α* ^ 2} and redetermine {*ω* ^ 1, *α* ^ 1} from **p** *<sup>f</sup>* <sup>1</sup> as in equation (49) and equation (50) above.

Iterate the previous two substeps until "practical convergence" is achieved then we can obtain {*ω* ^ *<sup>p</sup>*, *<sup>α</sup>* ^ *<sup>p</sup>*} *<sup>p</sup>*=1 <sup>2</sup> . Furthermore, by using *ω* ^ *<sup>p</sup>* <sup>=</sup> <sup>−</sup>2*πτ* ^ *<sup>p</sup>* / *NTs* the code delay *τ* ^ *<sup>p</sup>*, *p* =1, 2 can be obtained.

From the previous description one can find that the proposed algorithm can be implemented by simply FFT operation which leads to a less computation load.

The diagram of the novel code delay estimation algorithm is shown in Figure 5.

Efficient FFT-based Algorithms for Multipath Interference Mitigation in GNSS http://dx.doi.org/10.5772/59940 41

**Figure 5.** Diagram of the correlation domain code delay estimation algorithm

#### *3.2.2. Code delay estimation in the data domain via NLS*

Substitute equation (47) into equation (46) we have

40 Fourier Transform - Signal Processing and Physical Sciences

By minimizing the cost function *C*2, estimation of *ω*

and

*ω*

{*ω* ^ 1, *α* ^

{*ω* ^ *<sup>p</sup>*, *<sup>α</sup>* ^ *<sup>p</sup>*} *<sup>p</sup>*=1

**a***Η*(*ω*

^ *<sup>p</sup>*)(**<sup>R</sup>** *<sup>f</sup>* **<sup>p</sup>** *fp*) **R** *<sup>f</sup> <sup>F</sup>*

the following steps.

<sup>2</sup> .

w a

w

<sup>2</sup> <sup>2</sup>

 w

^ *<sup>p</sup>* and *<sup>α</sup>* ^

2

( )

ˆ

1} from equation (49) and equation (50).

2} and redetermine {*ω*

2} from **p** *<sup>f</sup>* <sup>2</sup> described above.

*<sup>p</sup>* / *NTs* the code delay *τ*

^ 1, *α* ^

^

1} from **p** *<sup>f</sup>* <sup>1</sup> as in equation

*<sup>p</sup>*, *p* =1, 2 can be obtained.

=

w w

*p p*

(48)

*<sup>p</sup>* can be obtained from

**a Rp** (49)

**<sup>R</sup>** (50)

, which can be efficiently computed by using FFT to the data

*<sup>p</sup>* is easily computed from the complex height of the peak of

^ *<sup>p</sup>* can be obtained as the location of the

2 1 ({ , } ) ( ) *C*

ˆ arg max ( ) H \* = *p p p f fp* w

H \*

**a Rp**

w

( ) <sup>ˆ</sup>

=

*p*

a

^

^ 1, *α* ^

by simply FFT operation which leads to a less computation load.

where *F* represents the Frobenius norm. Hence, *ω*

periodogram <sup>|</sup> **<sup>a</sup>***Η*(*ωp*)(**<sup>R</sup>** *<sup>f</sup>* **<sup>p</sup>** *fp*)| <sup>2</sup>

^ *<sup>p</sup>* with more accuracy.) Then *<sup>α</sup>*

Step 1. Assume *p* =1. Obtain {*ω*

(49) and equation (50) above.

1} obtained in Step (1). Obtain {*ω*

<sup>2</sup> . Furthermore, by using *ω*

Next, compute **p** *<sup>f</sup>* <sup>1</sup> with (47) by using {*ω*

*p pp*<sup>=</sup> = - **p Ra** *fp p f p* a

> w

( ) 2

sequence **R** *<sup>f</sup>* **p** *fp* padded with zeros. (Note the padding with zeros is necessary to determine

With the above simple preparations, we now proceed to present the correlation domain algorithm for the minimization of the NLS cost function. The proposed algorithm comprises

Step 2. Assume *p* =2 (the LOS signal and one reflected path). Compute **p** *<sup>f</sup>* <sup>2</sup> with (47) by using

Iterate the previous two substeps until "practical convergence" is achieved then we can obtain

From the previous description one can find that the proposed algorithm can be implemented

^ 2, *α* ^

> ^ 2, *α* ^

^ *<sup>p</sup>* <sup>=</sup> <sup>−</sup>2*πτ* ^

The diagram of the novel code delay estimation algorithm is shown in Figure 5.

*p f fp*

*<sup>f</sup> <sup>F</sup>*

Another DPE algorithm is proposed in this subsection. Different from the correlation domain algorithm, the unknown parameters {*αp*, *<sup>τ</sup>p*} *<sup>p</sup>*=1 <sup>2</sup> can also be estimated directly in the data domain, which can be described as follows.

Again, the signal after down-conversion is given by

$$\mathbf{y}(n) = \sum\_{\rho=1}^{2} \alpha\_{\rho} \mathbf{s}(n - \tau\_{\rho}) + \mathbf{e}(n) \tag{51}$$

To obtain the unknown parameters {*αp*, *τp*} *<sup>p</sup>*=1 <sup>2</sup> in equation (51), DFT is applied firstly

$$\mathbf{y}\_f(k) = \mathbf{s}\_f(k) \sum\_{p=1}^2 \alpha\_p e^{\beta \mathbf{e}\_p k} + \mathbf{e}\_f(k), -N/2 \le k \le N/2 - 1 \tag{52}$$

where *yf* (*k*), *sf* (*k*), *ef* (*k*) are the DFT of *y*(*n*), *s*(*n*) and *e*(*n*) respectively. *ω<sup>p</sup>* = −2*πτ<sup>p</sup> f <sup>s</sup>* / *N* with *f <sup>s</sup>* represents the sampling frequency. Till now, the code delay estimation is transformed to the angular frequency *ωp* estimation, after which *τp* can be obtained by the relation *τ<sup>p</sup>* = −*ωpN* / 2*π f <sup>s</sup>*.

To obtain (*α* ^ *<sup>p</sup>*, *ω* ^ *<sup>p</sup>*), a NLS cost function as follows is defined

$$\mathcal{Q}\_{i}\Big(\left\{\boldsymbol{\alpha}\_{p},\boldsymbol{\alpha}\_{p}\right\}\_{p=1}^{2}\Big) = \sum\_{k=-N/2}^{N/2-1} \left|\mathbf{y}\_{\boldsymbol{\cdot}}(k) - \mathbf{s}\_{\boldsymbol{\cdot}}(k)\sum\_{p=1}^{2} \boldsymbol{\alpha}\_{p}\boldsymbol{e}^{\boldsymbol{\cdot}\boldsymbol{\alpha}\_{p}k}\right|^{2}\tag{53}$$

Let

$$\mathbf{S}\_{f} = \text{diag}\left\{ \mathbf{s}\_{f}(-N/2), \mathbf{s}\_{f}(-N/2+1), \dots, \mathbf{s}\_{f}(N/2-1) \right\} \tag{54}$$

$$\mathbf{y}\_f = \left[ \mathbf{y}\_f(-N/2), \mathbf{y}\_f(-N/2+1), \dots, \mathbf{y}\_f(N/2-1) \right]^T \tag{55}$$

$$\mathbf{a}(\alpha o\_{\rho}) = \left[ e^{\beta o\_{\rho}(-N/2)}, e^{\beta o\_{\rho}(-N/2 + 1)}, \dots, e^{\beta o\_{\rho}(N/2 - 1)} \right]^T \tag{56}$$

Hence, minimizing the cost function in equation (53) is equivalent to minimizing the following cost function

$$\mathcal{Q}\_2\left(\left\{\boldsymbol{\alpha}\_{\rho}, \boldsymbol{\alpha}\_{\rho}\right\}\_{\rho=1}^2\right) = \left\|\mathbf{y}\_f - \sum\_{p=1}^2 \boldsymbol{\alpha}\_{\rho} \mathbf{S}\_f \mathbf{a}\left(\boldsymbol{\alpha}\_p\right)\right\|^2 \tag{57}$$

Assume {*α* ^ *<sup>q</sup>*, *ω* ^ *<sup>q</sup>*}*q*=1,*q*<sup>≠</sup> *<sup>p</sup>* <sup>2</sup> is known prior or has been estimated, then we have

$$\mathbf{y}\_{\mathcal{A}^p} = \mathbf{y}\_{\mathcal{f}} - \sum\_{q=1 \atop q \neq p}^2 \hat{\alpha}\_q \left[ \mathbf{S}\_{\mathcal{f}} \mathbf{a} \left( \hat{\alpha}\_p \right) \right] \tag{58}$$

Substitute equation (58) into equation (57)

$$\mathcal{Q}\_2\left(\boldsymbol{\alpha}\_{\boldsymbol{\rho}}, \boldsymbol{\alpha}\_{\boldsymbol{\rho}}\right) = \left\|\mathbf{y}\_{\boldsymbol{\rho}\boldsymbol{\rho}} - \boldsymbol{\alpha}\_{\boldsymbol{\rho}} \mathbf{S}\_f \mathbf{a}\left(\boldsymbol{\alpha}\_{\boldsymbol{\rho}}\right)\right\|^2 \tag{59}$$

Equation (59) gets the minimum while *y fp* =*αpS <sup>f</sup>* **a**(*ωp*), then the estimation results *ω* ^ *<sup>p</sup>* and *<sup>α</sup>* ^ *p* can be obtained in the following way

where *yf* (*k*), *sf* (*k*), *ef* (*k*) are the DFT of *y*(*n*), *s*(*n*) and *e*(*n*) respectively. *ω<sup>p</sup>* = −2*πτ<sup>p</sup> f <sup>s</sup>* / *N* with *f <sup>s</sup>* represents the sampling frequency. Till now, the code delay estimation is transformed to the angular frequency *ωp* estimation, after which *τp* can be obtained by the relation

({ } ) <sup>2</sup> /2 1 <sup>2</sup> <sup>2</sup>

<sup>=</sup> =- =

=- - + - é ù ( / 2), ( / 2 1),... ( / 2 1) ë û

( / 2) ( / 2 1) ( / 2 1) ( ) , ,..., - -+ - = é ù ë û *pp p <sup>T</sup> j N j N jN*

({ } ) ( )

1

( ) ( ) <sup>2</sup>

 a

= ¹ = - é ù å ë û *fp f qf p q q p* **y y Sa** a

**y Sa**

<sup>2</sup> is known prior or has been estimated, then we have

( ) <sup>2</sup>

ˆ ˆ

 w

> w

*p p fp p f p* , = - **y Sa** (59)

Hence, minimizing the cost function in equation (53) is equivalent to minimizing the following

<sup>2</sup> <sup>2</sup> <sup>2</sup>

a

 w

*<sup>p</sup> ee e* ww

**a**

<sup>2</sup> <sup>1</sup> <sup>1</sup> , <sup>=</sup> <sup>=</sup> *p p <sup>p</sup>* = -*f pf p* å*<sup>p</sup>*

*p p ff p <sup>p</sup> k N <sup>p</sup> Q yk sk e*

<sup>=</sup> å å - *<sup>p</sup> <sup>N</sup> j k*

w

*T*

(53)

(57)

(58)

a

**S** *ff f f* = - -+ - *diag s N s N s N* { ( / 2), ( / 2 1),..., ( / 2 1)} (54)

*ff f <sup>f</sup>* **y** *y N y N yN* (55)

(56)

 w

^ *<sup>p</sup>*), a NLS cost function as follows is defined

a w

42 Fourier Transform - Signal Processing and Physical Sciences

w

*Q* a w

Substitute equation (58) into equation (57)

2 *Q*

aw

<sup>1</sup> <sup>1</sup> / 2 <sup>1</sup> , () () -

*τ<sup>p</sup>* = −*ωpN* / 2*π f <sup>s</sup>*.

^ *<sup>p</sup>*, *ω*

To obtain (*α*

cost function

Assume {*α*

^ *<sup>q</sup>*, *ω* ^ *<sup>q</sup>*}*q*=1,*q*<sup>≠</sup> *<sup>p</sup>*

Let

$$\hat{\boldsymbol{\alpha}}\_{\rho} = \arg\max\_{\boldsymbol{\alpha}\_{\rho}} \left| \mathbf{a}^{\text{H}} (\boldsymbol{\alpha}\_{\rho}) \left( \mathbf{S}\_{\prime}^{\ast} \mathbf{y}\_{\prime \rho} \right) \right|^{2} \tag{60}$$

$$\hat{\alpha}\_{\rho} = \frac{\mathbf{a}^{\mathrm{H}} \left( \alpha\_{\rho} \right) \left( \mathbf{S}\_{\prime} \right) \mathbf{s}\_{\ \rho \ \rho}}{\left\| \mathbf{S}\_{\prime} \right\|\_{F}^{2}} \Bigg|\_{\alpha\_{\rho} - \hat{\alpha}\_{\rho}} \tag{61}$$

From equation (60), *ω* ^ *<sup>p</sup>* can be obtained as the location of the periodogram <sup>|</sup> **<sup>a</sup>***<sup>H</sup>* (*ω*)(*<sup>S</sup> <sup>f</sup>* \* *<sup>y</sup> fp*)| <sup>2</sup> , which can be efficiently computed by using the FFT with the data sequence *S <sup>f</sup>* \* *<sup>y</sup> fp* padded with zeros. (Note the padding with zeros is necessary to determine *ω* ^ *<sup>p</sup>* with more accuracy.) Then *α* ^ *p* is easily computed from the complex height of the peak of **a***<sup>H</sup>* (*ωp*)(*<sup>S</sup> <sup>f</sup>* \* *<sup>y</sup> fp*) *S <sup>f</sup> <sup>F</sup>* <sup>2</sup> .

The diagram of the data domain code delay estimation algorithm is shown as in Figure 6.

**Figure 6.** Diagram of the data domain code delay estimation algorithm

With the above simple preparations, we now proceed to present the relaxation algorithm for the minimization of the nonlinear least-squares cost function. The data domain WRELAX algorithm comprises the following steps.

Step 1. Assume *p* =1. Obtain {*ω* ^ 1, *α* ^ 1} from equation (61) and equation (60).

Step 2. Assume *p* =2 (the LOS signal and one reflected path). Compute *y <sup>f</sup>* <sup>2</sup> with (47) by using {*ω* ^ 1, *α* ^ 1} obtained in Step (1). Obtain {*ω* ^ 2, *α* ^ 2} from *y <sup>f</sup>* <sup>2</sup> described above.

Next, compute *y <sup>f</sup>* <sup>1</sup> with equation (47) by using {*ω* ^ 2, *α* ^ 2} and redetermine {*ω* ^ 1, *α* ^ 1} from *y <sup>f</sup>* <sup>1</sup> as in equation (49) and equation (50) above.

Iterate the previous two substeps until "practical convergence" is achieved then we can obtain {*ω* ^ *<sup>p</sup>*, *<sup>α</sup>* ^ *<sup>p</sup>*} *<sup>p</sup>*=1 <sup>2</sup> . Furthermore, by using *ω* ^ *<sup>p</sup>* <sup>=</sup> <sup>−</sup>2*πτ* ^ *<sup>p</sup>* / *NTs* the code delay *τ* ^ *<sup>p</sup>*, *p* =1, 2 can be obtained.

From the previous description one can find that the proposed algorithm can be implemented by simply FFT operation which deserves a less computation load.

#### *3.2.3. Comparison of the above two algorithms*

To further analyze the two proposed algorithms, the cost functions of them are further discussed in this subsection. Following the discussion in section 3.2.1 and section 3.2.2, we define

$$\mathbf{s}\_f = \left[ \mathbf{s}\_f \left( -N \,/\, \mathbf{2} \right), \mathbf{s}\_f \left( -N \,/\, \mathbf{2} + 1 \right), \dots, \mathbf{s}\_f \left( N \,/\, \mathbf{2} - 1 \right) \right]^\top \tag{62}$$

Then

$$\mathbf{R}\_{\boldsymbol{f}} = \mathbf{s}\_{\boldsymbol{f}} \mathbf{s}\_{\boldsymbol{f}}^H \tag{63}$$

In the noise free situation we have

$$\mathbf{p}\_{\rho\rho} = \alpha\_{\rho} \mathbf{R}\_{\rho} \mathbf{a} \begin{pmatrix} \alpha\_{\rho} \\ \end{pmatrix} \tag{64}$$

Then the correlation domain cost function in equation (49) can be further decomposed as

$$\begin{split} \hat{\boldsymbol{\alpha}}\_{\rho} &= \arg \max\_{\boldsymbol{\alpha}\_{\rho}} \left| \mathbf{a}^{\mathrm{H}} (\boldsymbol{\alpha}\_{\rho}) \mathbf{R}\_{\boldsymbol{\beta}} \, \mathbf{\dot{p}}\_{\rho} \right|^{2} \\ &= \arg \max\_{\boldsymbol{\alpha}\_{\rho}} \left| \mathbf{a}^{\mathrm{H}} (\boldsymbol{\alpha}\_{\rho}) \mathbf{R}\_{\boldsymbol{\beta}} \, \mathbf{\dot{s}} \left( \boldsymbol{\alpha}\_{\rho} \mathbf{R}\_{\boldsymbol{\beta}} \mathbf{a} (\boldsymbol{\alpha}\_{\rho}) \right) \right|^{2} \\ &= \arg \max\_{\boldsymbol{\alpha}\_{\rho}} \left| \mathbf{a}^{\mathrm{H}} (\boldsymbol{\alpha}\_{\rho}) \left( \mathbf{s}\_{f} \times \mathbf{s}\_{f} \, \mathbf{u} \right) \right| \left( \boldsymbol{\alpha}\_{\rho} \left( \mathbf{s}\_{f} \times \mathbf{s}\_{f} \, \mathbf{u} \right) \mathbf{a} (\boldsymbol{\alpha}\_{\rho}) \right) \right|^{2} \\ &= \arg \max\_{\boldsymbol{\alpha}\_{\rho}} \left| \mathbf{a}^{\mathrm{H}} (\boldsymbol{\alpha}\_{\rho}) \right| \left( \mathbf{s}\_{f}^{\mathrm{T}} \times \mathbf{s}\_{f} \, \mathbf{s} \right) \times \left( \mathbf{s}\_{f} \times \mathbf{s}\_{f}^{\mathrm{H}} \right) \left| \left( \boldsymbol{\alpha}\_{\rho} \mathbf{a} (\boldsymbol{\alpha}\_{\rho}) \right) \right|^{2} \end{split} \tag{65}$$

And in the same situation the data domain cost function in equation (60) can be deployed in the following way

With the above simple preparations, we now proceed to present the relaxation algorithm for the minimization of the nonlinear least-squares cost function. The data domain WRELAX

Step 2. Assume *p* =2 (the LOS signal and one reflected path). Compute *y <sup>f</sup>* <sup>2</sup> with (47) by using

Iterate the previous two substeps until "practical convergence" is achieved then we can obtain

From the previous description one can find that the proposed algorithm can be implemented

To further analyze the two proposed algorithms, the cost functions of them are further discussed in this subsection. Following the discussion in section 3.2.1 and section 3.2.2, we

( ) ( ) ( ) <sup>T</sup>

**p Ra** *fp p f p* = a

Then the correlation domain cost function in equation (49) can be further decomposed as

2

arg max ( ) ( )

**a ss ss a**

**a ss ss a**

*p fs p fs p*

 aw

arg max ( ) ( )

arg max ( ) ( )

\* <sup>H</sup>

**a R Ra**

= ´´

H \*

<sup>=</sup> é ù ´ ´´ ë û

H \*

**a Rp**

H \*

w

w

w  w

( )

( ) ( ( ) )

*p f f pf f p*

*p ff ff pp*

2

<sup>2</sup> H H

aw

a w

( ) ( ) ( )

<sup>2</sup> T H

^ 2, *α* ^

^ 2, *α* ^

^ *<sup>p</sup>* <sup>=</sup> <sup>−</sup>2*πτ* ^

by simply FFT operation which deserves a less computation load.

1} from equation (61) and equation (60).

2} from *y <sup>f</sup>* <sup>2</sup> described above.

*<sup>p</sup>* / *NTs* the code delay *τ*

=- - + - é ù / 2 , / 2 1 ,..., / 2 1 *ff f* ë û *<sup>f</sup>* **<sup>s</sup>** *s N s N sN* (62)

2} and redetermine {*ω*

^

<sup>H</sup> **R ss** *f ff* = (63)

( ) (64)

^ 1, *α* ^

*<sup>p</sup>*, *p* =1, 2 can be obtained.

1} from *y <sup>f</sup>* <sup>1</sup> as

(65)

algorithm comprises the following steps.

44 Fourier Transform - Signal Processing and Physical Sciences

1} obtained in Step (1). Obtain {*ω*

<sup>2</sup> . Furthermore, by using *ω*

*3.2.3. Comparison of the above two algorithms*

In the noise free situation we have

ˆ arg max ( )

w

*p*

*p p fs fp*

 w

*p*

w

*p*

w

*p*

w

=

w

=

in equation (49) and equation (50) above.

Next, compute *y <sup>f</sup>* <sup>1</sup> with equation (47) by using {*ω*

^ 1, *α* ^

Step 1. Assume *p* =1. Obtain {*ω*

{*ω* ^ 1, *α* ^

{*ω* ^ *<sup>p</sup>*, *<sup>α</sup>* ^ *<sup>p</sup>*} *<sup>p</sup>*=1

define

Then

( ) ( ) ( ( )) ( )( )( ( )) <sup>2</sup> H \* <sup>2</sup> H \* <sup>2</sup> H \*T ˆ arg max arg max arg max = = = *p p p p p f fp p f pf p p ff p p* w w w w w waw w aw **a Sy a S Sa a ss a** (66)

It has been discussed in previous subsection that *ω* ^ *<sup>p</sup>* can be obtained by using the FFT with the data sequence **R** *<sup>f</sup>* **p** *fp* and *S <sup>f</sup> <sup>H</sup> y fp*, respectively. Note (*αp***a**(*ωp*)) is an impulse signal, and FFT(*s<sup>f</sup>* \* *sf* T) convolutions with the impulse signal is still FFT(*s<sup>f</sup>* \* *sf* T) but with a code delay, which is the parameter to be estimated. From equation (65) and (66) we realized that the cost function are the FFT of *αp***a**(*ωp*) weighted by different window function, where window function in the correlation domain is *w***1**=FFT((*s<sup>f</sup>* T*sf* \*)(*s<sup>f</sup> <sup>s</sup><sup>f</sup>* H)), but the window function in data domain is *w*2=FFT(*s<sup>f</sup>* \* *sf* T). It's clear that the weighted window *w*1, *<sup>w</sup>*2 will broaden the cost function, and the impact of *w*<sup>1</sup> is more serious than that of *w*2, which can also be seen from the comparison of the cost function given in Figure 7. From this point of view, the proposed data domain algorithm would deserve a more accurate estimation result.

**Figure 7.** Comparison of the correlation domain and the data domain cost function

Correlation of the received data and the locally generated signal and the auto-correlation of the locally generated signal should be calculated in correlation domain WRELAX algorithm. Suppose the iteration times of the two proposed algorithms are the same, the correlation domain WRELAX algorithm is more complex than the data domain WRELAX algorithm.

#### *3.2.4. Numerical results*

To investigate the performance of the proposed algorithm in the presence of multipath, a simulation experiment was performed. The code delay estimation error can be represented as the function of *τ*, and as explained in [l] and [7], maximum and minimum errors occur when the multipath signal is in-phase *Δφ* =0 °or out-of-phase *Δφ* = ± 180 °with the LOS signal. The curve of the maximum and minimum errors is regarded as the error envelopes, which is used in the following numerical experiment to evaluate the performance. Chapter 1 reflected signal is 2 1 0.5 . The code delay difference 2 1 is varied from 0 to 1.5 chips. The error is calculated at the maximum points when the multipath signal is at 0°in phase with respect to the LOS signal. The results are shown in Figure 8(a).

In Figure 8, the curve 'classical DLL' denotes the error envelope of the classical DLL, and the curve 'narrow correlator' denotes the error envelope of the narrow correlator with spacing 0.1 *d* chip. Both the 'classical DLL' and the 'narrow correlator' were calculated by taking the DLL equations and solving for the tracking error. The curve 'MEDLL', 'Correlation Domain' and 'Data Domain' denote the code delay estimation error of MEDLL, proposed correlation domain algorithm and data domain algorithm, respectively. The results indicate that in the presence of multipath, the classical DLL technique has a large tracking error, which fails in multipath mitigation. The tracking error is reduced

Figure 8 Comparison of the error envelopes **Figure 8.** Comparison of the error envelopes

Suppose there is only one reflected multipath signal combined with the LOS signal. The noise and the receiver bandwidth are not considered in the experiment. The relative amplitude of the direct and the reflected signal is *α*<sup>2</sup> / *α*1=0.5. The code delay difference *Δτ* =*τ*2−*τ*1 is varied from 0 to 1.5 chips. The error is calculated at the maximum points when the multipath signal is at 0°in phase with respect to the LOS signal. The results are shown in Figure 8(a).

In Figure 8, the curve 'classical DLL' denotes the error envelope of the classical DLL, and the curve 'narrow correlator' denotes the error envelope of the narrow correlator with spacing *d* =0.1 chip. Both the 'classical DLL' and the 'narrow correlator' were calculated by taking the DLL equations and solving for the tracking error. The curve 'MEDLL', 'Correlation Domain' and 'Data Domain' denote the code delay estimation error of MEDLL, proposed correlation domain algorithm and data domain algorithm, respectively. The results indicate that in the presence of multipath, the classical DLL technique has a large tracking error, which fails in multipath mitigation. The tracking error is reduced in the narrow correlator. However, based on the previous conclusion, the tracking error is nearly constant when the correlation spacing is less than 0.1 chip. The tracking accuracy of the algorithms in this section and MEDLL show a considerable improvement over the narrow correlator technique and the classical DLL technique. The improved performance can be of great help in critical GPS applications where the multipath errors associated with using conventional receivers can easily exceed the accuracy requirements. To further compare the performance of the two proposed algorithms and MEDLL, we magnify the results in Figure 8(a) to obtain Figure 8(b). It is clear that the proposed two algorithms have a more favorable performance than MEDLL. Furthermore, the proposed data domain algorithm is superior to the correlation domain algorithm, which is consistent with the conclusion in section 3.2.3.

## **4. Conclusion**

Suppose the iteration times of the two proposed algorithms are the same, the correlation domain WRELAX algorithm is more complex than the data domain WRELAX algorithm.

To investigate the performance of the proposed algorithm in the presence of multipath, a simulation experiment was performed. The code delay estimation error can be represented as the function of *τ*, and as explained in [l] and [7], maximum and minimum errors occur when the multipath signal is in-phase *Δφ* =0 °or out-of-phase *Δφ* = ± 180 °with the LOS signal. The curve of the maximum and minimum errors is regarded as the error envelopes, which is used

0.5 . The code delay difference 2 1

Chapter 1

The error is calculated at the maximum points when the multipath signal is at 0°in phase with

(b). details of (a)

Figure 8 Comparison of the error envelopes

In Figure 8, the curve 'classical DLL' denotes the error envelope of the classical DLL, and the curve 'narrow correlator' denotes the error envelope of the narrow correlator with spacing 0.1 *d* chip. Both the 'classical DLL' and the 'narrow correlator' were calculated by taking the DLL equations and solving for the tracking error. The curve 'MEDLL', 'Correlation Domain' and 'Data Domain' denote the code delay estimation error of MEDLL, proposed correlation domain algorithm and data domain algorithm, respectively. The results indicate that in the presence of multipath, the classical DLL technique has a large tracking error, which fails in multipath mitigation. The tracking error is reduced

 

(a). error envelopes

is varied from 0 to 1.5 chips.

in the following numerical experiment to evaluate the performance.

respect to the LOS signal. The results are shown in Figure 8(a).

*3.2.4. Numerical results*

reflected signal is 2 1

**Figure 8.** Comparison of the error envelopes

 

46 Fourier Transform - Signal Processing and Physical Sciences

A novel acquisition algorithm is firstly proposed in this chapter. In the proposed algorithm, the Doppler frequency is obtained by utilizing FFT of the squared data. Then, the PRN index and initial code phases of the satellite are obtained based on the NLS criterion. It can be seen that the proposed algorithm can not only reduce the computation complexity but also attain a comparable performance to the conventional acquisition algorithm.

After that, two algorithms are presented to suppress the multipath interference by estimating the code delay. In the proposed algorithms, the code delay was obtained by solving a NLS equation, which can be further realized by FFT operation. Compared with the conventional estimation algorithm, the proposed two algorithms perform better in multipath propagating environments and bear lower computation burden.

## **Acknowledgements**

The work of this chapter is supported by the Project of the National Natural Science Foundation of China (Grant No. 61179064, 61172112, 61271404 and 61471363), and the Fundamental Research Funds for the Central Universities (Grant 3122014D008 and 3122014B001).

## **Author details**

Renbiao Wu1 , Qiongqiong Jia1 , Wenyi Wang1 and Jie Li2

1 Tianjin Key Lab for Advanced Signal Processing, Civil Aviation University of China, Tian‐ jin, P. R. China

2 Department of Electronic and Information Engineering, Zhonghuan Information College Tianjin University of Technology, Tianjin, P. R. China

## **References**


[10] Chengjun Li, Mingquan Lu, Zhenming Feng, Qi Zhang. Study on GPS Acquisition Algorithm and Performance Analysis. Journal of Electronics & Information Technol‐ ogy 2010; 32(2) 296-300.

**Author details**

, Qiongqiong Jia1

48 Fourier Transform - Signal Processing and Physical Sciences

2006; 42(2) 548-561.

57(1): 33-46

37(1) 183-195.

2999

Tianjin University of Technology, Tianjin, P. R. China

15-17 September 2010, Zadar, Croatia. 399-402.

PROPAGATION 2014; 61(5) 2775-2782.

tion Society 2010; 52(1) 104-113.

, Wenyi Wang1

and Jie Li2

1 Tianjin Key Lab for Advanced Signal Processing, Civil Aviation University of China, Tian‐

2 Department of Electronic and Information Engineering, Zhonghuan Information College

[1] Kalyanaraman S K, Braasch M S, Kelly J M. Code Tracking Architecture Influence on GPS Carrier Multipath. IEEE Transactions on Aerospace and Electronic Systems

[2] Kos T, Markezic I, Pokrajcic J. Effects of Multipath Reception on GPS Positioning Per‐ formance. 52nd International Symposium ELMAR-2010: conference proceedings,

[3] Scire-Scappuzzo F, Makarov S N. A Low-Multipath Wideband GPS Antenna With Cutoff or Non-Cutoff Corrugated Ground Plane. Antennas and Propagation 2009;

[4] Maqsood M, Gao S, Brown T, et al. A Compact Multipath Mitigating Ground Plane for Multiband GNSS Antennas. IEEE TRANSACTIONS ON ANTENNAS AND

[5] Ray J K, Cannon M E, Fenton P. GPS Code and Carrier Multipath Mitigation Using a Multiantenna System. IEEE Transactions on Aerospace and Electronic Systems 2001;

[6] Alfred R L. GPS Landing System Reference Antenna. IEEE Antennas and Propaga‐

[7] Michael S B. GPS Multipath Model Validation. IEEE Position Location and Naviga‐ tion Symposium 1996: conference proceedings, April 1996, Atlanta, GA. 672-678.

[8] Cannon M E, Lachapelle G, Qiu W, et al. Performance Analysis of a Narrow Correla‐ tor Spacing Receiver for Precise Static GPS Positioning. IEEE Position Location and Navigation Symposium 1994: conference proceedings, April 1994, Las Vegas. 2994—

[9] Tung Hai Ta, Son Hong Ngo. A Novel Signal Acquisition Method for GPS Dual-Fre‐ quency Ll C/A and L2C Receivers. International Conference on Advanced Technolo‐ gies for Communications (ATC 2011) : conference proceedings, Aug 2011. 231-234.

Renbiao Wu1

jin, P. R. China

**References**

