2.1. Signal model

BeiDou open service signals consist of the B1I and B2I signals. This paper mainly discusses the B1I and B2I signals broadcast by MEO/IGSO satellites which are modulated with NH code. These signals are composed of carrier, ranging code, NH code and navigation message. The ranging code period is one millisecond and the duration of one navigation message bit is 20 ms. The bit duration of NH code (0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0) is 1 ms with a length of 20 bits and a rate of 1 kbps. It is modulated on the ranging code synchronously with the navigation message bit (China Satellite Navigation Office, 2016).

The input signal of a receiver is obtained by down-converting, sampling and quantisation and is generally an Intermediate Frequency (IF) signal. The input IF signal can be written as follows:

(1)$$s_{IF}(n) = AD(nT_s)H(nT_s)c(nT_s-\tau _0)\exp [j\cdot 2\pi ( f_{IF} + f_d)nT_s + \phi _0]$$

where A, D and H represent the carrier amplitude, the navigation data and the NH code, T _s is sampling time, c(·) stands for the ranging code sequence, τ₀, ϕ ₀, f _IF and f _d are code propagation delay, initial carrier phase, the carrier IF and Doppler shift.

Signal acquisition is a correlation operation between the locally generated IF signals and the input IF signals. The locally generated signal sequence is:

(2)$$s_L(n) = c(nT_s-\hat{\tau }_L)\exp [j\cdot 2\pi (f_{IF} + \hat{f}_{d,L})nT_s]$$

Assuming the CC integration time is T _c, the output of the k-th integration is:

(3)$$Y(k) = \displaystyle{1 \over {N_c}}\sum\limits_{n = kN_c}^{(k + 1)N_c-1} {s_{IF}} (n)\cdot s_L(n) = 0\cdot 5AD_kH_kR(\delta \tau ){\rm sinc}[\delta f_d(k)T_c]\exp [j\cdot \phi _k + j\pi \delta f_d(k)T_c]$$

where k represents the index of the integration interval, N _c = T _c/T _s is the number of samples in the integration time T _c, D _k = ±1 is navigation data, H _k = ±1 is NH code, R(·) is the autocorrelation function of the ranging code, δτ and δ f _d are the errors of code delay estimation and the Doppler shift and ϕ _k is the initial carrier phase error at the k-th integration interval.

2.2. Effect analysis of NH code

In conventional 1 ms coherent integration, when the integration process crosses data bit reversions, the integration results are attenuated.

Suppose that bit reversions occur in N _v (0 ≤ N _v ≤ N _c), namely:

(4)$$D_kH_k\left\{ {\matrix{ 1 \hfill & {for\quad 0 < k < N_v} \hfill \cr {-1} \hfill & {for\quad N_v < k < N_c-1} \hfill \cr } } \right.$$

The CC integration in Equation (3) can be expressed as:

(5)$$\eqalign{Y(k) = &\displaystyle{1 \over {N_c}}\sum\limits_{n = 0}^{N_c-1} {s_{IF}} (n)\cdot s_L(n) \cr =& \displaystyle{1 \over {N_c}}\left[ {\sum\limits_{n = 0}^{N_v-1} {s_{IF}} (n)\cdot s_L(n) + \sum\limits_{n = N_v}^{N_c-1} {s_{IF}} (n)\cdot s_L(n)} \right] \cr =& \;0\cdot 5AR(\delta \tau )\{ {\rm sinc}[\pi \delta f_d(k)N_v]\exp [ j\cdot \phi _k + j\pi \delta f_d(k)(N_v-1)] \cr & -{\rm sinc}[\pi \delta f_d(k)(N_c-N_v)]\exp [ j\cdot \phi _k + j\pi \delta f_d(k)(N_c + N_v-1)]\} }$$

When N _v = N _c/2, Y(k) can be rewritten as:

(6)$$\eqalign{Y(k) = &\; 0\cdot 5AR(\delta \tau ){\rm sinc}[0\cdot 5\pi \delta f_d(k)N_c] \cr & \times \{ \exp [ j\cdot \phi _k + j\pi \delta f_d(k)(0\cdot 5N_c-1)]-\exp [ j\cdot \phi _k + j\pi \delta f_d(k)(1\cdot 5N_c-1)]\} \cr =&\; 0\cdot 5AR(\delta \tau ){\rm sinc}[0\cdot 5\pi \delta f_d(k)N_c] \cr & \times \exp (j\cdot \phi _k)\exp [ j\pi \delta f_d(k)(0\cdot 5N_c-1)]\{ 1-\exp [ j\pi \delta f_d(k)N_c]\} }$$

Then, the decision variable of CC integration is determined as:

(7)$$\eqalign{\vert Y(k)\vert ^2 =& \;\{ 0\cdot 5AR(\delta \tau ){\rm sinc}[0\cdot 5\pi \delta f_d(k)N_c]\exp (j\cdot \phi _k)\exp [ j\pi \delta f_d(k)(0\cdot 5N_c-1)]\} ^2 \cr & \times \{ 1-\exp [ j\pi \delta f_d(k)N_c]\} ^2 \cr = &\; 0\cdot 25A^2\{ R(\delta \tau ){\rm sinc}[0\cdot 5\pi \delta f_d(k)N_c]\} ^2 \cr & \times \{ \exp (j\cdot \phi _k)\exp [ j\pi \delta f_d(k)(0\cdot 5N_c-1)]\} ^2\cdot 2\{ 1-\cos (\pi \delta f_d(k)N_c)\} }$$

When δ f _d is equal to zero, {1 − cos(πδf _d(k) N _c)} is equal to zero too. When δ f _d is equal to 1/2 N _c, which is 500 Hz, {1 − cos(πδf _d(k) N _c)} will reach its maximum. Therefore, the peak in the frequency axis splits into two components and is no longer consistent with the right Doppler bin.

Figure 1 presents the accumulation results of a 1 ms coherent combination. The frequency bin width is 500 Hz and the bin width of the code dimension is 0·5 chips. Figure 1 shows that the peak in the frequency axis splits into two components when the integration process crosses bit transitions. The corresponding Doppler bin of the peak is 500 Hz and indicates a failed acquisition. Therefore, for BeiDou signals, the data bit transitions will bring in a large frequency error because of NH codes.

Figure 1. Accumulation results of 1 ms coherent combination.

2.3. Conventional acquisition algorithm of BeiDou signal

In every 6 ms segment of the NH code (0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0), there must be two consecutive bits which have the same sign. Therefore, the 1 ms CC method is repeated consecutively six times and the maximum peak magnitude of these six consecutive results can be detected to estimate acquisition parameters.

The zero-padding method is another effective strategy to eliminate the effect of NH code. This method involves three steps:

1. (1) An incoming signal of 2 ms is correlated with a local carrier in the frequency dimension with the frequency bin of 500 Hz for Doppler removal.

2. (2) The 1 ms locally generated code sequence is appended by 1 ms zeros to form the extended replica.

3. (3) The incoming signal after Doppler removal is correlated with the replica. Fast-Fourier Transformation (FFT) and Inverse FFT (IFFT) can be used here. The first 1 ms of the correlation result is used for peak detection. The location along the time axis is the millisecond boundary and the location along the Doppler axis is an estimate of the Doppler frequency.

In these two methods above, the effective integration time is only 1 ms and the acquisition success rate of 1 ms ranging code is very low for weak signals. Therefore, a long integration time is needed to increase the acquisition sensitivity.

2.4. A long NH-PRN code acquisition algorithm

To increase the acquisition sensitivity, first a long NH-PRN code is modified to extend the integration time. The long NH-PRN modulated code is formed by multiplying NH code with the locally generated PRN code. The locally generated signal sequence is:

(8)$$s_{{\rm NH - PRN}}(n) = H(nT_s-\hat{\psi }_L)c(nT_s-\hat{\tau }_L)\exp [j\cdot 2\pi (f_{IF} + \hat{f}_{d,L})nT_s]$$

Here, the equivalent CC integration time T _coh is 20 ms, the integration result is:

(9)$$Y_{NH{\rm 0}PRN}(k) = 0\cdot 5AD_kH(\delta \varphi )R(\delta \tau ){\rm sinc}[\delta f_d(k)T_{coh}]\exp [j\cdot \phi _k + j\pi \delta f_d(k)T_{coh}]$$

where δφ is the error of NH code delay estimation.

Here, sinc[δf _d(k) T _coh] is the amplitude attenuation of integration because of the frequency deviation and T _coh is 20 ms. Figure 2 shows that the amplitude shall be attenuated to zero with a frequency deviation of an integral multiple of 50 Hz. When the frequency deviation is 25 Hz, the amplitude attenuation 10*log 10(sinc[δf _d(k) T _coh]) shall be −1·9 dB. Therefore, the frequency bin is chosen to be 25 Hz and then the search time is 400 times for a frequency searching range of 10 kHz (Han et al., Reference Han, Liu, Li, Zeng and Wang2017).

Figure 2. Amplitude attenuation of integration due to frequency deviation. The integration time T _coh is 20 ms.

3.1. Differential coherent algorithm

The differential coherent algorithm is implemented on the coherent integration output. The present coherent integration result is multiplied with the delay conjugate of the last one. The mathematical representation is:

(10)$$\eqalign{Q(k) =\; &Y(k)\cdot Y(k-1)^* = [0\cdot 5AR(\delta \tau )]^2\exp (j\cdot \pi \delta f_dT_c) \cr & \times D_kD_{k-1}H_kH_{k-1}{\rm sinc}[\delta f_d(k)T_c]{\rm sinc}[\delta f_d(k-1)T_c]}$$

where ( · ) * denotes the conjugate operation, πδf _dT _c is the carrier phase deviation between Y(k) and Y(k − 1) and the coherent integration time T _c is 1 ms. The sign of Q(k) changes due to bit reversions.

Assuming the accumulation time is 20 ms, the accumulation result of differential coherent values is defined as:

(11)$$\eqalign{ \sum\limits_{k = 1}^u Q (k) = &\sum\limits_{k = 1}^u Y (k)\cdot Y(k-1)^* = [0\cdot 5AR(\delta \tau )]^2\exp (j\cdot 2\pi \delta f_dT_c) \cr & \times \sum\limits_{k = 1} {D_k} D_{k-1}H_kH_{k-1}{\rm sinc}[\pi \delta f_d(k)T_c]{\rm sinc}[\pi \delta f_d(k-1)T_c]}$$

where u is the accumulation time.

For GPS legacy signals, the sign of the product can only change at each bit boundary, so there is at most one change in 20 ms and it corresponds with the real bit boundary. However, for BeiDou signals, the data bit rate increases to 1 kbps and 10 bit reversions occur in every navigation data bit, so the differential coherent process of two consecutive correlation outputs is not suitable.

3.2. Zero-padding and differential coherent acquisition algorithm

The sign of Q(k) in Equation (10) depends on D _kD _k−1H _kH _k−1. The NH code is as long as the duration of a navigation message bit, and H _k = H _k−20m, where m is an integer. If the delay time of the differential coherent algorithm is 20 ms, H _k H _k−20 will equal to one and the sign of the differential coherent result will only depend on transitions of the navigation data and have nothing to do with NH code transitions.

In the proposed algorithm the differential delay time is modified to 20 ms to weaken the effect of NH code. The steps of the proposed algorithm are as follows:

Step 1. Continuous zero-padding.

1. (1) Incoming signal samples of 2 ms are used each time and half of the samples will be replaced in a first-in first-out fashion every 1 ms. The incoming signals are correlated with the local carrier in the frequency dimension with the frequency bin of 500 Hz for Doppler removal. Then, the FFT of the incoming signal samples is taken.

2. (2) The 1 ms locally generated code sequence is appended by 1 ms zeros to form the extended replica. FFT and complex conjugate of the extended replica are taken.

3. (3) The incoming signal after Doppler removal is multiplied with the replica for correlation operation. Then, IFFT is operated on the products. Instead of peak detection, the first 1 ms of the correlation result is stored.

Step 2. Differential coherent algorithm with long delay time.

The differential coherent algorithm is operated on the stored 1 ms correlation result. The delay time of differential coherent algorithm is modified to 20 ms and H _kH _k−20 = 1. The mathematical representation is:

(12)$$\eqalign{ Z(k) = &\;Y(k)\cdot Y(k-20)^* = [0\cdot 5AR(\delta \tau )]^2\exp ( j\cdot \pi \delta f_dT_c) \cr & \times D_kD_{k-20}H_kH_{k-20} {\rm sinc}[\delta f_d(k)T_c] {\rm sinc}[\delta f_d(k-20)T_c] \cr =& \;[0\cdot 5AR(\delta \tau )]^2\exp ( j\cdot \pi \delta f_dT_c) \cr & \times D_kD_{k-20} {\rm sinc}[\delta f_d(k)T_c] {\rm sinc}[\delta f_d(k-20)T_c]}$$

where πδf _dT _c is the carrier phase deviation between Y(k) and Y(k − 20).

The accumulation time is set to 10 ms and the accumulation result of differential coherent algorithm values is:

(13)$$\eqalign{\sum\limits_{k = 1}^{10} Z (k) =& \sum\limits_{k = 1}^{10} Y (k)\cdot Y(k-20)^* = [0\cdot 5AR(\delta \tau )]^2\exp (j\cdot 2\pi \delta f_dT_c) \cr & \times \sum\limits_{k = 1}^{10} {D_k} D_{k-20} {\rm sinc}[\pi \delta f_d(k)T_c] {\rm sinc}[\pi \delta f_d(k-20)T_c]}$$

here, $\sum\nolimits_{k=1}^{10}\, \hbox{sinc}\lsqb \delta f_{d}\lpar k\rpar T_{c}\rsqb \, \hbox{sinc}\lsqb \delta f_{d}\lpar k-20\rpar T_{c}\rsqb $ represents the effect of frequency deviation.

Step 3. Acquisition detection.

As the signs of the modified differential coherent result only depend on transitions of the navigation data, two continuous 10 ms accumulations of differential coherent values are needed, one of which crosses bit transition. Then the maximum peak magnitude of these two results can be detected to estimate acquisition parameters. The location along the time axis is the millisecond boundary and the location along the Doppler axis is an estimate of the Doppler frequency. Then, based on the estimated acquisition parameters, tracking loops can be carried out to obtain measurements (Xie et al., Reference Xie, Sun, Kang, Qian, Zhao and Zhang2017).

3.3. Performance analysis

In Equation (5), $\sum\nolimits_{k=1}^{10}\hbox{sinc}\lsqb \delta f_{d}\lpar k\rpar T_{c}\rsqb \, \hbox{sinc}\lsqb \delta f_{d}\lpar k-20\rpar T_{c}\rsqb $ is the effect of frequency deviation in the proposed algorithm. In Equation (9), sin c[δf _d(k) T _coh] represents the frequency deviation effect of the acquisition algorithm based on a long NH-PRN code. In order to compare the frequency deviation effects of these two algorithms, approximations should be made to $\sum\nolimits_{k=1}^{10}\hbox{sinc}\lsqb \delta f_{d}\lpar k\rpar T_{c}\rsqb \, \hbox{sinc}\lsqb \delta f_{d}\lpar k-20\rpar T_{c}\rsqb $. As δ f _d changes slightly in a single search, δf _d(k) can be approximately equal to δf _d(k − 20) and sinc[δf _d(k) T _c] sinc[δf _d(k − 20) T _c] can be approximately equal to sinc[δf _d(k) T _c] ² as:

(14)$$\sum\limits_{k = 1}^{10} {\rm sinc}[\delta f_d(k)T_c] {\rm sinc}[\delta f_d(k-20)T_c]\approx \sum\limits_{k = 1}^{10} {{\{ {\rm sinc}[\delta f_d(k)T_c]\} }^2} \approx 10\{ {\rm sinc}[\delta f_d(k)T_c]\} ^2$$

(15)$$10\{ {\rm sinc}[\delta f_d(k)T_c]\} ^2 = 10\{ {\rm sinc}[0\cdot 001\delta f_d(k)]\} ^2\ll {\rm sinc}[\delta f_d(k)T_{coh}] = {\rm sinc}[0\cdot 02\delta f_d(k)]$$

where T _c = 1 ms, T _coh = 20 ms.

The effect of frequency deviation in the proposed algorithm is much smaller than the acquisition algorithm based on long NH-PRN code. Because T _c = 1 ms, the frequency bin can be 500 Hz and then the number of searches is 20 for a frequency searching range of 10 kHz.

4.1. Monte Carlo simulation test

The simulation platform is shown in Figure 3, which consisted of BeiDou B1I IF signal generation and software receiver baseband processes. The BeiDou IF signal was generated according to the trace and ephemeris parameters and additional white Gaussian noise was generated according to the C/N0 settings. The parameters used in Monte Carlo simulations are shown in Table 1.

Figure 3. The diagram of the platform for Monte Carlo simulations.

Table 1. Parameters used in Monte Carlo simulations.

To prove the reliability of the proposed algorithm, Monte Carlo simulations were carried out to evaluate the sensitivity with different C/N0s. Additional white Gaussian noise was generated for each trial, and 1,000 trials were used for each probability. Monte Carlo simulations of different algorithms were conducted for comprehensive evaluation.

First, the performance of two traditional BeiDou signal acquisition methods were tested. The detection probabilities are plotted as a function of the C/N0s in Figure 4. The blue line indicates the performance of the 6 ms repeated search acquisition method and the red line indicates the performance of the zero-padding acquisition method. Using the repeated search acquisition method, the detection probability was 0·83 at C/N0 of 32 dB-Hz. However, the detection probabilities were 0·56 at C/N0 of 32 dB-Hz and 0·83 at 33 dB-Hz using the zero-padding method. The sensitivity of the repeated search acquisition method was about 1 dB better than the zero-padding method.

Figure 4. Detection probabilities of two traditional BeiDou signal acquisition methods under C/N0s of 27~37 dB-Hz.

Secondly, the performance of the proposed algorithm was tested. The black line in Figure 5 indicates the detection probability of the proposed algorithm. The detection probability was 0·94 at C/N0 of 22 dB-Hz using the proposed algorithm and the acquisition sensitivity was about 10 dB better than the traditional BeiDou signal acquisition methods.

Figure 5. Detection probabilities of the proposed algorithm under the C/N0s of 20 ~ 38 dB-Hz.

Finally, a Non-Coherent (NC) combination algorithm based on zero-padding was made to compare with the proposed algorithm. The non-coherent algorithm was operated on the zero-padding result and the mathematical representation is:

(16)$$\sum_{k=1}^uY(k)^2=\sum_{k=1}^u \vert0{\cdot}5AD_kH_kR(\delta\tau)\,\hbox{sinc}[\delta f_d(k)T_c]\exp[\,j\cdot\phi_k+j\pi\delta f_d(k)T_c]\vert^2$$

where u is the integration time.

The performance of the NC combination algorithm was tested. The detection probabilities are plotted as a function of the C/N0s in Figure 6. The green line indicates the performance of the NC algorithm with an integration time of 20 ms and the pink line indicates the NC algorithm with an integration time of 10 ms. Using the NC algorithm, the detection probability was only 0·79 at C/N0 of 27 dB-Hz with an integration time of 10 ms and was 0·64 at C/N0 of 25 dB-Hz with an integration time of 20 ms. This is because the noise was also squared and averaged toward a non-zero value. This value is referred to as the squaring loss which attenuates the SNR gain of the NC method significantly.

Figure 6. Detection probabilities of non-coherent combination algorithm under the C/N0s of 20 ~ 38 dB-Hz.

It is clear that the proposed algorithm outperforms the NC combination algorithm and the detection probability of the proposed algorithm is always higher than the NC combination algorithm.

The results show that the detection probabilities of the traditional BeiDou signal acquisition methods are very low in a weak signal environment. In contrast, the proposed algorithm outperforms the NC combination algorithm and the acquisition sensitivity of the proposed algorithm is about 10 dB better than the traditional BeiDou signal acquisition methods.

4.2. Real data tests

In order to confirm the feasibility of the proposed algorithm, real data tests were carried out. The BeiDou B1I receiver used was developed by the Navigation Research Center, Nanjing University of Aeronautics and Astronautics (NRC, NUAA) and the hardware structure was developed by Shanghai Yuzhi. As shown in Figure 7, the receiver test platform consisted of a DSP-TMS320C6713B and an FPGA-EP4CE115F23, and an AD8347 was used as the quadrature down-conversion mixer. The sampling frequency was 62 MHz and the IF carrier frequency was 4·098 MHz. The proposed algorithm and the traditional 6 ms repeated search algorithm were processed in the receiver at the same time.

Figure 7. The BeiDou B1I receiver test platform.

First, the antenna was placed in an open area. The number of visible satellites was ten using the traditional 6 ms repeated search algorithm and was 11 using the proposed algorithm. Satellite {#}9 is used as an example and the acquisition result is shown in Figure 8. The peak amplitude in the proposed algorithm was much higher than for the traditional algorithm.

Figure 8. Accumulation results of PRN-9 in an open area: Top traditional 6 ms repeated search algorithm; Bottom the proposed algorithm.

Next, the antenna was placed under dense foliage. The number of visible satellites was three using the traditional 6 ms repeated search algorithm and was 11 using the proposed algorithm. Satellite {#}9 is used as an example and the acquisition result is shown in Figure 9. Under dense foliage, the peak of traditional algorithm was not obvious and acquisition of satellite {#}9 failed. However, in the proposed algorithm, the peak is still obvious. Therefore, the acquisition sensitivity is higher than the traditional algorithm.

Figure 9. Accumulation results of PRN-9 under dense foliage: Top traditional 6 ms repeated search algorithm; Bottom the proposed algorithm.