Recent Tools

Herbgrind FP Bench

Conclusion

References IEEE-754 floating point arithmetic (1985,2008))

Discretisation (toy system) and precision

Normal floating point:

x = (-1) s • m • 2 e = ± 1.x1x2 . . . xp-1 p bits of mantissa ×2 e Precision: 2u = 1 + -1 = 2 -p 0 1/4 1/2 1 2 4 7.5
Rounding, correct rounding and unit roundoff

•(x) = x for x ∈ F, else •(x) = x(1 + e) with |e| ≤ u/2 (or u)
Correct rounding: best accuracy for +, -, ×, /, √

f : |x||f (x)| |f (x)| , [x||J(x)| |f (x)|
Motto: depends both on problem f and data x Example for summation:

cond( n x i ) = n |x i |/| n x

Recent Tools

Herbgrind FP Bench

Conclusion

Interval Arithmetic (1966) IA at a glance

•

Pros & Cons

• FPCore for fair comparison

• Small size cases, numerically safe case (worst 30% cases error = 5-6 bits)

Others benchmarks: SPEC FPU, Hamming's book,. . .

Conclusion

Numerical accuracy stuff: large and old subject, large literature, many tools but free space for human expertise up to the ideal tools Our Motto = hard issue to automatic tools Herbgrind: a gap in recent developments?

Corsika: tuning to low precision FP formats → full benefit of SIMD speedup e.g. . AVX512 = 16 × binary32

  binary32: u ≈ 5 • 10 -8 , p = 24, e ∈ {-126 . . . 127} binary64: u ≈ 10 -16 , p = 53, e ∈ {-1022 . . . 1023} Accuracy vs. Precision: The Rule of Thumb (RoT) RoT: Accuracy Condition Number ×u Forward error condition × backward error Backward stable in precision u: relative backward error ≈ u Condition number lim δ→0 sup |∆x|≤δ |∆y | |y | / |∆x| |x| , with y + ∆y = f (x + ∆x) and y = f (x).

Differentiable

  Data range or FP arithmetic → intervals + interval operation A sure (•) but too conservative (•) propagation of absolute errors (•) • Dependency problem, wrapping effect, variable decorrelation, conservative inclusion of convex set; intervals containing zero width([x] -[x]) = 2 width((x]), tight function range: tight interval [F ([x])]• Best computing flow driven convex set?endpoint pair, center+radius, subdivisions, Taylor expansions, affine arithmetic, zonotope, . . .[x k+1 ] = R([x k ] for R = R(0, π/4) and x0 = [-ε, ε]Interval RoT[2] width(f (X )) ≤ λ F (X )width(X ), where λ f : Lipschitz-constant of f .Interval sum: log 10 width(sn) vs. log 10 cond(sn) for sn = 1 [3]

  OpenTelemac2D simulation of Malpasset dam break (1959)A five year old dam break: 433 dead people and huge damage Triangular mesh: 26000 elements and 53000 nodes Water flow simulation →35min. after break, 2sec. time step

	1.5e-10 1.5e-10
	1e-10 1e-10	Exact value Horner algorithm
		Exact value
	Prerequisite
	5e-11 5e-11	FPA for dummies				p(x) = (x -2) 9 around x = 2 in IEEE binary64 p(x) = (x -2) 9 around x = 2 in IEEE binary64
	0 0	Errors and Measures			• expanded form • expanded form
	-5e-11 -5e-11	• developed polynomial + Horner algorithm Accuracy vs. Precision: The Rule of Thumb
	-1e-10 -1e-10 Tools
		Old Folks
	-1.5e-10 1.92 -1.5e-10 1.92	1.94 1.94 Interval arithmetic 1.96 1.98 2 1.96 1.98 2	2.02 2.02	2.04 2.04	2.06 2.06	2.08 2.08
		CADNA, verrou
	Interesting example!
		Problem? No problem: exact data!
	One problem + one algorithm + one precision Reproducible simulation? Accurate simulation?
		but different accuracy for different data velocity U Algorithms: The sequential run 0.4029747E-02 0.7570773E-02 0.3500122E-01 velocity V depth H
		one 64 procs run	0.4935279E-02 0.3422730E-02 0.2748817E-01
		one 128 procs run	0.4512116E-02 0.7545233E-02 0.1327634E-01

the rich vs. the poor the good vs. the ugly: summation Example: Industrial case

  i | arbitrarily larger than 1/u when catastrophic cancellation in n x i

	Accuracy Condition number ×u
		Accuracy of the polynomial evaluation [n=50]
	1
	0.01
	Prerequisite
	1e-04 FPA for dummies
	1e-12 1e-10 1e-08 1e-06 Errors and Measures ! 2n cond How to recover the numerical reproducibility of parallel FP computation? Accuracy vs. Precision: The Rule of Thumb Tools Old Folks relative forward error Interval arithmetic Reproducible algorithms, libraries vs. hand-made corrections CADNA, verrou
	1e-14
	1e-16				u
	1e-18		1/u
	100000	1e+10	1e+15	1e+20	1e+25
			condition number

Reproducible enough (i.e. modulo validation) vs. bitwise identical At least to debug parallel vs. sequential, also to validate for production step, to certify for legal process

Floating Point Arithmetic is Error Prone

Counter intuitive FPA Add is not associative • Absorption:

Error Free Transformations (EFT) for +, ×: Herbgrind (2018)

• Dynamic analysis, binaries (Valgring)

• Large programs, different languages, libraries

• Numerical tricks detection: compensation, EFT

• Open platform: front-end to "small sized oriented tools", ...


Sylvain Cluzel
email: sylvain.cluzel@tesa.prd.fr

Mathieu Dervin
email: mathieu.baudoin@thalesaleniaspace.com

José Radzik
email: jose.radzik@isae.fr

Sonia Cazalens
email: sonia.cazalens@cnes.fr

Cédric Baudoin

Daniela Dragomirescu
email: daniela@laas.fr

PHYSICAL LAYER ABSTRACTION FOR PERFORMANCE EVALUATION OF LEO SATELLITE SYSTEMS FOR IOT USING TIME-FREQUENCY ALOHA SCHEME

Keywords: Satellite, NB-IoT, LPWAN, Time-Frequency Aloha, Mutual information

One of the main issues in using a Low Earth Orbit (LEO) satellite constellation to extend a Low-Powered Wide Area Network is the frequency synchronization. Using a link based on random access solves this concern, but also prevents delivery guarantees, and implies less predictable performance. This paper concerns the estimation of Bit Error Rate (BER) and Packet Error Rate (PER) using physical layer abstractions under a time and frequency random scheme, namely Time and Frequency Aloha. We first derive a BER calculation for noncoded QPSK transmission with one collision. Then, we use the 3GPP LTE NB-IoT coding scheme. We analyze the interference that could be induced by repetition coding scheme and propose an efficient summation to improve the decoder performance. Finally, to estimate a PER for any collided scenario, we propose a physical layer abstraction, which relies on an equivalent Signal-to-Noise Ratio (SNR) calculation based on Mutual Information.

INTRODUCTION

Low Powered Wide Area Networks (LPWAN) [START_REF] Raza | Low Power Wide Area Networks: An overview[END_REF] for Internet of Things (IoT) communications is a growing market. LPWAN network operators try to outperform competition by having the widest coverage. Using a Low-Earth Orbit (LEO) satellite system to extend a LPWAN coverage is an attractive option, and has already been proposed [START_REF] Cluzel | 3GPP NB-IoT coverage extension using LEO satellites[END_REF][START_REF] Anteur | Modeling and performance analysis of ultra narrow band system for M2M[END_REF]. These solutions rely on random access protocols, both in the time and the frequency domains. The time randomness is inherent in state of the art LPWAN solutions, because leading market technologies (Sigfox, LoRa) use plain Aloha [START_REF] Abramson | The aloha system: Another alternative for computer communications[END_REF]. Additionally, when using a LEO satellite system, the channel suffers from Doppler shift. As carrier frequency offset compensation is too complex to be used with low-cost terminals, we consider the transmission frequency as random, so the compensation is not needed.

This leads to model the channel access as random both in the time and frequency domains. This subject has been studied under the name of Time-Frequency Aloha (TFA), or Time and Frequency Asynchronous Aloha [START_REF] Goursaud | Random unslotted time-frequency ALOHA: Theory and application to IoT UNB networks[END_REF][START_REF] Almonacid | Throughput performance of time-and frequency-asynchronous ALOHA[END_REF].

This paper focuses on the use of a Low Earth Orbit (LEO) satellite constellation to extend 3GPP LTE NB-IoT. More information about the system can be found in [START_REF] Cluzel | 3GPP NB-IoT coverage extension using LEO satellites[END_REF] and [START_REF] Anteur | Modeling and performance analysis of ultra narrow band system for M2M[END_REF], but this study can be adapted to terrestrial systems. In this system, we encounter a TFA scheme. For that matter, [START_REF] Gaudenzi | Random access schemes for satellite networks: from VSAT to M2M -A survey[END_REF] refers to the exploitation of this colliding scheme as a challenge for next generation systems.

In order to estimate the performance of such highly interfered system, the derivation of an expression of the Bit Error Rate (BER) and the Packet Error Rate (PER) as a function of the local interference scenario is an attractive alternative to computationally expensive physical layer simulation.

The collision probability depending on the traffic load has been well studied during the last years [START_REF] Goursaud | Random unslotted time-frequency ALOHA: Theory and application to IoT UNB networks[END_REF][START_REF] Almonacid | Throughput performance of time-and frequency-asynchronous ALOHA[END_REF]. The authors specify that with no Forward Error Control (FEC), the transmissions are considered as failed when a collision occurs.

Nowadays, in most systems, FEC techniques are used to improve the system performance. In [START_REF] Almonacid | An asynchronous highthroughput random access protocol for low power wide area networks[END_REF], a TFA scheme is studied under several low-rate coding rates, improving the system throughput. Though, the authors do not provide an abstraction method to derive PER estimations: the results are computed via simulations.

In [START_REF] Li | 2D time-frequency interference modelling using stochastic geometry for performance evaluation in Low-Power Wide-Area Networks[END_REF], an interesting and simple abstraction is proposed. The paper aims to model first-order performance for any TFA system. The packet error probability is computed independently of the physical layer implementation, and relies on an estimation of an equivalent Signal-to-Interference-plus-Noise Ratio (SINR), depending on all packets lengths and times of arrival. Then, if a relation linking SINR and PER is known, first-order performance can be derived. The authors' method assumes that collisions cause an equivalent white noise, proportionally to the overlapping areas. This implies a symmetric receiver reaction to frequency and time interference.

In [START_REF] Gaudenzi | Asynchronous Contention Resolution Diversity ALOHA: Making CRDSA Truly Asynchronous[END_REF], a time-only asynchronous Aloha scheme with FEC is assessed, so any collision causes a full frequency overlap. In this case, the interference can be considered as spectrally white. The authors models the behavior of the intra-system interference as an Additive White Gaussian Noise (AWGN), specifying that they believe this approximation to be loose when the number of colliding packets is low.

This paper proposes abstraction methods to compute the BER of a noncoded QPSK transmission, and the PER of a coded QPSK transmission under a TFA access scheme, for a low number of collisions. As described in [START_REF] Cluzel | 3GPP NB-IoT coverage extension using LEO satellites[END_REF], the transmission follows LTE NB-IoT specifications: both a turbocode and a repetition code spreading are used, as specified by 3GPP in LTE (Release 13) standard [11,12]. As the interference also follows the same transmission scheme, some common approximations, as whiteness, cannot be used.

We derive a BER calculation method for a single interferer scenario in section 2, after detailing the parameters used to understand the encountered phenomena.

In section 3, we study the impact of repetition code spreading in TFA schemes, and propose a decoding method.

Finally, we propose an estimation of the PER as a function of any interfering configuration in section 4. This method is adapted from the Mutual Information Effective SNR Mapping (MIESM), mainly used in LTE systems to estimate the PER using several orthogonal subcarriers. Section 5 draws some concluding remarks.

WORST CASE BER FOR NON-CODED TRANSMISSION

In this section we estimate the BER of noncoded transmission, of a TFA system. When the number of interfering transmissions is sufficiently high, the resulting interference is well modeled with a normal distribution. When the number of collisions is low, this approximation does not hold anymore. To the best of our knowledge, the BER of noncoded transmission with few intra-system collisions has not been studied yet. In this section, we present an approximation (c) Representation of all the repetitions of one symbol, and its summation. Fig. 1: Representation of a packet of interest (plain orange) and an interferer (dashed blue) value of a worst-case BER for a collided transmission plus thermal AWGN, in a single-interferer scenario.

We first consider a Packet of Interest (PoI). For the sake of simplification, we suppose the PoI central frequency to be constantly equal to zero. Let's now consider an interfering transmission. Fig. 1a and 1b present the overlapping parameters linking a PoI (plain orange) and an interferer (dashed blue). In Fig. 1a, τt is the normalized time shift (between -1 and 1), τ f is the normalized frequency shift, and ∆f is the difference of Doppler rate between the two transmissions. In Fig. 1b, δt is the normalized time shift between the optimal sampling times of the PoI and the interferer (between 0 and 0.5).

We first present an approximation value of the best-and worstcase BER for a collided transmission plus thermal AWGN. Then, we mix these models to obtain an approximation for the random case.

As presented in Fig. 1a, when the collision appears, only a part of the PoI is interfered by the collision ; we focus on this part of the packet. For the sake of simplification, we consider τt = 0.

The k-th symbol of the received signal r can be modeled as:

r(k) = APoI ePoI(k) + ατ f (k) einterf(k) + n.
(1) APoI is the received amplitude of the PoI, ePoI and einterf are the emitted QPSK symbols of the PoI and the interferer respectively. Then, r is the received symbol, n is the AWGN and ατ f is a complex function, representing the interference as a function of τf. The packet amplitude is APoI = √ PPoI, with PPoI the received power of the PoI. Knowing the behavior of ατ f as a function of τf leads to an estimation of the BER. When δt is close to 0.5, the optimal sampling times are separated, we model the resulting interference as an additive white Gaussian equivalent noise. Only a part of the interferer power, noted Pinterf, goes through the matched filter. The equivalent SNR is calculated as:

SNReq, 0.5 = PPoI/N0 1 + (1 -τf)Pinterf/N0 (2)
Then, by noting ΓBER the function that estimates the BER as a function of the Signal-to-Noise Ratio (SNR), we obtain the equivalent BEReq, 0.5 = ΓBER(SNReq, 0.5).

When δt is close to 0, the reception filter gets matched for the interferer: we observe the BER to be the highest.

We model ατ f to be a random variable, with U the uniform distribution and N the normal distribution, defined as:

ατ f (k) = ρ(k) e jφ(k) , such as ∀k, φ(k) ∼ U [0,2π] and ρ(k) ∼ N (µα, σ 2 α ) (3)
In our model, the power of the interferer is proportional to the overlapping area, leading to µα = Pinterf(1 -|τf|). We use simulations to approximate the value of σ 2 α as a function of τf. This function mainly depends on the used modulation. Now that the interfering signal is modeled, we compute the probability of getting an error in the decoder. This probability is calculated as the expected value of the BER over all the possible received symbols. We first calculate, for every possible received symbol r, the distance d(r) to the threshold, as illustrated in Fig. 2. We use a polar coordinate system r(θ, µα + u), centered in APoIePoI. Fig. 2a illustrates the probabilistic received symbol as being the sum of a PoI symbol (orange dot), and an interfering symbol (mean of the normal distribution in dotted blue). Fig. 2b represents the coordinate system used in the calculation, and the distance d to the threshold is:

d(θ, u) = PPoI 2 -(µα + u) • cos(θ) (4)
The probability of r is now described. So, noting E(•) the expected value and using the tail distribution of the normal distribution Q, we compute the equivalent BER as the expected value of the BER of r,

BEReq, 0 = E [Q(d(θ, u))]. The equivalent SNR is calculated as: SNReq, 0 = Γ -1 BER (BEReq, 0), (5) Γ -1
BER being the inverse fonction of ΓBER, defined earlier. When δt is a random value, uniformly distributed between 0 and 0.5 and using ( 2) and ( 5), we heuristically model the BER as being:

BEReq = ΓBER SNReq, 0.5 + SNReq, 0 2 (6)
In Fig. 3, we present the abstraction performance, with SNRPoI = SNRint = 0 dB, for values of τt between 0 and 1. When τf is close to 0, the interfering transmission is totally overlapping the PoI with the same power, leading to a high BER. Three types of simulations are presented. In blue, the parameter δt is fixed to 0.5. In red, the parameter is fixed to 0. As expected, the BER is the highest for δt = 0, but the difference can only be perceived when τf < 0.3. The green curves represent the simulation and its abstraction when δt is uniformly distributed between 0 and 0.5. The simulation seems to validate the proposed abstractions.

INTERFERED CODED TRANSMISSION DECODING

Forward Error Control (FEC) is a common way to improve the global throughput of a system. Both a turbocode and a repetition code are used in the system presented in [START_REF] Cluzel | 3GPP NB-IoT coverage extension using LEO satellites[END_REF]. In this section, we are going to study for this specific satellite system how the repetitions impact the decoding performance.

Interference analysis

In this system, the repetition coding helps the decoder by adding more diversity to demodulate the messages, but also by lengthening them, in order to simplify the detection and the synchronization.

Once the receiver is synchronized with a PoI, the coherent sum of the symbols is transmitted as an input to the turbodecoder. However, any colliding transmission from the same system is also composed of repetitions. Consequently, the same symbols get summed together, leading hypothetically to an increased interference. Because of the frequency shift, the interferer transmission undergoes a rotation: this phenomena is illustrated in Fig. 1c: one PoI and one interferer symbol are represented over all their repetitions, the PoI symbol having a smaller amplitude. Once the summation is performed (displayed as filled dots), we observe that the signal-to-interferer ratio is not only related to the original amplitudes, but also to the frequency rotation between the same symbol of several repetitions. As a result, the summed interferer could be amplified, or reduced, depending on the frequency shift between the PoI and the interferer.

In a terrestrial use of 3GPP standards, a scrambling should prevent this phenomena. The scrambling is performed by multiplying every repetition by a gold code, parameterized for each repetition by the frame number, the first slot number, a ID number of the cell, and an ID number of the terminal (see section 10.1.3.1 of [11]). Because of the frame numbering, the scrambling parameters will loop every 320 ms. By using 3.75 kHz subcarrier, 16 repetitions and 2 Resource Units by repetition as in [START_REF] Cluzel | 3GPP NB-IoT coverage extension using LEO satellites[END_REF], the scrambling is reinitialized to the same sequence every 5 repetitions, as displayed in Table 1. In this table, the scrambling parameters are presented as a function of the repetition number nrep. A repetition lasts 64 ms; a frame lasts 10 ms, and contains 5 slots, numbered from 0 to 4, as 3.75 kHz subcarriers are used. The frame number is calculated as n f = n f,init + (nrep -1) × 6.4 . The slot number is calculated as ns = 2×(n f,init +10 frac((nrep -1)×6.4)) mod 5, where frac(•) is the fractional part. The table is filled with initial values for frame and slot number equal to 0. We can observe that the sixth repetition follows the same scrambling as the first one. So the phenomena that we presented earlier is reduced by 3GPP scrambling as it is, but not fully prevented.

Decoding method

In order to decode a transmission, a sufficiently good channel estimation is needed. However the interfering power varies during the PoI reception, making the estimation more complex. As the demodulation is performed after the synchronization in [START_REF] Cluzel | 3GPP NB-IoT coverage extension using LEO satellites[END_REF], we suppose that the receiver is synchronized with the PoI. The noise power and PoI power are known. In order to improve the signal turbodecoding, we want to estimate the impact of the interfering transmissions on the PoI after the summation. The input of the turbodecoder is the likelihood ratio of the summed symbols. No scrambling is used in this section.

In order to sum the repetitions, a classic method is to use a Maximum Ratio Combining (MRC) [START_REF] Lo | Maximum ratio transmission[END_REF]. In this method, we multiply each repetition m by a weighting factor ρm, in proportion to the signal-to-noise ratio, before the summation. Evaluated in [START_REF] Sakai | Soft-decision Viterbi decoding with diversity combining[END_REF], this method proposes the best result for diversity combining for a collision free Viterbi decoding. However, as presented in the previous section, knowing the power of an interfering transmissions is not sufficient enough to estimate the signal to interferer power ratio.

We compare several summation methods to decode the interfered transmissions.

1. Unweighted summation. All the weighting coefficients are equal, and depend on the estimated power of the PoI. We present PER performance for perfectly estimated power (O.U., for Oracle Unweighted), and a realistic estimation of the power (R.U., for Realistic Unweighted), using the pilot symbols.

2. Maximum Ratio Combining summation, without knowledge on the interfering summation.

3. Optimization solving summation. We propose to solve an optimization problem, which consists in minimizing the power that is not the PoI signal. Let ρ = [ρm] be the the weighting vector, and APoI the amplitude of the received PoI. The total power of the interfered transmission Ptot(ρ) is easily estimated. The PoI power is:

PPoI, sum(ρ) =   Nrep m=1 ρmAPoI   2 (7)
We solve the following optimization problem:

ρopti = argmin ρ Ptot(ρ) -   Nrep m=1 ρm APoI   2 (8)
We use simulations to compare those three methods. We first generate a set of random configurations, with one interferer. The transmission parameters, i.e. power, Doppler shift, drift of the PoI and the interferer, follow the same distributions as presented in [START_REF] Cluzel | 3GPP NB-IoT coverage extension using LEO satellites[END_REF], and correspond to NB-IoT devices received by a LEO satellite. The terminals send their transmission over Nrep = 4, or 8 repetitions, and using a transmitting power PT X = 15, or 23 dBm. For each configuration, we measure a PER for every summation method. The average PER is presented in Table 2. The MRC summation does not provide the best result, because it tends to drop the interfered parts, leading to a lack of information to decode the transmission. The realistic summation method that provides the best result is the optimization-based proposed method. This solution is also the most computationally demanding.

PER ESTIMATION FOR CODED TRANSMISSION

In [START_REF] Li | 2D time-frequency interference modelling using stochastic geometry for performance evaluation in Low-Power Wide-Area Networks[END_REF], the authors propose an abstraction model of the lower layer in order to speed up a whole system simulation. In this section, we propose a novel abstraction method to estimate the Packet Error Rate (PER) of a given collided scenario, and compare it to abstraction proposed in [START_REF] Li | 2D time-frequency interference modelling using stochastic geometry for performance evaluation in Low-Power Wide-Area Networks[END_REF] for TFA schemes. The proposed method is not symmetrical in the time and frequency domains, and relates on a the Mutual Information of each temporal part of the transmission.

Proposed abstraction method

The abstraction is inspired from Mutual Information Effective SNR Mapping (MIESM), as described in [START_REF] Brueninghaus | Link performance models for system level simulations of broadband radio access systems[END_REF].

In a specific scenario, we know the parameters of the colliding transmissions. We know the power P i interf , and the parameters τ i f , and τ i t of the interfering transmission i, as previously defined, and represented in Fig. 1. We also know the power of the PoI PPoI, and the AWGN power N0. Let P be the number of areas defined by P -1 interferers. In Fig. 4, we represent these areas for P = 6. For each area p, we compute a relative duration tp of the area (such as p tp = 1), and we compute a SINR:

SINRp = PPoI N0 + N interf i=1 δ(p, i)(1 -τ i f )P i interf ( 9
)
The function δ(p, i) is equal to 1 if the interferer i is received during area p, and 0 otherwise. For most real world collided scenarios, and with a sufficient scrambling between repetition (at least as defined in the previous section), we believe that considering the intra-system interference as a Gaussian noise is a sufficient approximation to estimate the PER.

Then, the proposed method relies on mutual information. Instead of considering the channel as a single interfered channel, we consider the transmission to be distributed over P channels, with different SINR, and we combine the channels by using the following formula:

SNReq = β1I -1 Nrep P p=1 tp I SINRp β2 (10)
I is the mutual information for a QPSK constellation, and I -1 is its inverse function. This function can be computed easily and only depends on the constellation. However, the calibration parameters β1 and β2 are fixed by simulations without any intra-system interference, and depend on the coding rate and the number of repetitions.

Using interference-free simulations, we can interpolate the PER curve as a function of the SNR for a AWGN channel, Γ(SNR). Then, the β parameters are calibrated in order to minimize the error between the PER of the simulation and the estimated PLR without interference given by Γ(SNReq).

Results

In order to estimate the performance of the estimation of the PER, we define a score, that measures the difference between the simulated PER and the estimated PER. We aim to model first-order system performance, so the whole system transmission loss target is around 10 %, as proposed in [START_REF] Cluzel | 3GPP NB-IoT coverage extension using LEO satellites[END_REF]. The score has been chosen to favor correct coarse estimations of the PER. For every configuration, the score s is defined as s = |PERestim -PERsimu|, taking values between 0 and 1. A score close to 1 means that the packet demodulation probability has been badly estimated. We display in Fig. 5 the Cumulative Density Function (CDF) of this score. The best estimations show a CDF close to 1 for a low score, meaning that most of the collision configurations have been well estimated.

We present the CDF for two methods: the SNR-based abstraction based from [START_REF] Li | 2D time-frequency interference modelling using stochastic geometry for performance evaluation in Low-Power Wide-Area Networks[END_REF] (blue), and the proposed MIESM-based abstraction (dashed red). The estimations are compared with different number of repetitions, in order to validate these abstractions for different ranges of SINR. The number of collisions is fixed to 1, 5 and 10.

On the one hand, as expected, the two abstractions propose a good estimation of the PER, especially when the number of collisions is low. When the number of repetitions is higher, the expected PER is close to 0, this is why the two abstractions show good scores.

On the other hand, the proposed abstraction shows slightly better performance than the SNR-based abstraction, especially when the number of repetitions is low. This is due to the importance of knowing the impact of each interferer on every repetition: when the number of repetitions decreases, the importance of each repetition increases in relation to the total number of repetitions. We believe mutual information based methods could lead to more accurate results when the number of collision is low.

Having a good estimation for these moderately loaded scenarios with higher PER is a key feature to propose a relevant sizing of the whole system. Refinements such as taking the Doppler rate into account, should improve the estimation and should be investigated.

CONCLUSION

In this paper, we discuss the estimation of BER and PER for transmission following a TFA scheme for coded and noncoded transmission. The focus is made on low number of collisions. We derive a worst-case BER under the case of a single collision for noncoded transmissions. Then we discuss the impact of collisions when using repetitions, such as proposed by 3GPP NB-IoT, on a TFA scheme, and we propose a decoding method that takes into account these repetitions and minimizes the PER. Then, we propose an abstraction of this PER based on mutual information, and compare it to a stateof-the-art SNR-based one; our estimation shows better performance when the number of collisions is low and similar results when the number of collisions is high. This abstraction will be used to estimate the massive number of device the system should accommodate.

  Representation of a noiseless eye diagram of two transmissions. Dotted lines are the optimal sampling times.

  Distance to threshold Fig. 2: Representation of the parameters used in (4).

5 Fig. 3 :

 53 Fig. 3: BER abstraction for equal power PoI and interferer

Fig. 4 :

 4 Fig. 4: Representation of the distinct interfered parts (in grey) of a collided PoI (in orange).

10 (b) Nrep = 8 Fig. 5 :

 1085 Fig. 5: Performances of the abstraction methods.

Table 1 :

 1 Example of gold code parameters for a transmission.

	Repetition number nrep	1 2	3	4	5	6	7	8
	Frame number n f	0 6 12 19 25 30 36 42
	n f mod 2	0	0	0	1	1	0	0	0
	First slot number ns	0 3	1	4	2	0	3	1
	ns/2	0	1	0	2	1	0	1	0

Table 2 :

 2 Decoding performance of summation methods.

	Simulation parameters		Performance (PER)
	N interf	Nrep	P T X	R.U	O.U	MRC Optim.
	1	4	23 dBm	0.01 0.01	0.13	0.01
	1	4	15 dBm	0.86 0.86	0.93	0.86
	1	8	15 dBm	0.62 0.57	0.70	0.60


Sylvain Cluzel
email: sylvain.cluzel@tesa.prd.fr

Laurent Franck

José Radzik
email: jose.radzik@isae.fr

Sonia Cazalens
email: sonia.cazalens@cnes.fr

Mathieu Dervin
email: mathieu.baudoin@thalesaleniaspace.com

Cédric Baudoin

Daniela Dragomirescu
email: daniela@laas.fr

3GPP NB-IoT coverage extension using LEO satellites

Machine-Type Communications are meeting a growing interest on the consumer market. Dedicated technologies arise to support more robust communications involving a massive number of low cost, low energy-consuming devices This paper discusses the coverage extension of a Low-Powered Wide Area Network using a Low Earth Orbit satellite constellation, benefiting from the improved performance of a recent standard. The transmission complies with the user equipment specifications standardized as NB-IoT by 3GPP in Release 13. This radio technology is an update on LTE standard with enhanced performances: the supported path loss can be 20 dB higher than with legacy LTE. This improvement makes satellite-compatible the small and energy-constrained devices. A specific unidirectional system is defined, and a link budget is derived. Also, a receiver architecture is presented, that takes into consideration satellite channel specific impairments.

I. INTRODUCTION

The Internet of Things (IoT) relates to all technologies related to the automated communications between any devices [1], such as control or sensing. Machine-to-Machine (M2M) communications ensure these devices are able to connect and communicate together in a specific network.

M2M market is expected to grow in the forthcoming years. The networks are now designed to handle a massive number of devices. This trend, motivated by industrial needs, will lead to the possibility of interconnection of more and more devices together, on different scales. Wireless networking is one of the growing topics related to IoT.

Low Power Wide Area Networks (LPWAN) are dedicated to connect devices with very specific constraints. The transmission can be long-range (over 10 km), while the device autonomy is typically about 10 years with two AA batteries, constraining the power consumption. Moreover, the system sizing usually relies on sparse resource allocation to the users, so as to handle a massive number of connected devices: several hours, or days, can be allowed between successive messages, the data rate being usually very low, typically below 500 bps. Non-critical data communications are not driven by the transmission delay, which can reach up to a few hours.

Several technologies are emerging to support such transmissions [START_REF] Raza | Low power wide area networks: An overview[END_REF]. On the one hand, narrowband technologies use small bandwidth (less than 25 kHz, down to 100 Hz) to improve the link budget. On the other hand, spread spectrum technologies are well suited to manage intra-system interference.

Terrestrial networks involve costly ground infrastructure deployment, leading the service providers to focus on the most populated areas.

The use of satellites is an attractive option to extend LPWAN service coverage, and all the more if it does not imply specific user terminals. In [START_REF] Anteur | Modeling and performance analysis of ultra narrow band system for M2M[END_REF], the authors study a M2M satellite system using ultra narrowband.

In this paper, we explore the opportunity to address 3GPP LTE NB-IoT compatible devices [START_REF][END_REF], [5], in a satellite network based on a Low Earth Orbit (LEO) constellation. In terrestrial NB-IoT, the medium access is controlled by the eNB (evolved Node B), which ensures that transmissions are not interfered through a combination of supervised timing synchronization and orthogonal channel allocation. In the proposed system, the user terminals are operated regardless of the space segment. We consider a unidirectional return link from the terminal to the satellite. There is no clock synchronization between the devices and the satellites. The devices transmit the data on a specific carrier, but on a randomly chosen subcarrier. In the proposed architecture, the satellites continuously and transparently digitize the received signal.

In absence of synchronization and scheduling, the received signal at satellite input will be impaired by intra-system interference. On the one hand, the high speed motion of LEO satellites induce a Doppler shift, resulting in intercarrier interference. On the other hand, signal collisions will occur in this unscheduled transmission scheme. A similar twodimension interfering scheme has already been studied for Time-Frequency Aloha [START_REF] Goursaud | Random unslotted time-frequency ALOHA: Theory and application to IoT UNB networks[END_REF]. In is shown in [START_REF] Almonacid | Throughput performance of time-and frequency-asynchronous ALOHA[END_REF] that thanks to forward error correction (FEC), a high global system throughput may still be reached in such interfered system.

With no FEC, the expected throughput is low. But it is known that adding a FEC on the receiver side leads to a better global performance. NB-IoT standard allows the use of both turbocode spreading and repetition code spreading. Using a code makes the transmission more robust to interference. The performance improvement with a simple coding scheme has been assessed in [START_REF]An asynchronous high-throughput random access protocol for low power wide area networks[END_REF]. However, to the authors' knowledge, the performance of a repetition code spreading under a Time-Frequency Aloha interference scheme has not been studied yet. Repetitions are being used for a coherent summation, calling for a specific demodulation model. Thus, in order to improve the demodulation of these possibly highly interfered transmissions, a Successive Interference Cancellation (SIC) technique is used. In order to ease the satellite payload complexity, the digitized signal is assumed to be processed at the ground station.

Intra-system interference penalty is all the more impacting as the channel load grows. In this paper, we focus on the performance of a single end-to-end satellite link. The intrasystem interference cancellation is not addressed at this stage where the study purpose is to assess the satellite link feasibility and the related demodulation performance. Future work will address global throughput assessments for highly loaded systems, involving interference management.

The system model is presented in section II. A brief presentation of NB-IoT and a link budget is proposed. In section III, we define the performance metric and the load of such a system. Then, we propose a signal processing strategy that takes into account the satellite channel in section IV. Our algorithm is evaluated without interference in section V, in order to derive the link performance in terms of packet error rate (PER). In section VI, we draw some concluding remarks and propose some future work.

II. SYSTEM OVERVIEW

A. System presentation

A store and forward satellite system based on a LEO constellation is considered in this study, as illustrated on Figure 1. Satellites are used to digitize, and store the whole user bandwidth. Then, the data is transmitted to dedicated stations using a specific high throughput link. The satellite payload could be dedicated to the system, or hosted by an other constellation. In the ground station, the received signal is processed to extract the user messages from the aggregated signal. For this study, we consider that the satellite component allows to extend the terrestrial NB-IoT system, based on a dedicated bandwidth. The terminals are assumed to target the satellite system as soon as they are out of the terrestrial system coverage. The terminals transmit the data to the satellite in line of sight, in a dedicated bandwidth, but over a random subcarrier, without any previous time synchronization or coordination. No information on the user terminal locations is available in the satellite system. At the ground station level, the received signal is processed without a priori knowledge on the precise frequency, or packet time of arrival. Under these assumptions, to simplify the detection and the demodulation, we consider the same waveform parameters for all the terminals. More specifically, the modulation order, code rate, and message duration are identical for all terminals. As LEO satellite communications are considered, Doppler shift and rate must be taken into consideration in the transmission channel modeling. On the one hand, the received signal frequency localization bandwidth is extended in the Doppler shift range induced by the satellite move. On the other hand, the received signals are possibly impacted by nonnegligible Doppler drift rates, i.e. the time varying Doppler shift, with a subsequent degradation of the detection performance.

The studied system is characterized in Table I.

B. NB-IoT description

For this study, all transmissions are assumed to be based on NarrowBand-IoT (NB-IoT) air interface. Since June 2013, M2M communications are targeted by 3GPP. New releases have been published to progressively support low-cost, lowpowered devices, making the choice of using narrowband communications.

In August 2016, the LTE release 13 has been updated to handle a new category of terminal, providing low-cost, high autonomy devices, under the name of NB-IoT, with an improved coverage, allowing the devices to support a 20 dB-higher path loss than legacy LTE for low throughput terminals. This improvement is especially due to the possibility of using a smaller bandwidth, down to one single subcarrier transmission (3.75 or 15 kHz). If the 3.75 kHz subcarrier is used, then the symbol rate is four times lower than LTE as defined by former releases. In this paper, we consider only the narrowest available subcarrier size. The modulation order, the transport block size index I TBS , the resource units index I RU and the number of repetitions are fixed for every terminal The waveform settings considered in this study are summarized in Table II. The frame structure is illustrated in Figure 2. 1 These parameters, extracted from [5], design the transmissions. These values represent a transmission using QPSK, containing 104 bits payload, and of duration 1.024 s. Note that there is no delay between repetitions.

C. Link budget

In this paragraph, we aim to express the distribution of SNR, and Doppler drift rate of the transmissions. The SNR distribution, denoted f SNR , is derived analytically, and the Doppler rate is obtained by simulations. We model the satellite channel from the terminal to the satellite with an additive white Gaussian noise. Receiver signal to noise ratio (SNR) notably depends on the relative position of the terminal and the satellite. The parameter α refers to the angle between a fixed terminal position and the satellite nadir, and d to the distance between the terminal and the satellite, as represented in Figure 1. Let's α max be the maximum value taken by α, i.e. when a terminal is at the edge of coverage, and the elevation is minimal. Then, denoting k B the Boltzmann constant, and L other the other losses than the free space losses [START_REF] Maral | Satellite communications systems: systems, techniques and technology[END_REF]:

SNR(α) = EIRP • (G/T ) k B • BW • 4πf 0 c 2 • L other • d 2 (1)
We suppose that the user location probability is uniform in the satellite coverage. We use d u , the user spatial density per km with relation to α, defined as:

d u (α) = sin(α) (1 -cos(α max )) (2)
f SNR (x) = 1 SNR (SNR -1 (x)) d u (SNR -1 (x)) (3)
We use the formula of change of variable in the probability function, recalled in (3), to obtain f SNR represented in Figure 3a. The lower received signal-to-noise ratio (i.e. when the terminals are the edge of coverage) is around -5 dB.

In [START_REF] Adhikary | Performance evaluation of NB-IoT coverage[END_REF], SNR values about -12 dB leads to a block error rate lower than 10%. On Figure 3a, we can see that the worst SNR met in the proposed system is higher than the one supported by terrestrial NB-IoT. Although this is not evaluated here, it can be inferred that this margin is likely to be useful in case of intra-system interference.

Figure 3b displays the probability density function of the Doppler drift rate, obtained by simulations.

III. SYSTEM METRICS

The channel load is defined as the average number of simultaneously transmitting terminals, divided by the number of subcarriers. Due to Doppler spread, the signal bandwidth is widened at the satellite input with a favorable impact over the interference issue. We modify the load in relation to this spreading. The maximum frequency shift, with relation to our system, is:

δ f,max = f 0 • V sat,max /c ≈ 23 kHz (4)
where c refers to the speed of light, and V sat,max the satellite speed that implies the most Doppler shift. Since the Doppler shift range is significantly larger than the subcarrier width, we consider a uniform spectral distribution of the received carriers at satellite level. At satellite input, the aggregated received signal bandwidth is extended to BW sat = 2δ f,max + BW = 226 kHz. If the terminal spatial density and the message transmission rate are known, we define the channel load at satellite level ρ sat = ρ g •BW/BW sat , as a function of the system load ρ g .

The proposed figure of merit to assess the system performance is the satellite load that will ensure a PER of 10 -1 . In this paper, we only study interference free transmissions. This assumption is applicable to light load scenarios.

The channel load depends on the number of terminals, and the message transmission rate. The satellite cover area is S sat = π(α max R earth ) 2 , with R earth the Earth radius. Then, with the terminal density d t , expressed in km -2 , the number of terminals in the satellite coverage is N t = d t S sat . Every terminal is using the satellite system. If the devices transmit N m messages per day, each message lasting T message in seconds, the channel occupation ratio τ of one terminal is:

τ = N m • T message 3600 • 24 (5)
The average number of simultaneously transmitting terminals under the satellite coverage is λ n = τ N t , randomly spread over 48 subcarriers. The system load is then ρ g = λ n /48. (approximately half a million terminals in the satellite coverage), and N m = 1, the satellite load is ρ sat = 0.11. For this set of parameters, the channel load is considered as low. With this channel load, most of the transmissions do not undergo intra-system interference.

IV. GROUND SIGNAL PROCESSING

A. Receiver description

In the studied system, the satellite digitizes and stores the received data, which are then transferred to dedicated ground stations. The detection and demodulation of the received transmissions are performed in the ground station. The receiver has no specific information on those transmissions, especially about the Doppler-shifted central frequency, the Doppler rate, or the message time of arrival. As the pilot sequence is known, the receiver will perform a data-aided detection and demodulation.

For the synchronization, a classical approach [START_REF] Mengali | Synchronization techniques for digital receivers[END_REF] has been developed, and is presented in Figure 4.

First, a coarse time and frequency domain signal detection is performed to have an approximation of the packet time of arrival and the central frequency of the received signal, based on a power detector. This step relies on a shorttime Fourier transform. Note that this step could be omitted; consequently the next step would have to be parallelized to detect transmissions over the whole bandwidth; a massive computation power would be needed.

Then, a finer time-frequency detection is performed using a filter-bank. As detailed in the next section, in order to compensate for the possibly significant Doppler rate impact over the message duration (see Figure 3b), the frequency is separately estimated at the beginning and at the end of the packet.

Then, a phase correction is performed by solving an optimization problem, as described in [START_REF] Morelli | Doppler-rate estimation for burst digital transmission[END_REF].

Finally, a classical demodulation method is implemented. A power and SNR estimation is performed on the synchronized signal. Using this information, we compute the log-likelihood ratio (LLR) of each symbol. These LLR are summed over each repetition, as proposed in [START_REF] Sakai | Softdecision Viterbi decoding with diversity combining[END_REF] to optimize the decoding performance. The result of the summation is provided to the turbodecoder.

B. Presentation of our detection method

The receiver should be able to detect a 3.75 kHz transmission within 226 kHz wide sampled signal. We first compute a 2D-power detection, using a short-time Fourier transform (STFT). In our implementation, the STFT resolution is similar in the time and in the frequency domain, to represent the transmission with as many STFT samples as possible. Then, Then, a classical fine frequency and time estimation structure is considered, using banks of Wiener filter (FB), as described in [START_REF] Mengali | Synchronization techniques for digital receivers[END_REF]. Each filter is a frequency-shifted crosscorrelation between the pilot symbol sequence, and the received signal. The global maximum of the cross-correlation values gives an estimated time and frequency of the transmission. The number of required filters is deduced from the residual error after coarse detection, and from the frequency resolution of the filter bank.

∆ 2 D -MFB ∆ 1 D -MFB f beg δD -MFB ∆ 1 D -MFB ∆ 2 D -MFB
As displayed in Table II, the duration of the transmission is approx. 1 s; in order to cope with possibly high Doppler rate (up to 200 Hz/s), the cross-correlation is split into two shorter parts. We perform each cross-correlation over N p pilot symbols. The first and last symbols are respectively used to estimate the carrier frequency at the beginning and at the end of the packet. The value of N p has been heuristically set to N p = 150, that being approximately one third of the pilot symbol sequence.

To improve the detection performance in presence of Doppler rate, the filter banks are modified to include a frequency drift over the cross-correlation duration. In this purpose, the reference pilot sequence is multiplied by a quadratic complex chirp. The chirp frequency drift rate is denoted as ∆ D , and the MFB is mentioned ∆ D -MFB. If the difference between the chirp rate and the signal rate is too significant, the pilot sequence is not detected. So, to counterbalance the high range of Doppler rate, n d values of ∆ D are used, so that the difference between ∆ D and the effective frequency rate of the transmission is minimized. Once the beginning and ending frequencies, f beg and f end , are estimated by the MFB with the highest likelihood, we obtain an approximation of transmission's frequency rate δ D . This approximation is noted δD . Then we use the MFB a second time, using ∆ D = δD to estimate the f beg and f end with less error. This finer detection is represented in Figure 5 for n d = 2.

The frequency estimation at the beginning and at the end of the packet must be very accurate to ensure the packet demodulation. In our simulation scenarios, the residual error is found inferior to 10 Hz.

A joint phase and frequency estimation is finally performed as described in [START_REF] Morelli | Doppler-rate estimation for burst digital transmission[END_REF], by minimizing the mean square error between the known pilot sequence and the received data multiplied by a chirp, the parameters of which are estimated to fit the residual error after coarse and fine frequency detection.

V. SIMULATION RESULTS AND ANALYSIS

This section is dedicated to the performance assessment of the proposed receiver in absence of interference.

A NB-IoT standard compliant signal is generated on a random subcarrier. The transmission is affected by the channel, especially by the Doppler effect. For this transmission, we perform the detection, the synchronization, and then the demodulation.

If the receiver manages to detect the transmission and estimate its frequency parameters, we consider the synchronization as successful. In such case, the turbodecoder decodes the transmission and the Cyclic Redundancy Check (CRC) returns no error. However, if the receiver does not estimate the transmission parameters with enough accuracy, the synchronization is a failure: the CRC detects an error, and the packet is lost. As the turbo-decoder would guarantee quasierror free transmissions over an AWGN channel in this SNR range, all packet losses in our simulations can be imputed to the detection or to the synchronization.

In Figure 6, we compare the performance of our detection method to a reference case, without any frequency drift, using unmodified FB (dashed blue curve). These filters are designed to estimate the frequency shift. Then, we introduce the Doppler drift, without changing the method (crossed red curve). We can see that the Doppler drift degrades the performance. This is due to the nondetection of the transmissions that are affected by a high frequency drift rate, by the unmodified FB.

We present the performance of the MFB, using first n d = 1 in yellow; the PER is lower, as expected. So the modifications of the filter bank imply a better resilience to the frequency drift. Increasing the value of n d decreases the PER. Then, n d = 2 and n d = 3 are tested, and the performances are close to the case without frequency drift.

The very low PLR transmissions have a negligible impact on the system performance, which is mainly affected by the noisiest signals. A mean PER is estimated for the system, as a weighted average of SNR-specific PER estimations, using the SNR distribution illustrated. in Figure 3.

With this set of parameters, the system PER is found close to 3 %. As a conclusion, the goal of PER< 10 -1 is reached for very low loads.

VI. CONCLUSIONS

In this paper we describe a new system, designed to extend 3GPP NB-IoT terrestrial coverage to worldwide, using a constellation of LEO satellites. Our system is based on a unidirectional and unsynchronized link from the terminals to the satellite. Each satellite of the constellation digitizes and stores the targeted spectrum until it can be forwarded to a ground station, where the detection and demodulation is processed by the receiver. We derive an equivalent load for our satellite system, and we estimate an addressable number of terminals in the satellite coverage. The proposed detection algorithm is designed to minimize the impact of the satellite channel and more particularly the Doppler drift, in the demodulation performance. The global performance of our system is assessed under the assumption of low channel load.

Future research work will focus on the impact of intrasystem interference to extend the study towards high load channel scenarios.

Fig. 1 :

 1 Fig. 1: Illustration of the system

Fig. 3 :

 3 Fig. 3: SNR (3a) and a Doppler rate (3b) distribution for a 3.75 kHz subcarrier.

Fig. 4 :

 4 Fig. 4: Block diagram of the proposed synchronization strategy With the previously chosen parameters, d t = 0.05 km -2(approximately half a million terminals in the satellite coverage), and N m = 1, the satellite load is ρ sat = 0.11. For this set of parameters, the channel load is considered as low. With this channel load, most of the transmissions do not undergo intra-system interference.

Fig. 5 :

 5 Fig. 5: Block diagram of the detection method, for n d = 2

2 MFB, with n d = 3 Fig. 6 :

 236 Fig. 6: Performance of different detection methods.

TABLE I :

 I Some satellite parameters used by our system

	Uplink central frequency, f 0	950 MHz
	Satellite altitude	800 km
	Satellite G/T	-24.1 dBK -1
	Satellite antenna HPBW 1	120 •
	Equivalent minimal elevation	12.9 •
	Equivalent beam coverage radius 1900 km

2 HPBW = Half Power Beam Width.

TABLE II :

 II Reference NB-IoT waveform parameters

	Emitted power	23 dBm
	Bandwidth, BW	180 kHz
	Subcarrier	3.75 kHz
	Modulation	QPSK
	Modulation format	single-tone SC-FDMA
	Pilot symbol spreading one every 7 symbols
	Coding scheme	3GPP 1/3-Turbo code
	Number of repetitions	16
	I TBS	1	3
	I RU	1	1