Comprehensive Analysis of Acquisition Time for a Multi-Constellation and Multi-Frequency GNSS Receiver at GEO Altitude

3.2.1 Cell Level

The acquisition system compares the final correlation value, z(l_c, l_f), to a predefined threshold γ to determine the presence of the satellite signal. Therefore, the acquisition system can be examined through a hypothesis test on z(l_c, l_f). The null hypothesis (H₀) states that the navigation signal from a specific satellite is absent and the correlation function is created solely by noise. Consequently, the correlated value follows a central χ² distribution with 2MK degrees of freedom:

$P\left(z\mid H_0\right)\sim {\chi }^2\left(2MK\right)$ 8

Here, M represents the number of data or pilot components used to form a correlation within a single cell. In other words, when only one of the data or pilot components is used for acquisition, M = 1, whereas if both components are used, resulting in a combined acquisition, M = 2. The degrees of freedom of the χ² distribution are determined by the number of independent Gaussian variables involved in the distribution. As shown in Equation (5), y_X is a complex variable comprising two Gaussian components in its real and imaginary parts. Furthermore, as shown in Equation (7), each y_X is summed in proportion to the number of components used in the correlation (i.e., M) and then summed over the noncoherent accumulation count (i.e., K), resulting in 2MK degrees of freedom.

The alternative hypothesis (H₁) states that the navigation signal from a specific satellite is present. The correlation between the signals shifts the distribution away from the center, creating a non-central χ² distribution with the same degrees of freedom as H₀.H₁ is distributed as follows:

$P\left(z\mid H_1\right)\sim {\chi }^2(2MK,\beta )$ 9

with the non-centrality parameter β = 2KTC/N_{0, eff}. Here, C/N_{0, eff} represents the effective C/N₀, computed as follows:

$C/N_{0,\text{eff}}=\frac{C·\Delta C_y}{N_0}=\frac{\left(C_d+C_p\right)·\Delta C_y}{N_0}$ 10

This term accounts for power losses incurred during the acquisition process relative to the target C/N₀. The factor ΔC_y accounts for signal power losses resulting from mismatches in code and Doppler estimates owing to bin resolution limitations during the cross-ambiguity function computation (Equation (5)), expressed as a ratio relative to unity. Additionally, power losses owing to the utilization of only one component (data or pilot) from signals in which power is split between components are also considered.

Because the values of z(l_c, l_f) exceeding the threshold γ are considered to indicate the presence of a signal within the corresponding cell, the cell-level detection and false alarm probabilities can be expressed as integrals of the respective probability distributions (Equations (8) and (9)) over the range from γ as follows:

$P_d\left(\gamma \right)={\int }_{\gamma }^{\infty }P\left(z\mid H_1\right)dz$ 11

$P_{\text{fa}}\left(\gamma \right)={\int }_{\gamma }^{\infty }P\left(z\mid H_0\right)dz$ 12

3.2.2 Decision Level

Although the probabilities at the cell level are fundamental in determining the performance of the acquisition system, decisions are actually made over the entire search space (Borio et al., 2008; Musumeci et al., 2014), which has slightly different characteristics than the single cell level. Therefore, probability models should be extended to the decision level to realistically assess the acquisition system.

Decision-level probabilities depend on the searching strategy, which defines a method for visiting cells and making decisions. In this analysis, we assume that the GNSS receiver employs the maximum search strategy, which generates a cross-ambiguity function for the entire search space and makes a decision based on the maximum value. In other words, the bin with the maximum value in the search space is identified and compared with the threshold. If the satellite is determined to be visible, the estimates corresponding to that bin are assigned as the final estimates $\hat{\tau }$ and ${\hat{f}}_D$ as follows:

$\begin{array}{cc}\left(l_{c,\max },l_{f,\max }\right)& =\underset{\left(l_c,l_f\right)}{\arg \max }z\left(l_c,l_f\right)\\ z_{\max }& =z\left(l_{c,\max },l_{f,\max }\right)\\ \left(\hat{\tau },{\hat{f}}_D\right)& =\left({\hat{\tau }}_{l_{c,\max }},{\hat{f}}_{D,l_{f,\max }}\right), \text{if} z_{\max }\ge \gamma \end{array}$ 13

where l_{c, max} and l_{f, max} denote the indices of the code and Doppler bins yielding the maximum cross-ambiguity function value and z_max represents the corresponding correlation amplitude.

The hypotheses described in Equations (8) and (9) are confined to a single cell within the entire search space, implying that any cell with incorrect estimates ${\hat{\tau }}_{l_c}$ and ${\hat{f}}_{D,l_f}$ is assumed to contain only noise. This approach presumes that all signal power is concentrated in the cell corresponding to the exact estimates. However, as highlighted by Geiger and Vogel (2013), this assumption may not hold when the code and Doppler bin sizes are small, as the peak of the cross-ambiguity function may extend over multiple bins. Nevertheless, because this paper assumes the use of the maximum threshold crossing method, the impact of neighboring bins containing partial signal power is anticipated to be negligible. Moreover, this assumption facilitates a simplified probabilistic analysis of the acquisition system at the decision level.

Extending the assessment to the decision level introduces additional statistical conditions, which are depicted in Figure 3. The vertical line at the center delineates the presence or absence of a signal within the search space, whereas the horizontal line at the center distinguishes whether the maximum correlation value, z_max, exceeds the threshold γ. For cases in which the signal is absent (left of the vertical line), if z_max surpasses γ, a false alarm (absent) occurs. Conversely, when the signal is present (right of the vertical line) but z_max fails to exceed γ, this constitutes a missed detection. When the signal is present and z_max exceeds γ, the final code and Doppler estimation errors, $\Delta \tau =\tau -\hat{\tau }$ and $\Delta f_D=f_D-{\hat{f}}_D$ , are compared against the tolerable maximum error thresholds, Δτ_th and Δf_D,th, respectively, to distinguish between two scenarios:

1. Detection: If both |Δτ|≤Δτ_th and |Δf_D|≤Δf_D,th, the estimated values are considered sufficiently accurate, corresponding to a correct detection of the signal.

2. False alarm (present): If either |Δτ| > Δτ_th or |Δf_D| > Δf_D,th, the estimates are deemed inaccurate, resulting in a misclassification as a false alarm under the presence of a signal.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Statistical conditions for signal acquisition at the decision level

The thresholds Δτ_th and Δf_D,th are determined based on the pull-in range of the delay-locked and frequency-locked loops, respectively, to ensure that the tracking loops can maintain lock following acquisition.

Borio et al. (2008) developed probability models for decision-level acquisition events, which are utilized in this paper. Detection at the decision level using the maximum search strategy indicates that the maximum value exceeds the threshold while simultaneously yielding correct final estimates, resulting in the following probability:

$\begin{array}{cc}{P}_D\left(\gamma \right)& =P\left[\left(|\Delta \tau |\le \Delta {\tau }_{\text{th}}\right)\cap \left(\left|\Delta f_D\right|\le \Delta f_{D,\text{th}}\right)\mid z_{\max }\ge \gamma \right]\\ & ={\int }_{\gamma }^{\infty }{\left[1-P_{\text{fa}}\left(z\right)\right]}^{N_s-1}P\left(z\mid H_1\right)dz\end{array}$ 14

where N_s denotes the number of code-Doppler bins within the search space. Missed detection refers to a scenario in which the signal exists within the search space, yet no cell exceeds the threshold. The probability of this occurrence is given by the following:

$\begin{array}{cc}{P}_{\text{MD}}\left(\gamma \right)& =P\left[z_{\max }<\gamma \right]\\ & ={\left[1-P_{\text{fa}}\left(\gamma \right)\right]}^{N_s-1}{\int }_0^{\gamma }P\left(z\mid H_1\right)dz\end{array}$ 15

As demonstrated, a false alarm at the decision level can be classified based on whether a signal is present within the search space. A false alarm in the presence of a signal occurs when the maximum value exceeds the threshold, indicating satellite visibility, but the determined code and Doppler estimates are incorrect. As depicted in Figure 3, this event, along with detection and missed detection, constitutes a complete event when a signal is present. Therefore, the probability for this scenario is calculated as follows:

$\begin{array}{cc}{P}_{\text{FA}}^p\left(\gamma \right)& =P\left[\left(|\Delta \tau |>\Delta {\tau }_{\text{th}}\right)\cup \left(\left|\Delta f_D\right|>\Delta f_{D,\text{th}}\right)\mid z_{\max }\ge \gamma \right]\\ & =1-P_D\left(\gamma \right)-P_{\text{MD}}\left(\gamma \right)\end{array}$ 16

When employing the maximum search strategy, a false alarm in the absence of a signal occurs when the maximum value exceeds the threshold despite no signal being present. The probability of this occurrence is expressed as follows:

$P_{\text{FA}}^a\left(\gamma \right)=1-{\left[1-P_{\text{fa}}\left(\gamma \right)\right]}^{N_s}$ 17

3.3.1 Bin Size and Power Loss

The bin sizes are selected to be sufficiently small such that the estimation errors Δτ and Δf_D are smaller than the pull-in range of the discriminator and the signal power loss ΔC_y during the coherent correlation process does not exceed a certain threshold. However, if the bin sizes are too small, N_s increases, which reduces P_D(γ), increases $P_{\text{FA}}^a\left(\gamma \right)$ , and may also result in longer acquisition times. In this paper, the code bin size (b_c) is set to 0.1 chips, assuming the use of a narrow correlator to enhance code tracking sensitivity, except for the L5 band signals. As the subsequent results show, most L5 band signals use wideband signal structures, which would otherwise result in an excessively large number of cells to be searched in the code domain. Therefore, for L5 band signals, a slightly compromised code resolution of 0.5 chips is used to derive a reasonable MAT while preventing an impractically large search space. The Doppler bin size is determined in a similar manner as follows:

$b_f=\frac{1}{2T}$ 20

Because this paper focuses on the mean (i.e., expected) acquisition time, ΔC_y is calculated as the expected power loss of the cross-ambiguity function from Equation (5):

$\Delta C_y=E\left[\frac{R_X^2(\Delta \tau ){\mathrm{sinc}}^2\left(\pi T\Delta f_D\right)}{R_X^2\left(0\right){\mathrm{sinc}}^2\left(0\right)}\right]$ 21

where $E(·)$ denotes the expectation operator applied to the random variables Δτ and Δf_D, which are assumed to be uniformly distributed within their respective bins:

$\Delta \tau \sim U\left(-\frac{b_c}{2},\frac{b_c}{2}\right)$ 22

$\Delta f_D\sim U\left(-\frac{b_f}{2},\frac{b_f}{2}\right)$ 23

Because Δτ and Δf_D are independent and $R_X^2\left(0\right)={\mathrm{sinc}}^2\left(0\right)=1$ , Equation (21) can be expressed as follows:

$\Delta C_y=E\left[R_X^2(\Delta \tau )\right]E\left[{\mathrm{sinc}}^2\left(\pi T\Delta f_D\right)\right]$ 24

By analytically deriving ΔC_y for each signal using the given equations, the correlation loss corresponding to the bin size design can be analyzed.

The shape of R_X(Δτ) varies depending on the characteristics of the GNSS signal, such as modulation type and chipping rate. In this paper, the closed-form time-domain autocorrelation function expressions derived by Sousa and Nunes (2013) are utilized. For binary phase shift keying (BPSK) signals, a simple triangular model is employed:

$R_{X,\text{BPSK}}(\Delta \tau )=\{\begin{array}{cc}1-\frac{|\Delta \tau |}{T_{\text{chip}}},& |\Delta \tau |<T_{\text{chip}}\\ 0,& \text{otherwise}\end{array}$ 25

where T_chip represents the chip duration of the primary code. The remaining binary offset carrier (BOC)-family signals were approximated as sine-phase BOC (1,1) for simplification of the analysis, using the following expression:

$\begin{array}{cc}\kappa & =\lceil \frac{2|\Delta \tau |}{T_{\text{chip}}}\rceil \\ R_{X,\text{BOC}}(\Delta \tau )& ={(-1)}^{\kappa }\left[{\kappa }^2-3\kappa +1+(5-2\kappa )\frac{|\Delta \tau |}{T_{\text{chip}}}\right]\end{array}$ 26

The multiplexed BOC modulation used in civilian GNSS signals allocates most of the signal power to the sine-phase BOC (1,1), resulting in an autocorrelation function shape similar to that of BOC (1,1). Therefore, applying the above approximation is likely to have a negligible effect on the expected values.

3.3.2 Uncertainty

As shown in Equation (4), this paper utilizes the primary code lengthened by the secondary code to increase the coherent integration time T. However, directly applying this method would require simultaneous searching of the phases of the primary and secondary codes during the coherent integration stage, resulting in excessive code uncertainty and rendering the approach impractical. Therefore, in practice, the phase search processes for each code are assumed to be separated and sequentially conducted during signal acquisition. This assumption can be implemented using a straightforward exhaustive method, which will be described in detail in Section 3.5.4. Under this assumption, the code-domain uncertainty corresponds to the number of chips in the primary code period and can be expressed as follows:

$u_c=f_{\text{code}}{T}_{\text{code}}$ 27

where f_code is the chipping rate of the primary code and T_code is the period of the primary code.

The uncertainty in the Doppler domain (u_f) depends on the receiver starting option: cold start or warm start. During a cold start, the receiver lacks any a priori information and must search for all possible Doppler shifts. With the assumption of negligible receiver dynamics, which is a reasonable assumption for GEO applications, the uncertainty is expressed as follows:

$u_{f,\text{cold}}=2\left(f_{D,\text{max}}+f_{\text{osc}}\right)$ 28

where f_D,max is the expected maximum Doppler shift. f_osc is the frequency error induced by the local oscillator, which can be expressed as follows:

$f_{\text{osc}}=A_{\text{osc}}{f}_c$ 29

where A_osc is the accuracy of the oscillator and f_c is the carrier frequency.

To reduce the frequency search space, the receiver can utilize a warm start, which provides approximate position and time information to the receiver. By predicting the approximate Doppler shift, the uncertainty reduces to the following (Morton et al., 2020):

$u_{f,\text{warm}}=2\left(f_{\text{time}}+f_{\text{pos}}+f_{\text{osc}}\right)$ 30

where f_time and f_pos represent the Doppler prediction error caused by inaccurate time and position aiding information, respectively. f_time can be calculated from the maximum Doppler rate ${\dot{f}}_{D,\text{max}}$ and time aiding error δt as follows:

$f_{\text{time}}={\dot{f}}_{D,\text{max}}\delta t$ 31

f_pos has the following maximum value:

$f_{\text{pos}}\le \frac{v_{s,\max }\delta x}{{\lambda }_c{\rho }_{\min }}$ 32

where v_s,max is the maximum speed of the GNSS satellite, δx represents the aided receiver position error, λ_c is the carrier wavelength, and ρ_min is the minimum geometric range between the satellite and receiver.

3.5.1 Bit Inversion and Number of FFT Points

During the coherent correlation process, if a phase inversion occurs owing to the navigation message bits or the secondary code while the correlation is being performed, the correlation gain will be attenuated. If this inversion occurs in the middle of the correlation interval, it may completely nullify the correlation gain. Therefore, such effects must be carefully considered during the acquisition process. GNSS signals can be broadly classified as Type 1,2, or 3, each necessitating distinct algorithms tailored to their unique characteristics. Type 1 signals, such as GPS L1 C/A, are legacy signals in which multiple primary codes exist within a single navigation message bit. In contrast, Types 2 and 3 are modernized signals. Type 2 signals, such as GPS L2C, have primary code periods that match the length of navigation message bits. Type 3 signals, such as GPS L5, contain multiple primary codes within a single bit, but their code sequences are extended by a secondary code.

For Type 1 signals, the loss caused by phase inversion during coherent correlation can be mitigated by using the half-bit method proposed by Psiaki (2001). In other words, if the bit is split into two halves and coherent correlation is performed separately on each half, at least one of the regions will avoid experiencing phase inversion. However, for Type 2 and 3 signals, splitting the bit is not feasible, rendering this method inapplicable. Therefore, for such signals, the double-length zero-padding technique proposed by Yang (2001) is employed. With this technique, circular correlation is performed over a length corresponding to two periods of the primary code. The code replica is constructed from one period of the primary code concatenated with a zero vector of equal length. By discarding the last half-block of the correlation result, this method ensures that the correlation is conducted with an intact primary code within the two-period signal block, unaffected by bit transitions, thereby mitigating losses associated with bit inversions. However, as this approach doubles the required number of FFT points, it is computationally disadvantageous. Therefore, for Type 1 signals in which the half-bit method can be applied, we assume that the half-bit method is utilized instead. As a result, zero-padding is applied to the code replica generation only for Type 2 and 3 signals as follows:

${\tilde{c}}_X\left(k\right)=\{\begin{array}{cc}{c}_{1,X}\left(k{T}_s\right),& 0\le k<f_s{T}_{\text{code}}\\ 0,& f_s{T}_{\text{code}}\le k<2f_s{T}_{\text{code}}\end{array}\begin{array}{c}(\text{All} \text{Types})\\ (\text{Type} 2 \text{and} 3)\end{array}$ 35

To account for the change in the number of FFT points due to zero-padding, we define a variable N_zp such that N_zp = 2 when zero-padding is applied (as for Type 2 and 3 signals) and N_zp = 1 when zero-padding is not applied (as for Type 1 signals). Consequently, the duration for the FFT is given by the following:

$T_{\text{FFT}}=T_{\text{code}}{N}_{\text{zp}}$ 36

and the corresponding number of FFT points is calculated as follows:

$N_{\text{FFT}}=f_s{T}_{\text{FFT}}$ 37

3.5.2 Number of Forward FFTs

The code replica is generated during the receiver’s initialization phase and converted into the frequency domain with respect to the frequency sample index ω as follows:

${\tilde{C}}_X\left(\omega \right)={ℱ}^{*}\left[{\tilde{c}}_X\left(k\right)\right]$ 38

The code replica is then stored for subsequent reuse in each acquisition process. This code replica generation process must be performed only for the number of components, M, used in the acquisition. Therefore, the number of forward FFT operations required for code replica generation is constant and equal to M. In contrast, the carrier wipe-off for the n_code-th code sequence period and l_f-th Doppler bin, represented as:

$\tilde{r}\left(k,n_{\text{code}},l_f\right)=r\left(\frac{N_{\text{FFT}}}{N_{\text{zp}}}{n}_{\text{code}}+k\right)\exp \left[-j2\pi \left(f_{\text{IF}}+{\hat{f}}_{D,l_f}\right)k{T}_s\right]$ 39

and the corresponding FFT:

$\tilde{R}\left(\omega ,n_{\text{code}},l_f\right)=ℱ\left[\tilde{r}\left(k,n_{\text{code}},l_f\right)\right]$ 40

must be performed each time the target signal block for correlation changes or a new carrier is generated owing to a change in the Doppler bin. Here, n_code represents the period index of the primary code sequence.

Psiaki (2001) introduced a method to reduce the number of required FFT computations by employing circular shifts in the frequency domain. This method stores the carrier-wiped received signal block in the frequency domain and alters the Doppler value under test by applying a circular shift and phase correction as follows:

$\tilde{R}\left(\omega ,n_{\text{code}},l_f+f_{\text{res}}\right)=\tilde{R}\left(\omega -1,n_{\text{code}},l_f\right)\exp \left(-j2\pi f_{\text{res}}{T}_s\right)$ 41

which enables the results for other Doppler bins to be obtained without requiring additional FFT operations. Here, f_res represents the frequency resolution of the frequency domain representation obtained through the FFT and is calculated by dividing the two-sided sample bandwidth by the number of FFT points:

$f_{\text{res}}=\frac{f_s}{N_{\text{FFT}}}=\frac{1}{T_{\text{FFT}}}$ 42

Therefore, by applying a circular shift and phase correction to a single sample in the frequency domain, an effect equivalent to a Doppler value shift of f_res can be achieved. Thus, the number of forward FFTs that must be performed per received signal block, which would otherwise require N_f operations, is dramatically reduced.

In cases such as GPS L1 C/A, where T_code is short and N_zp = 1, the frequency resolution achieved by a circular shift, f_res, may be excessively coarse relative to the desired Doppler bin size, b_f. In such scenarios, when generating the initial $\tilde{R}$ , multiple carriers can be generated more densely to evenly divide f_res, and a circular shift is applied to each $\tilde{R}$ . The ratio of Doppler resolution is defined as follows:

$r_f=\frac{f_{\text{res}}}{b_f}=\frac{2T}{T_{\text{FFT}}}$ 43

Based on this value, the number of carriers that must be initially generated, N_carr, is determined. In this paper, the maximum number of carriers generated per f_res is limited to four:

$N_{\text{carr}}=\min \left(r_f,4\right)$ 44

Here, the upper limit of four carriers is imposed to ensure computational efficiency by restricting the number of carrier wipe-off operations and FFT calculations. Therefore, for each received signal block, N_carr instances of $\tilde{R}$ are initially generated simultaneously as follows:

$\begin{array}{c}\tilde{R}\left[l_f+f_{\text{res}}\left(\frac{0}{N_{\text{carr}}}\right)\right],\\ \tilde{R}\left[l_f+f_{\text{res}}\left(\frac{1}{N_{\text{carr}}}\right)\right], \dots ,\tilde{R}\left[l_f+f_{\text{res}}\left(\frac{N_{\text{carr}}-1}{N_{\text{carr}}}\right)\right]\end{array}$ 45

where the index information ω and n_code are temporarily omitted for simplicity in this example. Subsequently, circular shifts for the n_shift-th shift are applied to each instance:

$\begin{array}{c}\tilde{R}\left[l_f+f_{\text{res}}\left(n_{\text{shift}}+\frac{0}{N_{\text{carr}}}\right)\right],\tilde{R}\left[l_f+f_{\text{res}}\left(n_{\text{shift}}+\frac{1}{N_{\text{carr}}}\right)\right],\\ \dots ,\tilde{R}\left[l_f+f_{\text{res}}\left(n_{\text{shift}}+\frac{N_{\text{carr}}-1}{N_{\text{carr}}}\right)\right]\end{array}$ 46

enabling the Doppler bin search to be conducted efficiently. The number of required frequency circular shifts is determined as follows:

$N_{\text{shift}}=\lceil \frac{u_f}{f_{\text{res}}}\rceil$ 47

The number of forward FFT operations required for each received signal block is only N_carr, as a significant number of forward FFTs for the Doppler search are replaced by frequency circular shifts.

3.5.3 Number of IFFTs

The coherent correlation is completed by taking the IFFT of the product in the frequency domain, as follows:

${\tilde{y}}_X\left(k,n_{\text{code}},l_f\right)={ℱ}^{-1}\left[\tilde{R}\left(\omega ,n_{\text{code}},l_f\right)\odot {\tilde{C}}_X\left(\omega \right)\right]$ 48

By generating N_carr carriers, the achieved frequency resolution becomes f_res/N_carr. However, for Type 1 and 3 signals, where f_res is relatively large compared with the desired Doppler bin size b_f, this resolution may still be insufficient. In such cases, the post-correlation approximation method for nearby Doppler bins, as proposed by Psiaki (2001), can be effectively applied to resolve the problem. This method achieves an interpolation effect between sparse Doppler bin intervals by rotating the phase of the IFFT result, expressed by the following:

$\begin{array}{c}{\tilde{y}}_X\left(k,n_{\text{code}},l_f+n_{\text{interp}}{b}_f\right)\approx {\tilde{y}}_X\left(k,n_{\text{code}},l_f\right)\\ \exp \left[-j2\pi n_{\text{interp}}{b}_f{T}_{\text{FFT}}\left(n_{\text{code}}+0.5\right)\right]\end{array}$ 49

Thereby, the correlation results are approximated for the neighboring n_interp-th interpolated Doppler bin. Because the number of carriers, N_carr, is limited to a maximum of four, the required number of interpolation steps for each of the N_carr carriers can be derived as follows:

$N_{\text{interp}}=\lceil \frac{r_f}{4}\rceil$ 50

Equation (49) is applicable to Type 1 and 3 signals because they effectively form the correlation results for T by accumulating several correlation results over T_FFT. This process is further explained later in this section. In this context, if the Doppler bin spacing is set too sparsely, Doppler errors increase, causing phase rotation in the exponential term of Equation (5) for each correlation over the T_FFT block. Equation (49) corrects these phase rotations, ensuring that the phase remains continuous across each block, thus enabling coherent accumulation over multiple iterations, even with inaccurate Doppler values. Nevertheless, amplitude loss caused by the sinc function in Equation (5) will still occur, and this effect must be considered in the dwell time design. Therefore, for Type 1 and 3 signals, when interpolation is applied, Equation (24) can be reformulated as follows:

$\Delta C_y=E\left[R_X^2(\Delta \tau )\right]E\left[{\mathrm{sinc}}^2\left(\pi T_{\text{FFT}}\Delta f_{D,\text{coarse}}\right){\mathrm{sinc}}^2\left(\pi T\Delta f_{D,\text{fine}}\right)\right]$ 51

Here, Δf_D,coarse denotes the Doppler error resulting from the coarse Doppler bin interval and is assumed to be uniformly distributed within the bin width:

$\Delta f_{D,\text{coarse}}\sim U\left(-\frac{f_{\text{res}}}{2N_{\text{carr}}},\frac{f_{\text{res}}}{2N_{\text{carr}}}\right)$ 52

similar to Equation (23). Δf_D,fine represents the residual Doppler error after interpolation has been applied with a bin interval of b_f to Δf_D,coarse. This term can be calculated as follows:

$\Delta f_{D,\text{fine}}=\left| \mathrm{mod} \left(\Delta f_{D,\text{coarse}}-\frac{b_f}{2},b_f\right)-\frac{b_f}{2}\right|$ 53

Through this interpolation-based approximation method, the IFFT operation must be performed only for the coarsely spaced Doppler bins corresponding to each frequency circular shift. Because the operation for each component is independent and the frequency circular shift is performed N_shift times, the total number of IFFTs required per received signal block is given by N_carr MN_shift.

Additionally, evaluating Equation (49) requires N_FFT complex multiplications, which are repeated N_interp times for each of the N_carr carriers for each component. Thus, a total of N_FFTN_carrN_interp M complex multiplications are required per received signal block. In general-purpose processors, this operation may incur significant computational time. However, because all of these operations are mutually independent, they can be easily parallelized on dedicated hardware, assuming sufficient memory bandwidth. In this paper, we assume that such parallelization is implemented, and the time required for the computations in Equation (49) is not considered.

3.5.4 Number of Processing Blocks

As previously analyzed, the forward and inverse FFTs are repeatedly performed for each processing block. Therefore, to calculate the total number of FFTs, we must determine the number of processing blocks that fit within the dwell time. For Type 1 and 2 signals, which do not include a secondary code, the coherent correlation is completed by simply summing the correlation results for one primary code period N_coh times, as follows:

${\overline{y}}_X\left(k,l_f\right)=\sum _{n_{\text{code}}=0}^{N_{\text{coh}}-1}{\tilde{y}}_X\left(k,n_{\text{code}},l_f\right)$ 54

For Type 1 signals, the half-bit method is applied, necessitating that Equation (54) be calculated twice: once for the first half and once for the second half of each bit period. To generalize this characteristic for all signals, the processing block count incorporates a factor of 2/ N_zp.

However, for Type 3 signals, which utilize a secondary code, each primary code is modulated by different bits of the secondary code. Consequently, simply summing the results as in Equation (54) may result in correlation gain loss due to bit inversions. Relevant studies have proposed acquisition algorithms that account for the secondary code (Borio, 2011), and various acquisition algorithms have been analyzed in terms of both high sensitivity performance and hardware implementation aspects (Leclère et al., 2017). In this paper, we assume that the phase of the secondary code is determined via a relatively simple brute force approach. To mitigate the impact of the secondary code, we multiply the secondary code bits corresponding to each primary code period and coherently sum the result:

${\overline{y}}_{X,i_s}\left(k,l_f\right)=\sum _{n_{\text{code}}=0}^{L_s-1}{\tilde{y}}_X\left(k,n_{\text{code}},l_f\right){\tilde{c}}_{X,2}\left(n_{\text{code}}-i_s\right)$ 55

where ${\tilde{c}}_{X,2}(·)$ represents the secondary code sequence of component X expressed in bit units and the input bit index is taken modulo the secondary code period L_s, yielding values between 0 and L_s −1. i_s represents the initial bit offset of the secondary code to be multiplied, and for all possible phases of the secondary code, ${\overline{y}}_{X,i_s}\left(k,l_f\right)$ is generated for L_s instances. Finally, the offset that shows the maximum value among the L_s results is selected as the final secondary code offset, ${\hat{i}}_s$ , and the corresponding correlation result for this offset is used as follows:

$\begin{array}{cc}{\hat{i}}_s& =\underset{i_s\in \left[0, L_s\right)}{\text{arg} \text{max}}\left[\underset{\left(k,l_f\right)}{\text{max}} {\overline{y}}_{X,i_s}\left(k, l_f\right)\right]\\ {\overline{y}}_X\left(k, l_f\right)& ={\overline{y}}_{X,{\hat{i}}_s}\left(k, l_f\right)\end{array}$ 56

To generate one instance of ${\overline{y}}_{X,i_s}\left(k,l_f\right)$ , we require L_s instances of ${\tilde{y}}_X\left(k,n_{\text{code}},l_f\right)$ . Although L_s instances of ${\overline{y}}_{X,i_s}\left(k,l_f\right)$ must be generated, in practice, the L_s instances of ${\tilde{y}}_X\left(k,n_{\text{code}},l_f\right)$ can be generated and stored first. Subsequently, for each i_s, Equation (55) can be computed. Therefore, the number of correlation blocks required by this method is L_s.

In summary, the number of processing blocks required to complete the coherent correlation (i.e., the number of FFT–IFFT operations) is generalized for all signal types as 2N_cohL_s/N_zp. These operations must be repeated for the noncoherent accumulation count K; therefore, the final number of processing blocks is 2N_cohL_sK/N_zp.

3.5.5 Calculation Time

If the receiver employs the maximum threshold crossing method, which entails making a decision by comparing the maximum value with the decision threshold, the MAT model can be expressed as follows (Kassabian & Presti, 2012):

$t_{\text{MAT}}=t_{\text{cal}}\left[1+\frac{P_{\text{MD}}\left(\gamma \right)+P_{\text{FA}}^a\left(\gamma \right)}{P_D\left(\gamma \right)}\right]+\frac{t_p{P}_{\text{FA}}^a\left(\gamma \right)}{P_D\left(\gamma \right)}$ 57

where t_cal is the time consumed to form an entire cross-ambiguity function and t_p is the penalty time induced by the false lock.

Here, t_cal is approximated by considering the FFT-related operations, which are the most computationally demanding processes in the calculation of the cross-ambiguity function. As analyzed in the preceding subsections, the number of forward FFTs required to generate code replicas is M, whereas the numbers of forward FFTs and IFFTs required per received signal block are N_carr and N_carr MN_shift, respectively. Considering that the received signal block is processed 2N_cohL_sK/N_zp times and FFT and IFFT operations are computationally equivalent, the approximate value of t_cal is expressed as follows:

$t_{\text{cal}}\approx t_{\text{FFT}}\left[M+\frac{2}{N_{\text{zp}}}{N}_{\text{coh}}{L}_{s}K{N}_{\text{carr}}\left(1+M{N}_{\text{shift}}\right)\right]$ 58

Although the exact computational cost may vary depending on the specific algorithm used, for a basic radix-2 FFT algorithm, the numbers of real additions and multiplications required for a single FFT are approximately 3N_FFTlog₂N_FFT and 2N_FFTlog₂N_FFT, respectively (Heckbert, 1995). Based on this result, the time required for a single FFT can be approximated as follows:

$t_{\text{FFT}}\approx \frac{\alpha }{f_{\text{proc}}}\left(3N_{\text{add}}+2N_{\text{mul}}\right)N_{\text{FFT}} {\log }_2{N}_{\text{FFT}}$ 59

where N_add and N_mul are the numbers of processor cycles for addition and multiplication, respectively, and f_proc is the processor clock speed. α is a coefficient that reflects the unmodeled optimizations of the processor, including processing architectures (e.g., techniques such as pipelining) and the FFT algorithm. Modeling the multi-core and cache characteristics of modern processors, or the advanced FFT algorithm properties, with perfect accuracy is not feasible. Therefore, in this paper, a relatively simple basic model is used to predict t_FFT with appropriate precision, rather than attempting to derive it with full accuracy. However, as will be demonstrated in Section 6, even this simplified model is sufficiently accurate for evaluating the performance of various signals for a specific algorithm and can provide reliable estimates of t_cal and t_MAT.

4.2.1 Earth Blockage and Link Budget

The orbital characteristics of each system were incorporated by applying the nominal almanac parameters presented in Annex D of the Interoperable GNSS SSV booklet (United Nations Office for Outer Space Affairs, 2021). Based on the almanac parameters and time information, the position, velocity, and acceleration of each GNSS satellite were computed using the ephemeris determination algorithm detailed in the GPS interface specification document (Anthony & Kerns, 2022). If the target GNSS satellite was a GEO satellite from the BDS constellation, the satellite information was further corrected via the coordinate transformation algorithm described by the China Satellite Navigation Office (2019). Once the GNSS satellite positions were determined, the line-of-sight vector connecting the receiver (i.e., the GEO satellite) was evaluated to determine whether it intersected a spherical Earth model with a radius of 6,378 km, indicating whether the signal was blocked by the Earth. As noted by Estrada et al. (2024), signals passing through the ionosphere at the GEO induce pseudorange delays and velocity estimation errors, and removing such signals improves positioning accuracy. Hence, in this study, a 1,000-km ionospheric layer was added to the Earth’s radius, and signals passing through this layer were also considered to be blocked. The link budget was calculated for the GNSS satellites that were not blocked. The received signal power can be calculated as follows (Misra & Enge, 2012):

$C=\frac{P_T{G}_T\left({\varphi }_T\right)G_R\left({\varphi }_R\right)}{L_A{L}_R}{\left(\frac{{\lambda }_c}{4\pi \rho }\right)}^2$ 60

where P_T denotes the transmitted signal power and G_T (ϕ_T) and G_R (ϕ_R) represent the gains of the transmitting and receiving antennas relative to an isotropic antenna, respectively, as functions of the off-boresight angles of the transmitter, ϕ_T, and the receiver, ϕ_R. Here, the magnitude of P_T and shape of G_T (ϕ_T) can be determined based on the characteristics of the selected GNSS signal. L_A and L_R are the losses occurring in the atmosphere and at the receiver device, respectively, and ρ is the geometric distance between the satellite and receiver. In this study, L_A was assumed to be 0 dB, because the received signals did not traverse the atmosphere, whereas L_R was set to 2 dB to account for cable/filter and correlation losses.

4.2.2 Antenna Gain Patterns

To accurately reflect the actual environment, we used realistic transmitting and receiving antenna patterns during the link budget calculation, as shown in Figure 6. This figure shows the multi-band transmitting antenna patterns for GPS, Galileo, and QZSS, along with the receiving antenna pattern. All transmitting antenna patterns were obtained by averaging the antenna pattern data provided by the service providers for all satellites and azimuth angles. Consequently, the patterns became functions of the off-boresight angle of the transmitter, ϕ_T, which simplified the simulation and facilitated easier modifications of the pattern for signals from systems whose transmitting antenna patterns are not publicly available. The GPS antenna pattern was based on the L1/L2/L5 band directivity data for GPS Block III satellites, provided by Fischer (2022). The Galileo pattern utilized the recently released expected equivalent isotropic radiated power (EIRP) data for the Galileo E1/E5/E6 bands (Menzione et al., 2024). For QZSS, the directivity data for the L1/L2/L5 bands provided by the Cabinet Office (2023) were used. The receiving antenna pattern was designed to ensure optimal visibility from the GEO, using a pattern similar to the design result provided by Winkler et al. (2017).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 6

Transmitting antenna patterns of GPS, Galileo, and QZSS, along with the receiving antenna pattern used for the geometric simulation

For GPS and QZSS, directivity data were provided, and the antenna gain was obtained by adding a gain correction factor to the directivity (Steigenberger et al., 2018). Under the assumption of negligible antenna losses, the directivity was directly applied as the transmitting antenna gain, G_T (ϕ_T). In these cases, because detailed information about the transmitted signal power, P_T, was not available, a nominal value of 14.3 dBW was used (Misra & Enge, 2012). In contrast, because Galileo provides EIRP data, which combines transmitting power and gain, the pattern data were used as P_T G_T (ϕ_T).

For systems and signals without publicly available antenna patterns, the corresponding band pattern from the GPS antenna was modified in terms of power and beamwidth for use. Hence, the transmission characteristics of each signal, including the minimum received power at the GEO and the beamwidth of the mainlobe, are utilized. The nominal values, obtained from the Interoperable GNSS SSV booklet (United Nations Office for Outer Space Affairs, 2021), are listed in Table 3. The transmitting antenna patterns were modified according to the following procedure:

1. Determine the reference pattern: The GPS antenna pattern corresponding to the frequency band of the target signal is selected as the reference pattern. The nominal minimum received power and the mainlobe beamwidth of the simulation target signal are defined as P_{R, sim} and B_sim, respectively. The nominal characteristics of the GPS signal associated with the reference pattern are set as P_{R, ref} and B_ref, respectively.

2. Adjust the pattern height: The ratio of P_{R, sim} to P_{R, ref} (or their difference on a dB scale) is applied to adjust the height of the reference pattern, reflecting the nominal transmission power difference. Ideally, the transmission power adjustment should be reflected in P_T; however, to simultaneously account for the transmission characteristics within the antenna pattern modification process, we incorporate this adjustment into G_T (ϕ_T), yielding identical results.

3. Scale the beamwidth: The ratio of B_sim to B_ref is used to modify the indexing angle ϕ_T, stretching or compressing the entire pattern horizontally to match the nominal mainlobe beamwidth. This transformation applies to all ϕ_T values, proportionally adjusting the widths of the sidelobes in relation to the mainlobe.

Thus, the modified transmitting antenna gain can be expressed as follows:

$G_{T,\text{mod}}\left({\varphi }_T\right)=\frac{P_{R,\text{sim}}}{P_{R,\text{ref}}}{G}_T\left({\varphi }_T\frac{B_{\text{ref}}}{B_{\text{sim}}}\right)$ 61

By modifying the reference antenna pattern for signals without disclosed patterns, the simulation effectively incorporated the received signal power magnitude and the characteristics of the mainlobe and sidelobes. This approach enabled the prediction of realistic reception environment characteristics for these signals.

4.2.3 Visibility Threshold

The C/N₀ threshold for determining signal visibility in the simulation was set to 20 dB-Hz, which was selected based on prior studies analyzing visibility and navigation performance for GEO satellites (Ji et al., 2021; Lee et al., 2022). With this threshold, an average of 25 visible satellites was predicted for receivers utilizing GPS, GLONASS, Galileo, and BDS at the GEO altitude. Additionally, according to the United Nations Office for Outer Space Affairs (2021), using a 20-dB-Hz C/N₀ threshold ensures that at least one visible satellite is always available for GEO with interoperable GNSS. This result was derived considering only the mainlobe of the transmitting antenna pattern, representing a worst-case scenario compared with the actual reception environment. In actual conditions, where sidelobes are also utilized, the number of visible satellites increases significantly. Therefore, we can anticipate that sufficient visibility is secured with the 20-dB-Hz C/N₀ threshold, thereby reinforcing the validity of the selected threshold. Moreover, the GPS receiver onboard the GOES-R, using a similar threshold of 17 dB-Hz, could acquire and track an average of 11 visible satellites with standalone GPS in flight tests (Winkler et al., 2017). Even with a slight increase in threshold, as assumed in this study, the use of multi-constellation systems further increased the number of visible satellites, thereby confirming the appropriateness of the selected threshold.

4.2.4 Geometric Simulation Results

The parameters obtained through the geometric simulation are presented in Table 4. The values of each signal differed, even within the same system, because each signal used distinct antenna patterns and transmission characteristics, resulting in variations in the visible regions. Furthermore, the orbital characteristics of each system were revealed. As anticipated, inclined geosynchronous orbit (IGSO) and GEO satellites had a significantly larger ρ_min, in the range of 70,000 km, compared with MEO satellites, and they exhibited lower satellite dynamics. As introduced in Section 3.3.2, each parameter was used to determine the cold- and warm-start uncertainty of each signal, which enabled the orbital characteristics, antenna patterns, and transmitted power of each signal to be incorporated into the acquisition time analysis process.