ABBIA Myriam Foucras

Olivier Julien

Christophe Macabiau

Detailed Analysis of the Impact of the Code Doppler on the Acquisition Performance of New GNSS Signals

For several tens of years, when the main acquired signal was GPS L1 C/A and under normal circumstances, many results were derived assuming that there is no change in the received code frequency due to the incoming Doppler frequency. Nevertheless, [1] and [2] investigated the effect of the resulting code phase sliding during the correlation process and showed that the acquisition performance can significantly be degraded. Considering the new challenges (new environments and new modulations, signal structures, longer codes, faster code rate, and so on), this investigation needs to be reviewed. In this paper, the authors want to present some results dealing with the code Doppler effect on acquisition performances for several new GNSS signals.

INTRODUCTION

The number of transmitted GNSS (Global Navigation Satellite Systems) signals has significantly increased with the ongoing deployment of Galileo and Beidou constellations, the maintenance of GLONASS and the modernization of GPS. In order to share the same Lower band (L-band), the properties of these new GNSS signals were adapted. The main points of design and characterization of GNSS signals are:

-The carrier frequency -The spreading codes characterized by their length , code frequency also called chipping rate and which is the number of chips transmitted in 1 second and repetition period to improve [START_REF] Jiao | High Sensitivity GPS Acquisition Algorithm Based on Code Doppler Compensation[END_REF]. The acquisition process in these degraded conditions requires large total integration time (through coherent and/or non-coherent accumulations) to accumulate signal energy. The Doppler frequency, mainly caused by the satellite motion and the local oscillator [START_REF] Van Diggelen | A-GPS assisted GPS, GNSS, and SBAS[END_REF], impacts twice the emitted signal: by modifying the central frequency (change estimated by the acquisition process) and by modifying the code frequency, resulting in a code Doppler. If the incoming Doppler frequency is small, the code Doppler can be ignored. However, in high Doppler conditions, it can be significant. Indeed, the code Doppler leads to change the spreading code period (reduced or expanded). For instance, for an incoming Doppler frequency of 10 kHz, it translates into one spreading code chip slip in 154 ms for GPS L1 C/A and for the same Doppler frequency, the 1 chip slip occurs in only 12 ms for GPS L5. In this study, the Doppler frequency is assumed to be constant over the total integration time because the average rate of magnitude change of the Doppler frequency is less than 1Hz per second [START_REF] Van Diggelen | A-GPS assisted GPS, GNSS, and SBAS[END_REF] If the receiver does not take the code Doppler into account during the acquisition phase (which is the general case), this can result in degraded correlator outputs that will degrade the acquisition capability of the receiver. It is thus relevant to discuss the code Doppler effect on GNSS acquisition performance:

-Degradations on the acquisition performance in terms of probability of detection if code Doppler is not compensated by the receiver -Development of code Doppler compensation methods This paper specifically focuses on by providing a mathematical study and some innovative results for widely used new signals. The motivation behind this investigation is to give careful instructions concerning the acquisition of GNSS signals and in particular in low C/N 0 environments. The accurate knowledge of the code Doppler impact on acquisition allows adapting and improving acquisition methods. This paper is divided as follows:

-A first section serves as establishing the basis by presenting the considered GNSS signals, acquisition process and code Doppler problem leading to the problematic of the paper -A second section develops the mathematical model of the distorted autocorrelation when considering a received signal affected by code Doppler and a local one -Then, simulation results are presenting in the third section leading to the comparison of coherent/non-coherent accumulation methods and signals -A fourth section is dedicated to the state-of-theart of code Doppler compensation methods -At the end, the conclusion will summarize the main results and present future works

I. ACQUISITION AND CODE DOPPLER

GNSS Signals

Several GNSS signals are considered in this study. They encompass most of the GPS and Galileo civil signals: GPS L1 C/A, L1C, L5 and Galileo E1 OS and E5 whose frequency plans are presented in Figure 1.

Figure 1: GPS and Galileo frequency plans [START_REF]European GNSS (Galileo) Open Service Signal In Space Interface Control Document[END_REF] In order to develop the code Doppler study and bot being dependent from other facts degrading acquisition performance, some assumptions on the considered GNSS signals are taken. Because the code Doppler impact on the acquisition is similar on both components (for composite signals) and without loss of generalities, only results on one component are presented (data component). Moreover, the data sequence and secondary code are supposed to be constants to 1 to avoid degradations due to bit transitions [START_REF] Foucras | An Efficient Strategy for the Acquisition of Weak Galileo E1 OS Signals[END_REF]. Table I gives a brief overview of the data component signal features; [START_REF] Avila-Rodriguez | On Generalized Signal Waveforms for Satellite Navigation[END_REF] and the Interface Control Documents are reference documents for more details: [START_REF] Navstar | Navstar Global Positioning System GPS Space Segment / Navigation User Interfaces[END_REF], [START_REF] Navstar | Global Positioning Systems Directorate Systems Engineering & Integration Interface Specification IS-GPS-800 navstar GPS Space Segment/User Segment L1C Interface[END_REF], [START_REF] Navstar | Global Postioning Systems Directorate Systems Engineering & Integration Interface Specification IS-GPS-705 Navstar GPS[END_REF] and [START_REF]European GNSS (Galileo) Open Service Signal In Space Interface Control Document[END_REF] for GPS L1C/A, L1C, L5 and Galileo E1/E5 respectively. The data component of the received GNSS signal can be represented as follows over a time interval with constant Doppler:

( ) ( ) ( ) ( ) ( ) ( ( ) ) ( ) 
Where:

is the signal amplitude is the data sequence is the secondary code ( ) for GPS L5 -is the spreading code (PRN) periodic then ( ) ( ) for all is the incoming Doppler frequency whose uncertainty space is [ ] kHz is the incoming code delay is the initial phase of the incoming signal at the beginning of the time interval is the incoming noise which is assumed to be a centered Gaussian distribution white noise with a constant two-sided power spectral density equal to is the resulting subcarrier for BOC-modulated signals, if any, using:

( ) ( ) ( ( )) 
Table II focuses on the spreading code properties for the data component to complete the previous signal expression. The spreading code chips number is denoted by and the spreading code chip duration leading to a spreading code period equal to . Assuming the modulation of Galileo E5 a-I as a BPSK, only results on GPS L5 are presented because the Galileo E5a/b-I signal has the same characteristics as the GPS L5 signal so what is true for GPS L5 is also valid for Galileo E5a-I.

Acquisition

At the receiver level, through a correlation process as depicted in Figure 2, the acquisition process consists in giving a rough estimation of the incoming signal parameters. At the end, the in-phase ( ) and quadrature ( ) correlator outputs, resulting of the correlation process, are expressed by:

{ ( ) ( ) ( ) ( ( ( ) ) ) ( ) ( ) ( ) ( ) ( ( ( ) ) ) ( ) (3)
Where -| | is the code delay error -| ̂| is the Doppler frequency error and are the noise at the correlator output which follows a centered Gaussian distribution and their variances are [START_REF] Bitar | Advanced GPS signal processing techniques for LBS services[END_REF] -[ ( ) ] is the integration interval beginning at and lasting -( ) is the autocorrelation function of the spreading code defined by:

( ) ∑( ( ) ( ))( ( ) ( )) (4)
Figure 3 presents the autocorrelation functions for all of the considered GNSS signals. For BOC-modulated signals (Galileo E1 OS and GPS L1C), the autocorrelation function is sharpest than BPSK-modulated signals (GPS L1 C/A and GPS L5). A classical acquisition method, the serial search acquisition method is based on an acquisition grid which covers the frequency and time uncertainty search space. Then, the incoming signal is correlated with a locally generated one which the estimation parameters couple ( ̂ ) is read one by one in the acquisition grid. It is quite clear that a trade-off should be chosen between the search space discretization and the acquisition execution time (in general wanted as short as possible).

We decide, as it is done in [START_REF] Rtca | Assessment of Radio Frequency Interference Relevant to the GNSS L1 Frequency Band RTCA/DO-235B[END_REF], that the total maximal degradations due to the discretization acquisition grid should be lower than 3.4 dB (2.5 dB in the time space and 0.9 in the frequency space). This determines the width of a bin in the acquisition grid in the frequency and time dimension by assuming uniform distributions. If only 2.5 dB losses are accepted in the time domain (equivalent to ( ) ) then the bin width for the code delay should be lesser than:

for BPSK-modulated signals for GPS L1C -for Galileo E1 OS

In general, to evaluate the acquisition process, two criteria are used: the mean acquisition time and the probability of detection which is the probability to correctly acquire the signal under some conditions (sensitivity, integration time…). Herein, because we want to evaluate the degradations on the detection capability due to the code Doppler, we focus on the probability of detection.

To determine it, a detection hypothesis test is applied based on Neyman-Pearson's approach [START_REF] Bastide | Analysis of L5/E5 Acquisition, Tracking and Data Demodulation Thresholds[END_REF]:

-The null hypothesis assumes that the useful signal is not present (or its parameters are not well estimated), in other words, there is only noise -The alternative hypothesis assumes that the useful signal is present and its parameters are correctly recovered

-1 0 1 -0.5 0 0.5 1  (chip) R(  ) GPS L5-Galileo E5a/E5b -1 0 1 -0.5 0 0.5 1  (chip) R(  ) GPS L1 C/A 0.75 -1 0 1 -0.5 0 0.5 1  (chip) R(  ) GPS L1C -1 0 1 -0.5 0 0.5 1  (chip) R(  ) Galileo E1 OS
The acquisition detector is:

( ) ( ) ( ) ( ) ( ) ( ) ( ) (5)
The acquisition test becomes: against (6) Where is the acquisition threshold Knowing the probability of false alarm and fixed at , the threshold can be easily determined:

( ( ) √ √ ) ( ) ( √ ) (7)
The normalized acquisition detector under the null hypothesis is a chi-square distribution because it has the distribution of ( ( ) ( ))

. Two squared centered Gaussian distribution are summed then the number of degrees of freedom is 2. Then, the probability of detection, the probability that the acquisition detector under the alternative hypothesis , , exceeds the threshold is:

( ( ) √ √ ) ( ) ( √ ) (8)
Where is the non-centrality parameter equal to:

( ) ( ) (9)
The integration time is chosen to be a multiple of the spreading code period (

) and in general it is one. But in low SNR environment, signal acquisition cannot be successful in one spreading code period. Then, to reach higher probabilities of detection and acquisition sensitivity, it is possible to increase the total integration time (that means more than one spreading code period should be accumulated). Commonly used methods are the coherent integration and non-coherent summations. In both cases, the total integration time is spreading code periods. In general, the data bit transition is one of the main factors that limit the integration time but in the context of this paper, as previously said, we assume that there is no data.

1) Non-coherent integration Successive correlator outputs on the spreading code period are non-coherently summed; for non-coherent summations on the spreading code period , the acquisition detector (5) becomes:

∑( ( ) ( )) (10)
In this case, the number of degrees of freedom of the chisquare distribution under and is . 2) Coherent integration Only one integration is done but on . Then, the coherent integration interval as defined in (3) becomes:

[ ( ) ] (11)
The distribution of the acquisition detector is a chi-square with 2 degrees of freedom. To only study the acquisition degradation due to code Doppler, we assume that the incoming Doppler frequency and phase are well estimated. Moreover, the code delay is supposed to be null at the beginning of the correlation process; that means that the residual code delay error is only induced by code Doppler. To speed up the acquisition process, the Fourier transform is often used [START_REF] Lin | Comparison of acquisition methods for software GPS receiver[END_REF]. In this way, for a given Doppler frequency estimate ̂, the correlator outputs for all of the code delay estimates ̂ are computed in the same time. The results on acquisition performances obtained with Matlab use Fourier transform.

Code Doppler

The code Doppler, noted , is the code frequency derived from the incoming Doppler shift [START_REF] Tsui | Fundamentals of Global Positioning System Receivers: A Software Approach[END_REF]. As it is defined in [START_REF] Rtca | Assessment of Radio Frequency Interference Relevant to the GNSS L1 Frequency Band RTCA/DO-235B[END_REF], it depends on:

-The incoming Doppler frequency -The L-band central frequency -The code frequency [START_REF] Rtca | Assessment of Radio Frequency Interference Relevant to the GNSS L1 Frequency Band RTCA/DO-235B[END_REF] Where is the chip duration. is a priori known because it refers to the local spreading code whereas is unknown. The difference in code frequency leads to a change in the spreading code period as it can be seen in Figure 4 where 3 periods of a 4-chips spreading code are represented:

-A positive code Doppler frequency causes the spreading code duration to shrink ( ) -A negative Doppler shift causes the spreading code duration to expand ( )

Received code

Local code 1st period 2nd period 3rd period 1st period 2nd period 3rd period Mathematically, it can be modeled as two ways. The first one which is used to derive the mathematical model uses the rectangular function:

∑ ( ) ( ) (13)
Where -

( ) is the spreading code chip value - ( ) {
is the rectangular function The other way, which is used in [START_REF] Jiao | High Sensitivity GPS Acquisition Algorithm Based on Code Doppler Compensation[END_REF], [START_REF] O'driscoll | Performance analysis of the parallel acquisition of weak GPS signals[END_REF] and [START_REF] Psiaki | Block Acquisition of Weak GPS Signals in a Software Receiver[END_REF], consists in introducing the fractional perturbation caused by Doppler frequency and noted :

(( )( )) (14)
As it will be shown, the main effects of the code Doppler are degradation in the magnitude of the autocorrelation function and a shift of the true code delay. If the shift between the received spreading code and the local spreading code (without Doppler) exceeds 1 chip then the correlation process (summation) no longer makes sense because the power of the signal cannot be accumulated. This obvious limit of 1 chip represents the loss of matching between the two sequences. The maximum total integration time before the slip of 1 chip depending on the incoming Doppler frequency and GNSS signal, is given in Table III. For example, let us assume an incoming Doppler frequency of 10 kHz, for the GPS L1 C/A signal, this leads to a code Doppler of Hz. Then, a chip offset occurs after 154 ms; equivalent to 154 spreading code periods. However, for the L5 signals, a chip offset appears after 12 ms (12 spreading code periods). One can celarly understand that in 12 ms, it seems to be difficult to acquire weak signals because the total integration time is too short. Even for low Doppler frequencies (the incoming Doppler is lower for L5 signals rather than L1 signals), the total integration time is short. The slip in chip in one spreading code period is presented in Figure 5. This first result highlights the limits on the total integration time when several spreading code periods are required to acquire signals without Doppler compensation algorithms. Moreover, it is worth noting that the code Doppler impact is visible on the first spreading code due to the use of high sampling frequency. For example, for a sampling frequency of MHz let us simulate the first Galileo E1 OS spreading code (with a CBOC modulation) with no Doppler on the one hand and on the other hand with an incoming Doppler frequency of 10 kHz. Figure 6 presents the number of sample slips per chip for each chip in the first spreading code. For the last chips, around one half of the samples are slipped. Over the first spreading code, there are around 14000 differences over the 81920 samples (17%) even if the chip offset after one spreading code period is lesser than 1 sample ( chip corresponds to 0.52 sample). This section presented the problem that is addressed in this paper: the impact of the code Doppler on the acquisition performance. In the next sections, a simple mathematical model is given with some results.

II. MATHEMATICAL ANALYSIS

Mathematical model of distorted autocorrelation

To simplify and be able to develop the mathematical model of autocorrelation function, some assumptions need to be taken:

-There is no BOC subcarrier, then the values taken by the spreading code sequence are or -There is no incoming code delay, this means that ̂ and there is no a chip offset at the beginning of the autocorrelation process due to code Doppler We precise that if then [ ] ( ) and ∫ [ ] ( )

In this section, we focus on the autocorrelation function when the incoming spreading code is affected by the Doppler. Then, it can be written as: In the same way, the local spreading code is:

( ) ∑ ( ) ( ) (15)
( ) ∑ ( ) ( ) (16)
Then, on the integration interval [ ( ) ], the autocorrelation function affected by the code Doppler ̃( ) is given by [START_REF] Ziedan | GNSS receivers for weak signals[END_REF] which is obtained by dividing the total integration interval in small integration intervals where the local spreading code chip is constant. Then, as it is done in [START_REF] Foucras | A Novel Computationally Efficient Galileo E1 OS Acquisition Method for GNSS Software Receiver[END_REF], the product of the rectangular functions can be seen as an indicator function

[ (( ) ( ) ) (
)] . Applying [START_REF] Foucras | A Novel Computationally Efficient Galileo E1 OS Acquisition Method for GNSS Software Receiver[END_REF] to ( 19) gives [START_REF] Krasner | GPS Receiver and Method for Processing GPS Signals[END_REF] which can be simplified by [START_REF] Akopian | A Fast Satellite Acquisition Method[END_REF]. The double sum (on and on ) can be seen as the sum on of three terms. To be more precise, as it is depicted in Figure 7, the chip of the local sequence can be multiplied by the chips or of the received spreading code.

( ) ( )

[ ] ] corresponds to potential degradations (red zones) where the chips in the two spreading code sequences are not in the same position and their product result can be or with a probability of one half for each one The result can be easily verified by simulations but it seems to be difficult to extend in a general case (for CBOC-modulated signals, with initial delay…)

Resulting autocorrelation function approximation

Matlab simulations were run to verify that the resulting autocorrelation function can be approximated by the autocorrelation function for different signals, here only the GPS L1 C/A approximation is shown. On the one hand, the resulting autocorrelation function is evaluated at each spreading code period:

̃ ( ) ∑ ( ) ( ) (17)
Where:

corresponds to the integration period (for , we assume that the first chips of the received and local sequences are correctly synchronized) -is the number of samples in one local spreading code period -̃ ( ) represents the affected autocorrelation function in 0 after ( ) local spreading code periods On the other hand, the code delay is evaluated and the autocorrelation function taken in this point is calculated:

( ) | | (18) ̃( ) ∫ ( ) ( ) ( ) ∫ ∑ ( ) ( ) ∑ ( ) ( ) ( ) ∑ ( ) ∑ ( ) ∫ ∑ ( ) ( ) ( ) (19)
( ) ( ) { { ( ) ( ) [ (( ) ( ) ) ( )] ( ) (20)
̃( ) ∑ ( ) ∑ ( ) ∫ [ (( ) ( ) ) ( )] ( ) (21)
̃( ) ∑( ( ) ( )[ ( ) ] ( ) ( )[ ]) ( ( ) ( ) ( ) ∑ ( ) ( )) (22)
̃( ) ( ( ) ( ) ∑ [ ( ) ( ) [ ] ( ) ( ) [ ] ]) (23)

III. SIMULATION RESULTS

In this section, some simulation results are presented to illustrate the code Doppler effect on the acquisition performance. The results based on the normalized squared autocorrelation function differ with the summation technique (coherent or non-coherent). Firstly, the noncoherent results are observed.

Non-coherent summations

Based on the total integration time before the slip of 1 chip (Table IV), the autocorrelation function affected by the code Doppler is presented in Figure 9. The black curve is the squared autocorrelation function without Doppler, the reference one and the worst case is for a Doppler of 10 kHz in green.

For the BPSK-modulated signals (GPS L1 C/A and GPS L5), the autocorrelation function shape becomes rounded and offset compared to the reference triangular. The amplitude of the maximum value is reduced and the peak is shifted to the right for a negative Doppler. For the BOC-modulated signals (GPS L1C and Galileo E1 OS), the secondary peaks which characterize the autocorrelation function tend to disappear leading to a flat curve.

For the same maximum total integration time (154 ms for Galileo E1 OS, GPS L1 C/A and L1C), the impact on the resulting autocorrelation function shape is different mainly due to the modulation. Table IV provides the characteristics of the main peak of the normalized autocorrelation for an incoming Doppler of 10 kHz. The peak slips of 1 bin in the acquisition grid for BPSKmodulated signal and because the peak is not evident for BOC-modulated signals, there is a change of 4 or more bins.

For the GPS L1 signals, the probability of detection is equivalent even if there are more losses on the amplitude for the GPS L1C signals but it is compensated by longer coherent accumulation (on the spreading code period). For Galileo E1 OS, which has the same losses but a shorter spreading code period (4 ms instead of 10 ms), the probability of detection is close to 0.2 for high Doppler frequency (10 kHz). For GPS L5, the probability of detection is very low, this is explained by the very short integration time.

Coherent summations

Now, let us look at the coherent summation. The same incoming and local signals are processed but in a coherent way; this means that only one correlation on the total integration time is achieved.

As it was done in Figure 9, Figure 12 presents the autocorrelation function for the maximum integration time before a 1-chip offset for several incoming Doppler frequencies. For BPSK-modulated signals, the shape of the autocorrelation function seems to be similar as for noncoherent acquisition technique. However, for BOCmodulated signals, the flattening effect is very important leading to a quasi-null function for the worst case (in green). Indeed, Table V gives the amplitude of the peak of the resulting autocorrelation function. For the GPS L1C and Galileo E1 OS signals, the amplitude peak is around 0.035 which is lower than 0.2 for non-coherent acquisition technique (Table IV). The losses on the peak of the autocorrelation function for different coherent summations are presented in Figure 13 for the maximum total integration time. For the GPS L1 C/A and GPS L5 signals, the losses on the acquisition detector due to the code Doppler on the autocorrelation function (2.5 dB for 10 kHz) are similar that ones for noncoherent acquisition technique. As expected, the losses are important for the BOC-modulated signals (for 10 kHz more than 12 dB). Knowing the losses in dB on the main peak of the resulting autocorrelation function, the probability of detection can be evaluated depending on the Doppler frequency; results are shown in Figure 14.

( ) ( ( )) (( ̃ ) ) (26)
Figure 14: Probability of detection for weak signal and a coherent integration on the maximum allowed signal length

When not considering data sequence presence, it is preferable to coherently integrate rather than noncoherently sum for the same total integration time [START_REF] Bastide | Analysis of the Feasibility and Interests of Galileo E5a/E5b and GPS L5 for Use with Civil Aviation[END_REF]. Then, the probability of detection in presence of code Doppler is still better for coherent summation than noncoherent summation. Let us note for the GPS L1 C/A signal, the probability of detection seems to be not affected by code Doppler due to the long coherent integration time.

Discussion

To complete the coherent accumulation results, we present in Figure 15 the probability of detection in function of the received sensitivity for Galileo E1 OS on 80 ms (coherently accumulation). It is more degraded for low C/N 0 than high C/N 0 and it is really accentuated for high incoming Doppler frequency. Then, we will interest on the impact of the code Doppler for low C/N 0 . Without Doppler, it is interesting to assess the minimum integration time so that the receiver can acquire a signal with a target C/N0 and a target probability of detection. This is shown in Figure 16 in the case of a signal with a C/N0 of 27 dB•Hz -on the data component for composite signals (GPS L1C, Galileo E1 OS and GPS L5) or total signal for GPS L1 C/A. To complete Figure 16, Table VI provides the required time to acquire a signal at 27 dB•Hz with a high probability (90%) and without Doppler. This corresponds to the intersection of the red and blue curves with the black horizontal line materializing a probability of detection of 0.9. To study the impact of code Doppler on the acquisition performance; in the same conditions as previously presented, the probability of detection is computed considering Doppler and for the integration time described in Table VI. Let us remark that for GPS L5, the required total integration time to acquire at 27 dB•Hz exceeds the maximum total integration time before the slip of 1 chip. That means that if we want to acquire GPS L5, we suffer from the slip of 1 or 2 chips due to code Doppler.

When there is no Doppler, the probability of detection is around 0.9 but when an incoming Doppler is considered leading to a code Doppler, the probability of detection falls, particularly for non-coherent acquisition technique. The difference in acquisition performance between both acquisition accumulation techniques is mainly due to the integration time which is longer for non-coherent acquisition technique implying a higher offset between the incoming and the local spreading code sequence. The other argument is the better performance in terms of probability of detection for coherent accumulation. Of course, the higher the incoming Doppler frequency is, higher the degradations on the probability of detection are.

IV. CODE DOPPLER COMPENSATION METHODS

In the previous sections, the impact of the code Doppler on the acquisition performance has been demonstrated highlighting:

-The higher the incoming Doppler frequency is, the higher the losses are -The weaker the signal is, the higher the losses are So, now, knowing that, a state-of-the-art of the code Doppler compensation acquisition method is proposed.

Code Doppler compensation on the spreading code sequences 1) Local code generation based on the incoming Doppler

The first code compensation method is the simplest one and is proposed in [START_REF] Ziedan | GNSS receivers for weak signals[END_REF] for the Double Block Zero Padding (DBZP) which algorithm steps are presented in [START_REF] Foucras | A Novel Computationally Efficient Galileo E1 OS Acquisition Method for GNSS Software Receiver[END_REF]. It consists in generating for each possible incoming Doppler frequency a local version of the spreading code which code rate depends on the estimate of the Doppler frequency. Then, there are as many replica code versions as there are cells in the acquisition grid in the frequency search space. The th local code version is used for the calculation of the coherent integration of the th Doppler frequency estimate. Only the coherent result which is maximum (and corresponds to the right estimation of the incoming code Doppler) is kept. This computational analysis of this method brings to the conclusion that it is very expensive in terms of number of computations and storage space. 2) Interpolation This code Doppler compensation method is based on a modification of the local code and was developed in [START_REF] Psiaki | Block Acquisition of Weak GPS Signals in a Software Receiver[END_REF]. In this method, the local code is generated in much the same way as for classical acquisition method, in that the code frequency is varied with the Doppler frequency estimate. Then, for every Doppler frequency estimate, an interpolation of both the received signal and the local one is performed to ensure that there is exactly one full spreading code period of each present in the coherent integration interval. This method permits to provide code Doppler compensation within and between coherent integration intervals but requires large memory due to the storage of the sampled local codes.

3) Time shift theorem implementation

The time shift-theorem can be applied to achieve Doppler code compensation. Indeed, the time-shift frequency theorem is:

( ( ))( ) ( ) ( ) (27
) Where is the Fourier transform operator is a function is the Fourier transform of Then, for FFT-based acquisition method, the code delay induced by code Doppler can be compensated by multiplying in the frequency domain the FFT of the local spreading code by the complex exponential ( ̂ ). This technique was patented by Krasner [START_REF] Krasner | GPS Receiver and Method for Processing GPS Signals[END_REF] and also applied in [START_REF] Akopian | A Fast Satellite Acquisition Method[END_REF] which proposes an illustration given in Figure 19. The advantage of this method is that if the Fourier transform of the spreading code sequence has been built and stored in the memory, then the Fourier transform of the receded (or extended) spreading code can be transformed to the frequency domain in a simple way. Then the correct spreading code can be produced quickly. It can be seen as read more or less quick the stocked local spreading codes.

Code Doppler compensation on the correlator outputs

This method which is based on code Doppler picking points compensation is presented in [START_REF] Jiao | High Sensitivity GPS Acquisition Algorithm Based on Code Doppler Compensation[END_REF] and consists in picking points in coherent integration results to compensate phase sliding. It is illustrated by Figure 20. The estimation of the code Doppler (depending on the incoming Doppler frequency) permits to calculate the code phase shift moments to compensate phase sliding. The coherent summation results should compensate the chip offset to guarantee all of the correlation peaks appear in unique code phase. This method provides improvements in acquisition performance by still calculating the correct code phase even in more extreme environment. Based on fast Fourier transform, this method does not require important resource consumptions and storage.

CONCLUSION

In this paper, a study on the code Doppler impact on the acquisition performance of new GNSS signals was presented. To conclude it, we will summarize the main results. A first theoretic result gives the maximum total integration time before the slip between local and received sequences of 1 chip. To keep in mind is that for GPS L1 C/A, Galileo E1 OS and GPS L1C which share the same L-band, it is 154 ms for an incoming Doppler frequency of 10 kHz. Then, for GPS and Galileo signals in the L5 band and GPS L1C, one cannot sum more than 15 spreading codes for a 10 kHz incoming Doppler frequency if he does not want to suffer from 1-chip slip. The autocorrelation function affected by code Doppler is significantly degraded: rounded and offset. Due to BOC modulation, Galileo E1 OS and GPS L1C present more degradation on the peak amplitude than BPSK-modulated signals. The phenomenon is particularly visible for coherent accumulation rather than non-coherent accumulation. However, the coherent accumulation technique presents less acquisition degradations in terms of probability of detection. Finally, the degradations on the acquisition performance are considerable when weak signals should be acquired due to the long required integration time. Once again, the non-coherent accumulation technique suffers more than coherent accumulation when dealing with the probability of detection.

As it has been explained, for the new challenges (low C/N 0 environments, new GNSS signal structures, new applications…), one can hardly assume that there is no code Doppler. One possibility is to consider code Doppler compensation methods which are developed to limit acquisition performance losses when considering code Doppler. Maybe, some new code Doppler compensation methods can be explored to specific uses in terms of signal structures and sensitivity. Anyway, the residual losses when using these code Doppler compensation methods should be explained to permit to understand how these methods can be applied for the acquisition process.

Figure 2 :

 2 Figure 2: Acquisition process scheme

Figure 3 :

 3 Figure 3: Autocorrelation functions

Figure 4 :

 4 Figure 4: Code Doppler effect on the spreading code period

Figure 5 :

 5 Figure 5: Chip offset in one spreading code period

Figure 6 :

 6 Figure 6: Differences number of samples per chip

  -f s = 20.48 MHz -f d = 10 kHz -T I = 4 ms

Figure 7 :

 7 Figure 7: Correlation process with different code frequencies Then, by noting | | the absolute difference between the duration of the local and the received spreading code chips, the general expression of the resulting autocorrelation function is given by (23). It is composed of two terms: -( )

Figure 8

 8 Figure 8 depicts the verification result for the GPS L1 C/A spreading code.

Figure 8 :

 8 Figure 8: Approximation of the affected autocorrelation function For all of the incoming Doppler frequencies, the two curves are superposed which means that ̃ ( ) ( )(24)

Figure 9 :

 9 Figure 9: Affected autocorrelation function shape (Non-coherent)

Figure 10 provides

 10 Figure10provides the amplitude losses taken in the main peak versus the incoming Doppler frequency.

Figure 10 :Figure 11 :

 1011 Figure 10: Losses on the autocorrelation function due to code Doppler (Non-coherent)Due to the flattening of the autocorrelation function curve for BOC-modulated signals, losses in the main peak are higher (7 dB in the worst incoming Doppler frequency) than BPSK-modulated signals(2.4 dB). This means that for weak signals, the acquisition becomes very hard as it can be seen in Figure11with the probability of detection whose mathematical expression is given by (25) at 27 dB-Hz.

Figure 12 :

 12 Figure 12: Affected autocorrelation function shape (Coherent)

Figure 13 :

 13 Figure 13: Losses on the autocorrelation function due to code Doppler (Coherent)

Figure 15 :

 15 Figure 15: Probability of detection for different Doppler versus the sensitivity

Figure 16 :

 16 Figure 16: Required acquisition time to reach a high probability of detection without considering Doppler for weak signals (27 dB•Hz)

Figure 19 :

 19 Figure 19: Krasner' technique[START_REF] Akopian | A Fast Satellite Acquisition Method[END_REF]

Figure 20 :

 20 Figure 20: Code Doppler picking points compensation method [3]

Table I :

 I Data component signal features

				Central	Code
	Signal		Modulation	frequency	frequency
				(MHz)	(MHz)
	GPS L1 C/A		BPSK	1575.42	1.023
	GPS L1C-D	BOC	1575.42	1.023
	GPS L5-D	BPSK	1176.45	10.23
	Galileo E1 OS-B		CBOC (6,1,1/11)	1575.42	1.023
	Galileo	a-I	BPSK	1176.45	10.23
	E5	b-I	BPSK	1207.14	10.23

Table II :

 II PRN codes properties on the data component

			PRN	PRN
	Signal		length	period		Subcarrier
			(chips)	(ms)
	GPS L1 C/A	1023	1		None
	GPS L1C-D	10230	10		( ) ( )
	GPS L5-D	10230	1		None
	Galileo E1 OS B	4092	4	√	( ) ( ) √	( ) ( )
	Galileo	a-I	10230	1		None
	E5	b-I	10230	1		None

Table III :

 III Offset of 1 chip depending on the incoming Doppler frequency

			Incoming Doppler frequency
			Offset of 1 chip (ms)
		Offset of 1 chip (number of spreading code periods)
	Signal	2 kHz	4 kHz 6 kHz 8 kHz 10 kHz
	GPS L1 C/A	770 770	385 385	257 257	193 193	154 154
	GPS L1C-D	770 77	385 38.5	257 25.7	193 19.3	154 15.4
	GPS L5-D	58 58	29 29	20 20	15 15	12 12
	Galileo	770	385	257	193	154
	E1 OS B	192.5	96.25	64.25	48.25	38.5
	Galileo	58	29	20	15	12
	E5a-I	58	29	20	15	12

Table IV :

 IV Main autocorrelation function peak for a 10 kHz incoming Doppler

		GPS L1 C/A	GPS L1C	Galileo E1 OS	GPS L5
		154	15	38	12
	Maximum amplitude	0.58	0.20	0.20	0.59
	Argmax (chip)	0.54	0.67	0.75	0.56

Table V :

 V Main autocorrelation function peak for a 10 kHz incoming Doppler

		GPS L1 C/A	GPS L1C	Galileo E1 OS	GPS L5
	(ms)	154	154	154	12
	Maximum amplitude	0.56	0.06	0.08	0.55
	Argmax (chip)	0.5	0.46	0.46	0.54

Table VI :

 VI Required integration time for a probability of detection greater than 0.9 and a C/N0 of 27 dB•Hz (without Doppler)

		GPS L1 C/A	GPS L1C	Galileo E1 OS
	Non-	106 sum.	4 sum.	12 sum.
	coherent	106 ms	40 ms	48 ms
	Coherent	24 ms	30 ms	24 ms

Her work is funded by ABBIA GNSS Technologies, in Toulouse, France working on the development of a GPS/Galileo software receiver. is the head of the Signal Processing and Navigation (SIGNAV) research

As previously explained in section I, the serial search acquisition is based on an acquisition grid. Herein, it is assumed that the acquisition is achieved which means that the bin corresponding to the right estimation of the incoming code delay and Doppler frequency is found. The point of interest is the incoming Doppler frequency for which the chip offset due to code Doppler leads to a change in the grid acquisition bin. In Figure 18, the position of the main peak of the resulting autocorrelation function is given for each integration on and for several incoming Doppler frequencies. This permits to understand the relationship between the Doppler frequency and the number of successive summation. For instance, for the GPS L5 signal, the 3 first autocorrelations functions have a code shift lower than 0.25 (threshold given in Figure 3) for an incoming Doppler of 10 kHz; for 4 kHz, the 8 first ones. To complete this result, we are interested on the impact on the acquisition performance versus the acquisition accumulation technique (coherent or non-coherent). Table VII provides maximum Doppler frequency for which there is no acquisition grid cell change and this for all of the signals, acquisition techniques and integration. Of course, an average process occurs leading to higher Doppler for coherent or non-coherent accumulation on than the integration. This means that for an incoming Doppler of 5 kHz (or less), the peak of the resulting non-coherent autocorrelation function is for a chip offset lower than 0.25 when the maximum integration time is processed (154 ms). It is clear that for BOC-modulated signals, the incoming Doppler frequency limit is low because of the strong code Doppler impact. We concluded that the coherent accumulation provides better acquisition performances (in terms of probability of detection) even if the losses on the autocorrelation peak amplitude are higher.


Nicolas Delfosse
email: nicolas.delfosse@usherbrooke.ca

Jean-Pierre Tillich
email: jean-pierre.tillich@inria.fr

A decoding algorithm for CSS codes using the X/Z correlations

We propose a simple decoding algorithm for CSS codes taking into account the correlations between the X part and the Z part of the error. Applying this idea to surface codes, we derive an improved version of the perfect matching decoding algorithm which uses these X/Z correlations.

Introduction

Low Density Parity-Check (LDPC) codes are linear codes defined by low weight parity-check equations. It is one of the most satisfying construction of error-correcting codes since they are both capacity approaching and endowed with an efficient decoding algorithm. It is therefore natural to investigate their quantum generalization.

Besides their use for quantum communication, quantum LDPC codes could play a central role in quantum computing. A striking difference between classical and quantum information is the fact that every manipulation of quantum bits (qubits) is very noisy. Quantum gates must therefore be implemented in a fault-tolerant way. This is realized by applying operations on qubits encoded by a quantum error-correcting code. These qubits can then be regularly corrected. Some recent work of Gottesman [START_REF] Gottesman | What is the overhead required for fault-tolerant quantum computation[END_REF] has shown that quantum LDPC codes are well-suited for fault-tolerance. These codes, which are defined by low weight constraints on qubits, naturally limit the propagation of errors.

The first difficulty in the generalization of LDPC codes is that most of the constructions have bounded distance (see [START_REF] Tillich | Quantum LDPC codes with positive rate and minimum distance proportional to n 1/2[END_REF] and references therein). The rare families of quantum LDPC codes equipped with a growing distance are derived from Kitaev's construction. Kitaev's toric code is defined by local interactions between qubits placed on a square tiling of the torus. Similar constructions were proposed, based on tilings of surfaces [START_REF] Bombin | Homological error correction: Classical and quantum codes[END_REF][START_REF] Zémor | On Cayley graphs, surface codes, and the limits of homological coding for quantum error correction[END_REF], 3-colored tilings [START_REF] Bombin | Topological quantum distillation[END_REF][START_REF] Delfosse | Tradeoffs for reliable quantum information storage in surface codes and color codes[END_REF], Cayley graphs [START_REF] Couvreur | A construction of quantum LDPC codes from Cayley graphs[END_REF] or other geometrical objects [START_REF] Tillich | Quantum LDPC codes with positive rate and minimum distance proportional to n 1/2[END_REF][START_REF] Kovalev | Quantum Kronecker sum-product low-density paritycheck codes with finite rate[END_REF][START_REF] Freedman | Quantum systems on non-k-hyperfinite complexes: A generalization of classical statistical mechanics on expander graphs[END_REF][START_REF] Audoux | An application of Khovanov homology to quantum codes[END_REF].

The belief propagation decoding algorithm is an essential ingredient of the success of LDPC codes. Unfortunately, it is much less effective in the quantum setting due to two facts (i) the unavoidable presence of 4-cyles in the Tanner graph [START_REF] Camara | Constructions and performance of classes of quantum LDPC codes[END_REF] and (ii) the low weight generators can be considered as low-weight errors which are not detected by the belief propagation decoder but which are harmful for its convergence [START_REF] Poulin | On the iterative decoding of sparse quantum codes[END_REF]. To circumvent this obstacle, some techniques originating from classical coding theory were imported in quantum information recently [START_REF] Kasai | Sakaniwa Quantum error correction beyond the bounded distance decoding limit[END_REF][START_REF] Andriyanova | Spatially coupled quantum ldpc codes[END_REF]. Another direction to avoid the 4-cycles, is to consider the error, which is a quaternary vector, as a pair of binary vectors. These two binary vectors can then be decoded separately. The main problem of this point of view is that it does not consider the correlations between the two binary components of the error. In this work, we present a simple and general strategy to take into account these correlations. To illustrate this idea, we focus on surface codes equipped with the perfect matching decoding algorithm. This algorithm is usually unable to consider the correlations. Applying our method to a family of surface codes constructed from triangular tilings of a torus, we observe a clear improvement of the performance of the decoding algorithm. The depolarizing error threshold of these triangular codes is approximately 13.3% while it is close to 9.9% without considering the correlations.

This article is organised as follows. The definition of surface codes and the geometrical description of errors and syndrome over these codes are recalled in Section 2. Section 3 explains how decoding can be performed by using the aforementioned correlations. Section 4 is devoted to the description of the perfect matching decoding algorithm and its improvement to take into account the correlations.

Definitions and basic properties

Error model. We deal here with the depolarizing channel model which is one of the most natural quantum error model and the quantum analog of the binary symmetric channel. Over the depolarizing channel of probability p, each qubit is subjected, independently, to an error X, Y or Z with probability p/3 or is left unchanged with probability 1p where X, Y and Z denote the usual Pauli matrices. An error E over n qubits is therefore a tensor product ⊗ n i=1 E i where E i ∈ {I, X, Y, Z}. Errors are considered up to the phase {±1, ±i}, since quantum states are defined up to a phase.

Stabilizer and CSS codes. A quantum code is a subspace of dimension 2 k of (C 2 ) ⊗n . This code encodes k qubits into n qubits. A very useful way of constructing such codes is through the stabilizer code construction [START_REF] Gottesman | Stabilizer Codes and Quantum Error Correction[END_REF] where the code is described by the set of fixed states of a family of commuting Pauli operators {S 1 , . . . , S r }. In other words, the S i 's are generators of the stabilizer group of the quantum code. A particular case of this construction is the CSS construction due to Calderbank, Shor [START_REF] Calderbank | Good quantum error-correcting codes exist[END_REF] and Steane [START_REF] Steane | Multiple-particle interference and quantum error correction[END_REF]. It consists in choosing some of these Pauli operators in {I, X} ⊗n and the rest of them in {I, Z} ⊗n . This brings several benefits, first it simplifies the commutation relations and helps in constructing such codes and second decoding of such codes can be achieved by decoding two binary codes as will be explained in the next paragraph.

Syndrome measurement and decoding of CSS codes. For a stabilizer code with stabilizer generators {S 1 , . . . , S r } subjected to a Pauli error E, it is possible to perform a measurement which reveals the vector s(E) def = (E S i ) 1≤i≤r where E S i is equal to 0 if E commutes with S i and is equal to 1 otherwise. In the case of a CSS code, the syndrome splits into two parts, one corresponding to the commutation with the generators belonging to {I, X} ⊗n and the other one corresponding to the commutation relations with the generators in {I, Z} ⊗n . Moreover, if we decompose the error E as E = E X E Z where E X ∈ {I, X} ⊗n and E Z ∈ {I, Z} ⊗n and if we let S 1 , . . . , S r X be the generators which are in {I, X} ⊗n and S r X +1 , . . . , S r be the generators which are in {I, Z} ⊗n , then the syndrome part s X which corresponds to the generators in {I, X} ⊗n verifies s

X def = (E S i ) 1≤i≤r X = (E Z S i ) 1≤i≤r X whereas the syndrome part s Z which corresponds to the generators in {I, Z} ⊗n verifies s Z def = (E S i ) r X +1≤i≤r = (E X S i ) r X +1≤i≤r .
Notice that if we bring in the binary matrices H X and H Z whose rows are formed for H X , respectively H Z , by the generating elements belonging to {I, X} ⊗n , respectively {I, Z} ⊗n (and replacing I by 0 and X by 1, respectively replacing I with 0 and Z with 1), then s X is nothing but the syndrome H X e T Z of the binary error e Z (obtained from E Z by replacing I by 0 and Z by 1), whereas s Z is nothing but the syndrome H X e T Z of the binary error e Z (obtained from E Z by replacing I by 0 and Z by 1). In other words decoding a CSS code amounts to decode two binary codes. This is how decoding a CSS code is usually performed. We call this decoding technique the standard CSS decoder.

Tiling of a surface. A surface code is a CSS code associated with a tiling of surface. Let us recall the definition of a tiling of surface. A tiling of surface is defined to be a cellular embedding of a graph G = (V, E) in a 2-manifold, that is, a surface. Without loss of generality, we can assume that the surface is smooth. We assume that the graph G contains neither loops nor multiple edges. This embedding defines a set of faces F. Each face is described by the set of edges on its boundary. This tiling of surface is denoted

G = (V, E, F). The dual graph of G is the graph G * = (V * , E * ) of vertex set V * = F such
that two vertices are linked by an edge if and only if the two corresponding faces of G share an edge. There is a clear bijection between the edges of G and the edges of G * . This graph G * is endowed with a structure of tiling of surface and its faces correspond to the vertices of G: these faces are induced by the set of edges of G incident to a vertex v ∈ V.

Surface Codes. Surface codes are a special case of CSS codes. They have been introduced by Kitaev [START_REF] Kitaev | Fault-tolerant quantum computation by anyons[END_REF]. Assume that qubits are placed on the edges of a tiling of surface G = (V, E, F). The space of the system is e∈E H e , with H e = C 2 for every edge e ∈ E. The Pauli operators acting on this space are the tensor products ⊗ e∈E P e such that P e ∈ {I, X, Y, Z}. For every edge i ∈ E, denote by X i = ⊗ e P e the Pauli operator which is the identity on every edge except on edge i, where P i = X. The operators Z i are defined similarly for all i ∈ E. The site operators X v and the plaquette operators Z f are the Pauli operators defined by

X v = v∈e X e and Z f = e∈f Z e ,
for every vertex v ∈ V and for every face f ∈ F. Then, the surface code associated with the tiling of surface G is the CSS code fixed by the site operators and the plaquette operators. The commutation between these operators follows from the structure of the tiling. Note that H X is in this case the incidence matrix of the graph G and H Z the incidence matrix of its dual G * . The site operators and the plaquette operators of Kitaev's toric codes are represented in Figure 1 (a).

Syndrome of a surface code. In the case of a surface code, the syndrome has a graphical interpretation that we recall now. Consider the surface code associated with a tiling of surface G = (V, E, F). Assume that an error E Z acts on a path γ ⊂ E of G. In other words, we have E Z = e∈γ Z e . Then, the syndrome s X = (E Z X v ) v∈V of this error is indexed by the vertices of the graph and is non-trivial if and only if the vertex v is an end-point of the path γ. This follows from the fact that E Z commutes with all the operators X v , except the two operators centered on the end-points of γ. More generally, the support of the error E Z can be decomposed as a union of disjoint paths and its syndrome indicates the end-points of the (a) To obtain an analogous description of the error E X and its syndrome, replace the graph by its dual. Indeed, this transformation exchanges the roles of X and Z in the definition of the code.

X X X X Z Z Z Z (b) X X X XX X Z Z Z
The following well known lemma summarizes the graphical description of the error.

Lemma 1. Let G = (V, E, F) be a tiling of surface and let G * = (V * , E * , F * ) be its dual. An error acting on the surface code associated with G corresponds to a pair (E X , E Z ) such that E X ⊂ E * and E Z ⊂ E, and its syndrome is the pair

(s X , s Z ) such that s X ⊂ V is the set ∂(E Z ) of end-points of E Z and s Z ⊂ V * is the set ∂(E X ) of end-points of E X .
3 Decoding by using correlations between errors in X and Z Virtually all decoders of CSS codes try to recover the E X and E Z part of the error independently by decoding two binary codes as explained in Section 2. There is however some correlation between the X part of the error and the Z part as shown by the following conditional probabilities computed for a single error E = E X E Z generated by the depolarizing channel of depolarizing probability p:

P(E Z = I|E X = X) = 1/2 (1)
P(E Z = Z|E X = X) = 1/2 (2)
whereas

P(E Z = I|E X = I) = 1 -p 1 -2p/3
(3)

P(E Z = Z|E X = I) = p/3 1 -2p/3 • (4)
When E X = X, we recognize an erasure channel and in the second case this corresponds to a binary symmetric channel of probability p def = p/3 1-2p/3 . This can be exploited by the following strategy for decoding. First, decode the X component of the error. Then, erase the coefficients of E Z corresponding to the X errors. Finally, decode the Z component of the error, which is subjected to a combination of errors and erasures. We call such a decoder a CSS decoder using X/Z correlations.

It is insightful to calculate the capacity of the two classical channels that both decoders face. The X decoder has to work for a binary symmetric channel of crossover probability p def = 2p/3 whereas the Z decoder has to work for a binary error and erasure channel, where a bit gets erased with erased with probability p and corrupted with probability (1p )p . The capacity of the first channel is equal to 1h(p ) whereas the capacity of the second channel is equal to (1p )(1h(p )). It can be readily observed that the second capacity is always larger than the first one. This suggests two things (i) if the two binary codes have the same rate (that is if the number of X generators is the same as the number of Z generators), then we may expect that the second decoder behaves much better than the first decoder and that the probability of the whole decoding is essentially the probability that the first decoder fails instead of being essentially twice this probability as is usually the case for the standard CSS decoder described in the previous section.

(ii) In order to fully use this decoder, the best strategy for choosing the CSS code (without using the possible degeneracy of the code) is to choose an asymmetric CSS code where the number of Z generators of the CSS code is chosen such that the binary code associated to H Z has rate slightly below 1-h(p ) whereas the X generators are chosen such that the rate of the binary code associated to H X has rate slightly below (1p )(1h(p )). This strategy of decoding is able to reach the hashing bound, which is equal to 1 + p log p 3 + (1p) log(1p) for a depolarizing channel as explained by the following theorem.

Theorem 2. For any > 0, there exists a family of CSS codes of quantum rate ≤ 1 + p log p 3 + (1p) log(1p)for which the error probability after decoding with the CSS decoder using X/Z correlations goes to 0 as the length goes to infinity. This theorem is proved by random coding techniques and will be given in the full version of this paper. Notice that the hashing bound is significantly bigger than 1 -2h(p ) which is the biggest quantum rate that random CSS codes may have in order to be decoded succesfully by the standard CSS decoder.

Improvement of the Perfect Matching Decoding

In this section, we recall the perfect matching decoding algorithm [START_REF] Dennis | Topological quantum memory[END_REF] for surface codes, we discuss about its two main weakness and improve its performance by using the strategy outlined in Section 3.

The Perfect Matching Decoding Algorithm

We consider that a surface code is subjected to a random error (E X , E Z ) generated by a depolarizing channel of probability p. The goal of this algorithm is to determine a most likely error E Z which corresponds here to an error of minimum weight (since in general we are in a situation where p ≤ 1/2) which has syndrome s X . The component E X is decoded similarly in the dual graph.

To determine an error E Z ⊂ E of minimum weight, given its end-points s X = ∂(E Z ) ⊂ V, we are looking for a set of paths whose end-points are exactly s X and whose size is minimum.

Algorithm 1 computes such a set using Edmonds' minimum weight perfect matching algorithm [START_REF] Edmonds | Maximum matching and a polyhedron with 0-1 vertices[END_REF][START_REF] Edmonds | Path, trees, and flowers[END_REF][START_REF] Kolmogorov | Blossom V: a new implementation of a minimum cost perfect matching algorithm[END_REF]. This decoding algorithm first computes the distance graph associated with a syndrome s ⊂ V. It is the weighted complete graph K(s), with vertex set s = {s 1 , s 2 , . . . , s m }, such that the weight of the edge {s i , s j } is the distance d(s i , s j ) in G. The second step of the algorithm is the determination of a minimum weight perfect matching M in K(s). Recall that a perfect matching in a graph H is a set of edges of H meeting all the vertices of H exactly once. With each edge {v i , v j } ∈ M , we associate a geodesic of G joining the vertices v i and v j . Denote by G(v i , v j ) this geodesic. The algorithm returns the symmetric difference of all the geodesics corresponding to the edges of M . It is the support of a most likely error of syndrome s.

Algorithm 1 Perfect Matching Decoding

Require: A graph G = (V, E), a subset s ⊂ V of odd size. Ensure: A subset x ⊂ E of minimum size with end-points ∂(x) = s.

1: Construct the distance graph K(s) associated with s.

2: Determine a minimum weight perfect matching M 3: return the symmetric difference of all the geodesics G(v i , v j ) for {v i , v j } ∈ M .

Degeneracy and Correlations

We now discuss of two cases of failure of the perfect matching decoding algorithm and their effect on the performance. First, by definition, surface codes are fixed by the plaquette operators of the tiling. Thus two errors which differ in a sum of plaquettes (or faces) have exactly the same effect on the quantum code. We should thus look for the most likely error coset up the sums of faces instead of the most likely error. This phenomenon is called degeneracy. The threshold of the toric code obtained by taking account optimally of the degeneracy has been estimated using an Ising model interpretation of the decoding problem [START_REF] Dennis | Topological quantum memory[END_REF]. This threshold is close to p = 0.163 whereas the perfect matching algorithm reaches its threshold at approximately p = 0.155. Note that the renormalization group approach of [START_REF] Duclos-Cianci | Fast decoders for topological quantum codes[END_REF] is one of the rare decoding algorithm of the toric code which is able to make use of the degeneracy of the code.

The second possibility of improvement of the decoding algorithm is the most important potential gain in the performance. It is the correlation between the 2 components, E X and E Z , of the error and consists in using the decoding strategy explained in Section 3. The threshold of the toric code using the X/Z correlations has been estimated close to 0.189 with the Ising model correspondence [START_REF] Bombin | Strong resilience of topological codes to depolarization[END_REF] and is approximately 0.185 with the non-efficient Metropolis decoding algorithm [START_REF] Wootton | High threshold error correction for the surface code[END_REF].

These two remarks are generally true for all surface codes.

A Correlated Perfect Matching Algorithm

To implement the decoding strategy of Section 3 we need to be able to correct errors and erasures when decoding the E Z part. The correction of combinations of errors and erasures for topological codes has been considered by Stace, Barrett, and Doherty [START_REF] Stace | Thresholds for topological codes in the presence of loss[END_REF]. We choose here to adapt Algorithm 1 to find a most likely error for this error model. Denote by E e Z the restriction of E Z to the erased positions, that is the positions (or edges) such that E X = X, and denote by E ē Z its restriction to the non-erased positions. To find the error E Z of syndrome s X such that the weight of E X E Z is minimum, we just have to modify the distance function in Algorithm 1. We introduce the e-distance d e , associated with an erasure e. The usual distance between two vertices of a graph G is the length of a shortest path joining these two vertices. The distance d e is defined similarly but the length of a path is its number of non-erased edges. An e-geodesic between two vertices u and v of G is a path of G of length d e (u, v) joining these two vertices. This provides us a version of the perfect matching algorithm to correct combinations of errors and erasures. It is presented in Algorithm 2. The distance graph based on the e-distance d e is denoted K e (s). The notation G e (u, v) represents a e-geodesic between u and v.

Algorithm 2 Perfect Matching Decoding for errors and erasures

Require: A graph G = (V, E), a subset s ⊂ V of odd size, a set of erased edges e ⊂ E. Ensure: A subset x ⊂ E with end-points ∂(x) = s such that the cardinality of x\e is minimum. 1: Construct the e-distance graph K e (s) associated with s. 2: Determine a minimum weight perfect matching M ⊂ E(K e (s)). 3: return the symmetric difference of all the e-geodesics G e (v i , v j ) for {v i , v j } ∈ M .

Combining Algorithm 1 and Algorithm 2, we obtain Algorithm 3, which is an improved version of the Perfect Matching Decoding algorithm taking partially account of the X/Z correlations.

An example of error over the toric code that can not be corrected with the usual perfect matching decoding but that is corrected with Algorithm 3 is represented in Figure 2.

For Kitaev's toric codes, we obtain a slight improvement of the decoding performance

Algorithm 3 Correlated Perfect Matching Decoding

Require: A tiling G = (V, E, F), a syndrome (s X , s Z ) ⊂ V × V * . Ensure: An error (E X , E Z ) ⊂ E * × E of syndrome (s Z , s X ), such that |E X | is minimum and |E X E Z | is minimum given E X . 1:
Compute E X by applying Algorithm 1 to s Z in the dual graph G * . 2: Compute E Z by applying Algorithm 2 to s X and e = E X in the graph G. 3: return (E X , E Z ). using Algorithm 3, but we cannot overcome the usual threshold since the E X part of the error is decoded using a standard perfect matching algorithm. Nevertheless, as explained in Section 3, the use of the X/Z correlations is well appropriate to asymmetric CSS codes. To define asymmetric surface codes, it suffices to consider non-self dual tilings.

A natural construction of asymmetric surface codes is the family derived from triangular lattices of the torus. For example, the Cayley graph of the group Z/mZ × Z/mZ and the generating set {±(1, 0), ±(0, 1), ±(1, -1)}, described in Figure 1 (b), clearly defines a triangular tiling of the torus. Using Algorithm 1, we remark a threshold for the correction of phase errors at p = 0.066 in Figure 3, whereas the bit-flip error threshold, observed in the dual graph (which is a hexagonal lattice), is very high (more than p = 0.14 for this family of tiling). This implies a depolarizing error threshold at p = 3p /2 = 0.099 for the standard perfect matching algorithm, while Algorithm 3 leads to a depolarizing error threshold at approximately p = 0.133. This good performance is explained by the fact that, while the phase error threshold is low, the error correction in the dual graph exhibits a very good performance and the bit-flip error threshold is high. This allows Algorithm 3 to take into account the X/Z correlations.

Concluding remarks

We proposed a decoding algorithm for CSS codes partially taking into account the correlations between the X component and the Z component of the error for a depolarizing channel. Applied to triangular toric codes, this algorithm exhibits a good performance and clearly improves the threshold. It could be applied to other classes of codes, for instance for color codes, where the decoding algorithm by projection onto 3 surface codes can be adapted to take into account the correlations between the 3 surface codes [START_REF] Delfosse | Decoding color codes by projection onto surface codes[END_REF].

Figure 1 :

 1 Figure 1: A plaquette operator and a site operator acting on a square tiling of the torus in Fig. (a) and on a triangular tiling of the torus in Fig. (b). The opposite boundaries are identified.

Figure 2 :

 2 Figure 2: An example of error correction using Algorithm 3. (a) An error for Kitaev's toric code. (b) The component E X computed at Step 1. of Algorithm 3. (c) The syndrome s X of E Z is given by the vertices marked with '1'. The dashed edges form the erasure defined from E X . (d) The Z component estimated in Step 2. of Algorithm 3. It is the an error of syndrome s X which has minimum weight on the non-erased qubits.

Figure 3 :

 3 Figure 3: Phase decoding performance of Algorithm 1 and depolarizing decoding performance of Algorithm 3 for triangular toric codes of length 3.2 2m

5

  • 10 -2 6 • 10 -2 7 • 10 -2 8 • 10 -2 9 • 10 -2

		Phase error threshold in triangular lattices	Error threshold in triangular lattices using correlation
		10 0		10 0
	block error probability after decoding	10 -1	m=2 m=3 m=4 m=5	block error probability after decoding				m=2 m=3 m=4
			0.1	0.1 10 -1	0.11	0.12	0.13	0.14	0.15
		phase error probability				depolarizing probability

Acknowledgments

The authors wish to thank David Poulin for useful discussions. The work of both authors was supported in part by the French PEPS ICQ2013 program (TOCQ project) and Nicolas Delfosse was supported by the Lockheed Martin Corporation.