An Analytical Model for BDS B1 Spreading Code Self-Interference Evaluation Considering NH Code Effects

The short spreading code used by the BeiDou Navigation Satellite System (BDS) B1-I or GPS Coarse/Acquistiion (C/A) can cause aggregately undesirable cross-correlation between signals within each single constellation. This GPS-to-GPS or BDS-to-BDS correlation is referred to as self-interference. A GPS C/A code self-interference model is extended to propose a self-interference model for BDS B1, taking into account the unique feature of the B1-I signal transmitted by BDS medium Earth orbit (MEO) and inclined geosynchronous orbit (IGSO) satellites—an extra Neumann-Hoffmann (NH) code. Currently there is no analytical model for BDS self-interference and a simple three parameter analytical model is proposed. The model is developed by calculating the spectral separation coefficient (SSC), converting SSC to equivalent white noise power level, and then using this to calculate effective carrier-to-noise density ratio. Cyclostationarity embedded in the signal offers the proposed model additional accuracy in predicting B1-I self-interference. Hardware simulator data are used to validate the model. Software simulator data are used to show the impact of self-interference on a typical BDS receiver including the finding that self-interference effect is most significant when the differential Doppler between desired and undesired signal is zero. Simulation results show the aggregate noise caused by just two undesirable spreading codes on a single desirable signal could lift the receiver noise floor by 3.83 dB under extreme C/N0 (carrier to noise density ratio) conditions (around 20 dB-Hz). This aggregate noise has the potential to increase code tracking standard deviation by 11.65 m under low C/N0 (15–19 dB-Hz) conditions and should therefore, be avoided for high-sensitivity applications. Although the findings refer to Beidou system, the principle weakness of the short codes illuminated here are valid for other satellite navigation systems.


Introduction
The BeiDou Navigation Satellite System (BDS) is a relatively new member of the Global Navigation Satellite System (GNSS) family. Its global role has recently been recognized by formal acceptance by the International Maritime Organization (IMO) as a key part of its World-Wide Radio-Navigation System (WWRNS). This is in addition to the adoption of 1575.42 MHz for the future B1c component (instead of 1561.098 MHz currently used by B1 signal) for interoperability with other GNSS's. For both mass-market and professional GNSS users, BDS will provide improved accuracy, integrity, continuity and availability. There are several factors that determine a GNSS receiver's tracking sensitivity, one of Figure 1 shows a generic correlator model, whose input is the undesired (interfering) signal, x 1 (k) (t − ∆ k )e j2πf , and the desired (victim) signal is x 2 (t). There is a differential Doppler frequency,

Signal and Correlator Model
i.e., difference between Doppler frequencies, denoted as f, between the undesired and desired signals. These two signals are correlated by mixing and integration over an arbitrary navigation data bit period, [kT b , (k + 1)T b ), where T b = 20 ms and k = 0, 1, 2, .... and carrier tracking as self-interference. The accuracy of the model is ensured by using a twoparameter ACF to characterize a B1-I code either with or without navigation data. This type of cyclostationary analytical model will be validated as an accurate method to evaluate self-interference for short codes, producing excellent agreement with actual receiver observables, using hardware simulator intermediate frequency (IF) data and a software receiver. The Section "Signal and Correlator Model" provides basic assumptions and signal models to be used in the development of the proposed analytical model. Section "SSC and Equivalent White Noise Level" derives analytical models for the variance at the correlator output for three situations: the data bits of the desired and undesired signals are (1) aligned or (2) misaligned by an integer number of code chips or (3) misaligned by a fraction of a code chip. These three models are the basis of the proposed analytical model. Section "Proposed Model for BDS MEO/IGSO B1-I Self-Interference" presents the proposed model based on SSCs derived in the previous section. Section "Model Validation" demonstrates the validity of the model using live BDS data record and a modified software receiver. Section "Results" analyzes typical behavior of B1-I self-interference and its impact on receiver range observables using the developed model. Section "Conclusions" summarizes the expected contributions and presents a framework for further efforts. Figure 1 shows a generic correlator model, whose input is the undesired (interfering) signal, x1 (k) (t − Δk)e j2πf , and the desired (victim) signal is x2(t). There is a differential Doppler frequency, i.e., difference between Doppler frequencies, denoted as f, between the undesired and desired signals. These two signals are correlated by mixing and integration over an arbitrary navigation data bit period, [kTb, (k + 1)Tb), where Tb = 20 ms and k = 0, 1, 2, .... To simplify the problem without loss of generality, we focus on an individual undesired signal including a specific spreading code identified by PRN (Pseudo Random Noise) number, x1(t). This undesired signal can be any one of the undesired signals x1 (k) (t) in the input to the correlator in Figure 1. Therefore, the undesired signal without Doppler or time delay is modelled as:

Signal and Correlator Model
where pTc(t) is a rectangular pulse of width Tc: (2) and cv ∈ {−1, 1} is spreading code chips with a chipping rate 1/Tc = 2.046 MHz. The chip width is Tc; dk is navigation data bits with bit rate 1/Tb = 50 bps, and the bit width is Tb. αu is the NH code adopted for B1-I signal and takes a fixed pattern of (0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0), which spans one navigation data bit. The spreading code here is treated as a random periodic sequence of length 2046 chips, with chips independent from each other. The navigation data bit stream is also treated as a random periodic sequence, but with an infinite length. To simplify the problem without loss of generality, we focus on an individual undesired signal including a specific spreading code identified by PRN (Pseudo Random Noise) number, x 1 (t). This undesired signal can be any one of the undesired signals x 1 (k) (t) in the input to the correlator in Figure 1.
Therefore, the undesired signal without Doppler or time delay is modelled as: where p Tc (t) is a rectangular pulse of width T c : and c v ∈ {−1, 1} is spreading code chips with a chipping rate 1/T c = 2.046 MHz. The chip width is T c ; d k is navigation data bits with bit rate 1/T b = 50 bps, and the bit width is T b . α u is the NH code adopted for B1-I signal and takes a fixed pattern of (0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0), which spans one navigation data bit. The spreading code here is treated as a random periodic sequence of length 2046 chips, with chips independent from each other. The navigation data bit stream is also treated as a random periodic sequence, but with an infinite length. In order to facilitate the development in the appendices, spreading code period, T, is introduced, and there is an obvious relationship among T c , T, and T b :T b = 20T = 40,920T c . This relationship is illustrated in Figure 2. The receiver local code replica or desired signal is modelled similarly as: using the same notation as Equation (1). The only difference between the desired and undesired code is that the undesired one has a time delay of ∆ seconds, which is normally within the range of ±20 ms for BDS users on or near the horizon of the Earth. Both signals are cyclostationary, not WSS, and therefore, their characteristics are accurately represented using a two-parameter ACF [16]. Since the period of cyclostationarity is 20 ms, the model derived in the following sections only takes into account time delays ranging from 0 to 20 ms, i.e., 0 ≤ ∆ ≤ 20 ms, which is equal to the time delay modulo 20 ms. In order to facilitate the development in the appendices, spreading code period, T, is introduced, and there is an obvious relationship among Tc, T, and Tb:Tb = 20T = 40,920Tc. This relationship is illustrated in Figure 2. The receiver local code replica or desired signal is modelled similarly as: using the same notation as Equation (1). The only difference between the desired and undesired code is that the undesired one has a time delay of ∆ seconds, which is normally within the range of ±20 ms for BDS users on or near the horizon of the Earth. Both signals are cyclostationary, not WSS, and therefore, their characteristics are accurately represented using a two-parameter ACF [16]. Since the period of cyclostationarity is 20 ms, the model derived in the following sections only takes into account time delays ranging from 0 to 20 ms, i.e., 0 ≤ ∆ ≤ 20 ms, which is equal to the time delay modulo 20 ms. In a real receiver, both the incoming signal and local replica have a Doppler shift. However, for simplicity and without loss of generality, we assume here that only the incoming signal has its Doppler shift, f, while the local replica has zero Doppler. This is simply because the receiver carrier tracking loop is after all, tracking the differential Doppler between the incoming and the local replica, not the absolute Doppler of each signal.

SSC and Equivalent White Noise Level
The extent of the undesired signal's impact on the desired signal is characterized by an equivalent white noise level, which is derived from SSC, whose analytical forms are derived in this section. The k-th correlator output, yk, is a random variable. Its mean value is: (4) and its variance is: ( 1) ( 1) 2 2 1 1 2 2 ( 1) ( 1) 2 2 11 22 , ,  (5) where the dot within the curly brace means "dot product". The development of Equation (5) uses the definition of auto-correlation function of two complex signals. The two-parameter ACF's of both the undesired signal and receiver local replica are derived in Appendix A. The result is as follows: In a real receiver, both the incoming signal and local replica have a Doppler shift. However, for simplicity and without loss of generality, we assume here that only the incoming signal has its Doppler shift, f, while the local replica has zero Doppler. This is simply because the receiver carrier tracking loop is after all, tracking the differential Doppler between the incoming and the local replica, not the absolute Doppler of each signal.

SSC and Equivalent White Noise Level
The extent of the undesired signal's impact on the desired signal is characterized by an equivalent white noise level, which is derived from SSC, whose analytical forms are derived in this section. The k-th correlator output, y k , is a random variable. Its mean value is: and its variance is: where the dot within the curly brace means "dot product". The development of Equation (5) uses the definition of auto-correlation function of two complex signals. The two-parameter ACF's of both the undesired signal and receiver local replica are derived in Appendix A. The result is as follows: where R 22 (s, t) is the two-parameter ACF of the local replica, the integer v is selected as: and: where R 11 (s, t) is the two-parameter ACF of the undesired signal, ∆ is the time delay of the undesired signal and the integer k and v are selected as: where function mod(a, m) returns the remainder after division of a by m where a is the dividend and m is the divisor. Function floor(x) rounds x to the nearest integer less than or equal to x. The development of Equations (8)-(10) has taken into account the special structure of the B1-I signal which has an NH code modulated on it: in subsection "Desired signal ACF" of Appendix A it can be seen that since the NH code is a fixed pattern finite-length code, it is treated as a deterministic variable here. Figure 3 illustrates an example of ACF for the desired signal R 22 (s, t) with s ∈ [0, T c ). For every fixed value of the first parameter, s, this two parameter ACF is composed of 20 pulses of width T c with unity amplitude, over arbitrary 20-ms correlation interval, which is equal to one data bit width. For example, if s falls within the time duration of the repetition of the first chip of the spreading code sequence, the ACF takes on a value of unity during each repetition of this chip and zero elsewhere. The same ACF results if s fell within repetition of the first chip, e.g., s ∈ [7T, 7T+T c ). The ACF for the undesired signal can be obtained by shifting the desired signal's ACF by ∆.
where R22(s, t) is the two-parameter ACF of the local replica, the integer v is selected as: and: where R11(s, t) is the two-parameter ACF of the undesired signal, ∆ is the time delay of the undesired signal and the integer k and v are selected as: where function mod(a, m) returns the remainder after division of a by m where a is the dividend and m is the divisor. Function floor(x) rounds x to the nearest integer less than or equal to x. The development of Equations (8)-(10) has taken into account the special structure of the B1-I signal which has an NH code modulated on it: in subsection "Desired signal ACF" of Appendix A it can be seen that since the NH code is a fixed pattern finite-length code, it is treated as a deterministic variable here.   signal are aligned, (ii) data bits of these two signals are not aligned by an integer multiples of a code chip, and (iii) data bits of these two signals are misaligned by an arbitrary delay. This variance is then converted to the equivalent white noise level using SSC, which is the metric to measure how the self-interference affects receiver performance [17].

Data Bits Aligned
When the data bits modulated on the undesired and desired signal are aligned, i.e., ∆ = 0, the situation is illustrated in Figure 4. converted to the equivalent white noise level using SSC, which is the metric to measure how the selfinterference affects receiver performance [17].

Data Bits Aligned
When the data bits modulated on the undesired and desired signal are aligned, i.e., ∆ = 0, the situation is illustrated in Figure 4. In this case, ACF's of the undesired and desired signals become identical over any arbitrary integration interval. The variance of the correlator output could be represented as: R s t e e dsdt (11) the last line of which comes from the fact that the ACF's of the incoming and local code are identical: both of them are unity pulse trains. As derived in Appendix B, substituting Equation (6) into Equation (11), we have: (12) where: is the sampling function. On the other hand, assume that the correlator shown in Figure 1 was driven by additive Gaussian white noise with power spectral density I0, then its output should have zero mean and a variance of: which is the product of the power spectral density and equivalent bandwidth 1/Tb. Comparing Equations (12) and (14), and an equivalent white noise with power spectral density I0 is obtained defined as: In this case, ACF's of the undesired and desired signals become identical over any arbitrary integration interval. The variance of the correlator output could be represented as: the last line of which comes from the fact that the ACF's of the incoming and local code are identical: both of them are unity pulse trains. As derived in Appendix B, substituting Equation (6) into Equation (11), we have: where: is the sampling function. On the other hand, assume that the correlator shown in Figure 1 was driven by additive Gaussian white noise with power spectral density I 0 , then its output should have zero mean and a variance of: which is the product of the power spectral density and equivalent bandwidth 1/T b . Comparing Equations (12) and (14), and an equivalent white noise with power spectral density I 0 is obtained defined as: Sensors 2017, 17, 663 7 of 22 because this equivalent white noise generates the same variance values as a true interference does at the output of the correlator. Based on this definition, we can further define the spectral separation coefficient [17] for B1-I as SSC B1-I : where, P R is inserted for the more general case where an undesired signal has an arbitrary power level in watts.

Data Bits Misaligned by an Integer Number of Code Chips
When the data bits modulated on the undesired and desired signals are aligned, i.e., ∆ is not zero, the situation is illustrated in Figure 5. because this equivalent white noise generates the same variance values as a true interference does at the output of the correlator. Based on this definition, we can further define the spectral separation coefficient [17] for B1-I as SSCB1-I: where, PR is inserted for the more general case where an undesired signal has an arbitrary power level in watts.

Data Bits Misaligned by an Integer Number of Code Chips
When the data bits modulated on the undesired and desired signals are aligned, i.e., ∆ is not zero, the situation is illustrated in Figure 5. In this case, the ACF of the incoming/undesired signal is: where k and v is set up using Equations (7), (9) and (10). As derived in Appendix B, assuming that the data bit misalignment is an integer multiple of the spreading code period T, i.e., ∆ = KT for 0 ≤ K < 20, where K is an integer, results in: which, for users near or at the Earth's horizon, could be approximated as: More generally, assuming that the data bit misalignment is an integer multiple of the spreading code chip width Tc, i.e., ∆ = KT + CTc, where 0 ≤ K < 20, 0 < C < 2046, K and C is an integer, yields:  (20) where SSCB1-I,K is the SSC in Equation (18) for any users or SSC in Equation (19) for users near or at the Earth's horizon, when ∆ = KT. In this case, the ACF of the incoming/undesired signal is: where k and v is set up using Equations (7), (9) and (10). As derived in Appendix B, assuming that the data bit misalignment is an integer multiple of the spreading code period T, i.e., ∆ = KT for 0 ≤ K < 20, where K is an integer, results in: which, for users near or at the Earth's horizon, could be approximated as: More generally, assuming that the data bit misalignment is an integer multiple of the spreading code chip width T c , i.e., ∆ = KT + CT c , where 0 ≤ K < 20, 0 < C < 2046, K and C is an integer, yields: where SSC B1-I,K is the SSC in Equation (18)

Data Bits Misaligned by a Fraction of a Code Chip
Section 3.2 "Data Bits Misaligned by an Integer Number of Code Chips" summarizes SSC expressions that is valid when the differential time delay between undesired and desired signals, ∆ is multiples of a single code chip T c . With this constraint, R 11 (s − ∆, t − ∆)R 22 (s,t) is a pulse train of width T c , given by Equation (17). Now we take a step further to examine what the SSC should look like when ∆ is only a fraction of a code chip. In this case, R 11 (s − ∆, t − ∆)R 22 (s,t) is still a pulse train of uniformed width, but the pulse width is no longer equal to T c . As illustrated in Figures 6 and 7, for fixed s, R 11 Sensors 2017, 17, 663 8 of 23

Data Bits Misaligned by a Fraction of a Code Chip
Section 3.2 "Data Bits Misaligned by an Integer Number of Code Chips" summarizes SSC expressions that is valid when the differential time delay between undesired and desired signals, ∆ is multiples of a single code chip Tc. With this constraint, R11(s − ∆, t − ∆)R22(s,t) is a pulse train of width Tc, given by Equation (17). Now we take a step further to examine what the SSC should look like when ∆ is only a fraction of a code chip. In this case, R11(s − ∆, t − ∆)R22(s,t) is still a pulse train of uniformed width, but the pulse width is no longer equal to Tc. As illustrated in Figures 6 and 7 is either a train of pulses of width ∆ or Tc − ∆.  When the differential time delay takes on arbitrary value, we proceed by using τ = mod(∆, Tc), 0 ≤ τ < Tc. Therefore, R11(s − ∆, t − ∆)R22(s,t) will be a train of pulses of width τ for τ/Tc percent of Tc, or a train of pulses of width Tc -τ for (Tc -τ)/Tc percent of Tc. The total effect of such arbitrary time delay, for the differential Doppler, is that a factor of Tc in the SSCB1-I expressions in the Section 3.2 "Data Bits Misaligned by an Integer Number of Code Chips" is replaced by [12]: which implies the expectation of receiving such a time difference. It can be seen from Equation (21) that spreading code self-interference is maximized when τ = 0, i.e., when the misalignment between the desired and undesired signals happens to be multiples of a code chip; it can be seen also that self-interference is minimized then τ = Tc/2. Moreover, the average of Equation (21) can be computed as:

Data Bits Misaligned by a Fraction of a Code Chip
Section 3.2 "Data Bits Misaligned by an Integer Number of Code Chips" summarizes SSC expressions that is valid when the differential time delay between undesired and desired signals, ∆ is multiples of a single code chip Tc. With this constraint, R11(s − ∆, t − ∆)R22(s,t) is a pulse train of width Tc, given by Equation (17). Now we take a step further to examine what the SSC should look like when ∆ is only a fraction of a code chip. In this case, R11(s − ∆, t − ∆)R22(s,t) is still a pulse train of uniformed width, but the pulse width is no longer equal to Tc. As illustrated in Figures 6 and 7 is either a train of pulses of width ∆ or Tc − ∆.  When the differential time delay takes on arbitrary value, we proceed by using τ = mod(∆, Tc), 0 ≤ τ < Tc. Therefore, R11(s − ∆, t − ∆)R22(s,t) will be a train of pulses of width τ for τ/Tc percent of Tc, or a train of pulses of width Tc -τ for (Tc -τ)/Tc percent of Tc. The total effect of such arbitrary time delay, for the differential Doppler, is that a factor of Tc in the SSCB1-I expressions in the Section 3.2 "Data Bits Misaligned by an Integer Number of Code Chips" is replaced by [12]: which implies the expectation of receiving such a time difference. It can be seen from Equation (21) that spreading code self-interference is maximized when τ = 0, i.e., when the misalignment between the desired and undesired signals happens to be multiples of a code chip; it can be seen also that self-interference is minimized then τ = Tc/2. Moreover, the average of Equation (21) can be computed as: When the differential time delay takes on arbitrary value, we proceed by using τ = mod(∆, T c ), 0 ≤ τ < T c . Therefore, R 11 (s − ∆, t − ∆)R 22 (s,t) will be a train of pulses of width τ for τ/T c percent of T c , or a train of pulses of width T cτ for (T cτ)/T c percent of T c . The total effect of such arbitrary time delay, for the differential Doppler, is that a factor of T c in the SSC B1-I expressions in the Section 3.2 "Data Bits Misaligned by an Integer Number of Code Chips" is replaced by [12]: which implies the expectation of receiving such a time difference. It can be seen from Equation (21) that spreading code self-interference is maximized when τ = 0, i.e., when the misalignment between the desired and undesired signals happens to be multiples of a code chip; it can be seen also that self-interference is minimized then τ = T c /2. Moreover, the average of Equation (21) can be computed as: which, when substituted for T c in Equation (19) yields: which is used to compute SSC in the following sections.

Proposed Model for BDS MEO/IGSO B1-I Self-Interference
Based on Equations (20) and (23), the proposed model for BDS MEO/IGSO B1-I self-interference, or the aggregate equivalent white noise due to self-interference can be expressed as: where I 0,m is the aggregate equivalent white noise imposed on the m-th signal by self-interference, p m is the received power at the antenna, (SSC B1-I ) i m is the SSC between the i-th and m-th signals, which is determined by differential Doppler f i m and differential time delay K i m T + C i m T c . Table 1 shows an example of how to use the developed analytical model, i.e., Equations (20) and (23) to predict SSC B1-I and equivalent white noise I 0 . Note that the last line of Table 1 comes from the sum of I 0 of PRN 8 and 9. At a specific time, there are 8 satellites visible to the receiver. In Table 1, the first column is the PRN number. The second column is the relative received power, relative to −160 dBW. The third column is the true (geometric) range from satellite to the receiver. The fourth column is the true Doppler. These four columns are the parameters used to compute the fifth and sixth columns, i.e., SSC B1-I and equivalent white noise I 0 , respectively. Since here it is assumed that PRN 6 is the desired signal, it is unnecessary to compute its own SSC and I 0 . Table 2 presents the procedure to compute SSC and I 0 . The procedure described in Table 2 is simply a re-iteration of the development of the model Equation (24). 1: while GPS measurement do 2: Use true range r to obtain transit time t by t = r/c, where c is the speed of light in vacuum. 3: Use relative received power p_rel to obtain actual received power p by p = p_rel + p_ref, where p_ref is the reference power set in the software simulator. 4: Obtain differential Doppler f by differencing the first satellite with the rest. 5: Obtain differential transit time ∆ by differencing the first satellite with the rest. 6: Obtain the number of 1ms in ∆, K by K = floor (absolute value of ∆ × 10 3 ). 7: Obtain the number of code chips in ∆, C by C = round( (∆ − K) × chip_rate × 10 −3 ). 8: if C is equal to code period (2046) then 9: reset C. 10: increment K. 11: end if 12: Use Equations (20) and (23) to obtain SSC B1-I . 13: I 0 = p + SSC B1-I . 14: end while

Model Validation
Live BDS signal records from a hardware signal simulator are collected to test if the proposed model is correct within a certain boundary, set by established benchmark work.
The validity of the model is not directly proved because the proposed model includes the NH code effect and there is no current work covering this unique signal component on B1 self-interference evaluation, though there does exist effective ways to overcome NH code impact on various stages of signal processing [15]. Since no other formula is currently available, an indirect method is applied to justify the assumptions and the proposed model. Self-interference effect on receiver noise is intangible or indirect to see or to feel but can appropriately be shown and visualized by its impact in the measurement domain (i.e., DLL output or a step further, code range), and in this regard, we choose to apply Betz's well-established expressions linking C/N 0 with discriminator output statistics. Since no validation of Betz's model was provided in Betz's paper [18] or his later publications, and therefore in order to use Betz's formula, we had first to prove it (or at least, show its effectiveness using real measurements), which is presented in this section and does form an integral part of the completeness and soundness of the validation logic.

Criterion
Model validity is checked by comparing the predicted and measured tracking standard deviation, using the expression [18]: where K is code tracking loop gain and C s the measured signal power. var{e(ε) | τ s k } is measured discriminator output variance. var{τ u k | τ s k } is variance of the unsmoothed code delay estimate τ u k , which is predicted by using Equation (26) [19]: which is the analytical form of code tracking variance. The term '(C s /N 0 ) eff ' is effective carrier to noise ratio, and is equal to C s − I 0 , where C s is the measured signal power and I 0 the equivalent white noise as defined in Equation (16) and transformed into dB. Here I 0 includes self-interference effects, to be evaluated using the proposed model (Equations (16) and (20)). Both C s and I 0 are in dB and (C s /N 0 ) eff should be converted back to non-logarithmic unit when it is used in Equation (26). If the model is correct, the squared root of the two sides of Equation (25) should agree within a small error tolerance for a specific channel in a software receiver: where λ is a tolerance value and we set it as follows. We compare our prediction accuracy with accuracy of most recent analytical methods for self-interference, such as the work by Cerruti et al. [4], where the accuracy of prediction, i.e., the difference between predicted and measured carrier tracking standard deviation is 0.5 degrees. This translates into a prediction accuracy of 0.0264 cm (0.5/360 × 19) in range by phase measurements, which is in turn, translated into code tracking prediction accuracy of 6.6 cm by a scale factor of 250, since code tracking loop is 100 (with multipath) to 250 (without multipath; this is the case with current validation settings described in Table 3, and is also the case with the benchmark work by Cerruti et al. [3]) times noisier than carrier tracking loop in terms of tracking standard deviation [20]. Therefore, we set λ = 0.066 m.

Test Setup
The test setup is illustrated in Figure 8. Since we want to compare predicted and measured tracking standard deviations, we have to collect live signal records to measure the code tracking standard deviation using a software receiver. Without loss of generality, we choose two IGSO satellites, PRN 13 and 14 instead of the whole current constellation to test the developed model. This is because: (i) In validation, only two satellites (desired and undesired signal) are necessary; (ii) PRN 13 & 14 are the two of the four currently operational IGSO satellites in BDS, and since MEO and IGSO in BDS use the same signal modulation scheme; IGSO signal will suffice to represent a combination of MEO and IGSO signals. Therefore, first we collect IF data containing both PRN 13 and PRN 14 signals and then we collect IF data containing only the PRN 13 signal. Finally, we obtain the measured code tracking standard deviation of PRN 13 using the first dataset, which includes the interference effect imposed by PRN 14. This interference effect is predicted using the estimated power (C s ) from the second dataset (free of interfering signal PRN 14) and our proposed model. As a result, a hardware signal simulator and not live signals from satellites, must be used as the signal source. as defined in Equation (16) and transformed into dB. Here I0 includes self-interference effects, to be evaluated using the proposed model (Equations (16) and (20)). Both Cs and I0 are in dB and (Cs/N0)eff should be converted back to non-logarithmic unit when it is used in Equation (26). If the model is correct, the squared root of the two sides of Equation (25) should agree within a small error tolerance for a specific channel in a software receiver: where λ is a tolerance value and we set it as follows. We compare our prediction accuracy with accuracy of most recent analytical methods for self-interference, such as the work by Cerruti et al. [4], where the accuracy of prediction, i.e., the difference between predicted and measured carrier tracking standard deviation is 0.5 degrees. This translates into a prediction accuracy of 0.0264 cm (0.5/360 × 19) in range by phase measurements, which is in turn, translated into code tracking prediction accuracy of 6.6 cm by a scale factor of 250, since code tracking loop is 100 (with multipath) to 250 (without multipath; this is the case with current validation settings described in Table 3, and is also the case with the benchmark work by Cerruti et al. [3]) times noisier than carrier tracking loop in terms of tracking standard deviation [20]. Therefore, we set λ = 0.066 m.

Test Setup
The test setup is illustrated in Figure 8. Since we want to compare predicted and measured tracking standard deviations, we have to collect live signal records to measure the code tracking standard deviation using a software receiver. Without loss of generality, we choose two IGSO satellites, PRN 13 and 14 instead of the whole current constellation to test the developed model. This is because: (i) In validation, only two satellites (desired and undesired signal) are necessary; (ii) PRN 13 & 14 are the two of the four currently operational IGSO satellites in BDS, and since MEO and IGSO in BDS use the same signal modulation scheme; IGSO signal will suffice to represent a combination of MEO and IGSO signals. Therefore, first we collect IF data containing both PRN 13 and PRN 14 signals and then we collect IF data containing only the PRN 13 signal. Finally, we obtain the measured code tracking standard deviation of PRN 13 using the first dataset, which includes the interference effect imposed by PRN 14. This interference effect is predicted using the estimated power (Cs) from the second dataset (free of interfering signal PRN 14) and our proposed model. As a result, a hardware signal simulator and not live signals from satellites, must be used as the signal source. We use a HWA-RNSS-7400, which is a multi-constellation GNSS signal simulator by HWA Create Co., Ltd (Beijing, China) [21]. It features BDS B1, B2, and B3 signal simulation and is widely used within BDS industries. We also use a SAS6862A, which is a wide-band, dual channel (GPS L1+BDS B1 and BDS B3) radio frequency front end (RFFE) with a highly stable (better than 0.005 ppm) oven controlled crystal oscillator (OCXO) [22]. We use a multi-constellation software receiver called UNIversal COmmunications and Radio Navigation receiver (UNICORN), developed by the authors at the Centre for Transport Studies (CTS), Imperial College London and used to obtain the ESA's We use a HWA-RNSS-7400, which is a multi-constellation GNSS signal simulator by HWA Create Co., Ltd. (Beijing, China) [21]. It features BDS B1, B2, and B3 signal simulation and is widely used within BDS industries. We also use a SAS6862A, which is a wide-band, dual channel (GPS L1+BDS B1 and BDS B3) radio frequency front end (RFFE) with a highly stable (better than 0.005 ppm) oven controlled crystal oscillator (OCXO) [22]. We use a multi-constellation software receiver called UNIversal COmmunications and Radio Navigation receiver (UNICORN), developed by the authors at the Centre for Transport Studies (CTS), Imperial College London and used to obtain the ESA's certificate for the "First 50 Successful Galileo Position Fixes". The receiver has been adapted to process the B1 signal.
Parameters used to configure UNICORN and drive the analytical model are the same and are provided in Table 3. The simulation starts at UTC 00:10:00, 27 Feb 27, 2014 and two 37-s IF data was recorded using SAS6862A. The receiver is assumed to be located at 31 • N, 121.5 • E, at an altitude of 5 m in Shanghai, China. Table 3. Parameters used in model validation.

Sections Parameters Values
Satellite UNICORN has to produce two measurements for the analytical model: C/N 0 and code loop gain, which is indispensable for successful evaluation of tracking error variance through the model. Within UNICORN, several lines of code were inserted to calculate tracking error standard deviation.

Results
Finally, the results from UNICORN and the analytical model are compared to see if the predicted values agree with measured values according to Equation (27). Figure 9 shows the results. certificate for the "First 50 Successful Galileo Position Fixes". The receiver has been adapted to process the B1 signal. Parameters used to configure UNICORN and drive the analytical model are the same and are provided in Table 3. The simulation starts at UTC 00:10:00, 27 Feb 27, 2014 and two 37-s IF data was recorded using SAS6862A. The receiver is assumed to be located at 31° N, 121.5° E, at an altitude of 5 m in Shanghai, China. UNICORN has to produce two measurements for the analytical model: C/N0 and code loop gain, which is indispensable for successful evaluation of tracking error variance through the model. Within UNICORN, several lines of code were inserted to calculate tracking error standard deviation.

Results
Finally, the results from UNICORN and the analytical model are compared to see if the predicted values agree with measured values according to Equation (27). Figure 9 shows the results.  In Figure 9, it can be observed that the discrepancy between measured and predicted standard deviations are bounded by 0.0002 B1 primary code chips, which translates into a prediction accuracy of 0.0293 m, less than half of the threshold value of 0.66 m. This proves the validity of the analytical model.

Results
Based on the validity of the model, we now apply this model to show some insights into the typical behaviour of self-interference in BDS B1-I signals. Firstly, we try to find when the equivalent noise introduced by the self-interference approaches its maximum value.
The SSC predicted from Equation (23) is illustrated in Figures 10 and 11. Here the differential time delay between the desired and undesired signals is an integer multiple of spreading code periods. Since the SSC values for K > 10 is equal to SSC values with 20 − K, only cases with K ≤ 10 are shown here. To facilitate viewing, the results are shown separately when K is odd or even. From Figures 10 and 11 we can conclude that when the differential Doppler between desired and undesired signals is bounded by 50 Hz, self-interference will introduce the most serious noise into receiver processing chain. In Figure 9, it can be observed that the discrepancy between measured and predicted standard deviations are bounded by 0.0002 B1 primary code chips, which translates into a prediction accuracy of 0.0293 m, less than half of the threshold value of 0.66 m. This proves the validity of the analytical model.

Results
Based on the validity of the model, we now apply this model to show some insights into the typical behaviour of self-interference in BDS B1-I signals. Firstly, we try to find when the equivalent noise introduced by the self-interference approaches its maximum value.
The SSC predicted from Equation (23) is illustrated in Figures 10 and 11. Here the differential time delay between the desired and undesired signals is an integer multiple of spreading code periods. Since the SSC values for K > 10 is equal to SSC values with 20 − K, only cases with K ≤ 10 are shown here. To facilitate viewing, the results are shown separately when K is odd or even. From Figures 10 and 11 we can conclude that when the differential Doppler between desired and undesired signals is bounded by 50 Hz, self-interference will introduce the most serious noise into receiver processing chain. Figure 10. SSCB1-I variation as a function of differential Doppler (between desired and undesired signals) and differential time delay. K is even. Figure 11. SSCB1-I variation as a function of differential Doppler (between desired and undesired signals) and differential time delay. K is odd. Figure 10. SSC B1-I variation as a function of differential Doppler (between desired and undesired signals) and differential time delay. K is even. In Figure 9, it can be observed that the discrepancy between measured and predicted standard deviations are bounded by 0.0002 B1 primary code chips, which translates into a prediction accuracy of 0.0293 m, less than half of the threshold value of 0.66 m. This proves the validity of the analytical model.

Results
Based on the validity of the model, we now apply this model to show some insights into the typical behaviour of self-interference in BDS B1-I signals. Firstly, we try to find when the equivalent noise introduced by the self-interference approaches its maximum value.
The SSC predicted from Equation (23) is illustrated in Figures 10 and 11. Here the differential time delay between the desired and undesired signals is an integer multiple of spreading code periods. Since the SSC values for K > 10 is equal to SSC values with 20 − K, only cases with K ≤ 10 are shown here. To facilitate viewing, the results are shown separately when K is odd or even. From Figures 10 and 11 we can conclude that when the differential Doppler between desired and undesired signals is bounded by 50 Hz, self-interference will introduce the most serious noise into receiver processing chain. Figure 10. SSCB1-I variation as a function of differential Doppler (between desired and undesired signals) and differential time delay. K is even. Figure 11. SSCB1-I variation as a function of differential Doppler (between desired and undesired signals) and differential time delay. K is odd. Figure 11. SSC B1-I variation as a function of differential Doppler (between desired and undesired signals) and differential time delay. K is odd.
Secondly, in order to analyze the self-interference effects in current BDS configuration, a simulation case is setup. The orbit propagation, receiver acquisition and tracking are performed by the Spirent SimGEN simulator software suite. SimGEN is capable of full and versatile case generation. Users are allowed to control multiple GNSS and regional satellite constellations including BDS, signal propagation, multipath and obstruction effects, antenna patterns, vehicle trajectory and various error models, thus vastly facilitating our simulation. For the current simulation, three inputs are needed from SimGEN, namely true range, received power, and carrier Doppler shift. However, due to the version of SimGEN, it is not able to simulate BDS constellation. The problem is bypassed by using Two Line Element (TLE) file retrieved from NORAD [23] and feeding it to the Motion control panel for GPS in SimGEN. The simulation length is 20 minutes, starting from UTCG1927, 26 April 2015, with a measurement rate of 1 Pulse Per Second (PPS). The receiver is placed in Shanghai, China. The result is shown in Figure 12. The constellation information is shown in Appendix C. Secondly, in order to analyze the self-interference effects in current BDS configuration, a simulation case is setup. The orbit propagation, receiver acquisition and tracking are performed by the Spirent SimGEN simulator software suite. SimGEN is capable of full and versatile case generation. Users are allowed to control multiple GNSS and regional satellite constellations including BDS, signal propagation, multipath and obstruction effects, antenna patterns, vehicle trajectory and various error models, thus vastly facilitating our simulation. For the current simulation, three inputs are needed from SimGEN, namely true range, received power, and carrier Doppler shift. However, due to the version of SimGEN, it is not able to simulate BDS constellation. The problem is bypassed by using Two Line Element (TLE) file retrieved from NORAD [23] and feeding it to the Motion control panel for GPS in SimGEN. The simulation length is 20 minutes, starting from UTCG1927, 26 April 2015, with a measurement rate of 1 Pulse Per Second (PPS). The receiver is placed in Shanghai, China. The result is shown in Figure 12. The constellation information is shown in Appendix C. As illustrated in Figure 12, we can obtain such observations: (1) Within most of the simulation interval, the individual equivalent noise and aggregate noise is below −215 dBW/Hz, which is much lower than a typical receiver noise floor of about −201.5 dBW/Hz, and therefore would not cause much trouble for most applications. Here trouble means worsening of such receiver capabilities as acquisition sensitivity and tracking jitter. (2) During 200 and 240 s, the equivalent noise of PRN 8 and the aggregate noise reaches −200 dBW/Hz. For a typical value of receiver noise floor (−201.5 dBW/Hz), this noise floor will be transferred to −197.67 dBW/Hz. In this case, PRN 6's effective carrier to noise density ratio (Cs/N0)eff is 3.83 dB lower than the original carrier to noise density ratio, should the selfinterference caused by PRN 8 be absent. This could pose a serious problem to high-sensitivity processing during acquisition, tracking or navigation bit modulation. The potential hazard caused by this (C/N0)eff decrease could be best exemplified by Figure 13, which shows the code tracking standard deviation versus (Cs/N0)eff , using established analytical method by Zhang and Zhan [19]. The B1-I receiver used in Figure 13 is assumed to have a precorrelation bandwidth of 4 MHz, a sampling frequency of 16 MHz, and 2-bit quantization. The coherent integration interval is 20 ms and a Dot Product discriminator is used. When (C/N0)eff decreases from 19 to 15 dB-Hz, tracking standard deviation increases from 11.07 m to 22.72 m. In this regard, B1-I selfinterference could pose a potential hazard on pseudorange measurement for high-sensitivity receivers. However, it must be underlined that for most receivers (from mass market to As illustrated in Figure 12, we can obtain such observations: (1) Within most of the simulation interval, the individual equivalent noise and aggregate noise is below −215 dBW/Hz, which is much lower than a typical receiver noise floor of about −201.5 dBW/Hz, and therefore would not cause much trouble for most applications. Here trouble means worsening of such receiver capabilities as acquisition sensitivity and tracking jitter. (2) During 200 and 240 s, the equivalent noise of PRN 8 and the aggregate noise reaches −200 dBW/Hz. For a typical value of receiver noise floor (−201.5 dBW/Hz), this noise floor will be transferred to −197.67 dBW/Hz. In this case, PRN 6's effective carrier to noise density ratio (C s /N 0 ) eff is 3.83 dB lower than the original carrier to noise density ratio, should the self-interference caused by PRN 8 be absent. This could pose a serious problem to high-sensitivity processing during acquisition, tracking or navigation bit modulation. The potential hazard caused by this (C/N 0 ) eff decrease could be best exemplified by Figure 13, which shows the code tracking standard deviation versus (C s /N 0 ) eff , using established analytical method by Zhang and Zhan [19]. The B1-I receiver used in Figure 13 is assumed to have a precorrelation bandwidth of 4 MHz, a sampling frequency of 16 MHz, and 2-bit quantization. The coherent integration interval is 20 ms and a Dot Product discriminator is used. When (C/N 0 ) eff decreases from 19 to 15 dB-Hz, tracking standard deviation increases from 11.07 m to 22.72 m. In this regard, B1-I self-interference could pose a potential hazard on pseudorange measurement for high-sensitivity receivers. However, it must be underlined that for most receivers (from mass market to medium-grade receiver), a (C/N 0 ) eff value under 20 dB-Hz, though may be possible, is an exceptional case. When (C/N 0 ) eff decreases from 24 to 20 dB-Hz, tracking standard deviation increases from 5 m to 9 m, which does not impose significant impact on final positioning accuracy for most low-to medium-grade receivers.
(3) Around 240 s for PRN 8 and around 1200 s for PRN 9, the self-interference effect on the desired signal PRN 6 reaches its maximal value for each undesired signal respectively. This corresponds to the time when these two PRN's Doppler shift (modulo 1 kHz) are aligned (zero differential Doppler shift) with the desired signal, PRN 6. The same situation has also been found in GPS C/A-to-C/A self-interference [4].
Sensors 2017, 17, 663 15 of 23 medium-grade receiver), a (C/N0)eff value under 20 dB-Hz, though may be possible, is an exceptional case. When (C/N0)eff decreases from 24 to 20 dB-Hz, tracking standard deviation increases from 5 m to 9 m, which does not impose significant impact on final positioning accuracy for most low-to medium-grade receivers. (3) Around 240 s for PRN 8 and around 1200 s for PRN 9, the self-interference effect on the desired signal PRN 6 reaches its maximal value for each undesired signal respectively. This corresponds to the time when these two PRN's Doppler shift (modulo 1 kHz) are aligned (zero differential Doppler shift) with the desired signal, PRN 6. The same situation has also been found in GPS C/A-to-C/A self-interference [4]. In this regard, the proposed model could be used in B1 receiver as a code-tracking performance indicator, since code tracking standard deviation is directly determined by (Cs/N0)eff., the effective carrier-to-noise power density ratio, which is in turn, directly affected by self-interference effects, as indicated in Figure 13. In the simplest case, the common satellite selection step (e.g., depending on DOP, Dilution of Precision) should be complemented by a new step to narrow down the list of satellites participating in navigation solution, based on the result of self-interference calculation according to Equations (20) and (23), as suggested in Figure 14.
In Figure 14, a typical receiver baseband (correlator + measurements) plus navigation module is shown, aided by a satellite selection module, whose action depends on a Signal Quality Monitoring (SQM) module [24] on top of its normal operations, monitoring the equivalent noise spectral density imposed by self-interference,. This may be an indication to receiver manufacturers for their design of strategies for satellite selection in Position, Velocity and Time (PVT) computation and signal acquisition in case of very low C/N0, which could be estimated stably and accurately through novel estimators [25].  In this regard, the proposed model could be used in B1 receiver as a code-tracking performance indicator, since code tracking standard deviation is directly determined by (C s /N 0 ) eff ., the effective carrier-to-noise power density ratio, which is in turn, directly affected by self-interference effects, as indicated in Figure 13. In the simplest case, the common satellite selection step (e.g., depending on DOP, Dilution of Precision) should be complemented by a new step to narrow down the list of satellites participating in navigation solution, based on the result of self-interference calculation according to Equations (20) and (23), as suggested in Figure 14.
In Figure 14, a typical receiver baseband (correlator + measurements) plus navigation module is shown, aided by a satellite selection module, whose action depends on a Signal Quality Monitoring (SQM) module [24] on top of its normal operations, monitoring the equivalent noise spectral density imposed by self-interference,. This may be an indication to receiver manufacturers for their design of strategies for satellite selection in Position, Velocity and Time (PVT) computation and signal acquisition in case of very low C/N 0 , which could be estimated stably and accurately through novel estimators [25]. medium-grade receiver), a (C/N0)eff value under 20 dB-Hz, though may be possible, is an exceptional case. When (C/N0)eff decreases from 24 to 20 dB-Hz, tracking standard deviation increases from 5 m to 9 m, which does not impose significant impact on final positioning accuracy for most low-to medium-grade receivers. (3) Around 240 s for PRN 8 and around 1200 s for PRN 9, the self-interference effect on the desired signal PRN 6 reaches its maximal value for each undesired signal respectively. This corresponds to the time when these two PRN's Doppler shift (modulo 1 kHz) are aligned (zero differential Doppler shift) with the desired signal, PRN 6. The same situation has also been found in GPS C/A-to-C/A self-interference [4]. In this regard, the proposed model could be used in B1 receiver as a code-tracking performance indicator, since code tracking standard deviation is directly determined by (Cs/N0)eff., the effective carrier-to-noise power density ratio, which is in turn, directly affected by self-interference effects, as indicated in Figure 13. In the simplest case, the common satellite selection step (e.g., depending on DOP, Dilution of Precision) should be complemented by a new step to narrow down the list of satellites participating in navigation solution, based on the result of self-interference calculation according to Equations (20) and (23), as suggested in Figure 14.
In Figure 14, a typical receiver baseband (correlator + measurements) plus navigation module is shown, aided by a satellite selection module, whose action depends on a Signal Quality Monitoring (SQM) module [24] on top of its normal operations, monitoring the equivalent noise spectral density imposed by self-interference,. This may be an indication to receiver manufacturers for their design of strategies for satellite selection in Position, Velocity and Time (PVT) computation and signal acquisition in case of very low C/N0, which could be estimated stably and accurately through novel estimators [25].

Conclusions
A simple analytical model is developed to predict the impact of B1-I code self-interference on BDS receiver functions that are dependent primarily on the sum of the correlations (e.g., carrier phase tracking and data demodulation). The model is based on the fact that the code is actually a cyclostationary process rather than WSS, as previously assumed in SSC evaluation. A two-parameter ACF, in contrast to the usual one-paramter ACF, is applied to characterize B1-I code either with or without navigation data. The assumption of cyclostationarity is essential in evaluating self-interference effects for short codes, producing excellent agreement with actual receiver observables. The model is found to be able to accurately predict self-interference effects on receiver code tracking standard deviation, when the difference of time delay and Doppler between desired and undesired signals varies.
Even during the days when BDS is evolving into its Phase III stage, there is a large user group still using the legacy BDS signal. The proposed model could be used in B1 receivers as a code-tracking performance indicator and an aid to satellite selection during navigation solution.
Finally, we combine Equation (A17) through Equation (A20) and will arrive at: sin 2 (π f T) +C · sin 2 (π f (K+1)T)+sin 2 (π f (19−K)T) This completes the development of the variance of correlator output when the data bits are misaligned by an integer multiple of primary code chips. Specifically, when C = 0, i.e., the differential time delay between undesired and desired signals are integer multiples of a code period, we arrive at: E |y k | 2 = T · T c [sin c( f T c )] 2 sin 2 (π f KT)+sin 2 (π f (20−K)T) which is a special case and we put it here for ease of reference.

Appendix C. Simulation Case For Section "Results"
The receiver is placed at 31 • 13.3 N, 121 • 27.12 E, Shanghai, China. The TLE retrieved for the simulation results is provided in Table A1.
For ease of reference, TLE elements used in simulation are extracted and titled.