GNSS data demodulation over fading environments: Antipodal and M ‐ ary CSK modulations

This article investigates new strategies to compute accurate low ‐ complexity log ‐ likelihood ratio (LLR) values based on the Bayesian formulation under uncorrelated fading channels for both antipodal and code shift keying modulations when no channel state information (CSI) is available at the receiver. These LLR values are then used as input to modern error ‐ correcting schemes used in the data decoding process of last ‐ generation Global Navigation Satellite System signals. Theoretical analysis based on the maximum achievable rate is presented for the different methods in order to evaluate the performance degradation with respect to the optimal CSI channel. Finally, the frame error rate simulation results are shown, validating the appropriate performance of the proposed LLR approximation methods.


| INTRODUCTION
Reliable and precise position, navigation and timing information is fundamental in safety-critical applications such as intelligent transportation systems, automated aircraft landing or autonomous unmanned ground/air vehicles (robots/ drones), to name a few. The main source of positioning information is provided by Global Navigation Satellite Systems (GNSS) [1], a technology which has attracted great interest in the last decade. Although most of the research has been in the signal processing aspects [2], in order to overcome the system limitations under non-nominal conditions [3,4], and data fusion strategies with alternative ranging technologies [5][6][7][8], a key part of GNSS receivers is the data demodulation stage which allows recovery of essential information. The latter has been long disregarded but may be a critical point in harsh environments, which is the main aim of this contribution.
In the last generation of GNSS signals, modern errorcorrecting codes (i.e. such as low-density parity check [LDPC] or convolutional codes) were considered in the GNSS signal design in order to enhance the data demodulation performance, especially over harsh scenarios [9,10]. The inputs to the corresponding soft decoding algorithms are the so-called loglikelihood ratio (LLR) values [11,12], which represent a statistical test to compare the goodness-of-fit between probabilities of receiving a positive or negative logic bit. These LLR values can be shown to be sufficient statistics for the decoding and detection processes [13]. Typically, in order to compute the LLRs, the entire knowledge of the propagation channel behaviour referred to as perfect channel state information (CSI) is considered. However, this assumption does not hold necessarily true in real-life applications since the CSI might not be fully available at the receiver [14], yielding a possible decoding loss due to the incorrect information at the decoding input. This situation can be further aggravated under urban environments where effects such as shadowing or multipath reduce the channel capacity at the receiver.
In this work, we focus on the uncorrelated fading channel [11] which is commonly used to model phenomena such as shadowing or multipath. This channel is modelled by a fading gain h n and an additive Gaussian noise with variance σ 2 . Note that if these parameters are perfectly known at the receiver (i.e. perfect CSI), the LLR has a well-known closed-form expression [11]. Otherwise, the LLR expression is unknown and LLR approximations are required. Indeed, when no CSI is available at the GNSS receiver, only one method to compute such LLR approximation is available in the literature [11], which has several limitations: (1) this method can only be used with antipodal modulations, and it is not valid for M-ary modulations, such as the code shift keying (CSK) modulation [15], and (2) it requires a high decoder complexity in order to compute the LLR approximation, which may not be useful in practice.
Bearing in mind the lack of practical solutions in the literature, the goal of this contribution is to provide new strategies to compute low-complexity closed-form LLR approximation expressions for both antipodal and CSK modulations, considering the uncorrelated fading channel and with no CSI available at the receiver. Thus, the article focuses on the following cases: � Antipodal modulation and GNSS pilot signal (3): Most of the new-generation GNSS signals are composed of a data and a pilot component; therefore, two LLR approximation methods are presented first considering the GNSS pilot signal case, which allows estimation of some channel parameters at the receiver. The first method was proposed in Ref. [11], based on Ref. [16], which maximizes the mutual information between the transmitted symbol and the LLR, and is given for completeness. This method provides a good LLR approximation at the expense of a high complexity burden. To reduce the computational complexity, and resorting to a Bayesian formulation, a second method based on our previous work [17] is presented. This method not only reduces the complexity but also provides similar LLR values. In order to evaluate the performance with respect to the perfect and statistical CSI solutions, a study of the capacity is also provided. � Antipodal modulation and GNSS data signal (4): Legacy GNSS signals may not have a pilot component and/or simple receiver implementations may only track the data component.
In that context, and with respect to the pilot signal case in Section 3, new alternatives must be studied. We propose two Bayesian approximations considering an uncorrelated fading channel, no CSI and a data signal component. Again, we also provide the channel capacity performance analysis of the proposed methods, with respect to the perfect CSI case and the Bayesian approximation with a pilot signal component. � M-ary CSK modulation and GNSS pilot signal (6): The CSK is a M-ary modulation that can increase the data rate without losing synchronization performance, thus being a suitable signal candidate for future GNSS applications [15]. We derive a new low-complexity LLR approximation expression that can be used over uncorrelated fading channels when no CSI is available at the receiver. Moreover, since this modern signal is expected to be transmitted as a data component [12,18] along with a GNSS pilot component, the latter can be exploited to infer some of the channel parameters. Finally, we also provide the channel capacity performance analysis of the proposed methods.
The article is organized as follows: the system model for the antipodal modulation and background on LLR expressions considering both perfect CSI and statistical CSI [17,19] are provided in Section 2; LLR approximations for an antipodal modulation without CSI considering a GNSS pilot signal are provided in Section 3, and without a pilot signal in Section 4; the system model for the CSK modulation and background on LLR expressions considering perfect CSI are provided in Section 5; the LLR approximations for the M-ary CSK modulation without CSI considering a GNSS pilot signal are provided in Section 6; the results are analysed for two uncorrelated fading channels in Section 7 and conclusions are drawn in Section 8.

| System model and LLR with perfect CSI
Current GNSS signals transmit binary data information; we assume the transmitted message as a binary vector u = [u 1 , …, u K ] of K bits. This message is encoded into a codeword c = [c 1 , …, c N ] of length N > K and mapped to antipodal symbols (e.g. binary phase-shift keying) x n = μ(c n ) ∈{−1, 1}, where we impose μ(0) = 1 and μ(1) = −1. Each symbol x n is then spread using a pseudo-random noise (PRN) sequence that can be expressed in vector form as p n ∈ R L , where L corresponds to the number of chips of the PRN sequence. Then, the transmitted symbol per coded bit is given as where convention vectors are defined as column vectors. Then, chip-level rectangular pulse shaping is used before transmission. Considering the data demodulation stage of a GNSS receiver, a key task is to obtain the posterior probability of a transmitted code symbol c n given the observed sample y n . The received signal models for two relevant (open sky and fading) scenarios are discussed in this subsection, after relevant LLR concepts are reviewed. The information in y n is used to compute the LLR value, defined for the nth symbol as This LLR can be used to feed the input of a Soft Input-Soft Output (SISO) error correction decoder [13]. Assuming that c n are identically and uniformly distributed (i.u.d.) [13] ∀n = 1, …, N, Equation (2) can also be written as where equiprobable symbols are assumed. Note that p(y n |x n ) represents the likelihood distribution given a transmitted symbol, which directly depends on the transmission channel.

| Open sky environment
Standard techniques typically assume an additive white Gaussian noise (AWGN) propagation model. Assuming perfect time and frequency synchronization, the received baseband symbol sequence at the chip-level can be written as where w n ∼ N ð0; L 2 ⋅ σ 2 I L Þ with I L being the identity matrix of size L. Thus, we denote the normalized output of the matched filter as y n ¼ y ⊤ n p n L ∈ R. Then, the normalized received symbol sequence is where w n ∼ N ð0; σ 2 Þ and σ 2 are known. In a standard GNSS receiver, this is the symbol model at the output of the prompt correlator, for which the channel is considered to be static over a symbol period. Note that all symbols y n are affected by the same noise statistics, then the Gaussian likelihood p(y n |x n ) is pðy n |x n Þ ¼ 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2πσ 2 and considering equiprobable symbols, the LLR can be computed as

| Fading environment
If we consider now a GNSS environment characterized by effects such as shadowing or multipath, the detection function typically used in this context corresponds to an uncorrelated fading channel with additive real-valued AWGN. Again, assuming perfect time and frequency synchronization, the received baseband symbol sequence at the chip-level can be written as where h n denotes the fading gain per chip and w n ∼ N ð0; L 2 ⋅ σ 2 I L Þ, with I L being the identity matrix of size L. Thus, we denote the normalized output of the matched filter as y n ¼ y ⊤ n p n L ∈ R. Then, the normalized received symbol sequence is where both w n and h n are independent random processes. w n are i.i.d. centred Gaussian random variables with variance σ 2 , that is, w n ∼ N ð0; σ 2 Þ. h n are also i.i.d. random variables with an associated probability density function (pdf) given by p(h), that is, h n ∼ p(h). It is assumed that h n ≥ 0 and change independently from one sample to another. All the symbols are again affected by the same noise statistics and the LLR simplifies to [16] L n ¼ 2 which explicitly implies perfect CSI, that is, h n and the variance σ 2 are assumed known. In practice, this assumption does not hold and even if σ 2 can be precisely estimated, the fading gain remains unknown in most of the situations.

| LLR with statistical CSI
A relaxation of the perfect CSI situation is to consider that full statistical CSI is available at the receiver, that is, σ 2 is assumed known or accurately estimated and h n is an unknown random quantity whose pdf and parameters are well characterized. Additionally, we consider a binary input memoryless channel where symbol x n is unknown. From a Bayesian perspective [20], since both x n and h n are unknown quantities, it is approriate to consider them as random variables. All the statistically relevant information about these variables is contained in their joint posterior distribution p(x n , h n |y n ).
Assuming that x n and h n are independent, pðx n ; h n |y n Þ ∝ pðy n |x n ; h n Þpðx n ; h n Þ ¼ pðy n |x n ; h n Þpðx n Þpðh n Þ; ð11Þ where the first term corresponds to the likelihood of observations given the unknowns and the second term represents the a priori knowledge on x n and h n . From Equation (9), the likelihood distribution is with mean μ = h n ⋅ x n and known variance σ 2 . According to the LLR definition (2), we are interested in the marginal distribution ORTEGA ET AL.
-115 which leads to ð15Þ since x n are equiprobable. In order to model a GNSS urban environment, it is common practice to use the two-state Prieto model [21]. Nevertheless, this model does not have a closedform expression for the channel gain pdf p(h). An alternative proposed in Ref. [12] is to consider a Rice distribution. However, as noted in Ref. [12], the LLR has no closed form, and in practice, its use is computationally too complex. In addition, the statistical CSI required to compute the LLR may not be available; therefore, different alternatives must be accounted for, which is the objective of the rest of the article.

| ANTIPODAL GNSS DATA DEMODULATION WITHOUT CSI FOR A PILOT SIGNAL COMPONENT
In this section, we present two methods to compute the LLR approximation when no CSI is available at the receiver, but some channel parameters can be inferred by means of a GNSS pilot component, which is typically available in new-generation GNSS signals. For instance, the GPS L1C signal is composed of two different components: a data component used to transmit the C/NAV-2 message [22] and a pilot component which transmits a secondary known code [22]. This secondary code can be used as a learning/training sequence as described in this section.

| Best linear LLR approximation
For completeness we introduce the method proposed in Ref. [16] and used for a GNSS data component in Ref. [12]. Motivated by Equation (10), this method computes the coefficient α that provides the best linear approximation of the LLR as The scaling factor α is obtained by maximizing the mutual information I L ; X ð Þ between the transmitted symbol X and the detector input L , both being random variables whose realizations x n and L n are observed at the receiver: where the mutual information is defined as with H X ð Þ and H X|L ð Þ the entropy of X and the conditional entropy of X given L , respectively. When considering a memoryless binary input symmetric output channel and consistent LLR values, this expression can be expressed as a function of the LLR pdf at the receiver input [23], considering {X = + 1}: This approximate criterion is derived from the capacity associated with a binary input memoryless channel, , for which the conditional pdf of the true LLRs has been replaced by the conditional pdf of the approximated ones. When considering C ¼ I L ; X ð Þ, it can be shown that the conditional pdfs given the true LLRs are both symmetric and consistent (see [16]).
The latter condition is not necessarily fulfilled for p c L |X � � and the quantity b I c L ; X � � cannot be interpreted as a true mutual information quantity. However, this quantity can be used as a good approximate measure, as proved in [16], where this quantity is maximized for the best linear LLR approximation (BLA). Notice that the proposed optimization method [16] assumes the knowledge of the linearly approximated LLRs conditional pdf. However, in real scenarios, this is unknown at the receiver and one has to resort to a numerical resolution (which is often computationally demanding) in order to estimate the approximated LLRs conditional pdf. To overcome this limitation, one can resort to the corresponding empirical mean estimator as done in Ref. [11,12] through the time average estimation proposed in Ref. [24]: where K is the number of samples used to estimate b I c L ; X � � . Notice that this method needs a learning sequence (i.e. known values x k ). This information can be directly obtained through the symbols of the GNSS pilot component. Finally, in order to compute α, one can apply a one-dimensional search method [25] based on the objective function (20). We underline that this method does not require knowledge of σ 2 ; therefore, no CSI is required.

| Bayesian LLR linear approximation
Recall from Section 2.2 that the problem of computing the LLR values involves solving the integrals in Equation (14), for which we have to select a prior distribution for h n . In [11], the pdf p(h) was selected to be a Rice distribution, leading to a complex LLR expression for practical applications. In Bayesian inference, a common approach is to select a prior distribution to be a conjugate of the likelihood distribution, which results in a posterior distribution that is of the same family as the a priori, where only the parameters need to be updated [26]. This idea was exploited in [17] to provide a simple low-complexity closed-form LLR approximation for M-ary modulations over uncorrelated fading channels. Given that the likelihood (12) is Gaussian, the conjugate prior distribution for h n in Equation (11) is also Gaussian [20]: where the pdf parameters (i.e. μ h and σ 2 h ) need to be adjusted according to the unknown parameters' uncertainty. The marginal distribution in Equation (14) is obtained by solving the integral which can be shown to be another Gaussian (refer to [17]): and the corresponding LLR are where similar to Equation (16), the resulting LLR approximation is a linear function of y n . Additionally, we underline that under the AWGN channel (i.e. μ h = 1 and σ 2 h ¼ 0) the LLR in Equation (24) reduces to the Gaussian LLR solution in Equation (7). Notice that Equation (24) requires the knowledge of σ 2 , as well as μ h and σ 2 h . The latter values can be obtained by resorting to the maximum likelihood (ML) estimates [17]: where K is the number of samples used to estimate μ h and σ 2 h , x k is the kth symbol of a binary learning sequence (i.e. the known secondary code from the GNSS pilot component) and y k is the kth received pilot symbol. Finally, the Bayesian LLR linear approximation is given as Figure 1 summarizes the different approaches in Section 3.

| Performance of the LLR approximation methods for the antipodal modulation and a pilot signal
In this section, we evaluate the LLR approximations over wellknown uncorrelated fading channels in order to compare the channel capacity performance with the perfect CSI and statistical CSI cases.

| Normalized Rayleigh fading channel
The normalized Rayleigh fading channel [27] is typically used to describe phenomena such as shadowing or multipath, and it is well-known to provide a closed-form solution of the LLR (14), that is, when full statistical CSI is considered. As a consequence, we can obtain an upper bound for the LLR approximation performances. Following [27], the LLR expression (14) when considering a normalized Rayleigh fading channel is . ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi . ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where ΦðzÞ ¼ 1 þ ffi ffi ffi π p ze z 2 erfcð−zÞ and erfc(⋅) represent the complementary error function. We evaluate the upper bound of the maximum achievable rate R 0 [28] of the LLR linear approximations as well as the LLR expression with full statistical CSI (27) and perfect CSI (7), (10). To compute such maximum achievable rate, we follow the methodology proposed in [24] based on the extrinsic information transfer chart analysis. Rayleigh channel, the R 0 upper bound is shown for (iii) perfect CSI (10), (iv) full statistical CSI (27), (v) Bayesian LLR linear approximation (26), and (vi) the LLR approximation (16), considering that the mutual information I L ; X ð Þ in Equation (17) is computed from Equation (19) when 12 s of pilot symbols are retrieved. We can observe that: (i) Channel capacity loss caused by the fading effect: Considering an ideal coding scheme of rate R = 1/2, the channel capacity loss between the AWGN channel and the normalized Rayleigh channel (when perfect CSI solutions are assumed) is around 1.6 dB. Notice that this loss can be reduced when applying lower rate channel coding schemes.
As an example, the channel capacity difference with a channel coding scheme of rate R = 1/3 is around 0.9 dB. (ii) Channel capacity loss due to channel uncertainty: Considering an ideal coding scheme of rate R = 1/2, the channel capacity loss is around 0.8 dB between the full statistical CSI solution and the perfect CSI solution. Moreover, the best and Bayesian LLR linear approximations provide the same channel capacity than the full statistical CSI solution, proving that when no perfect CSI is available, only the first-and second-order moments of the fading distribution are required to achieve an optimal solution. Finally, when the transmission channel is characterized by an AWGN, the Bayesian solution (26) converges to the perfect CSI solution (7).

| Two-state Prieto channel
In a second scenario, we propose to evaluate the performance in a more realistic GNSS urban scenario. We consider a twostate Prieto channel model [21] for a vehicle speed of 50 km/h and a satellite elevation angle of 40°. In this example, we consider the data component of the signal GPS L1C which is characterized by a symbol rate of 100 symbols/s and a PRN code of length 10,230 chips. The results are shown in Figure 2 (bottom plot), where we illustrate again the upper bound of the maximum achievable rate. In this case, considering perfect CSI and a coding scheme of rate 1/2, the channel capacity loss is around 6 dB with respect to the AWGN scenario. Notice that for this type of scenarios, errorcorrecting schemes with lower rates are highly recommended. On the other hand, we verify that the best and Bayesian linear approximations, Equations (16) and (24)) almost converge to the perfect CSI solution (loss around 0.8 dB), proving the validity of such approximations. Finally, we underline that no full statistical CSI expression is available since p L |X ¼ þ1 ð Þ has no closed form and is unknown at the receiver.
To conclude, we provide a brief comment on the complexity of the LLR linear approximation methods (refer to Figure 1): (1) the first method requires an online estimation technique which needs to resort to an iterative one-dimensional search method, based on a cost function involving log/ exp function evaluations. Note that the complexity of this method directly depends on the number of samples K to estimate b I ð c L ; XÞ; (ii) the second method avoids to compute the one-dimensional search, but instead the first-and second-order moments of the fading distribution have to be estimated. Note from Equation (25) that only simple arithmetical operations are required. Again, the complexity of the method depends on the number of samples K used to estimate μ h and σ 2 h , but for equal number of samples, the Bayesian solution is computationally less expensive than the best linear one.

| ANTIPODAL GNSS DATA DEMODULATION WITHOUT CSI FOR A DATA SIGNAL COMPONENT
The previous section focused on LLR approximations when a pilot component (i.e. training sequence) is available. However, legacy GNSS signals may not have a pilot component; therefore, different alternatives must be accounted for. In the sequel, we introduce data demodulation strategies considering an uncorrelated fading channel, no CSI and a data signal component.

| Bayesian LLR approximation without training data: MLE of μ h and σ 2 h through expectation-maximization
In contrast with Section 3.2, where the underlying pilot signal assumption allowed to compute the ML estimates in (25) which are in turn used to compute the LLR approximation (26), in this case, we do not have access to such training sequence. Therefore, we propose a method to derive the μ h and σ 2 h ML estimates when no learning sequence is available. The marginal likelihood p(y n ) = p(y n |x n = 1)p(x n = 1) + p (y n |x n = −1)p(x n = −1) is a mixture of two Gaussian distributions: with and obtain the ML estimates of μ h and σ 2 a as the roots of the partial derivatives with respect to μ h and σ 2 a . The partial derivative with respect to μ h is and we have to solve for dlog Λ ð Þ dμ h ¼ 0, which has no analytical solution. However, conditional on a specific realization of the latent variables x n we could use the μ h estimate from Section 3.2. We first obtain the (discrete) posterior distribution of x n given the observations: and Then, Equation (30) can be expressed as and we can estimate μ h as Similarly, we can proceed to estimate σ 2 a : and obtain σ 2 a as (i.e. dlog Λ ð Þ dσ 2 Notice that in the previous equations, it was not taken into account that γ x 1 and γ x −1 depend on the unknown parameters. An iterative expectation-maximization (EM) algorithm can be used in order to estimate b μ h and b σ 2 a . The resulting EM algorithm is as follows: 4. Evaluate the log-likelihoods with the updated parameter estimates. If the log-likelihood change is below a given small threshold ϵ, stop. Otherwise, go back to the E-step.
Finally, the resulting expression for the LLR values with no CSI and a data channel signal is

| Bayesian LLR approximation without training data: rough estimation of μ h and σ 2 h
To avoid the iterations in the previous EM-based solution, a simpler low-complexity approach to derive μ h and σ 2 h ML estimates is proposed. The ML estimates in Equation (25) can be approximated as where y n is the nth received symbols, and N is the number of received symbols used to estimate μ h and σ 2 h . The corresponding LLR expression without CSI and a data signal is Notice that the approximated ML estimates in Equation (38) are biased, in contrast to the true ones in Equation (25). As expected, this bias is made apparent only at low signal-to-noise ratio (SNR) values. As shown in Section 4.3, this has a minor impact on the maximum achievable rate R 0 with respect to the perfect CSI case, but a clear impact is present at low SNR for the fading distribution first and second order estimation performance (refer to Figure 4). The fact that such bias has almost no impact on the maximum achievable rate translates into a marginal impact in the frame error rate (FER) analysis shown in Section 7.1.

| Performance of the LLR approximation methods for the antipodal modulation and a data signal
As described in Section 3.3 for the pilot signal scenario, we can assess the upper bound of the maximum achievable rate R 0 for approximations (37) and (39); the results are shown in Figure 3. The two approximations are compared to the perfect CSI case for both AWGN and Rayleigh channel, as well as to the Bayesian LLR linear approximation (26) considering a pilot signal.
Notice that all the LLR approximation methods converge to the same solution. However, it is well-known that the finite length codes used in real coding schemes are not optimal, thus it is likely that in real scenarios, these methods will not be equivalent. Then, in order to evaluate this issue, we illustrate in Figure 4 the estimation accuracy of the firstand second-order moments of the fading distribution when 12 s of signal are retrieved (i.e. obtained from 2000 Monte Carlo runs). Note from these results that when a pilot is available at the receiver, an accurate estimation of the firstand second-order moments is achieved independently of the SNR. On the other hand, when no pilot is available, the EMbased method does not correctly estimate the parameters in low SNR regimes. In addition, the estimation accuracy using the rough approximation (38) is significantly degraded and only performs well for high SNR. Thus, it is expected that the error-correcting performance of real channel coding schemes when considering both LLR approximations (37) and (39) will be degraded with respect to the Bayesian approximation exploiting a pilot signal (26). Because the EMbased solution performs better than (38), we expect the error correction performance using the former (37) to be better than with the latter (39).

| SYSTEM MODEL FOR THE M-ary CSK AND LLR WITH PERFECT CSI
The CSK modulation [29] is a M-ary orthogonal modulation which was first proposed as a GNSS signal candidate in [15]. Each symbol CSK x ℓ corresponds to a different circular shift of a unique PRN sequence c. Let S ℓ = {ℓ, 1 ≤ ℓ ≤ 2 Q = M} be the set of data symbols, with Q the number of bits to be transmitted, then the PRN sequence c ℓ associated to the symbol x ℓ , ℓ ∈ S ℓ , satisfies the following rule:  where i represents the PRN chip, m ℓ is the integer number corresponding to the ℓth symbol shift, L is the number of chips in the PRN sequence and mod(x, y) is the modulus operation. As an example, in Figure 5, it is illustrated the PRN sequences associated to the 4-ary CSK modulation with a number of chips L = 10,230.
At the transmitter, the information bits are usually encoded by an error correction code, generating a codeword of length N bits. Then, the N codeword bits are grouped in N/Q CSK symbols of Q bits. Finally, the M-ary CSK associates each CSK symbol x ℓ to a PRN sequence c ℓ by right shifting the fundamental PRN sequence c.

| CSK data demodulation in open sky environments
As presented in Section 1, an open sky environment can be modelled by an AWGN channel. Then, assuming an AWGN channel and perfect time and frequency synchronization, the nth received sequence y n corresponding to the transmitted CSK PRN sequence c n,ℓ associated to the vector sequence [b 1 , b 2 , …, b Q ] and the CSK symbol x ℓ can be expressed as y n;i ¼ c n;i;ℓ þ n n;i ; where i represents the PRN chip, n n;i ∼ N ð0; σ 2 Þ are zeromean i.i.d. Gaussian noise samples with variance σ 2 = N 0 /2. Note that coherent reception is a valid assumption since a GNSS receiver capable to demodulate the CSK signal is also tracking in parallel the pilot signal component, which may provide a precise phase estimation. Let us now define x j , 1 ≤ j ≤ 2 Q−1 the transmitted symbol if b q = 1, and x t , 1 ≤ t ≤ 2 Q−1 the transmitted symbol if b q = 0. Therefore, considering perfect synchronization and following the LLR derivation in [18], the LLR expression for the bit b q is given as where P(b j,z ) denotes the probability of b j,z which is zth bit of the transmitted symbol x j and P(b t,z ) denotes the probability of b t,z which is the zth bit of the transmitted symbol x t . Moreover, the term 1 L ∑ L i¼1 y n;i c n;i;x l corresponds to the normalized correlation between the nth transmitted and nth received PRN sequences. When a Bit Interleaver Coded Modulation (BICM) scheme [18,28] is implemented at the receiver and equiprobable transmission bits are assumed, the LLR can be simplified to

| CSK data demodulation in fading environments with perfect CSI
Assuming the uncorrelated fading channel defined in Section 2, the received sequence can be written as y n;i ¼ h n ⋅ c n;i;ℓ þ n n;i ; ð44Þ where n n;i ∼ N ð0; σ 2 Þ are zero-mean i.i.d. Gaussian noise samples with variance σ 2 = N 0 /2 and h n is the fading gain, which is assumed to be invariant within the symbol, and is also defined as an i.i.d. random variable with associated pdf give p (h) and h ≥ 0. Again, coherent reception and perfect synchronization are assumed. The LLR expression for the bit b q over the uncorrelated fading channel can be derived from Equation (42) as where h n is considered known at the receiver. Note that the previous equation can be simplified considering a BICM scheme:

F I G U R E 4
Fading distribution of first-order (top) and second-order (bottom) moment estimation, using the Bayesian log-likelihood ratio approximation methods with pilot and data signal components.
ORTEGA ET AL.
where perfect CSI is assumed in order to compute the LLR.

| A BAYESIAN APPROACH TO CSK DEMODULATION IN FADING ENVIRONMENTS WITHOUT CSI
In this sequel, we derive the LLR values for a CSK modulation considering that no CSI is available at the receiver. In that perspective, we adapt the Bayesian method in Section 3.2 to compute a closed-form LLR expression. From the definition of the LLR in Equation (3), we redefine the LLR expression for the bit b q over the uncorrelated fading channel as where pðy n |b q ¼ 1; We follow the approach in [17] and consider a conjugate prior distribution for h n : whereas in Equation (21), μ h and σ 2 h need to be adjusted according to the channel uncertainty. From Equations (48) and (49), we are interested in In order to compute the LLR expression, we only need to compute those terms which depend on c n;i;x j and c n;i;x t (refer to Appendix 9): Finally, Equation (47) is given by Considering that the qth bit of the symbol x j is always equal to 1 and the qth bit of the symbol x t is always equal to 0, then Equation (55) is written as Considering that a BICM scheme is implemented at the receiver, the previous expression can be written as Notice that to compute Equation (57), several exponential operations are required. A useful metric to reduce the computational complexity is the log-sum approximation [30], and the LLR approximation is The previous expression avoids the use of log/exp function evaluations, reducing the computational burden at the receiver. Notice that the first-order μ h and second-order σ 2 h moments of p (h) are assumed to be known, that is, partial statistical CSI is assumed. However, these parameters might not be available at the receiver and therefore they must be estimated online. Assuming a binary learning sequence (e.g. symbols of a pilot component), we can infer μ h and σ 2 h as in Equation (25). Moreover, this LLR approximation considers σ 2 known at the receiver, as typically done in the literature [17,27], a result which also holds true in those scenarios where σ 2 was precisely estimated before the fading effect. In any case, if one wants to avoid the knowledge of σ 2 and provide an LLR without CSI close-form expression, Equation (58) can be replaced by y n;i c n;i;xt þβ⋅ 1 L ∑ L i¼1 y n;i c n;i;xt where b μ h is estimated as in Equation (25), K is the number of pilot symbols used to estimate the channel fading parameters and β is a coefficient which weighs the second term in Equations (53) and (54). Notice that β = 0 involves neglecting the second term in Equations (53) and (54) and provides an LLR approximation based on the metric obtained in Equation (24), where the observed symbol is the output of the matched filter of the CSK demodulator 1 L ∑ L i¼1 y n;i c i;x l . On the other hand, high values of β can mask the information provided by the first term in Equations (53) and (54), generating an inaccurate LLR approximation. Based on simulations, we have found that β = 1 is the optimal value, yielding Equation (59) finally to

| Performance of the LLR approximation methods for the CSK modulation
As seen in Section 3.3, we compute the maximum achievable rate of the Bayesian LLR approximation method (60) for different CSK modulation orders Q = {2, 4, 6}, considering a BICM CSK demodulator, and for both uncorrelated normalized Rayleigh fading channel and two-state Prieto channel. These results are summarized in Figure 6. The new Bayesian CSK demodulation is compared to the perfect CSI cases given by Equations (43) and (46). Both fading scenarios are also compared to the LLR approximation case (58), where partial statistical CSI is available at the receiver (i.e. where μ h , σ 2 h and σ 2 are assumed to be known).

| Normalized Rayleigh fading channel
As seen in Section 3.3, two different effects can cause a channel capacity loss. In Figure 6 (top), we can first see the impact of the fading channel, which induces for an ideal coding scheme of rate R = 1/2 a loss of 1.2, 1.3 and 1.4 dB for modulation orders of Q = 2, Q = 4 and Q = 6, respectively, with respect to the AWGN case. Moreover, an additional 0.5 dB is lost due to channel uncertainty. Note that the channel capacity loss can be reduced by using coding scheme of lower rates. The latter is highly recommended for modulation orders greater than Q = 2, because transmitting more bits in one symbol increases the demodulation threshold.

| Two-state Prieto channel
Results for the maximum achievable rate considering a twostate Prieto channel model, for a vehicle speed of 50 km/h and a satellite elevation angle of 40°, are shown in Figure 6 (bottom). As an example, we consider the data component of the GPS L1C signal which is characterized by a symbol rate of 100 symbols/s and a PRN code of length 10,230 chips. With respect to the Rayleigh fading channel, we can see a channel capacity loss which is further degraded. For an ideal channel coding scheme of rate R = 1/2, a loss of around 6 dB is found for the CSK modulation scheme of order Q = 2, Q = 4 and Q = 6. Moreover, an additional 1 dB is lost due to the channel uncertainty. Note that in this limit fading channel scenario, high-order CSK modulations are not recommended without low-rate channel coding schemes. Furthermore, due to the bad channel quality, specific channel coding structures such as rate compatible Root-LDPC codes [9], which allows retrieving the entire diversity of the channel, are also highly recommended.

| RESULTS: FER PERFORMANCE FOR THE CED WITH THE GPS L1C SIGNAL AND LDPC CODES
In this section, we compare soft decoding performance corresponding to the different LLR approximations introduced in the previous sections. Particularly, as an example, we provide the FER (i.e. with respect to the carrier-to-noise density [C/ N 0 ]) performance for the clock and ephemeris data (CED) considering the GPS L1C signal [31] with irregular LDPC codes, decoded with a sum-product algorithm [23].

| Results for antipodal GNSS modulations
First, we consider a normalized Rayleigh fading channel and a channel coding scheme based on the standard irregular LDPC code of rate 1/2 used to encode the GPS L1C subframe 2 [31]. The FER results are summarized in Figure 7, where the performance of the LLR approximation in Section 3 is shown, (16) (best LLR linear approximation [BLA]) and (26) (Bayesian), when 12 s of pilot symbols are retrieved. For comparison, we also show the FER results corresponding to the perfect CSI LLR (10) (Perfect CSI) and the full statistical CSI LLR (27) (Stats CSI). From the FER results we can see that both BLA and Bayesian approximation methods show a similar data demodulation performance w.r.t the statistical CSI case, which proves that when no CSI is available only a good estimation of the first and second order moments of the fading distribution are required. Moreover, when the channel transmission is characterized by an AWGN channel, the LLR approximation in (26) converges to the perfect CSI LLR solution (7). Figure 8 shows the same comparison but for a data signal component. For comparison, we show the FER results corresponding to the LLR approximation with a pilot component (i.e., with a learning sequence and 12 s of pilot symbols retrieved) (26). Considering that no pilot sequence is available, we show the FER performance for the approximations introduced in Section 4, (37) (No Pilot EM) and (39) (No Pilot Approx). Notice first that the method which estimates the first and second order moments of the fading distribution through the EM algorithm provides a FER very similar to the one obtained with the pilot case (26), with only a 0.1 dB performances loss. On the other hand, the simpler method which estimates the corresponding first and second order moments of the fading distribution based on the absolute value of the received symbol reduces the data demodulation performance to around 0.4 dB. Note from Figure 4 that the LLR values from the methods in Section 4 were expected to provide lower error correction performance than the method obtained with the pilot case (26).
As a second scenario, we consider a two-state Prieto channel model for a vehicle speed of 50 km/h and a satellite elevation angle of 40°. We show the FER results again for the GPS L1C data within subframe 2, but in this case with a regular rate compatible Root LDPC codes of rate R = 1/3 (proposed in [9] to reduce the demodulation threshold) channel coding scheme. The comparison of the different methods with and without a pilot signal is illustrated in Figure 9. Notice that the EM-based method provides a FER very similar to the one obtained with the pilot case (26). In addition, a very low performance loss is observed with respect to the perfect CSI case. On the other hand, the simpler approximation (39) reduces the data demodulation performance around 1 dB with respect to the EM-based solution.

| Results for M-ary CSK modulations
As in the previous subsection, we first consider a normalized Rayleigh fading channel, with a channel coding scheme based on the standard irregular LDPC used to encode the GPS L1C subframe 2. The FER results for the CSK modulation with Q = {2, 4, 6, 8, 10} are summarized in Figure 10 of pilot symbols are retrieved (60). For comparison, we also show the FER results corresponding to the perfect CSI LLR (46) and the perfect CSI LLR values under AWGN channel (43). From the FER results, we can see that both partial statistics CSI (58) or no CSI (60) achieve a similar data demodulation performance, independently of the CSK modulation order. On the other hand, we can appreciate a FER performances loss in the order of 0.6-0.8 dB due to the channel uncertainty. The channel uncertainty impact seems to increase along with the CSK modulation order. Finally, we can see a significant FER performance loss due to the fading effect, that is, in the order of 2-4 dB. Again, this loss of performances is related to the modulation order. From these results, we suggest the use of channel coding schemes of a lower rate for higher order CSK modulations.
The results for a second scenario, considering a two-state Prieto channel model for a vehicle speed of 50 km/h and a satellite elevation angle of 40°, and a channel coding scheme based on rate compatible Root LDPC codes of rate R = 1/4, are summarized in Figure 11. We show the FER performance results for CSK modulations with Q = {2, 4, 6} * and the different approximations in Section 6: (i) partial statistical CSI (58), and (ii) when 12 s of pilot symbols are retrieved (60). For comparison, we also show the FER results corresponding to the perfect CSI LLR case (46). From the FER results, we can see that both (58) and (60) achieve almost a similar data demodulation performance independently of the CSK modulation order, as for the previous normalized Rayleigh fading channel scenario. On the other hand, we can appreciate an FER performance loss in the order of 1.5-2 dB due to the channel uncertainty. Again, these results suggest that a low-rate error correcting schemes should be considered in order to reduce the demodulation threshold. F I G U R E 7 Frame error rate of standard GPS L1C clock and ephemeris data over a normalized Rayleigh channel for antipodal modulations and a pilot signal component.

F I G U R E 8
Frame error rate of standard GPS L1C clock and ephemeris data over a normalized Rayleigh channel for antipodal modulations and a data signal component.

F I G U R E 9
Frame error rate GPS L1C clock and ephemeris data with a regular rate compatible Root-LDPC code of rate R = 1/3 over a two-state Prieto channel, for antipodal modulations.  . Since the product of two Gaussian distributions is in turn a Gaussian distribution, we proceed by finding the resulting mean (μ a ) and variance (σ 2 a ) as where κ n is an auxiliary constant, and after identifying terms on both sides of Equation (62) Then, the constant κ n can be computed as where the pdf definition is applied [17]. Note that in order to compute the LLR expression, we are only interested in those terms which depend on c n;i;x j . Then, we define the constant κ n,1 as those values of κ n which depends on c n;i;x j , κ n;1 ¼ − 2μ h 1 N ∑ N i¼1 y n;i c n;i;x j þ and the corresponding κ n,2 for c n;i;x t (see Equations (53) and (54)).
ORTEGA ET AL.