1. INTRODUCTION
The BeiDou Navigation Satellite System, which is a Chinese global satellite navigation system, has gradually been brought online in the last few years. The China Satellite Navigation Office (2013) has released the Interface Control Document (ICD) for open service signals B1 and B2. This indicates that BeiDou can provide users with a dual-frequency high precision positioning service. A satellite navigation receiver can receive the signals broadcast by the satellites in view and process them to provide location and navigation information. The key part is baseband signal processing, which is based on the acquisition algorithm. The operating time of the acquisition algorithm will directly affect the Time To First Fix (TTFF) of the receiver (Jan and Lin, Reference Jan and Lin2009). C B2I is the same as the Coarse/Acquisition (C/A) code of the Global Positioning System (GPS) in sequence generation, which is generated by two series of Gold code sequences. The difference is that the Gold code in C B2I is generated through 11 shift registers, while the Gold code in C/A code is generated from ten shift registers. The code rate of BeiDou C B2I code is twice that of GPS C/A code. In order to ensure acquisition accuracy, the sampling rate of the B2 signal needs to be increased accordingly. The number of sampling points to be processed by the acquisition algorithm will also increase. Using the traditional algorithm will extend the acquisition time. Therefore, it is necessary to study the fast acquisition algorithm (Simone et al., Reference Simone, Fittipaldi and Sanchez2011) to reduce the receiver delay (Zhang et al., Reference Zhang, Xu, Cai and Li2016; Jin et al., Reference Jin, Yang, Huang and Qin2015).
When the software receiver demodulates and spreads received signals, there are three important parameters: the number of satellites in view, the Doppler frequency offset and the code phase shift. The process of signal acquisition can be seen as a three-dimensional search for the above three parameters (Tsui, Reference Tsui2005). The Serial Search Acquisition (SSA) algorithm mainly takes advantage of the correlation property of the ranging code (that is the C/A code or C B2I code) and scans all the cells in the search space (Polydoros and Weber, Reference Polydoros and Weber1984). To speed up the acquisition process, Van Nee and Coenen (Reference Van Nee and Coenen1991) applied a Fast Fourier Transform (FFT) to the correlation of the ranging code and proposed the Parallel Code phase Acquisition (PCA) algorithm. This converted correlation in the time-domain to a multiplication in the frequency-domain and obtained correlation results at every sampling point at one time. Thus, only a two-dimensional search for the received signal was needed. In this way, the computational load was significantly reduced. However, the computational load is still large when the satellite signal is long or the sampling frequency is high, so many researchers have studied acquisition algorithms based on PCA algorithms. For instance, Starzyk and Zhu (Reference Starzyk and Zhu2001) combined an Averaging Correlations (AC) method with the PCA algorithm and proposed the AC-PCA algorithm. This divided the 5,000-point FFT to six 1,023-point FFTs in the correlation operation. Fantino et al. (Reference Fantino, Pini, Mulassano, Girau, Nicola and Nordio2008) analysed the performance of this algorithm thoroughly, and this work shows that the process speed can be improved without loss of energy. Chun (Reference Chun2005) applied segmented FFT to the acquisition of L2 signals. In this algorithm, L2 signals were divided into small segments and processed after arriving at the receiver. In this way, the latency, memory and computation requirements per segment were reduced.
Lin et al. (Reference Lin, Tsui and Howell1999) proposed the Double-Block Zero-Padding (DBZP) method based on a delay-and-multiply approach. This algorithm was used to capture the P(Y)-Code directly. It broke down the long data into M subsets and performed local correlation in the frequency domain, where M was the number of Doppler bins. In this way, the length of the FFT block was decreased and the calculation time could be reduced. Jan and Lin (Reference Jan and Lin2009) proposed a multi-C/A code acquisition method. It replaced the single local C/A code in the PCA algorithm to a sum of two or more C/A codes of different satellites. This can search for several satellites at one time and speed up the acquisition process. Leclère et al. (Reference Leclère, Botteron and Farine2012) combined the Fast Finite Impulse Response (FIR) Algorithms (FFAs) and PCA algorithm into a new algorithm. This algorithm decomposed the initial circular correlation into several smaller circular correlations and reduced the length of the FFT blocks. Patel and Shukla (Reference Patel and Shukla2011) exploited the symmetry of the power spectrum of the C/A code to reduce the number of points to perform correlation and FFT. The reduction in FFT/Inverse FFT (IFFT) points can decrease the computation time and hardware complexity. The above-mentioned algorithms all improved the PCA algorithm and reduced the computational load or resource consumption, but at the cost of increasing algorithm complexity or sometimes decreasing acquisition accuracy.
In this paper, we propose a new acquisition algorithm. It divides the signal into K parts to conduct FFT, and obtains the correlation results in the frequency domain with a fusion algorithm. This algorithm can shift K sampling points in one correlation operation. What is more, it can also be applied to capture any satellite signals modulated in Binary Phase Shift Keying (BPSK) and Quadrature Phase Shift Keying (QPSK) modes with a lower computational load. In this paper, this algorithm will be denoted as the K-PCA algorithm.
The rest of the paper is organised as follows. In Section 2, we briefly introduce a model of the B2 signal. A quick review of the SSA and PCA algorithms is presented in Section 3. Then the principle of the K-PCA algorithm is presented in Section 4. In Section 5, we evaluate the computational complexity of the K-PCA algorithm and the PCA algorithm. The method to choose the optimal K value is also described in this section. Section 6 shows the algorithm operation simulation results of the K-PCA algorithm and conclusions are drawn in Section 7.
2. MODEL OF THE B2 SIGNAL
The China Satellite Navigation Office (2013) has released the ICD for open service signal B2, describing its structure, basic parameters, ranging code features and navigation data format. Unlike the BPSK modulation of the GPS signal, BeiDou adopts the QPSK modulation mode (China Satellite Navigation Office, 2013). Compared with BPSK, the QPSK transmits In-phase/Quadrature-phase (I/Q) signals simultaneously and utilises the frequency band more efficiently. The B2 signal consists of the I branch (open) and the Q branch (authorisation), and can be expressed as:
$$S_{B2}^{\,j}(t) = A_{B2I}C_{B2I}^{\,j}(t)D_{B2I}^{\,j}(t)\cos (2\pi f_{B2}t + \varphi_{B2I}^{\,j}) + A_{B2Q}C_{B2Q}^{\,j}(t)D_{B2Q}^{\,j}(t)\sin (2\pi f_{B2}t + \varphi_{B2Q}^{\,j})$$In this equation, superior j represents the serial number of the satellite. A B2I and A B2Q indicate the signal intensity of the I branch and the Q branch. C B2I and C B2Q are ranging codes modulated in the I/Q branches. D B2I and D B2Q are navigation data. φB2I and φB2Q are carrier phases of the I/Q branches respectively. The carrier frequency f B2 is 1207.14 MHz. This paper only studies the I branch signal, as the Q branch signal needs authorisation and is not open to civil use. The C B2I code is a Gold sequence composed of 2,046 chips, and the code rate is 2.046 Mcps.
The satellite signal is obtained by modulating the C B2I code and navigation data with the carrier in the QPSK mode. Unlike the traditional QPSK mode, the B2 signal digital data (C B2I code and navigation data) do not go through a multiplexer but multiply with the carrier of the I (or Q) branch directly (Xie et al., Reference Xie, Liu, Li and Feng2016). The modulation process is illustrated in Figure 1.
Modulation mode of B2 signal.
The code rate of C B2I is twice that of the GPS C/A code. To ensure acquisition accuracy, the signal sampling frequency needs to be improved. Thus, the number of sampling points that need to be processed will increase too, which will affect the real-time quality of the receiver. The heavy computational load correspondingly will make the software receiver much more difficult to implement on an embedded platform with limited resources (Jin et al., Reference Jin, Ye, Lin and Hua2008).
3. PRINCIPLES OF CONVENTIONAL ACQUISITION ALGORITHMS
When the signal is transmitted from the satellite to the GNSS receiver, its carrier frequency and code phase will change due to the Doppler effect and the unknown distance. Also, it will be submerged in noise. Thus, it is necessary for the GNSS receiver to raise the signal from the noise by acquisition blocks. The acquisition algorithm mainly takes advantage of the correlation property of the ranging code. The process of the acquisition algorithm is illustrated in Figure 2(a). Once the local carrier and Pseudorandom Number (PRN) code are synchronised with the received signal, there is a peak and the signal can be captured (Akopian, Reference Akopian2005).
Principle diagrams of acquisition algorithms.
The PCA algorithm converted the correlation in the time-domain to multiplication in the frequency-domain by using FFT, while the correlation result at every sampling point can be obtained simultaneously by applying IFFT to the product. The process is shown in Figure 2(b). For each satellite, this algorithm only requires searching for the Doppler frequency offset, so the computational load can be significantly reduced. Though the PCA approach is computationally efficient, it must compute the correlation results at every sampling point due to the FFT and IFFT. Thus, its structure is less flexible than the SSA algorithm illustrated in Figure 2(a), which can shift several sampling points in one correlation operation. Provided that the signal sampling frequency is relatively high or the number of signal sampling points that need to be processed is large (e.g. weak signals and long period signals), this algorithm will need large FFT processing blocks to process the data, which will increase the load on the receiver (Chun et al., Reference Chun, Mikel and Thao2006). Therefore, if we can combine the advantages of the SSA and PCA algorithms, the length of FFT blocks and computational load can be significantly decreased. Based on this idea, we improve the structure of the PCA - more precisely the structure of FFT blocks in the correlation operation, and propose the following algorithm. It shifts K sampling points in every correlation operation, and the K value can be arbitrarily selected according to the situation; the algorithm is thus called the K-PCA algorithm.
4. THE CONCEPT OF K-PCA ALGORITHM
4.1. The principle of the 3-PCA algorithm (take K=3 as an example)
The correlation (Wang et al., Reference Wang, Hu, Xiong and Song2013) between mixed signal x(n) and local C B2I code h(n) is shown in Equation (2).
$$y(n) = \sum_{m = 0}^{N - 1} x(m) h(m-n)$$where x(n), h(n) and y(n) are finite sequences whose lengths are N(N is the number of sampling points in one period of the signal). x(n) means mixed signal x(n), h(n) means local C B2I, and the period of h(n) is 1 ms.
We then transform the correlation in the time domain to multiplication in the frequency domain by applying FFT to y(n). The transformed result can be expressed as Equation (3), where 
$W_{N}^{nk} = e^{-j2\pi kn/N}$ represents the twiddle factor, and H*(k) is the complex conjugate of H(k).
$$\eqalign{Y(k)&=\sum_{n=0}^{N-1}y(n)W_{N}^{nk}=\sum_{n=0}^{N-1}\left(\sum_{m=0}^{N-1}x(m)h(m-n)\right)W_{N}^{nk}=\sum_{n=0}^{N-1}\left(\sum_{m=0}^{N-1}x(m)h(m-n)\right)W_{N}^{mk}W_{N}^{-(m-n)k}\cr &=\underbrace{{\sum_{m=0}^{N-1}x(m)W_{N}^{mk}}}_{X(k)}\underbrace{\sum_{n=0}^{N-1}h(m-n)W_{N}^{-(m-n)k}}_{H^\ast(k)}=X(k)H^{\ast}(k)}$$Meanwhile, X(k) can be represented as follows:
$$\eqalign{X(k)&=\sum_{n=0}^{N-1}x(n)W_{N}^{nk}\,=\,\mathop{\sum_{n=0}^{N-1}x(n)W_N^{nk}}\limits_{{\rm mod}(n,3)=0}\,+\,\mathop{\sum_{n=0}^{N-1}x(n)W_{N}^{nk}}\limits_{{\rm mod}(n,3)=1}\,+\,\mathop{\sum_{n=0}^{N-1}x(n)W_{N}^{nk}}\limits_{{\rm mod}(n,3)=2}\cr &=\sum_{r=0}^{{N}/{3^{-1}}}x(3r)W_{N}^{3rk}\,+\,\sum_{r=0}^{{N}/{3^{-1}}}x(3r+1)W_{N}^{(3r+1)k}\,+\,\sum_{r=0}^{{N}/{3^{-1}}}x(3r+2)W_{N}^{(3r+2)k}\cr &=\underbrace{\sum_{r=0}^{{N}/{3^{-1}}}x_{0}(r)W_{N/3}^{rk}}_{X_0(k)}\,+\,W^k_N \underbrace{\sum_{r=0}^{{N}/{3^{-1}}}x_{1}(r)W_{N/3}^{rk}}_{X_1(k)}\,+\,W_N^{2k} \underbrace{\sum_{r=0}^{{N}/{3^{-1}}}x_{2}(r)W_{N/3}^{rk}}_{X_2(k)}\cr &=X_{0}(k)+W_{N}^{k}X_{1}(k)+W_{N}^{2k}X_{2}(k)}$$where
$$\eqalign{&\left\{\matrix{x_{0}(r)=x(3r)\hfill \cr x_{1}(r)=x(3r+1)\hfill \cr x_{2}(r)=x(3r+2)\hfill} \right. \hfill &r\in[1,2,\ldots,N/3-1] \cr &X_{i}(k) = FFT[x_{j}(r)]\quad i \in [0,2]}$$Similarly,
$$H^{\ast}(k)=H_{0}^{\ast}(k)+W_{N}^{-k}H_{1}^{\ast}(k)+W_{N}^{-2k}H_{2}^{\ast}(k);\quad Y(k)=Y_{0}(k)+W_{N}^{k}Y_{1}(k)+W_{N}^{2k}Y_{2}(k)$$where
$$h_{i}(r)=h(3r+i),H_{i}^{\ast}(k)=[FFT[h_{i}(r)]]^{\ast}\quad i\in[0,2],\quad r\in[1,2,\ldots,N/3-1]$$According to Equations (3), (4) and (6), we can obtain:
$$\eqalign{Y(k)&=[X_{0}(k)+W_{N}^{k}X_{1}(k)+W_{N}^{2k}X_{2}(k)]\times[H_{0}^{\ast}(k)+W_{N}^{-k}H_{1}^{\ast}(k)+W_{N}^{-2k}H_{2}^{\ast}(k)]\cr &=\underbrace{X_{0}(k)H_{0}^{\ast}(k)+X_{1}(k)H_{1}^{\ast}(k)+X_{2}(k)H_{2}^{\ast}(k)}_{Y_0(k)}+ \underbrace{W_{N}^{-2k}X_{0}(k)H_{2}^{\ast}(k)+W_{N}^{k}X_{1}(k)H_{0}^{\ast}(k)+W_{N}^{k}}_{W_N^kY_1(k)}\cr &\quad +\underbrace{W_{N}^{-k}X_{0}(k)H_{1}^{\ast}(k)+W_{N}^{-k}X_{1}(k)H_{2}^{\ast}(k)+W_{N}^{2k}X_{2}(k)H_{0}^{\ast}(k)}_{W^{2k}_NY_2(k)}}$$
$Y_{i}(k)(i\in[0,2])$ can be represented as:
$$\eqalign{Y_{0}(k)&=H_{0}^{\ast}(k)X_{0}(k)+H_{1}^{\ast}(k)X_{1}(k)+H_{2}^{\ast}(k)X_{2}(k)\cr Y_{1}(k)&=H_{0}^{\ast}(k)X_{1}(k)+H_{1}^{\ast}(k)X_{2}(k)+W_{N}^{-3k}H_{2}^{\ast}(k)X_{0}(k)\cr Y_{2}(k)&=H_{0}^{\ast}(k)X_{2}(k)+W_{N}^{-3k}H_{1}^{\ast}(k)X_{0}(k)+W_{N}^{-3k}H_{2}^{\ast}(k)X_{1}(k)}$$ The result in the time domain can be derived through IFFT after 
$Y_{i}(k)(i\in[0,2])$ is computed, and can be represented as:
$$y_{i}(n)=IFFT(Y_{i}(k))\quad i\in[0,2]$$ The final correlation result in the time domain can be obtained by combining 
$y_{i}(n) (i\in [0,2])$.
Next, we will take y 0(n) as an example to analyse the correlation process of the signal. The algorithm divides the mixed signal and local C B2I code into three parts (Jin et al., Reference Jin, Lu, Liu, Qin and Luo2013), which is shown in Figure 3(a), and the FFT operation is applied separately to each part. From Equation (9), we can see that, in the frequency domain, Y 0(k) is a sum of three parts—
$X_{0} (k) \times H_{0}^{\ast} (k)$, 
$X_{1} (k) \times H_{1}^{\ast}(k)$ and 
$X_{2} (k) \times H_{2}^{\ast} (k)$. In the time domain, y 0(n) is a sum of three parts—
$x_{0}(r)^{\ast}h_{0}(r)$, 
$x_{1}(r)^{\ast}h_{1}(r)$ and 
$x_{2}(r)^{\ast}h_{2}(r)$ where * means convolution. The computing process of y 0(0) and y 0(1) is shown in Figure 3(b). From Figure 3(a) we can see that when 
$h_{i}(r)(i \in [0,2])$ shifts by one bit to the right, the h(n) shifts by three bits to the right compared to the original sequence. So, in this algorithm, y 0(n) is still the cross-correlation result between x(n) and h(n). The difference is that h(n) shifts by three bits in every cross-correlation operation.
Computing process of y 0(n).
The computing process of y 1(n) and y 2(n) is similar to that of y 0(n), but in Equation (9), 
$H_{i}^{\ast}(k)$ is multiplied by the twiddle factor 
$W_{N/3}^{-k}$. In light of the circumference-shifting characteristic of an FFT (Cheng, Reference Cheng2012), we can see that this operation is equal to shifting h i(r) to the right by one bit first, and then applying the FFT and the conjugate operation to it. The computing process of y i(n) is shown in Figure 4. The final correlation result y′(n) can be obtained by combining 
$y_{i}(n)(i\in [0,2])$ together, as shown in Equation (11).
Computing process of y i(n).
$$y^{\prime}(3n) = y_{0}(n); y^{\prime}(3n + 1) = y_{1}(n); y^{\prime} (3n + 2) = y_{2}(n) \quad n\in[0,N/3-1]$$If the difference between local carrier frequency and received signal frequency is bigger than the search step (for example, 1 kHz or 500 Hz), there will be no peak in the correlation results between the mixed signal and local C B2I code. Similarly, there will be no value exceeding the threshold in y 0(n), y 1(n) and y 2(n). If the local carrier frequency is the same as the received signal frequency, there will be peaks in y(n). Also, the correlation results in the sampling points near the code phase will exceed the threshold. The peaks in y 0(n), y 1(n) and y 2(n) will also exceed the threshold. Hence, it can be judged whether the frequency of the local carrier is the same as the frequency of the received signal by the peak value of y 0(n). When a correlation peak in y 0(n) occurs that exceeds the threshold, it indicates that the search for frequency offset is accomplished. Otherwise, the frequency of the local carrier should be changed and the acquisition steps need to be repeated until a correlation peak appears. Afterwards, the value of y 1(n) and y 2(n) should be calculated. We then compare the peak values of y 0(n), y 1(n) and y 2(n). If the biggest peak value is located at the mth points in y i(n) (m is a constant), the final code phase is located at the (3m+i)th sampling point.
The acquisition process of the 3-PCA algorithm is demonstrated in Figure 5 and as follows.
Acquisition process of 3-PCA algorithm.
Step 1: First, strip off the carrier from the received signal by the process shown in Equation (12), where S IF represents the received signal and L(n) represents the local carrier. Next, divide the mixed signal and local C B2I code into three parts by the process shown in Equation (13).
$$x(n)=S_{IF}(n)L(n)$$
$$\matrix{\left\{\matrix{x_{0}(r)=x(3\times r) \hfill \cr x_{1}(r)=x(3\times r+1) \hfill \cr x_{2}(r)=x(3\times r+2) \hfill} \right. &\left\{\matrix{h_{0}(r)=h(3\times r) \hfill \cr h_{1}(r)=h(3\times r+1) \hfill \cr h_{2}(r)=h(3\times r+2) \hfill} \right. &r\in \left[0, \displaystyle{N \over 3}-1\right]}$$ Step 2: Apply the FFT to 
$x_{i}(r), h_{i}(r)i\in\{0,1,2\}$, then take the conjugate value of H i(k). The procedures are as follows:
$$\eqalign{X_{i}(k) &= FFT(x_{i}(r))\cr H_{i}(k) &=FFT(h_{i}(r))\cr H_{i}^{\ast}(k) &=conj(H_{i}(k)) \quad i \in \{0,1,2\}, \quad r\in [1,2, \ldots, N/3 - 1]}$$Step 3: Obtain Y 0(k) according to the fusion algorithm shown in Equation (15).
$$Y_{0}(k)=X_{0}(k)H_{0}^{\ast}(k) + X_{1} (k) H_{1}^{\ast} (k) + X_{2} (k) H_{2}^{\ast} (k)$$Step 4: Then obtain the correlation results y 0(n) in the time-domain by applying IFFT to Y 0(k). If the peak value of y 0(n) exceeds the threshold value, the value of Y 1(k) and Y 2(k) will be computed according to the fusion algorithm shown in Equation (9). Then obtain the correlation results y 1(n) and y 2(n) in the time-domain by applying IFFT to Y 1(k) and Y 2(k). The final code phase can be determined by comparing the peak values of y 0(n), y 1(n) and y 2(n). Otherwise the frequency of the local carrier should be changed and Steps 1 to 4 need to be repeated until a correlation peak appears.
The proposed algorithm shortens the length of signals for FFT (Kurz et al., Reference Kurz, Kappen, Coenen and Noll2010) operation by dividing them into three parts, which makes it easy to achieve on an embedded platform with limited resources (Humphreys et al., Reference Humphreys, Psiaki and Kintner2006). Moreover, in the PCA algorithm (Leclère et al., Reference Leclère, Botteron and Farine2013), the local C B2I code shifts by only one bit in every correlation operation, while it shifts by three bits every time in the 3-PCA algorithm. In addition, the number of sampling points involved in every correlation operation remains the same in the two algorithms. Hence, the proposed algorithm will not decrease the acquisition accuracy.
4.2. The principle of the K-PCA algorithm
The above 3-PCA algorithm divides the mixed signal and local C B2I code into three parts. Furthermore, the signal can be divided into K parts to perform the correlation operation. However, with the increase of the K value, the complexity of the algorithm will also be increased. The principle of the K-PCA algorithm is shown as follows.
First, the mixed signal is divided into K parts:
$$\matrix{\left\{\matrix{x(rK)=x_{0}(r) \hfill \cr x(rK+1)=x_{1}(r) \hfill \cr \vdots\cr x(rK+K-1)=x_{K-1}(r) \hfill}\right. &r\in \left[0, \displaystyle{N \over K} -1 \right] \hfill}$$Then a FFT is applied to x(n), and X(k) can be written as Equation (17):
$$\eqalign{X(k) &=\sum_{n=0}^{N-1} x(n)W_{N}^{nk}= \sum_{r=0}^{N/K-1} x(rK)W_{N}^{rKk} + \sum_{r=0}^{N/K-1} x(rK +1)W_{N}^{(rK+1)k}\cr &\quad +\cdots \sum_{r=0}^{N/K-1} x(rK +K -1)W_{N}^{(rK+K-1)k}}$$
$$\eqalign{&= \sum_{r=0}^{N/K-1} x_{0}(r)W_{N/K}^{rk}+W_{N}^{k} \sum_{r=0}^{N/K-1} x_{1}(r)W_{N/K}^{rk}\cr &\quad +\cdots W_{N}^{(K-1)k} \sum_{r=0}^{N/K-1} x_{K-1}(r)W_{N/K}^{rk}=\sum_{i=0}^{K-1} W_{N}^{ik}X_{i}(k)\cr X_{i}(k)&=FFT(x_{i}(n))}$$Similarly, H*(k) and Y(k) can be represented as:
$$H^{\ast}(k)=\sum_{i=0}^{K-1}W_{N}^{-ik}H_{i}^{\ast}(k),\quad Y(k)=\sum_{i=0}^{K-1}W_{N}^{ik}Y_{i}(k)$$The correlation result in the frequency domain is shown in Equation (20):
$$Y(k)=X(k)H^{\ast}(k)=\left(\sum_{i=0}^{K-1}W_{N}^{ik}X_{i}(k)\right)\times\left(\sum_{i=0}^{K-1}W_{N}^{-ik}H_{i}^{\ast}(k)\right)$$
according to Equations (19) and (20), 
$Y_{i}(k)(i\in[0,K-1])$ can be defined as follows
$$\eqalign{Y_{0}(k)&=H_{0}^{\ast}(k)X_{0}(k)+H_{1}^{\ast}(k)X_{1}(k)+\ldots+H_{K-2}^{\ast}(k)X_{K-2}(k)+H_{K-1}^{\ast}(k)X_{K-1}(k)\cr Y_{1}(k)&=H_{0}^{\ast}(k)X_{1}(k)+\ldots+H_{K-3}^{\ast}(k)X_{K-2}(k)+H_{K-2}^{\ast}(k)X_{K-1}(k)+W_{N/K}^{-k}H_{K-1}^{\ast}(k)X_{0}(k)\cr &\vdots \cr Y_{p-1}(k)&=H_{0}^{\ast}(k)X_{K-1}(k)+W_{N/K}^{-k}H_{1}^{\ast}(k)X_{0}(k)+W_{N/K}^{-k}H_{2}^{\ast}(k)X_{1}(k)+\ldots+W_{N/K}^{-k}H_{K-1}^{\ast}(k)X_{K-2}(k)}$$The correlation results in the time domain y i(n), therefore, can be obtained by applying an IFFT to Y i(k). The process of the correlation is the same as the process shown in Figure 5 where K=3. In this algorithm, h(n) shifts by K bits in every correlation operation. We can choose the value of K according to the situation. Meanwhile, the value of K should be smaller than the number of sampling points contained in half a C B2I code chip. This is because when the local C B2I code shifts more than half a chip every time, there will be no peaks in y i(n) even though the frequency of the received signal and local carrier is the same. In addition, the choice of K is also related to the number of sampling points N, which will be analysed in the next section.
4.3. Method to reduce the computational complexity
y 0(n) is the correlation result between h(n) and x(n), where h(n) shifts by K bits every time. In this acquisition process, if the peak value of y 0(n) exceeds the threshold, 
$y_{i}(n)(i\in[1,K-1])$ will be further computed to determine at which sampling point the code phase is located. This operation is not efficient enough as it involves a number of multiplications and IFFTs which are computation-intensive. Since the range code phase can be derived from the peak position of y 0(n), the final code phase can be asserted by comparing the correlation values in the sampling points around the peak position of y 0(n).
h(n) shifts by K bits every time it is correlated with x(n). If the coordinate of the peak of y 0(n) is (m, R), the range of the code phase should be (K×m−K~K×m+K). Then the final code phase can be determined by comparing the correlation results at the (K×m−K~K×m+K)th sampling point. The process is illustrated in Figure 6 (take K=3 as an example) and Equation (22). This method only needs 2K×N multiplications and 2K×(N−1) additions, and the procedure is simpler.
Process of the method to reduce the computational load.
$$y_{K}(n)=\sum_{m=0}^{N-1}x(m)h(m-n)\quad n\in[K\times m-K,K\times m+K]$$5. ANALYSIS OF COMPUTATIONAL COMPLEXITY
When implemented on an embedded platform with limited resources, the time complexity of an acquisition algorithm will affect the start-up speed of the software receiver. In a traditional hardware receiver, signal acquisition is achieved by an Application Specific Integrated Circuit (ASIC), whose computing speed is high. In a software receiver, signal acquisition is achieved by software running in the Central Processing Unit (CPU). It has a more flexible structure but a slower speed (Liu et al., Reference Liu, Jin and Qin2011). Therefore, a software receiver will become more useful than a hardware receiver if we can reduce the time complexity of the acquisition algorithm.
The K-PCA algorithm decomposes the N-point FFT to K parts. The first issue to consider is the number of points of the FFT when the acquisition algorithm is implemented on the hardware platform. If the algorithm uses the Cooley-Tukey radix-2 FFT structure, the FFT input vector needs to be fulfilled up to a power of two with zeros. If the split-radix FFT structure is chosen, this issue should not cause any concern (Tamazin et al., Reference Tamazin, Noureldin, Korenberg and Massound2016; Molino et al., Reference Molino, Girau, Nicola, Fantino and Pini2008). We will analyse the computational load in these two situations in the following parts.
5.1. Time complexity analysis with radix-2 FFT structure
The number of points the FFT needs to evaluate in the conventional PCA algorithm and K-PCA algorithm are shown in Equations (23) and (24) respectively:
$$\hbox{Length of FFT blocks in the PCA algorithm:} M = 2^{ceil(\log_{2}^{N})}$$
$$\hbox{Length of FFT blocks in the} K\hbox{-PCA algorithm:} F = 2^{ceil(\log_{2}^{N/K})}$$where N is the number of sampling points to be processed, and ceil represents the function that rounds the number up to the nearest integer.
In a common scenario, the Doppler frequency ranges between −10kHz to 10kHz. Assuming the scan step is 1 kHz, the acquisition algorithm needs to search for the carrier frequency offset 21 times. If the frequency of the two signals (received signal and local carrier) is different, the K-PCA algorithm only needs to compute the result of y 0(n). While the frequency of the two signals is the same, y i(n)(i∈[0,K−1]) will be calculated. The computational complexity of different algorithms is listed in Table 1.
Computational complexity comparison.
The information related to Table 1 is as follows:
Step 1: Strip the carrier off the received signal.
Step 2: Apply FFT to x(n) or x i(n)(i∈[0,K−1]). Due to the local C B2I code, h(n) is fixed during the algorithm operation process, the H*(k) and H i*(k)(i∈[0,K−1]) are stored in memory without being calculated.
Step 3: Multiply X(k) with H*(k) to obtain Y(k) or combine X i(k)(i∈[0,K−1]) with H i*(k)(i∈[0,K−1]) together to obtain Y i(k)(i∈[0,K−1]).
Step 4: Apply IFFT with the results obtained in Step 3.
MP: Multiplication operation.
AP: Addition operation.
With the increase of sampling points, the growing trend of computational load with different K values is shown in Figure 7.
Comparison of computational load.
If the K-PCA algorithm has the lowest computational load, we call the corresponding K the optimal K value. For example, if the computational load of 3-PCA is lowest, the optimal K is three. Provided that the sampling frequency and the number of sampling points are certain, we can choose the optimal K value to achieve the lowest computational complexity. In Figure 7, the optimal K value changes when the computational load has a step change. The parameter that can cause the step change of the computational load is F. Next, we will analyse its relationship with the optimal K value.
In Equation (24), 
$F=2^{ceil(\log_{2}^{N/K})}$, and by using t K to represent the exponential part ceil (log2N/K), then 
$F=2^{t_{K}}$. The change of F is related to the parameter t K, so we will analyse the step change of t K.
$$t_{K}=ceil(\log_{2}^{N/K})=ceil(\log_{2}^{N}-\log_{2}^{K})$$log2N and log2K can be presented as the forms in Equation (26), where D N and D K represent the decimal part of log2N and log2K respectively, and I N and I K represent the integer part.
$$\log_{2}^{N}= \underbrace{I_{N}}_{{{\rm integer}\atop{\rm part}}}\cdot \underbrace{D_{N}}_{{{\rm decimal}\atop{\rm part}}}\quad \log_{2}^{K}=\underbrace{I_{K}}_{{{\rm integral}\atop{\rm part}}}\cdot \underbrace{D_{K}}_{{{\rm decimal}\atop{\rm part}}}$$If D N>D K, t K=ceil(log2N−log2K)=I N−I K+1; otherwise, t K=ceil(log2N−log2K)=I N−I K.
When D N=D K, t K will change suddenly and F will have a step change. Therefore, the optimal K value is related to the values of D N and D K. The length of the FFT block in different situations is shown in Table 2.
Length of FFT block.
By analysing Table 2 and Figure 7, we can see that the optimal K value can be determined by the relationship between D N and D K. If 
$D_{K_{1}} \lt D_{N} \leq D_{K_{2}}$, the optimal K value is K 2. Since 
$D_{2K_{2}} = D_{K_{2}}$, 2K 2 is also the optimal K value. However, as K increases, the algorithm will become more and more complex. Therefore, the K value should be selected according to the actual situation.
In summary, the optimal K value selection method is: 1. Calculate the value of l, which is the number of sampling points contained in half a C B2I code chip. 2. Compare the values of D N and D K. If 
$D_{K_{1}} \lt D_{N}\leq D_{K_{2}}$, the optimal K value is 2kK 2 (
$k \in [0, \log_{2}^{l}], k$ is an integer). If there is more than one optimal value, the final optimal K value can be selected according to the structure complexity of the algorithm. As shown in Figure 8, the computational load can be reduced by as much as about 50% with an appropriate K value.
The reduction ratio of computational load.
5.2. Time complexity analysis with split-radix FFT structure
There is no need to consider the zero-padding problem with the use of a split-radix FFT structure, and the computational load will also change. The Doppler frequency range is still set to −10 kHz to 10 kHz and the scan step is 1 kHz. The acquisition algorithm needs to search for the carrier frequency offset 21 times. The acquisition process is the same as the process shown in Figure 5. So, in the K-PCA algorithm, y 0(n) should be computed 21 times, and y i(n)(i∈[1,K−1]) only needs to be computed once. In the PCA algorithm, y(n) should be computed 21 times. The computational complexities of the PCA algorithm and K-PCA algorithm are listed in Table 3. With the number of sampling points varying from 1×103 to 80×103, the trend of the computational load of different algorithms is shown in Figure 9. Compared with the PCA algorithm, the K-PCA algorithm exhibits a lower computational burden. Figure 10 shows the reduction ratio of computational load with different numbers of sampling points. Since the reduction ratio is proportional to the value of K, a larger K value can be chosen when acquisition accuracy is not affected.
Multiplication computational load of algorithms with different K values.
Reduction ratio of computational load of algorithms with different K values.
Computational complexity of two algorithms with split-radix FFT structure.
6. SIMULATION RESULTS
To validate the performance of the K-PCA algorithm, we applied the algorithm to capture the B2 signal generated by MATLAB, which has been filtered, down-converted to an Intermediate Frequency (IF) and digitised to a digital signal. The related parameters are shown in Table 4.
Parameters of BeiDou intermediate frequency signal.
From Table 4, we can see that the signal period is 1 ms and sampling frequency is 18·073 MHz, so the number of sampling points to be processed is 18,073. The number of sampling points contained in a half C B2I code chip is four or five, so K should not exceed four. As log218073=14·1415, log23=1·585 and D N<D 3, the optimal k value should be three. Hence, we chose K=3.
In this simulation process, a two-dimensional search for carrier frequency offset and code shift is carried out. The frequency offset range is set between −10 kHz and 10 kHz and is scanned with a 1 kHz step. The frequency offset is set to 2 kHz. The code phase is located at the 13,665th sampling point. If the frequencies of the received signal and local carrier are different, the result of y 0(n) is similar to Figure 11(a). Accordingly, the frequency of the local carrier should be changed and the acquisition process repeated until a correlation peak appears. If the peak value in y 0(n) is much larger than the others, as shown in Figure 11(b). This indicates that the search for carrier frequency offset and a rough estimation of code phase have been accomplished. At this stage, two methods can be used to find the final code phase. Method 1: As illustrated in Figure 5, the results of y 1(n) and y 2(n) are computed and their peak values are compared. If the largest peak value belongs to y i(n) and its coordinate is (m, R), the final code phase is located at the (3×m+i)th sampling point. Method 2: Use the methods shown in Figure 6 and Equation (22) to compute the correlation results in the sampling points near the peak position of y 0(n). The code phase can be determined according to the position of the largest correlation result. The results of y 1(n) and y 2(n) in Method 1 are displayed in Figure 12, where the peak value coordinates of y 0(n), y 1(n) and y 2(n) are (4555, 1·253×10−4), (4554, 1·182×10−4) and (4554, 1·242×10−4) respectively. So, the final code phase is located at the (4555×3+0=13663)th sampling point. Figure 13 shows the results of y K(n) obtained with Method 2. From the figure, we can see that the correlation result at the 13,665th sampling point is the largest. Accordingly, the code phase is located at the 13,665th sampling point. Obviously, the final code phase obtained by the two methods is the same as the pre-set value, so both methods can find the final code phase accurately. The Doppler frequency offset results are shown in Figure 14, and the frequency offset is 2 kHz which is the same as the pre-set value. In conclusion, the K-PCA algorithm can successfully estimate the Doppler shift and the code phase of the B2 signal.
Results of y 0(n) in different local carrier frequency situations.
Results of y 1(n) and y 2(n).
Results of y K(n).
Doppler frequency offset map.
From Section 5.1, we know that the value of K is related to the number of sampling points to be processed. When the signal sampling frequency is changed from 18·073 MHz to 15·491 MHz, the number of sampling points that need to be processed in one period is 15,491. The number of sampling points contained in half a C B2I code chip is three to four, so K should not exceed three. As log215491=13·9191, log23=1·713 and D 3<D 15491<D 2, (Define D2=1) the optimal K value should be two. Hence, we choose K=2. In this simulation, the C/N 0 of the signal is 40 dB·HZ. The code phase and the Doppler frequency are set at the 1,446th chip and 2 kHz. The acquisition results of 2-PCA are as shown in Figures 15–18.
Results of y 0(n) in different local carrier frequency situations.
Results of y 1(n).
Results of y K(n).
Doppler frequency offset map.
When the frequencies of the received signal and local carrier are different, the result of y 0(n) is as Figure 15(a). If the received signal and local carrier frequencies are the same, the result of y 0(n) is shown in Figure 15(b). The peak value in y 0(n) is much larger than the others. Then the result of y 1(n) and y K(n) should be calculated. The final code phase can be determined by comparing the peak values of y 0(n) and y 1(n), or by comparing the values of y K(n). Figures 16 and 17 show the results of y 1(n) and y K(n). The peak value coordinates of y 0(n), y 1(n) are (5471, 1·034×10−4), (5476, 9·7491×10−5) respectively. So, the final code phase is located at the (5471×2)/15491×2046≈1446th chip. From Figure 3, we can see that the peak value of y K(n) is located at the 10,952nd sampling point. So, the final code phase is located at the 10953/15491×2046≈1446th chip. The two methods have the same acquisition results. The Doppler frequency offset results are shown in Figure 18, and the frequency offset is 2 kHz, which is the same as the pre-set value.
For the SSA algorithm, when the code phase search step changes, its acquisition results still obey the Rician distribution. Compared with the PCA algorithm, the acquisition results of K-PCA are still the cross-correlation results between x(n) and h(n), and the difference is that h(n) shifts by three bits in every cross-correlation operation. The code phase search step changes, but the acquisition results of K-PCA should still obey the Rician distribution. The threshold values of PCA and K-PCA can be calculated according to Equation (27). In this equation, P fa is the probabilities of false alarm, σ is the root value of noise power. σ2 can be calculated according to Equation (28), where V(n) is the acquisition value in one search bin and M is the number of search bins.
$$V_{T} =\sigma \sqrt{-2\ln(P_{fa})}$$
$$\sigma^2 = {1 \over 2M} \sum_{n=1}^M V^{2}(n)$$To assess the performance of the K-PCA, Monte Carlo simulation is performed to the Receiver Operating Characteristic (ROC) curve (Yang et al., Reference Yang, Jin, Huang and Qin2014; Reference Yang, Jin, Huang and Qin2016). ROC represents the detection probability under different probabilities of false alarm. Figure 19 shows ROC curves of PCA and 3-PCA. From Figure 19 we can see that the detection probabilities under the same P fa of 3-PCA is slightly lower than that of PCA. This is because the code search steps of the two algorithms are different, and the peak value of the 3-PCA algorithm is slightly lower than that of the PCA algorithm.
ROC comparisons between PCA and 3-PCA.
7. CONCLUSION
A novel acquisition algorithm is proposed in this paper. Unlike the PCA algorithm, which obtains correlation results at every sampling point, our algorithm does not evaluate correlation results at every sampling point but at intervals of K points. The K value can be arbitrarily selected according to the sampling frequency and the number of sampling points to be processed. If the K value is optimal, the computational load can be reduced by about 50%; otherwise, its detection probability is slightly lower than that of the PCA algorithm, under the same probability of false alarm. Therefore, this algorithm can be used to achieve fast acquisition when the power of the signal is weak and the sampling frequency is high or the number of sampling points to be processed is large.
FINANCIAL SUPPORT
This work was supported by the National Natural Science Foundation of China (51505221).
