5G-NR

    5G New Radio
3GPP

    3rd Generation Partnership Project
AC

    address coding
ACF

    autocorrelation function
ACR

    autocorrelation receiver
ADC

    analog-to-digital converter
AIC

    Analog-to-Information Converter
AIC

    Akaike information criterion
ARIC

    asymmetric restricted isometry constant
ARIP

    asymmetric restricted isometry property
ARQ

    automatic repeat request
AUB

    asymptotic union bound
AWGN

    Additive White Gaussian Noise
AWGN

    additive white Gaussian noise
PSK

    asymmetric PSK
AWRICs

    asymmetric weak restricted isometry constants
AWRIP

    asymmetric weak restricted isometry property
BCH

    Bose, Chaudhuri, and Hocquenghem
BCHSC

    BCH based source coding
BEP

    bit error probability
BFC

    block fading channel
BG

    Bernoulli-Gaussian
BGG

    Bernoulli-Generalized Gaussian
BPAM

    binary pulse amplitude modulation
BPDN

    Basis Pursuit Denoising
BPPM

    binary pulse position modulation
BPSK

    binary phase shift keying
BPZF

    bandpass zonal filter
BU

    Bernoulli-uniform
BER

    bit error rate
BS

    base station
BC

    backscatter communication
CP

    Cyclic Prefix
CDF

    cumulative distribution function
CDF

    cumulative distribution function
CDF

    cumulative distribution function
CCDF

    complementary cumulative distribution function
CCDF

    complementary CDF
CCDF

    complementary cumulative distribution function
CD

    cooperative diversity
CDMA

    Code Division Multiple Access
ch.f.

    characteristic function
CIR

    channel impulse response
CoSaMP

    compressive sampling matching pursuit
CR

    cognitive radio
CS

    compressed sensing
CS

    Compressed sensing
CS

    compressed sensing
CSS

    chirp spread spectrum
CSI

    channel state information
CCSDS

    consultative committee for space data systems
CC

    convolutional coding
COVID-19

    Coronavirus disease
SSBC

    spectrum sharing backscatter communications
CW

    continuous wave
COT

    commercial off-the-shelf
DAA

    detect and avoid
DAB

    digital audio broadcasting
DCT

    discrete cosine transform
DDS

    direct digital synthesis
DFT

    discrete Fourier transform
DR

    distortion-rate
DS

    direct sequence
DS-SS

    direct-sequence spread-spectrum
DTR

    differential transmitted-reference
DVB-H

    digital video broadcasting – handheld
DVB-T

    digital video broadcasting – terrestrial
DL

    downlink
DSSS

    Direct Sequence Spread Spectrum
DFT-s-OFDM

    Discrete Fourier Transform-spread-Orthogonal Frequency Division Multiplexing
DAS

    distributed antenna system
DNA

    Deoxyribonucleic Acid
EC

    European Commission
EED

    exact eigenvalues distribution
EIRP

    Equivalent Isotropically Radiated Power
ELP

    equivalent low-pass
eMBB

    enhanced mobile broadband
EMF

    electric and magnetic fields
EU

    European union
EVT

    extreme value theorem
ETSI

    European Telecommunications Standards Institute
FC

    fusion center
FCC

    Federal Communications Commission
FD

    full-duplex
FEC

    forward error correction
FER

    frame error rate
FFT

    fast Fourier transform
FH

    frequency-hopping
FH-SS

    frequency-hopping spread-spectrum
FS

    Frame synchronization
FS

    frame synchronization
FDMA

    Frequency Division Multiple Access
GA

    Gaussian approximation
GF

    Galois field
GG

    Generalized-Gaussian
GIC

    generalized information criterion
GLRT

    generalized likelihood ratio test
GPS

    Global Positioning System
GMSK

    Gaussian minimum shift keying
GSMA

    Global System for Mobile communications Association
HAP

    high altitude platform
HD

    half-duplex
IDR

    information distortion-rate
IFFT

    inverse fast Fourier transform
IHT

    iterative hard thresholding
i.i.d.

    independent, identically distributed
IoT

    Internet of Things
IR

    impulse radio
LRIC

    lower restricted isometry constant
LRICt

    lower restricted isometry constant threshold
ISI

    intersymbol interference
ITU

    International Telecommunication Union
ICNIRP

    International Commission on Non-Ionizing Radiation Protection
IEEE

    Institute of Electrical and Electronics Engineers
ICES

    IEEE international committee on electromagnetic safety
IEC

    International Electrotechnical Commission
IARC

    International Agency on Research on Cancer
IS-95

    Interim Standard 95
LB

    LoRa backscatter
LEO

    low earth orbit
LF

    likelihood function
LLF

    log-likelihood function
LLR

    log-likelihood ratio
LLRT

    log-likelihood ratio test
LOS

    Line-of-Sight
LRT

    likelihood ratio test
LWRIC

    lower weak restricted isometry constant
LWRICt

    LWRIC threshold
LPWAN

    low power wide area networks
LoRaWAN

    low power long range wide area network
NLOS

    non-line-of-sight
MB

    multiband
MC

    multicarrier
MDS

    mixed distributed source
MF

    matched filter
m.g.f.

    moment generating function
MI

    mutual information
MIMO

    multiple-input multiple-output
MISO

    multiple-input single-output
MJSO

    maximum joint support cardinality
ML

    maximum likelihood
MMSE

    minimum mean-square error
MMV

    multiple measurement vectors
MOS

    model order selection
${M}$-PSK

     $M$-ary phase shift keying
${M}$-PSK

     $M$-ary asymmetric PSK
MTC

    machine type communication
MGF

    moment generating function
$M$-QAM

     $M$-ary quadrature amplitude modulation
MRC

    maximal ratio combiner
MSO

    maximum sparsity order
M2M

    machine to machine
MUI

    multi-user interference
mMTC

    massive machine type communications
mm-Wave

    millimeter-wave
MP

    mobile phone
MPE

    maximum permissible exposure
MAC

    media access control
NB

    narrowband
NBI

    narrowband interference
NLA

    nonlinear sparse approximation
NLOS

    Non-Line of Sight
NTIA

    National Telecommunications and Information Administration
NTP

    National Toxicology Program
NHS

    National Health Service
NB-IoT

    narrowband Internet of things
OC

    optimum combining
OC

    optimum combining
ODE

    operational distortion-energy
ODR

    operational distortion-rate
OFDM

    orthogonal frequency-division multiplexing
OOK

    ON-OFF Keying
OMP

    orthogonal matching pursuit
OSMP

    orthogonal subspace matching pursuit
OQAM

    offset quadrature amplitude modulation
OQPSK

    offset QPSK
OFDMA

    Orthogonal Frequency-division Multiple Access
OPEX

    Operating Expenditures
OQPSK/PM

    OQPSK with phase modulation
PAM

    pulse amplitude modulation
PAR

    peak-to-average ratio
PDF

    probability density function
PDF

    probability density function
PDF

    probability distribution function
PDP

    power dispersion profile
PMF

    probability mass function
PMF

    probability mass function
PN

    pseudo-noise
PPM

    pulse position modulation
PRake

    Partial Rake
PSD

    power spectral density
PSEP

    pairwise synchronization error probability
PSK

    phase shift keying
PD

    power density
$8$-PSK

     $8$-phase shift keying
PR

    primary receiver
PT

    primary trasmitter
FSK

    frequency shift keying
QAM

    Quadrature Amplitude Modulation
QPSK

    quadrature phase shift keying
OQPSK/PM

    OQPSK with phase modulator
RD

    raw data
RDL

    ”random data limit”
RFEH

    RF energy harvesting
RIC

    restricted isometry constant
RICt

    restricted isometry constant threshold
RIP

    restricted isometry property
ROC

    receiver operating characteristic
RQ

    Raleigh quotient
RS

    Reed-Solomon
RSSC

    RS based source coding
r.v.

    random variable
R.V.

    random vector
RMS

    root mean square
RFR

    radiofrequency radiation
RIS

    reconfigurable intelligent surface
RNA

    RiboNucleic Acid
Rx

    receiver
SA-Music

    subspace-augmented MUSIC with OSMP
SCBSES

    Source Compression Based Syndrome Encoding Scheme
SCM

    sample covariance matrix
SEP

    symbol error probability
SER

    symbol error rate
SG

    sparse-land Gaussian model
SIMO

    single-input multiple-output
SINR

    signal-to-interference plus noise ratio
SIR

    signal-to-interference ratio
SISO

    single-input single-output
SMV

    single measurement vector
SNR

    signal-to-noise ratio
SP

    subspace pursuit
SS

    spread spectrum
SW

    sync word
SAR

    specific absorption rate
SSB

    synchronization signal block
SR

    secondary receiver
ST

    secondary trasmitter
SF

    spreading factor
TH

    time-hopping
ToA

    time-of-arrival
TR

    transmitted-reference
TW

    Tracy-Widom
TWDT

    TW Distribution Tail
TCM

    trellis coded modulation
TDD

    time-division duplexing
TDMA

    time division multiple access
Tx

    transmitter
UAV

    unmanned aerial vehicle
URIC

    upper restricted isometry constant
URICt

    upper restricted isometry constant threshold
UWB

    ultrawide band
UWB

    Ultrawide band
URLLC

    Ultra Reliable Low Latency Communications
VCO

    voltage-controlled oscillator
UWRIC

    upper weak restricted isometry constant
UWRICt

    UWRIC threshold
UE

    user equipment
UL

    uplink
URLLC

    ultra reliable low latency communications
WiM

    weigh-in-motion
WLAN

    wireless local area network
WM

    Wishart matrix
WMAN

    wireless metropolitan area network
WPAN

    wireless personal area network
WRIC

    weak restricted isometry constant
WRICt

    weak restricted isometry constant thresholds
WRIP

    weak restricted isometry property
WSN

    wireless sensor network
WSS

    wide-sense stationary
WHO

    World Health Organization
Wi-Fi

    wireless fidelity
SpaSoSEnc

    sparse source syndrome encoding
VLC

    visible light communication
VPN

    virtual private network
RF

    radio frequency
FSO

    free space optics
IoST

    Internet of space things
GSM

    Global System for Mobile Communications
2G

    second-generation cellular networks
3G

    third-generation cellular networks
4G

    fourth-generation cellular networks
5G

    5th-generation cellular networks
gNB

    next generation node B base station
NR

    New Radio
UMTS

    Universal Mobile Telecommunications Service
LTE

    Long Term Evolution
QoS

    Quality of Service

# LoRa Backscatter Communications: Temporal, Spectral, and Error Performance
Analysis

Ganghui Lin, Ahmed Elzanaty, , and
Mohamed-Slim Alouini Ganghui Lin and Mohamed-Slim Alouini are with the
Division of Computer, Electrical and Mathematical Sciences and Engineering,
King Abdullah University of Science and Technology, Thuwal 23955-6900, Saudi
Arabia (e-mail<EMAIL_ADDRESS>slim.alouini@kaust.edu.sa).A.
Elzanaty is with the 5GIC & 6GIC, Institute for Communication Systems (ICS),
University of Surrey, Guildford, GU2 7XH, United Kingdom (e-mail:
a.elzanaty@surrey.ac.uk).The source codes can be accessed at
https://github.com/SlinGovie/LoRa-Backscatter-Performance-Analysis.

###### Abstract

LoRa backscatter (LB) communication systems can be considered as a potential
candidate for ultra low power wide area networks (LPWAN) because of their low
cost and low power consumption. In this paper, we comprehensively analyze LB
modulation from various aspects, i.e., temporal, spectral, and error
performance characteristics. First, we propose a signal model for LB signals
that accounts for the limited number of loads in the tag. Then, we investigate
the spectral properties of LB signals, obtaining a closed-form expression for
the power spectrum. Finally, we derived the symbol error rate (SER) of LB with
two decoders, i.e., the maximum likelihood (ML) and fast Fourier transform
(FFT) decoders, in both additive white Gaussian noise (AWGN) and double
Nakagami-m fading channels. The spectral analysis shows that out-of-band
emissions for LB satisfy the European Telecommunications Standards Institute
(ETSI) regulation only when considering a relatively large number of loads.
For the error performance, unlike conventional LoRa, the FFT decoder is not
optimal. Nevertheless, the ML decoder can achieve a performance similar to
conventional LoRa with a moderate number of loads.

###### Index Terms:

symbol error rate (SER); LoRa backscatter (LB); Internet of Things (IoT);
power spectral density

##  I Introduction

The goal of 5G and beyond networks is to realize three core services, i.e.,
enhanced mobile broadband (eMBB), ultra reliable low latency communications
(URLLC), and massive machine type communications (mMTC) [1, 2], to support
extreme network capacity, latency-sensitive critical missions, and a massive
number of connected devices, respectively. In mMTC, a colossal variety of
inexpensive devices is needed to communicate at low power consumption. With
nearly 30 billion devices expected to be connected by 2030 [3], the prospect
of mMTC makes it possible for IoT applications such as smart cities, homes,
and agriculture.

In this regard, many wireless technologies have recently emerged to realize
such low-power communications with low-priced devices. For instance, low power
long range wide area network (LoRaWAN), SigFoX, and narrowband Internet of
things (NB-IoT) are three prominent technologies supporting LPWAN [4, 5, 6,
7]. However, some application scenarios put forward more stringent
requirements on the power consumption that even technologies such as LPWAN
cannot meet.

In this context, backscatter communication (BC) is one of the most remarkable
technologies for ultra-low-power communications with its ability to operate
with sub-milliwatt power consumption. Typical BC systems contain backscatter
tags and transceivers. The tag reflects the incoming excitation signals while
modulating it with its data by changing the incident signal amplitude,
frequency, or phase [8]. The innovation of BC lies in its radio-passive
property that removes the power-consuming analog components such as RF
oscillators, decoupling capacitors, and crystals. Hence, BC devices can
operate with emerging flexible and printed batteries or even harvest
environmental energy[9, 10].

BC can be divided into three categories based on the adopted architecture: (i)
monostatic backscatter, where a tag modulates an excitation signal generated
by a transceiver, then reflects it back to the transceiver; (ii) bistatic
backscatter, where the transmitter (Tx) and receiver (Rx) are separated; (iii)
ambient backscatter, where the excitation signal comes from existing
surrounding RF sources in the environment [11, 12].

Despite the low-cost, low-power merits that BC owns, one of the most
detrimental defects of BC is its limited communication range. The
electromagnetic waves suffer from severe multipath fading and strong
attenuation in urban or indoor environments, significantly reducing the
coverage of BC [13]. In some medical or wearable applications, the human body
also shortens the communication range[14].

One possible solution to this problem is LB, where the backscatter tag
generates chirp spread spectrum (CSS) modulated signals to improve the
robustness of BC[15]. LoRa modulated signals have exceptionally high
sensitivity up to $-149$ dBm, compensating for the long-range path loss as
well as the multipath fading effects. Although LB is similar, to some extent,
to LoRa modulation in terms of high sensitivity, they differ in many aspects.
First, the LB modulated signals only have a finite number of phases because of
the limited number of antenna loads in the backscatter tags. Second,
conventional LoRa Txs are active devices that directly transmit signals to the
gateway, hopefully through a line of sight link. In contrast, in LB the
signals from the Tx pass through the channel between the Tx and tag before
backscattered to the Rx, i.e., double-fading channel.

The aforementioned aspects make a huge difference in analyzing the performance
of LB. More precisely, a double fading channel should be considered in LB,
resulting in a more complicated probability density function (PDF) for the
channel amplitude. The discrete phases also make the waveforms representing
different symbols non-orthogonal at the Nyquist sampling rate, rendering non-
equivalency between the ML and FFT decoders.

###  I-A Related work

In the following, we review the state-of-the-art LB schemes and the
performance analysis regarding BC and LoRa as well as LB.

####  I-A1 The State-of-the-art of LB

In order to increase the range for BC, a LB communication system with a
harmonic cancellation scheme is proposed in [15]. The proposed system consists
of a single-tone carrier Tx, a backscatter tag, and a Rx. The tag uses the
single tone to synthesize CSS modulated signals and reflect them to the Rx. To
avoid interference from the RF source, a frequency shift is introduced in the
tag by multiplying the incoming signals by an approximated cosine and sine
with discontinuous step transitions, resulting in high-frequency components,
i.e., the harmonics. The harmonic cancellation is achieved by adding more
voltage levels, i.e., more antenna loads in the tag, to better approximate a
pure sinusoid. The design presents a reliable wide area network coverage,
i.e., $475$ m from the Tx and Rx, provided by a low-cost device that consumes
$1000$x lower power than a normal LoRa Tx. However, the tag still needs a
battery as the power source rather than harvesting energy from the
environment. In [16], an RF energy harvesting (RFEH) LB scheme is proposed. It
is reported that the tag can self-start up while harvesting RF energy from the
excitation signal as low as $-22.5$ dBm and send the acquired data to a Rx
that is $381$ m away. There are two possible reasons for the shorter range
compared to [15]. One is that the RFEH LoRa tag only uses two antenna loads,
resulting in more harmonics. Another reason is that the RFEH LoRa tag harvest
energy from the RF source instead of having an embedded battery so that a
proportion of the energy is used to power the IC. To resolve the problem of
deployment difficulties of traditional half-duplex LB, the first full-duplex
LB architecture is proposed in [17]. Nevertheless, the full-duplex LB achieves
a shorter range due to a different sensitivity protocol and a reduced link
budget introduced by a hybrid coupler architecture used to reduce cost.

The above-mentioned LB systems all use single-tone carriers and the tag
modulates its bits into a subcarrier to generate LoRa packets. However, there
are other LB systems that use normal LoRa modulated waves as excitation
signals [18, 19, 20, 21, 22]. These ambient backscatter designs use simpler
modulation methods such as ON-OFF Keying (OOK), which leads to less
complicated tag design at a cost of either shorter communication range, more
complicated decoding algorithms, or lower throughput.

####  I-A2 Performance Analysis for BC, LoRa, and LB

For BC, several works have considered its performance analysis in terms of
detection schemes, error rate, and information rate. In [23], the authors
consider frequency shift keying (FSK) with a coherent receiver for the
bistatic backscatter radio channel. Nevertheless, coherent detection is
complex and requires channel state information (CSI) estimation. To avoid this
difficulty, non-coherent detection schemes for ambient BC have been considered
in [24, 25, 26]. With the different detection schemes mentioned above, the
exact bit error rate (BER) for the ambient backscatter system is derived in
[25]. Since it uses unknown ambient signals as carriers, the actual capacity
is of interest. In [26], the achievable rate of the ambient BC system is
analyzed under different channels. Although the performance analysis of
several BC systems with conventional modulation schemes (e.g., binary phase
shift keying (BPSK), FSK, OOK) have been analyzed in the literature, few works
consider long-range CSS modulation such as LB. In the following, we first
discuss the performance analysis of conventional LoRa then LB.

For conventional LoRa, the power spectrum characteristics of normal LoRa are
analyzed in terms of the Fresnel functions[4], in which the derived
expressions are not applicable for LB because of its limited number of phases.
In [27], the BER of LoRa in an AWGN channel is derived by using the Monte-
Carlo approximation, but the method requires lots of computing resources. The
exact SER of orthogonal signaling with non-coherent detection is presented in
[28, eq: 4.5-44], which is applicable for LoRa. However, the equation is only
derived in an AWGN scenario and the high-order combination involved in the
equation will introduce severe numerical problems. The exact SER of LoRa
modulation in various fading scenarios, i.e., Rayleigh, Rician, and Nakagami
fading channels, are analyzed in [29]. Nevertheless, the expressions are
hardly computable because of the remaining high-order combination. In [30],
simple asymptotic BER expressions of LoRa modulation over Rician and Rayleigh
channels are derived at cost of accuracy. The numerical results show poor
accuracy at low signal-to-noise ratio (SNR). In [31], an accurate closed-form
approximation of the BER in both AWGN and Rayleigh fading channels is derived.
However, the analysis for other fading models such as Rician and Nakagami-m
fading channels is not included. In [32], the author proposed a new approach
based on the Marcum function to estimate the BER of LoRa in different
propagation environments, namely, AWGN, Rayleigh, Rician, and Nakagami
channels. In [33], a modification of a union bound on the error probability is
proposed to calculate the BER with less complexity. The approximation is
applicable for both coherent and non-coherent detection. The theoretical
performance analysis is also extended to Hamming-coded LoRa systems. In [34],
a novel modulation scheme is introduced that employs both up-chirp and down-
chirp simultaneously, and its BER is derived in the AWGN channel. Numerical
results demonstrate that this scheme achieves a doubled throughput without a
significant increase in BER. Further scenarios have been studied such as
channel coding[35], interference channels [36, 37], and imperfect
orthogonality[38, 39] for LoRa. Although the performance analysis of
conventional LoRa modulation over AWGN and various fading channels have been
extensively studied in the literature, they do not consider the double fading
channel and the limited number of loads that characterize LB.

In [40], the error performance of LB is discussed in AWGN channels. However,
the essential characteristics of LB for the signal model are not considered.
For instance, the fact that LB modulated signals only have limited available
phases is not considered. Also, since the performance was only in an AWGN
channel, it does not account for the double-fading effects which are crucial
for LB.

###  I-B Contributions

Although some prototypes have proved the feasibility of LB, more theoretical
insights are needed in order to support the design and deployment of LB
systems. Nevertheless, a comprehensive performance analysis of LB
communication systems is still lacking in the literature. In fact, none of the
aforementioned work has considered the main features that characterize LB such
as finite phases and the double fading channel.

In this paper, we comprehensively analyze LB modulation from various aspects,
i.e., temporal, spectral, and error performance characteristics. More
precisely, we provide the first expression for the LB modulated signal that
accounts for the reduced complexity of the system compared to LoRa. The reason
is that the number of antenna loads is finite, limiting the number of
available phases for the signal.

Based on the provided signal expression, we analytically derive a closed-form
expression for the power spectrum of baseband LB signals. The derivation of
the power spectral density involves integration over the non-linear phase
quantization function. The spectral analysis provides a better understanding
of the adjacent channel interference for LB modulation and its relation to the
number of loads (i.e., system complexity) [41].

Also, we derive the optimal decoder for LB, i.e., ML, and compare it with the
conventional FFT decoder. For ML and FFT decoder, we conduct error performance
analysis in AWGN and fading channels. More precisely, we derive closed-form
approximations for the SER in both AWGN and double Nakagami-m fading channels
for fixed transmit power. Besides, given a constraint on the average transmit
power, we provide the optimal power allocation scheme using water-filling and
analyze the corresponding average SER.

The performance analysis of LB is quite involving compared to normal LoRa. For
instance, unlike normal LoRa, the LB waveforms representing different symbols
are not orthogonal when sampled at the Nyquist rate. This leads to a SER
expression involving a product of Marcum Q functions with different shape
parameters. Additionally, the PDF of the cascaded Nakagami-m fading channel
entails integration over a modified Bessel function of the second kind,
putting further the difficulty in numerical evaluation.

The contribution of this paper can be summarized as follows:

  * •

We provide a novel signal expression in the time-domain for LB signal with a
generalized number of loads.

  * •

We derive a closed-form expression for the power spectrum of LB signals and
investigate whether LB meets the regulation of ETSI for adjacent channel
interference.

  * •

We provide a ML-based decoder for LB.

  * •

We derive closed-form approximations for the SER of LB in both AWGN and double
Nakagami-m fading channels, considering both FFT and ML decoders.

  * •

We study the SER performance of LB for various power allocation schemes, i.e.,
fixed and optimized power allocation techniques.

###  I-C Notations

TABLE I: Table of Symbols Symbol | Meaning
---|---
$D$ | decoder type
$\sigma^{2}$ | noise power per dimension
$\gamma$ | SNR
$E_{\textrm{s}}$, $\overline{E_{\textrm{s}}}$ | instantaneous, average symbol energy
$a$, $i$ | transmit symbol index, decoder output bin index
$\mathcal{L}^{D}_{\left(a,i\right)}$ |  $i$-th decoder output bin value for symbol $a$
$\hat{\mathcal{L}}^{D}_{\left(a\right)}$ | $\mathop{\max}\limits_{i,i\neq a}\left\\{\mathcal{L}^{D}_{\left(a,i\right)}\right\\}$
$\xi^{D}_{\left(a,i\right)}$ | decoder outputs without noise
$\mathcal{W}^{D}_{\left(a,i\right)}$ | noise at the decoder
$\mu_{a}$, $\sigma_{a}^{2}$ | expectation and variance of $\mathcal{L}^{D}_{\left(a,a\right)}$
$d_{1}$, $d_{2}$ | link distances of the Tx-tag and tag-Tx links
$m$, $\Omega$ | shape and spread parameter of Nakagami distribution
$r$, $v$, $n$ |  $m/\Omega$, $m_{1}+m_{2}$, and $|m_{1}-m_{2}|$
$x_{t}$, $\omega_{t}$, $N_{\textrm{GH}}$ | points, weights, and number of points in GH quadrature
$x_{e}$, $\omega_{e}$, $N_{\textrm{GL}}$ | points, weights, and number of points in GL quadrature

Throughout this paper, we denote the Tx-tag and tag-Rx links with subscripts 1
and 2, respectively, vectors with bold small letters, random variables with
calligraphic letters, e.g., $\mathcal{L}$, the PDF of $\mathcal{L}$ with
$f_{\mathcal{L}}\left(l\right)$, the cumulative distribution function (CDF) of
$\mathcal{L}$ with $F_{\mathcal{L}}\left(l\right)$, the expectation and
variance of $\mathcal{L}$ with $\mathbb{E}\left[\mathcal{L}\right]$ and
$\mathbb{V}\left[\mathcal{L}\right]$, respectively, complex Gaussian
distribution that has mean $\mu$ and variance $\sigma^{2}$ with
$\mathcal{CN}\left(\mu,\sigma^{2}\right)$, complex conjugate of a complex
variable $x$ with $x^{*}$, and inner product of vectors $\boldsymbol{\alpha}$
and $\boldsymbol{\beta}$ with $<\boldsymbol{\alpha},\boldsymbol{\beta}>$. The
main symbols considered in this paper are listed in Table I.

###  I-D Structure of the paper

The remainder of this paper is structured as follows. First, the system model
is presented in Section II, and the spectral analysis of LB is provided in
Section III. In Section IV, the error performance in AWGN and fading channels
are investigated. Then, the numerical results are presented in Section V, and
Section VI concludes this paper.

##  II System Setup and Signal Model

In this section, we first explain the system setup and channel model for LB
communication systems. Then, we provide the analytical expressions for LB
signals, and we propose two decoders for LB.

###  II-A System Setup and Channel Model

Figure 1: Considered LB communication system.

We consider a LB communication system with a Tx generating a carrier wave, a
backscatter tag, and a Rx, as shown in Fig. 1. The Tx sends a single-tone
excitation wave to the tag with carrier frequency $f_{\textrm{c}}$. Then, the
tag modulates the excitation signal with its data using CSS modulation and
reflects it to the Rx. We consider a scheme where the tag introduces a
frequency offset $\Delta f$ to the incoming signal, shifting it to a channel
centered at $f_{\textrm{c}}+\Delta f$111$\Delta f$ is typically $3$ orders of
magnitude lower than $f_{\text{c}}$.to avoid interference from the direct link
between Tx and Rx [15]. LB networks can be scaled by considering different SFs
and/or carrier frequencies with negligible interference among users.

Also, we consider a flat double fading channel, where the complex channel
gains of the Tx-tag and tag-Rx links are represented by $h_{1}$ and $h_{2}$,
respectively. The received signal, sampled at the Nyquist rate, can be
expressed as

${r_{a}}\left[k\right]=h_{\textrm{1}}\,h_{\textrm{2}}\,\sqrt{E_{\textrm{s}}}\,{x_{a}}\left[k\right]+{\omega}\left[k\right],\,\forall
k\in\\{0,1,\cdots,M-1\\},$ (1)

where $E_{\textrm{s}}$ denotes the baseband symbol energy, $x_{a}[k]$ denotes
$k^{\textrm{th}}$ sample of the baseband LB signal with data symbol $a$, and
${\omega}\left[k\right]$ represents the thermal noise drawn from
$\mathcal{CN}\left(0,2\sigma^{2}\right)$. For simplicity, the multiplication
of two channel gains can be expressed in an overall channel gain $h$, i.e.,
$h\triangleq h_{1}h_{2}$. This channel model includes the case of monostatic
BC by considering $h_{1}=h_{2}=\sqrt{h}$.

In the following, we explain the backscatter process and provide analytical
expressions for LB baseband signals.

###  II-B Backscatter Process and Signal Model

The reflected LoRa signals are synthesized by the RF switch that toggles
between different antenna loads to create a set of reflection coefficients.
The set of coefficients is selected to change the amplitude and phase of the
incoming signal $S_{\textrm{in}}$. More precisely, let us define
$\Gamma_{\textrm{L}}$ and $\Gamma_{\textrm{A}}$ as the load impedance and
antenna impedance, respectively. The reflected signal $S_{\textrm{out}}$ can
be expressed as

$S_{\textrm{out}}\left(t\right)=\frac{\Gamma_{\textrm{L}}\left(t\right)-\Gamma_{\textrm{A}}}{\Gamma_{\textrm{L}}\left(t\right)+\Gamma_{\textrm{A}}}\,S_{\textrm{in}}\left(t\right)=\left|\Gamma_{\textrm{T}}\left(t\right)\right|e^{j\theta_{\textrm{T}}\left(t\right)}\,S_{\textrm{in}}\left(t\right),$
(2)

where the change of $\Gamma_{\textrm{L}}\left(t\right)$ is discrete so that
the possible number of phases $\theta_{\textrm{T}}\left(t\right)$ and
amplitudes $\left|\Gamma_{\textrm{T}}\left(t\right)\right|$ are limited.
However, to perfectly synthesize LoRa signals, we need an infinite number of
phases so that the backscattered signals are only approximations of the ideal
sinusoidal waves. For this reason, increasing the number of phases in LB can
enhance the system performance and approach traditional LoRa systems. On the
other hand, increasing the number of possible phases raises not only the
design complexity by adding more loads, but also the switching frequency,
reducing the tag service life.

(a)

(b)

(c)

Figure 2: IQ components comparison of backscattered waves (red) and sinusoidal
waves (blue). The first, second, and third columns represent two, four, and
eight phases, respectively: (a) phase diagram, (b) in-phase components, and
(c) quadrature components.

We first elaborate on some typical cases, i.e., using two, four, and eight
phases to approximate sine and cosine. Then, we generalize the model. As shown
in the first, second, and third columns in Fig. 2, we compare the IQ
components of approximated sinusoidal waves (marked with red) with pure
sinusoidal waves (marked with blue).

If we consider only two phases for LB, the in-phase component of the
backscattered wave is zero while the quadrature component is a square wave. To
avoid interference with the direct link, a frequency shift $\Delta f$ should
be introduced. In normal non-BC systems, the frequency shift is realized by
multiplying the original signals centered at $f_{c}$ by a sinusoid with
frequency $\Delta f$. However, due to the limited number of loads in BC tags,
it is introduced through multiplication with square waves with frequency
$\Delta f$, resulting in additional harmonics [16]. More precisely, while
moving the LoRa signals into the channel centered at $f_{c}+\Delta f$, a
mirror copy centered at $f_{c}-\Delta f$ is also created. Also, there are
harmonics centered at $f_{c}\pm 3\Delta f$, $f_{c}\pm 5\Delta f$, $f_{c}\pm
7\Delta f$, etc.

If the system adopts four loads, both in-phase and quadrature components are
square waves. This scheme cancels the mirror copies of the spectra.
Nevertheless, the harmonics are preserved at $f_{\textrm{c}}-3\Delta f$,
$f_{\textrm{c}}+5\Delta f$, $f_{\textrm{c}}-7\Delta f$, etc[14]. When
considering eight loads, two more voltage levels are added to the approximated
sinusoids. The staircase-like waveforms cancels at least the harmonics
centered at $f_{\textrm{c}}-3\Delta f$ and $f_{\textrm{c}}+5\Delta f$[15].
Also, higher-order harmonics can be canceled by adding more phases if
required.

In this regard, unlike conventional LoRa, LB has a finite number of phases in
addition to several harmonics. In the following, we provide the signal
expressions for the baseband signals accounting for the finite phases.

####  II-B1 Continuous-Time Description

The baseband LB signals are synthesized by switching between different loads
so that they only have a limited number of phases, making them an
approximation to the normal LoRa signals.

To generalize the model, let us consider the number of loads written in a form
of $2^{N}$, where $N\in\left\\{1,2,3,...\right\\}$. The phases are evenly
distributed and rotationally symmetric in the phase diagram at position
${\left(2n-1\right)\pi}/{2^{N}},\,\forall
n\in\left\\{1,2,\cdots,2^{N}\right\\}$ with some examples shown in Fig. 2a.

We introduce to the phase of conventional LoRa signal a mid-rise quantizer to
perfectly model the discrete phases of LB. The quantizer maps the infinite
phases of the LoRa modulation into a finite set of phases in LB. More
precisely, the mid-rise quantization function with $2^{N}$ levels for input
value $x\in\left(-\pi,\pi\right)$ is given by

$\widetilde{Q}_{N}\left(x\right)=\dfrac{\lfloor
2^{N-1}x/\pi\rfloor+\frac{1}{2}}{2^{N-1}/\pi},$ (3)

where $\lfloor\cdot\rfloor$ is the floor function[42]. Let us also specify the
periodicity of the quantization function as
$\widetilde{Q}_{N}\left(x+2\pi\right)\triangleq\widetilde{Q}_{N}\left(x\right)$
for any real input value $x$, making it consistent with the periodic change of
the phase.

Let $T_{\textrm{s}}\triangleq M/B$ be the symbol duration, the instantaneous
phase of LoRa modulation with symbol $a$ for interval
$\left[0,T_{\textrm{s}}\right]$ can given in[4] as

$\hat{\phi}_{a}\left(t\right)=2\pi
Bt\left[\frac{a}{M}-\frac{1}{2}+\frac{Bt}{2M}-u\left(t-\tau_{a}\right)\right],$
(4)

where $B$ denotes the bandwidth of LoRa, $M=2^{\text{SF}}$,
$\text{SF}\in\left\\{7,8,\cdots,12\right\\}$ denotes the spreading factor, and
$\tau_{a}=({M-a})/{B}$ denotes the time instant of the sudden frequency
change. Thus, the instantaneous phase of LB with symbol $a$ is calculated as

$\displaystyle\begin{aligned}
\phi_{a}\left(t\right)&=\widetilde{Q}_{N}\left[\hat{\phi}_{a}\left(t\right)\right].\end{aligned}$
(5)

Thus, the complex envelope for LB waveform of symbol $a$ with unit power is
calculated using (5) as

$\displaystyle x_{a}\left(t\right)$
$\displaystyle=\exp\left\\{j\phi_{a}\left(t\right)\right\\}$ (6)
$\displaystyle=\exp\left\\{j\widetilde{Q}_{N}\left[2\pi
Bt\left[\frac{a}{M}-\frac{1}{2}+\frac{Bt}{2M}-u\left(t-\tau_{a}\right)\right]\right]\right\\}.$

####  II-B2 Discrete-Time Description

Let us consider the simple receiver implementation that sample the received
signal at the chip rate, i.e., ${1}/{T_{\text{c}}}=M/{T_{\text{s}}}=B$, the
$M$ samples in the interval $\left[0,T_{\textrm{s}}\right]$, for
$k\in\left\\{0,1,\cdots,M-1\right\\}$, are

$\displaystyle\phi_{a}\left(kT_{\text{c}}\right)$
$\displaystyle=\widetilde{Q}_{N}\left[\frac{k\pi}{M}\left(2a-M+k\right)-2M\pi
u\left(\frac{k+a-M}{B}\right)\right]$
$\displaystyle=\widetilde{Q}_{N}\left[\frac{k\pi}{M}\left(2a-M+k\right)\right],$
(7)

where the last equality is due to the fact that $2M\pi u\left(\cdot\right)$ is
always an integer multiple of $2\pi$ for all $k$ and
$\widetilde{Q}_{N}\left(x\right)$ is periodic of $2\pi$. Normalizing the total
energy of the samples, we have the discrete-time expression of the baseband LB
signal with $2^{N}$ quantization phases

$\displaystyle{x_{a}}\left[k\right]=\sqrt{\dfrac{1}{M}}\exp\left\\{j\widetilde{Q}_{N}\left[\dfrac{k\pi}{M}\left(2a-M+k\right)\right]\right\\},$
$\displaystyle k\in\left\\{0,1,\cdots,M-1\right\\}.$ (8)

###  II-C LoRa backscatter Decoders

We study the error performance of LB under two decoders, i.e., the ML and FFT
decoders. We consider that the overall channel amplitude
$\left|h\right|=\left|h_{1}\right|\left|h_{2}\right|$ is known to the Rx, and
it performs non-coherent detection.222Non-coherent detectors for LB can
provide a good trade-off between the performance and complexity since the
performance improvement using coherent detectors is only 0.7 dB for
traditional LoRa [43].

####  II-C1 Maximum Likelihood Decoder

The optimum detector in an AWGN channel is the ML detector if all the symbols
are equiprobable[44]. The decoder performs cross-correlation on the received
signal with all possible waveforms and chooses the symbol that has the largest
absolute value of cross-correlation with the received signal, i.e.,

$\displaystyle\hat{a}$ $\displaystyle=\mathop{\arg\max}\limits_{0\leq i\leq
M-1}\left|\left<\boldsymbol{r_{a}},\boldsymbol{x^{*}_{i}}\right>\right|$ (9)
$\displaystyle=\mathop{\arg\max}\limits_{0\leq i\leq
M-1}\left|\sum_{k=0}^{M-1}{r_{a}}\left[k\right]\cdot{x_{i}^{*}}\left[k\right]\right|,$
(10)

where $\hat{a}$ is the symbol decision made by the Rx,
${x_{i}^{*}}\left[k\right]$ is the complex conjugate of
${x_{i}}\left[k\right]$. The ML decoder is the optimum decoder for LB as shown
in Appendix A. However, LB is a relatively high-order modulation that reaches
$\textrm{M}=4096$ for $\textrm{SF}=12$, resulting in a more complicated design
and implementation of the decoder.

####  II-C2 FFT Decoder

The process of FFT decoder follows three steps. First, The received signal is
firstly multiplied by a complex down-chirp signal
${x_{\textrm{d}}}\left[k\right]$

${\tilde{r}_{a}}\left[k\right]={r}_{a}\left[k\right]\,{x_{\textrm{d}}}\left[k\right]={r}_{a}\left[k\right]\,\sqrt{\frac{1}{M}}e^{-j2\pi\frac{k^{2}}{2M}+j\pi
k},$ (11)

where the dechirped signal ${\tilde{r}_{a}}\left[k\right]$ is called the
twisted signal. The symbol decision is made by choosing the maximum value
among M-bin discrete Fourier transform (DFT) output on
${\tilde{r}_{a}}\left[k\right]$

$\displaystyle\hat{a}$ $\displaystyle=\mathop{\arg\max}\limits_{0\leq i\leq
M-1}\left|\textrm{DFT}\left({\tilde{r}_{a}}\left[k\right]\right)\right|$ (12)
$\displaystyle=\mathop{\arg\max}\limits_{0\leq i\leq
M-1}\left|\sum_{k=0}^{M-1}{\tilde{r}_{a}}\left[k\right]\,e^{-j2\pi
ki/M}\right|.$ (13)

(a) LoRa backscatter ($N=2$)

(b) Normal LoRa

Figure 3: The cross-correlation of LB and LoRa symbols when $\textrm{SF}=7$.

It has been proven that the two decoders are equivalent for norma LoRa[4] but
this is not valid for LB. Also, the waveforms representing different symbols
sampled at the chip rate $B$ are non-orthogonal. Fig. 3 compares the cross-
correlation of LB waveforms (4 quantization phases) with that of normal LoRa
at $\textrm{SF}=7$. It shows that the cross-correlation of LB waveforms has
numerous non-zero values while the cross-correlation between different symbols
are all zero for normal LoRa333LoRa symbols are non-orthogonal, it can be
proven orthogonal only when sampled at Nyquist rate.. However, most of the
non-zero correlation values are relatively small, indicating that the symbols
are not strongly correlated.

In Table. II, we show the maximum cross-correlation value between different LB
waveforms for $N\in\left\\{2,3,4,5\right\\}$ and
$\text{SF}\in\left\\{7,8,\cdots,12\right\\}$. The cross-correlation generally
decreases as $N$ and SF increase.

TABLE II: Maximum cross-correlation value between different LB waveforms for
$N\in\left\\{2,3,4,5\right\\}$ and
$\text{SF}\in\left\\{7,8,\cdots,12\right\\}$.

| SF=7 | SF=8 | SF=9 | SF=10 | SF=11 | SF=12
---|---|---|---|---|---|---
N=2 | 0.250 | 0.156 | 0.156 | 0.117 | 0.086 | 0.067
N=3 | 0.125 | 0.082 | 0.107 | 0.064 | 0.071 | 0.050
N=4 | 0.000 | 0.000 | 0.053 | 0.032 | 0.043 | 0.024
N=5 | 0.000 | 0.000 | 0.000 | 0.000 | 0.022 | 0.013

##  III Spectral analysis

The power spectrum of the LB modulation is analytically determined in closed
form in this section. We consider a source that produces a random symbol
sequence where all the symbols are independent, identically distributed with
random variables $\mathcal{B}_{n}$. Assuming that the symbols are
equiprobable, we have

$P\left(\mathcal{B}_{n}=a\right)=\frac{1}{M},\quad
a\in\left\\{0,1,\cdots,M-1\right\\}.$ (14)

The time sequence of the LB waveforms can be expressed with a random process

$\mathcal{I}\left(t\right)=\sum_{n}x_{s}\left(t-nT_{\text{s}};\mathcal{B}_{n}\right)g_{T_{\text{s}}}\left(t-nT_{\text{s}}\right),$
(15)

where $g_{T_{\text{s}}}$ is a rectangular function with width
${T_{\text{s}}}$. The power spectrum density of the random process
$\mathcal{I}\left(t\right)$ can be derived as

$\mathcal{G}_{\mathcal{I}}\left(f\right)=\mathcal{G}_{\mathcal{I}}^{\left(c\right)}\left(f\right)+\mathcal{G}_{\mathcal{I}}^{\left(d\right)}\left(f\right).$
(16)

where $\mathcal{G}_{\mathcal{I}}^{\left(c\right)}\left(f\right)$ and
$\mathcal{G}_{\mathcal{I}}^{\left(d\right)}\left(f\right)$ are the continuous
and a discrete parts of the spectrum, respectively. The continuous and
discrete spectrum in (16) can be obtained by applying frequency domain
analysis of randomly modulated signals for the random process [45]

$\displaystyle\mathcal{G}_{\mathcal{I}}^{\left(c\right)}\left(f\right)$
$\displaystyle=\frac{1}{T_{\text{s}}}\left[\sum_{a=0}^{M-1}\frac{1}{M}\left|S_{a}\left(f\right)\right|^{2}-\left|\sum_{a=0}^{M-1}\frac{1}{M}S_{a}\left(f\right)\right|^{2}\right],$
(17) $\displaystyle\mathcal{G}_{\mathcal{I}}^{\left(d\right)}\left(f\right)$
$\displaystyle=\frac{1}{M^{2}T^{2}_{\text{s}}}\sum_{l=-\infty}^{\infty}\left|\sum_{a=0}^{M-1}S_{a}\left(\frac{lB}{M}\right)\right|^{2}\delta\left(f-\frac{lB}{M}\right),$
(18)

where $\left\\{S_{a}\left(f\right)\right\\}_{a=0}^{M-1}$ are the Fourier
transforms of the waveforms $\left\\{x_{a}\left(t\right)\right\\}_{a=0}^{M-1}$
from the modulator.

The Fourier transform of the complex envelope can be expressed as

$\displaystyle S_{a}\left(f\right)$
$\displaystyle=\int_{0}^{T_{\text{s}}}e^{j\widetilde{Q}_{N}\left[2\pi
Bt\left[\frac{a}{M}-\frac{1}{2}+\frac{Bt}{2M}-u\left(t-\frac{M-a}{B}\right)\right]\right]}e^{-j2\pi
ft}\mathrm{d}t$
$\displaystyle=\int_{0}^{\frac{M-a}{B}}e^{j\widetilde{Q}_{N}\left[2\pi
Bt\left(\frac{a}{M}-\frac{1}{2}+\frac{Bt}{2M}\right)\right]}e^{-j2\pi
ft}\mathrm{d}t+$
$\displaystyle\quad\quad\int_{\frac{M-a}{B}}^{\frac{B}{M}}e^{j\widetilde{Q}_{N}\left[2\pi
Bt\left(\frac{a}{M}-\frac{3}{2}+\frac{Bt}{2M}\right)\right]}e^{-j2\pi
ft}\mathrm{d}t.$ (19)

Leveraging the periodicity of $\widetilde{Q}_{N}\left(\cdot\right)$, we
rewrite the second term in (19) as

$\int_{\frac{M-a}{B}}^{\frac{B}{M}}e^{j\pi\widetilde{Q}_{N}\left[2\pi
Bt\left(\frac{a}{M}-\frac{3}{2}+\frac{Bt}{2M}\right)+2\pi\left(M-a\right)\right]}e^{-j2\pi
ft}\mathrm{d}t.$ (20)

Let us also define $f_{1}\left(t\right)\triangleq 2\pi
Bt\left(\frac{a}{M}-\frac{1}{2}+\frac{Bt}{2M}\right)$, and
$f_{2}\left(t\right)\triangleq 2\pi
Bt\left(\frac{a}{M}-\frac{3}{2}+\frac{Bt}{2M}\right)+2\pi\left(M-a\right)$ as
the two terms within $\widetilde{Q}_{N}\left(\cdot\right)$ of two integral
parts. Let us define $f\left(t\right)$ as

$f\left(t\right)\triangleq\left\\{\begin{aligned}
f_{1}\left(t\right)\quad&,\quad\quad\ 0\quad\leq t<\frac{M-a}{B}\\\
f_{2}\left(t\right)\quad&,\quad\frac{M-a}{B}\leq
t\leq\quad\frac{M}{B}\end{aligned}\right.$

It is worth noting that, $f\left(t\right)$ is a continuous function with
$f_{1}\left(\frac{M-a}{B}\right)=f_{2}\left(\frac{M-a}{B}\right)$. Thus, the
sum of two separated integrals in (19) can be replaced with one integral

$S_{a}\left(f\right)=\int_{0}^{\frac{B}{M}}e^{j\widetilde{Q}_{N}\left[f\left(t\right)\right]}e^{-j2\pi
ft}\mathrm{d}t.$ (21)

$t_{0}$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$t_{5}$$t_{6}$$t_{7}$$t_{\psi}$$-\pi$$-\frac{\pi}{2}$$0$$\frac{\pi}{2}$$\pi$$t$$\widetilde{Q}_{N}\left[f\left(t\right)\right]$
Figure 4: An example of $\widetilde{Q}_{N}\left[f\left(t\right)\right]$ as a
function of time at $N=2$.

It is challenging to directly compute the integral in (21) with the
quantization function inside. However, one possible solution is to divide the
integral interval into slots, as shown in Fig. 4. The value of
$\widetilde{Q}_{N}\left[f\left(t\right)\right]$ remains unchanged in each
intervals divided by a time set $\left\\{t_{m}\right\\}_{m=1}^{\psi}$. Details
about obtaining $\left\\{t_{m}\right\\}_{m=1}^{\psi}$ are provided in Appendix
B. The integral in each slot can be calculated as

$\displaystyle I_{m}\left(f\right)$
$\displaystyle=\int_{t_{m}}^{t_{m+1}}e^{j\widetilde{Q}_{n}\left[f\left(t\right)\right]}e^{-j2\pi
ft}\mathrm{d}t$ (22)
$\displaystyle=\frac{j}{2f}e^{j\pi\widetilde{Q}_{n}\left[f\left(\frac{t_{m}+t_{m+1}}{2}\right)\right]}\left(e^{-j2\pi
ft_{m+1}}-e^{-j2\pi ft_{m}}\right),$

where $({t_{m}+t_{m+1}})/{2}$ can be replaced by any real value between
$t_{m}$ and $t_{m+1}$.

Since we have $f\left(0\right)=f\left({M}/{B}\right)=0$, the overall integral
interval begins with $t_{0}$ and ends with $t_{\psi}$. Thus, the integral in
(19) is calculated as a sum of $I_{m}\left(f\right)$

$\displaystyle\begin{aligned} S_{a}\left(f\right)=\frac{j}{2\pi
f}\sum_{m=1}^{\psi-1}e^{j\pi\widetilde{Q}_{N}\left[f\left(\frac{t_{m}+t_{m+1}}{2}\right)\right]}\left(e^{-j2\pi
ft_{m+1}}-e^{-j2\pi ft_{m}}\right).\end{aligned}$ (23)

By substituting (23) into (17) and (18), we can calculate both the continuous
and discrete parts of the spectrum for baseband LB signals.

##  IV Performance Analysis

In this section, we derive the SER of the LB communication system in terms of
different decoders, channel models, and power strategies.

###  IV-A LoRa backscatter SER Performance over AWGN Channels

The SER in an AWGN channel is a conditional probability conditioned on the
overall channel gain $h$. Thus, we consider $h$ as a constant in this
subsection. We define the r.v. $\mathcal{L}^{D}_{\left(a,i\right)}$ as the
absolute value of $i$-th bin of
$D\in\mathbb{D}\triangleq\left\\{\textrm{ML},\textrm{FFT}\right\\}$ decoder
when symbol $a$ is transmitted. For ML decoder,
$\mathcal{L}^{\textrm{ML}}_{\left(a,i\right)}$ is calculated as

$\displaystyle\mathcal{L}^{\text{ML}}_{\left(a,i\right)}$
$\displaystyle=\left|\sum_{k=0}^{M-1}{r_{a}}\left[k\right]{x_{i}^{*}}\left[k\right]\right|$
(24)
$\displaystyle=\left|h\sqrt{E_{\text{s}}}{\xi_{\left(a,i\right)}^{\text{ML}}}+{\mathcal{W}^{\text{ML}}_{\left(a,i\right)}}\right|,$
(25)

where
$\xi_{\left(a,i\right)}^{\text{ML}}=\sum_{k=0}^{M-1}{x_{a}}\left[k\right]{x_{i}^{*}}\left[k\right]$
is the cross-correlation between the transmitted waveform and the $i$-th
reference signal and
$\mathcal{W}^{\text{\text{ML}}}_{\left(a,i\right)}=\sum_{k=0}^{M-1}{\omega}\left[k\right]{x_{i}^{*}}\left[k\right]$
is the noise projection of the reference signals, following complex Gaussian
distribution with variance $2\sigma^{2}$. For FFT decoder,
$\mathcal{L}^{\textrm{FFT}}_{\left(a,i\right)}$ is

$\displaystyle\mathcal{L}^{\text{FFT}}_{\left(a,i\right)}$
$\displaystyle=\left|\textrm{DFT}\left({r_{a}}\left[k\right]{x_{\textrm{d}}}\left[k\right]\right)\right|$
(26)
$\displaystyle=\left|h\sqrt{E_{\text{s}}}{\xi_{\left(a,i\right)}^{\text{FFT}}}+{\mathcal{W}^{\text{FFT}}_{\left(a,i\right)}}\right|,$
(27)

where
$\xi_{\left(a,i\right)}^{\text{FFT}}=\textrm{DFT}\left({x_{a}}\left[k\right]{x_{\textrm{d}}}\left[k\right]\right)$
denotes the i-th bin of the DFT for the dechirped reference signals and
$\mathcal{W}^{\text{FFT}}_{\left(a,i\right)}=\text{DFT}\left({\omega}\left[k\right]{x_{\textrm{d}}}\left[k\right]\right)$
is the $i$-th bin of the DFT for the dechirped complex white Gaussian noise,
which also follows complex Gaussian distribution with variance $2\sigma^{2}$.
Since (25) and (27) have the same pattern, i.e., the absolute value of a sum
of a complex constant and a complex Gaussian r.v.,

$\mathcal{L}^{D}_{\left(a,i\right)}=\left|h\sqrt{E_{\textrm{s}}}\xi^{D}_{\left(a,i\right)}+\mathcal{W}^{D}_{\left(a,i\right)}\right|,$
(28)

which follows Rician distribution with shape parameter
$\kappa^{D}_{\left(a,i\right)}={\left|h\sqrt{E_{\textrm{s}}}\xi_{\left(a,i\right)}^{D}\right|^{2}}/{2\sigma^{2}}$,
as shown in Appendix C. Without losing generality, if
$\xi_{\left(a,i\right)}^{D}=0$, the shape parameter is zero so that the Rician
distribution becomes Rayleigh distribution.

It is worth noting that the r.v.s in the set
$\left\\{\mathcal{L}^{\text{ML}}_{\left(a,i\right)}\right\\}_{i=0}^{M-1}$ are
not independent because $\mathcal{W}^{\text{ML}}_{\left(a,i\right)}$ are the
noise projection of non-orthogonal waveforms, as shown in Fig. 3a. On the
other hand, the r.v.s in the set
$\left\\{\mathcal{L}^{\text{FFT}}_{\left(a,i\right)}\right\\}$ are independent
since the noise is projected to the directions of orthogonal basis in the DFT
process. The dependence of
$\left\\{\mathcal{L}^{\text{ML}}_{\left(a,i\right)}\right\\}_{i=0}^{M-1}$ make
it challenging to find their joint distribution. However, they are
approximately independent for the following three reasons. First, as shown in
Fig. 3a, most of the cross-correlation values are zero, namely, most of the
r.v.s within
$\left\\{\mathcal{L}^{\text{ML}}_{\left(a,i\right)}\right\\}_{i=0}^{M-1}$ are
independent. Second, the non-zero cross-correlation values are relatively
small so that the dependent r.v.s are not strongly correlated. Also, the non-
zero cross-correlation values decrease as the number of quantization phases or
the spreading factor increases. Nevertheless, ignoring the correlation makes
the derived expression an upper bound to the exact SER[44, Section 4.5-4]. In
the following, we assuming that
$\left\\{\mathcal{L}^{\text{ML}}_{\left(a,i\right)}\right\\}_{i=0}^{M-1}$ are
independent so that the derivation applies to both decoders.

To obtain an expression for the SER, let us define the instantaneous SNR as

$\gamma=\dfrac{\left|h\right|^{2}E_{\textrm{s}}/T_{\textrm{s}}}{2\sigma^{2}B}=\dfrac{\left|h\right|^{2}E_{\textrm{s}}}{2\sigma^{2}M}.$
(29)

The PDF and CDF of $\mathcal{L}^{D}_{\left(a,i\right)}$ conditioned on the
overall channel gain $h$ can be written for $D\in\mathbb{D}$ as

$\displaystyle f_{\mathcal{L}^{D}_{\left(a,i\right)}|h}\left(l\right)$
$\displaystyle=\frac{l}{\sigma^{2}}\exp\left[\dfrac{-\left(l^{2}+C^{2}\right)}{2\sigma^{2}}\right]I_{0}\left(\dfrac{lC}{\sigma^{2}}\right),$
(30) $\displaystyle F_{\mathcal{L}^{D}_{\left(a,i\right)}|h}\left(l\right)$
$\displaystyle=1-Q_{1}\left(\dfrac{C}{\sigma},\dfrac{l}{\sigma}\right),$ (31)

where $C=\left|h\sqrt{E_{\textrm{s}}}\xi_{\left(a,i\right)}^{D}\right|$,
$Q_{1}\left(\cdot\right)$ is the Marcum Q-function of order one [46, eq: 7],
$I_{0}\left(\cdot\right)$ is the modified Bessel function of the first kind
and order zero[47, eq: 9.6.1].

The detection error happens when $\hat{a}\neq a$, i.e.,
$\mathcal{L}^{D}_{\left(a,a\right)}$ is not maximum among the set
$\left\\{\mathcal{L}^{D}_{\left(a,i\right)}\right\\}_{i=0}^{M-1}$. In other
words, the symbol error probability (SEP) can be calculated by comparing the
correct bin, namely, $\mathcal{L}^{D}_{\left(a,a\right)}$, with the maximum of
the noisy bins. Thus, we define
$\hat{\mathcal{L}}^{D}_{\left(a\right)}\triangleq\mathop{\max}\limits_{i,i\neq
a}\left\\{\mathcal{L}^{D}_{\left(a,i\right)}\right\\}$ whose CDF of
$\hat{\mathcal{L}}^{D}_{\left(a\right)}$ can be computed using order
statistics as

$\displaystyle F_{\hat{\mathcal{L}}^{D}_{\left(a\right)}|h}\left(l\right)$
$\displaystyle=\prod_{i=0,i\neq
a}^{M-1}F_{\mathcal{L}^{D}_{\left(a,i\right)}|h}\left(l\right)$
$\displaystyle=\prod_{i=0,i\neq
a}^{M-1}\left[1-Q_{1}\left(\dfrac{\left|h\sqrt{E_{\textrm{s}}}\xi^{D}_{\left(a,i\right)}\right|}{\sigma},\dfrac{l}{\sigma}\right)\right].$
(32)

The conditional SEP using $D\in\mathbb{D}$ decoder given a transmit symbol $a$
is calculated as

$\displaystyle P^{D}_{\textrm{e}|a}$
$\displaystyle=P\left(\hat{\mathcal{L}}^{D}_{\left(a\right)}>\mathcal{L}^{D}_{\left(a,a\right)}\big{|}h\right)$
$\displaystyle=\int_{0}^{\infty}\left[1-F_{\hat{\mathcal{L}}^{D}_{\left(a\right)}|h}\left(l\right)\right]\,f_{\mathcal{L}^{D}_{\left(a,a\right)}|h}\left(l\right)\mathrm{d}l,$
(33)

where the shape parameter of $\mathcal{L}^{D}_{\left(a,a\right)}$ is

$\displaystyle\kappa^{D}_{\left(a,a\right)}$
$\displaystyle=\left|h\sqrt{E_{\text{s}}}\xi^{D}_{\left(a,a\right)}\right|^{2}/2\sigma^{2}$
(34) $\displaystyle=\left|\xi^{D}_{\left(a,a\right)}\right|^{2}\,M\,\gamma.$
(35)

The Rician distributed r.v. $\mathcal{L}^{D}_{\left(a,a\right)}$ can be
approximated as a Gaussian r.v. if the shape parameter
$\kappa^{D}_{\left(a,a\right)}$ is larger than 2 [48]. In our setting, as
shown in (35), $\xi^{D}_{\left(a,a\right)}\approx 1$ and $M=2^{\textrm{SF}}$
ranging from $2^{7}=128$ to $2^{12}=4096$ is a large number so that
$\mathcal{L}^{D}_{\left(a,a\right)}$ can be approximated as a Gaussian r.v. in
an acceptable SNR range. For example, if we consider $\textrm{SF}=8$ and ML
decoder, we have $\xi^{\textrm{ML}}_{\left(a,a\right)}=1$ and $M=256$. The
proper SNR range for the approximation is $\gamma\geq
1/128\approx-21.1~{}\textrm{dB}$, which covers most of the SER range of
interest. Hence, we have

$f_{\mathcal{L}^{D}_{\left(a,a\right)}|h}\left(l\right)\approx\dfrac{1}{\sqrt{2\pi\sigma_{a}^{2}}}\exp\left[-\frac{\left(l-\mu_{a}\right)^{2}}{2\sigma_{a}^{2}}\right],$
(36)

where $\mu_{a}$ and $\sigma_{a}$ are the mean and variance of
$\mathcal{L}^{D}_{\left(a,a\right)}$, respectively, i.e.,

$\displaystyle\mu_{a}$
$\displaystyle=\mathbb{E}\left[\mathcal{L}^{D}_{\left(a,a\right)}\right]=\sigma\sqrt{\pi/2}\cdot
L_{1/2}\left(-\kappa^{D}_{\left(a,a\right)}\right),$ (37)
$\displaystyle\sigma_{a}^{2}$
$\displaystyle=\mathbb{V}\left[\mathcal{L}^{D}_{\left(a,a\right)}\right]=2\sigma^{2}\left(1+\kappa^{D}_{\left(a,a\right)}\right)-\mu_{a}^{2},$
(38)

where $L_{q}\left(\cdot\right)$ denotes a Laguerre polynomial[47, eq:
13.6.27], and for the case $q=1/2$, we have

$L_{1/2}\left(x\right)=e^{x/2}\left[\left(1-x\right)I_{0}\left(-\dfrac{x}{2}\right)-xI_{1}\left(-\dfrac{x}{2}\right)\right].$
(39)

TABLE III: Typical values of $\kappa^{D}_{\left(a,a\right)}$, $\mu_{a}$, and
$\sigma_{a}^{2}$ for $N\in\left\\{2,3,4,5\right\\}$ at $\text{SF}=7$,
$\gamma=-10~{}\textrm{dB}$, and $a=h=E_{s}=1$.

$N$ | $\kappa^{D}_{\left(a,a\right)}$ | $\mu_{a}$ | $\sigma_{a}^{2}$
---|---|---|---
ML | FFT | ML | FFT | ML | FFT
2 | 12.80 | 10.51 | 26.11 | 23.75 | 25.08 | 24.96
3 | 12.80 | 12.27 | 26.11 | 25.58 | 25.08 | 25.05
4 | 12.80 | 12.71 | 26.11 | 26.01 | 25.08 | 25.07
5 | 12.80 | 12.78 | 26.11 | 26.09 | 25.08 | 25.08

In Table. III, we show some typical values of $\kappa^{D}_{\left(a,a\right)}$,
$\mu_{a}$, and $\sigma_{a}^{2}$ as a function of $N$ for both ML and FFT
decoders.

The complexity of evaluating integral (33) mainly comes from the product of
$M-1$ non-identical Rician CDFs. As a possible solution, we can leverage an
asymptotic approximation of the order statistics of a sequence of independent
and non-identically distributed (i.n.i.d.) Rician r.v.s [49] using the extreme
value theorem (EVT). Nevertheless, it shall result in an exponential-in-
exponential expression, rendering difficulties in the integration.

Alternatively, we consider Gauss-Hermite quadrature to evaluate the integral
numerically. In this procedure, the integral is converted to a weighted sum of
function values using the Gauss-Hermite quadrature [47]. Since the Rician PDF
(30) is zero for $l<0$, the lower limit of integral (33) can be substituted
with $-\infty$. Thus, we get

$\displaystyle\begin{aligned}
P^{D}_{\textrm{e}|a}&\overset{}{\approx}\int_{-\infty}^{\infty}\dfrac{\left[1-F_{\hat{\mathcal{L}}^{D}_{\left(a\right)}|h}\left(l\right)\right]}{\sqrt{2\pi\sigma_{a}^{2}}}\,\exp\left[-\frac{\left(l-\mu_{a}\right)^{2}}{2\sigma_{a}^{2}}\right]\mathrm{d}l.\end{aligned}$
(40)

Let us define $\hat{l}\triangleq\left(l-\mu_{a}\right){\sqrt{2}\sigma_{a}}$
and substitute $l$ with $\sqrt{2}\sigma_{a}\hat{l}+\mu_{a}$ so that the
integral becomes

$\displaystyle P^{D}_{\textrm{e}|a}$
$\displaystyle\overset{}{\approx}\int_{-\infty}^{\infty}{\left[1-F_{\hat{\mathcal{L}}^{D}_{\left(a\right)}|h}\left(\sqrt{2}\sigma_{a}\hat{l}+\mu_{a}\right)\right]}\frac{\exp\left[-\hat{l}^{2}\right]}{\sqrt{\pi}}\mathrm{d}\hat{l}$
(41)
$\displaystyle\overset{}{\approx}\dfrac{1}{\sqrt{\pi}}\sum_{t=1}^{N_{\textrm{GH}}}\omega_{t}\left[1-F_{\hat{\mathcal{L}}^{D}_{\left(a\right)}|h}\left(\sqrt{2}\sigma_{a}x_{t}+\mu_{a}\right)\right],$
(42)

where $N_{\textrm{GH}}$ is the number of function samples used to approximate
the integral, $x_{t}$ are the roots of the physicists’ version of the Hermite
polynomial $H_{N_{\textrm{GH}}}\left(x\right)$, and $\omega_{t}$ are the
corresponding weights. The integration formula and corresponding weights can
be found in [47, eq: 25.4.46].

We assume that the transmit symbols are equiprobable so that the average SER
of LB in the AWGN channel can be expressed as

$\displaystyle\overline{P^{D}_{\textrm{e}|h}}$
$\displaystyle=\frac{1}{M}\sum_{a=0}^{M-1}P^{D}_{\textrm{e}|a}$
$\displaystyle=\frac{1}{M\sqrt{\pi}}\sum_{a=0}^{M-1}\sum_{t=1}^{N_{\textrm{GH}}}\omega_{t}\left[1-F_{\hat{\mathcal{L}}^{D}_{\left(a\right)}|h}\left(\sqrt{2}\sigma_{a}x_{t}+\mu_{a}\right)\right].$
(43)

###  IV-B LoRa backscatter SER Performance over Double Nakagami-m Fading
Channels

We consider a double Nakagami-m fading scenario with different shape
parameters $\left\\{m_{1},m_{2}\right\\}$ and spread parameters
$\left\\{\Omega_{1},\Omega_{2}\right\\}$. The spread parameters are related to
the link distance, i.e., $\Omega_{1}\propto 1/d_{1}^{2}$ and
$\Omega_{2}\propto 1/d_{2}^{2}$, where $\left\\{d_{1},d_{2}\right\\}$ are the
distances of the two links. Let us define r.v.s $\left|\mathcal{H}_{1}\right|$
and $\left|\mathcal{H}_{2}\right|$ as the channel amplitudes of Tx-tag and
tag-Rx links, respectively. The PDFs of $\left|\mathcal{H}_{1}\right|$ and
$\left|\mathcal{H}_{2}\right|$ can be expressed as

$\displaystyle
f_{\left|\mathcal{H}_{1}\right|}\left(\left|h_{1}\right|\right)$
$\displaystyle=\dfrac{2m_{1}^{m_{1}}}{\Gamma\left(m_{1}\right)\Omega_{1}^{m_{1}}}|h_{1}|^{2m_{1}-1}\exp\left[-\dfrac{m_{1}}{\Omega_{1}}|h_{1}|^{2}\right],$
(44) $\displaystyle
f_{\left|\mathcal{H}_{2}\right|}\left(\left|h_{2}\right|\right)$
$\displaystyle=\dfrac{2m_{2}^{m_{2}}}{\Gamma\left(m_{2}\right)\Omega_{2}^{m_{2}}}|h_{2}|^{2m_{2}-1}\exp\left[-\dfrac{m_{2}}{\Omega_{1}}|h_{2}|^{2}\right].$
(45)

Let us also define the r.v.
$\left|\mathcal{H}\right|\triangleq\left|\mathcal{H}_{1}\right|\left|\mathcal{H}_{2}\right|$
as the overall channel amplitude, whose PDF is governed by [50]

$f_{\left|\mathcal{H}\right|}\left(\left|h\right|\right)=\dfrac{4\left(r_{1}r_{2}\right)^{v/2}}{\Gamma\left(m_{1}\right)\Gamma\left(m_{2}\right)}|h|^{v-1}K_{n}\left(2\sqrt{r_{1}r_{2}}|h|\right),$
(46)

where $\Gamma\left(\cdot\right)$ denotes the Gamma function[47, eq: 6.1.1],
$K_{n}\left(\cdot\right)$ denotes the modified Bessel function of the second
kind with order $n$[47, eq: 9.6.1], $r_{1}=m_{1}/\Omega_{1}$,
$r_{2}=m_{2}/\Omega_{2}$, $v=m_{1}+m_{2}$, and $n=|m_{1}-m_{2}|$.

In the following, we analyze two cases with different power allocation
schemes, i.e., fixed transmit power and varying transmit power with an average
power constraint.

####  IV-B1 Fixed Transmit Power

Considering the fixed symbol energy, the average SER in fading channels is
calculated as

$\displaystyle\overline{P^{D}_{\textrm{e}}}=\int_{0}^{\infty}\overline{P^{D}_{\textrm{e}|h}}\,f_{\left|\mathcal{H}\right|}\left(\left|h\right|\right)\mathrm{d}|h|.$
(47)

Substituting (43) and (46) into (47), and substituting $\left|h\right|$ with
${\tilde{h}}/{2\sqrt{r_{1}r_{2}}}$, we obtain

$\displaystyle\begin{aligned}
\overline{P^{D}_{\textrm{e}}}&=\dfrac{2^{2-v}\sqrt{1/\pi}}{M\Gamma\left(m_{1}\right)\Gamma\left(m_{2}\right)}\sum_{a=0}^{M-1}\sum_{t=1}^{N_{\textrm{GH}}}\omega_{t}I_{1},\end{aligned}$
(48)

where

$\displaystyle\begin{aligned}
I_{1}={\int_{0}^{\infty}{\left[1-F_{\hat{\mathcal{L}}^{D}_{\left(a\right)}|\tilde{h}}\left(\sqrt{2}\sigma_{a}x_{t}+\mu_{a}\right)\right]}\tilde{h}^{v-1}K_{n}\left(\tilde{h}\right)\mathrm{d}\tilde{h}}.\end{aligned}$
(49)

In the following, we compute $I_{1}$ depending on the value of $n$ since for
certain cases the Bessel function in (49) can be expressed in closed-form
while for other cases it is challenging to find a closed-form expression.

#####  $n$ is half-integer

We consider the order of the Bessel function as a half-integer, i.e.,
$n=u+\frac{1}{2}$ for all $u\in\left\\{0,1,2,\cdots\right\\}$. Thus, the
modified Bessel function of the second kind with order $n$ can be expressed in
closed form[51, eq: 8.468]

$K_{u+\frac{1}{2}}\left(z\right)=\sqrt{\dfrac{\pi}{2z}}\sum_{k=0}^{u}\dfrac{\left(u+k\right)!}{k!\left(u-k\right)!\left(2z\right)^{k}}\,e^{-z}.$
(50)

Substituting (50) into $I_{1}$ yields

$\displaystyle I_{1}$
$\displaystyle=\int_{0}^{\infty}\sqrt{\frac{\pi}{2\tilde{h}}}\,{\left[1-F_{\hat{\mathcal{L}}^{D}_{\left(a\right)}|\tilde{h}}\left(\sqrt{2}\sigma_{a}x_{t}+\mu_{a}\right)\right]}$
$\displaystyle\quad\quad\quad\quad\quad\times\sum_{k=0}^{\left\lfloor
n\right\rfloor}\dfrac{\left(\left\lfloor
n\right\rfloor+k\right)!}{k!\left(\left\lfloor
n\right\rfloor-k\right)!\left(2\tilde{h}\right)^{k}}\tilde{h}^{v-1}e^{-\tilde{h}}\mathrm{d}\tilde{h}.$
(51)

where an a power and an exponential terms $\tilde{h}^{v-1}e^{-\tilde{h}}$ are
included in the integrand. Thus, unlike using Gauss-Hermite quadrature in (42)
that integrates on an integrand with an exponential term $e^{-\hat{l}^{2}}$,
we consider Gauss-Laguerre quadrature to evaluate the integral numerically

$\displaystyle
I_{1}\overset{}{\approx}\sqrt{\frac{\pi}{2}}\,\sum_{e=1}^{N_{\textrm{GL}}}\omega_{e}\dfrac{G_{t}\left(x_{e}\right)}{\sqrt{x_{e}}}\sum_{k=0}^{\left\lfloor
n\right\rfloor}\dfrac{\left(\left\lfloor
n\right\rfloor+k\right)!}{k!\left(\left\lfloor
n\right\rfloor-k\right)!\left(2x_{e}\right)^{k}},$ (52)

where $G_{t}\left(x_{e}\right)$ is given by

$\displaystyle G_{t}\left(x_{e}\right)=1-\prod_{i=0,i\neq
a}^{M-1}\left[1-Q_{1}\left(\dfrac{\left|x_{e}\sqrt{E_{\textrm{s}}}\xi^{D}_{\left(a,i\right)}\right|}{2\sqrt{r_{1}r_{2}}\sigma},\dfrac{\sqrt{2}\sigma_{a}x_{t}+\mu_{a}}{\sigma}\right)\right],$

$N_{\textrm{GL}}$ is the number of function samples, $x_{e}$ are the roots of
the generalized Laguerre polynomial
$L^{\left(v-1\right)}_{N_{\textrm{GL}}}\left(x\right)$, and $\omega_{e}$ are
the corresponding weights[52].

Substituting (52) into (48), we obtain the average SER in double Nakagami-m
fading channel for the half-integer order $n$

$\displaystyle\overline{P^{D}_{\textrm{e}}}=\dfrac{2^{2/3-v}}{M\Gamma\left(m_{1}\right)\Gamma\left(m_{2}\right)}\sum_{a=0}^{M-1}\sum_{t=1}^{N_{\textrm{GH}}}\sum_{e=1}^{N_{\textrm{GL}}}\omega_{t}\omega_{e}\dfrac{G_{t}\left(x_{e}\right)}{\sqrt{x_{e}}}$
$\displaystyle\times\sum_{k=0}^{\left\lfloor
n\right\rfloor}\dfrac{\left(\left\lfloor
n\right\rfloor+k\right)!}{k!\left(\left\lfloor
n\right\rfloor-k\right)!\left(2x_{e}\right)^{k}}.$ (53)

#####  $n$ is not half-integer

If $n$ is not half-integer, it is challenging to obtain a closed-form
expression for the modified Bessel function of the second kind with an order
$n$. Thus, we provide an approximation for the Bessel function using the two
Bessel functions with adjacent half-integer orders. The approximation is tight
according to numerical results. We also assume that $n$, namely, the
difference between the shape parameters of two channels, is larger than $1/2$.
This is because the tag is typically placed closer to the Tx to send signals
to a Rx farther away. The link with a shorter distance, i.e., the Tx-tag link,
is more likely to have a Line-of-Sight (LOS) path, while the tag-Rx link is
more likely to experience a more severe multipath fading effect.

For all $n>{1}/{2}$, there exists only one integer $u$ satisfying the
condition

$u-\frac{1}{2}<n<u+\frac{1}{2}.$ (54)

Both $K_{u-1/2}\left(x\right)$ and $K_{u+1/2}\left(x\right)$ can be written in
form of elementary functions with an exponential term, as shown in (50). Let
us define a closed-form function $K_{u,n}\left(x\right)$ as

$\displaystyle\begin{aligned}
K_{u,n}\left(x\right)=\left[\dfrac{K_{u-1/2}\left(x\right)}{K_{u+1/2}\left(x\right)}\right]^{u-n}\sqrt{K_{u-1/2}\left(x\right)K_{u+1/2}\left(x\right)},\end{aligned}$
(55)

The approximation is provided in [53] as

$\displaystyle K_{n}\left(x\right)$ $\displaystyle\approx
C_{u,n}\,K_{u,n}\left(x\right),$ (56)

where $C_{u,n}$ is a constant given by

$C_{u,n}=\dfrac{\left(u-{1}/{2}\right)^{u-n+{1}/{2}}\Gamma\left(n\right)}{\Gamma\left(u+{1}/{2}\right)}<1.$
(57)

Let us also define the summation term in $K_{u-1/2}\left(x\right)$ and
$K_{u+1/2}\left(x\right)$ as

$\displaystyle\Psi_{u-1/2}\left(x\right)$
$\displaystyle=\sum_{k=0}^{u-1}\dfrac{\left(u-1+k\right)!}{k!\left(u-1-k\right)!\left(2x\right)^{k}},$
(58) $\displaystyle\Psi_{u+1/2}\left(x\right)$
$\displaystyle=\sum_{k=0}^{u}\dfrac{\left(u+k\right)!}{k!\left(u-k\right)!\left(2x\right)^{k}},$
(59)

and $\Psi_{u,n}\left(x\right)$ as

$\displaystyle\begin{aligned}
\Psi_{u,n}\left(x\right)=\left[\dfrac{\Psi_{u-1/2}\left(x\right)}{\Psi_{u+1/2}\left(x\right)}\right]^{u-n}\sqrt{{\Psi_{u+1/2}\left(x\right)\Psi_{u-1/2}\left(x\right)}}.\end{aligned}$
(60)

Substituting (60) into (56) results in

$\displaystyle K_{n}\left(x\right)$
$\displaystyle\approx\sqrt{\frac{\pi}{2x}}C_{u,n}\,\Psi_{u,n}\left(x\right)\,e^{-x},$
(61)

which has an exponential term so that the generalized Gauss-Laguerre
quadrature can be applied to compute $I_{1}$. Following similar steps, the
approximation for the average SER under double Nakagami-m fading channels is
calculated as

$\displaystyle\overline{P^{D}_{\textrm{e}}}$
$\displaystyle=\dfrac{2^{2/3-v}C_{u,n}}{M\Gamma\left(m_{1}\right)\Gamma\left(m_{2}\right)}\sum_{a=0}^{M-1}\sum_{t=1}^{N_{\textrm{GH}}}\sum_{e=1}^{N_{\textrm{GL}}}\frac{\omega_{t}{\omega_{\textrm{e}}}G_{t}\left(x_{\textrm{e}}\right)}{\sqrt{x_{e}}/\Psi_{u,n}\left(x_{e}\right)},$
(62)

####  IV-B2 Limited Average Power

We consider a scenario where the transmit power is under certain constraints.
Also, the symbol energy $E_{\textrm{s}}\left(\gamma\right)$ can be adjusted
according to the channel condition, or similarly, instantaneous SNR. Let us
define the instantaneous SNR as

$\gamma=\dfrac{\overline{E_{\textrm{s}}}|h|^{2}}{2\sigma^{2}M}=\tilde{\gamma}\left|h\right|^{2},$
(63)

where $\overline{E_{\textrm{s}}}$ denotes the average symbol energy and
$\tilde{\gamma}\triangleq{\overline{E_{\textrm{s}}}}/{2\sigma^{2}M}$. The
energy constraint is governed by

$\int_{0}^{\infty}E_{\textrm{s}}\left(\gamma\right)f_{\Gamma}\left(\gamma\right)\mathrm{d}\gamma\leq\overline{E_{\textrm{s}}},$
(64)

where $f_{\Gamma}\left(\gamma\right)$ is the PDF of the instantaneous
SNR.444Note that the randomness in the SNR permits analyzing energy harvesting
schemes where the harvested energy can be considered as a r.v.. It is
determined only by the channel power gain $\left|h\right|^{2}$. Thus, we have

$p_{\Gamma}\left(\gamma\right)=\dfrac{2\left(r_{1}r_{2}/\tilde{\gamma}\right)^{v/2}}{\Gamma\left(m_{1}\right)\Gamma\left(m_{2}\right)}\gamma^{v/2-1}K_{n}\left(2\sqrt{r_{1}r_{2}\gamma/\tilde{\gamma}}\right).$
(65)

We adopt water-filling scheme, the optimal power allocation scheme, to adjust
the transmit power [54]. The scheme is governed by

$\displaystyle\dfrac{E_{\textrm{s}}\left(\gamma\right)}{\overline{E_{\textrm{s}}}}=\left\\{\begin{aligned}
\frac{1}{\gamma_{0}}-\frac{1}{\gamma},\quad\quad&\gamma>\gamma_{0}\\\
\quad\quad 0\quad,\quad\quad&\gamma\leq\gamma_{0}\end{aligned}\right.$ (66)

where the outage SNR $\gamma_{0}$ is found by numerically solving

$\int_{\gamma_{0}}^{\infty}\left(\frac{1}{\gamma_{0}}-\frac{1}{\gamma}\right)p\left(\gamma\right)\mathrm{d}\gamma=1.$
(67)

Substituting (66) into (1) yields the expression of the received signal with
power allocation

${r_{a}}\left[k\right]=\sqrt{\dfrac{\overline{E_{\textrm{s}}}|h|^{2}}{\gamma_{0}}-2\sigma^{2}M}\,{x_{a}}\left[k\right]+{\omega}\left[k\right].$
(68)

Following similar derivation, the average SER in double fading channels with
water-filling scheme is calculated as

$\displaystyle\overline{P^{D}_{\textrm{e}}}=\frac{1}{{\int_{h_{0}}^{\infty}f_{\left|\mathcal{H}\right|}\left(\left|h\right|\right)\mathrm{d}|h|}}\int_{h_{0}}^{\infty}\int_{0}^{\infty}\left[1-F_{\hat{\mathcal{L}}^{D}_{\left(a\right)}|h}\left(l\right)\right]$
$\displaystyle\times
f_{\mathcal{L}^{D}_{\left(a,a\right)}|h}\left(l\right){f_{\left|\mathcal{H}\right|}\left(\left|h\right|\right)}\mathrm{d}l\
\mathrm{d}|h|,$ (69)

where $h_{0}$ is the amplitude of the outage channel amplitude defined as
$h_{0}=\sqrt{{\gamma_{0}}/{\tilde{\gamma}}}$. The integral is divided by
${\int_{h_{0}}^{\infty}f_{\left|\mathcal{H}\right|}\left(\left|h\right|\right)\mathrm{d}|h|}$
since it is a conditional probability that the SNR is higher than the outage.

##  V Numerical Results and Discussion

In this section, we show the performance of LB in terms of the power spectrum
and error performance in AWGN and fading channels.

###  V-A Power Spectrum

$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$5$$0$$-5$$-10$$-15$$-20$$-25$$-30$$f/B$Continuous
Power Spectrum (dB)$\text{SF}=7$, $N=2$$\text{SF}=9$, $N=2$$\text{SF}=9$,
$N=3$ (a)

$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$0$$-20$$-40$$-60$$-80$$-100$$-120$$f/B$Discrete
Power Spectrum (dB)$\text{SF}=7$, $N=2$$\text{SF}=9$, $N=2$$\text{SF}=9$,
$N=3$ (b)

Figure 5: Continuous and discrete spectrum of baseband LoRa backscatter
signals with $\text{SF}\in\left\\{7,9\right\\}$ and
$N\in\left\\{2,3\right\\}$: (a) continuous part of the spectrum, (b) discrete
spectrum.

In Fig. 5, we present the double-sided power spectrum of the baseband LB
signals as a function of normalized frequency $f/B$ with different spreading
factors and the number of phases, i.e., $\text{SF}\in\left\\{7,9\right\\}$ and
$N\in\left\\{2,3\right\\}$. We only show
$\mathcal{G}_{\mathcal{I}}\left(f\right)$ for $f\geq 0$ because
$\mathcal{G}_{\mathcal{I}}\left(f\right)=\mathcal{G}_{\mathcal{I}}\left(-f\right)$.
We show both the normalized power spectral density,
$10\log\left(\mathcal{G}_{\mathcal{I}}^{\left(c\right)}\left(f\right)*B\right)$,
and the discrete part of the spectrum. For the discrete spectrum, we show the
power
$\left|\sum_{a=0}^{M-1}S_{a}\left(\frac{lB}{M}\right)\right|^{2}/{M^{2}T^{2}_{\text{s}}}$
at frequency $lB/M$. We see from the continuous part of the spectrum that it
has a staircase shape where the power spectrum does not always decrease as
frequency increases but maintains relatively stable over some frequency
ranges. Also, for a higher number of phases, i.e., loads, the spectrum drops
faster compared to that with a smaller number of quantization phases. We can
also see that as SF increases, the power spectrum gets increasingly condensed
so that more power of the complex envelope is contained between $-B/2$ and
$B/2$.

$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$10$$0$$-10$$-20$$-30$$-40$$-50$$f/B$Power
Spectrum (dB)Normal LoRaLB, $N=2$, TheoreticalLB, $N=3$, TheoreticalLB, $N=4$,
TheoreticalLB, $N=2$, MCLB, $N=4$, MC Figure 6: Power Spectrum of baseband
LoRa and LoRa backscatter modulated signals at $\text{SF}=9$ for
$N\in\left\\{2,3,4\right\\}$.

In Fig. 6, we compare the derived total power spectrum with both continuous
and discrete parts with that of normal LoRa at $\text{SF}=9$, for various
numbers of quantization phases, i.e., $N\in\left\\{2,3,4\right\\}$. The
frequency range of interest is split into multiple bins of $\Delta f=B/1024$
width. We verified the analytical results using Monte Carlo simulations
generated by estimating the power spectral density in (15) using Welch’s
method. Only the cases for $N\in\left\\{2,4\right\\}$ are included for space
limitation. We see from the figure that the power spectrum of LB modulated
signals has a strong similarity with that of normal LoRa in the range of
$f\in\left[0,B/2\right]$. However, the power spectrum of LB signals begins to
saturate for $f>B/2$ while that of normal LoRa decreases continuously. For an
increasing number of phases, LB matches normal LoRa for a larger frequency
range. For example, at $f/B=1$, the power spectrum of LB with $N=5$ is nearly
equivalent to normal LoRa, while for $N=2$, the gap is approximately over 25
dB.

$867.6$$867.7$$867.8$$867.9$$868$$868.1$$868.2$$868.3$$868.4$$868.5$$868.6$$868.7$$868.8$$868.9$$869$$0$$-10$$-20$$-30$$-40$$-50$$-60$$f$
(MHz)Power Spectrum (dBm/1kHz)Mask$\text{SF}=7$, $N=3$$\text{SF}=7$, $N=4$
Figure 7: One-sided power spectrum for LoRa backscatter passband signals
compared with the mask of the ETSI regulation in the G1 sub-band, for
$\text{SF}=7$, $N\in\left\\{3,4\right\\}$, $B=250~{}\text{kHz}$, RBW=1 kHz and
$P_{s}=14~{}\text{dBm}$.

Also, we investigate the LB spectrum in the ISM band with the ETSI
regulations[41]. In Fig. 7, we report the one-sided power spectrum calculated
with $\text{SF}=7$, the resolution bandwidth $\text{RBW}=1~{}\text{kHz}$ and
transmission power $P_{s}=\text{14~{}dBm}$, i.e., the maximum allowed power.
The spectrum exceeds the spectral mask for $N=3$ while it is within the mask
for $N=4$ at the maximum transmission power. Thus, to satisfy the ETSI
regulation on the ISM band, one should either increase the number of loads
used in the tag or decrease the transmission power. The same approach may be
used to investigate the compliance of LB for different ISM bands, spreading
factors, number of quantization phases, and bandwidths in accordance with
other regional requirements.

###  V-B Error Performance in AWGN Channels

$-18$$-17$$-16$$-15$$-14$$-13$$-12$$-11$$-10$$-9$$-8$$-7$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$SNRSER$\mathrm{SF}=7$
LoRa$\mathrm{SF}=7$ LB-2$\mathrm{SF}=7$ LB-3$\mathrm{SF}=7$
LB-4$\mathrm{SF}=8$ LoRa$\mathrm{SF}=8$ LB-2$\mathrm{SF}=8$
LB-3$\mathrm{SF}=8$ LB-4$\mathrm{SF}=9$ LoRa$\mathrm{SF}=9$
LB-2$\mathrm{SF}=9$ LB-3$\mathrm{SF}=9$ LB-4 Figure 8: The SER performance in
the AWGN channel using ML decoder at $\text{SF}\in\left\\{7,8,9\right\\}$ of
normal LoRa and LB for $N\in\left\\{2,3,4\right\\}$.

$-20$$-19$$-18$$-17$$-16$$-15$$-14$$-13$$-12$$-11$$-10$$-9$$-8$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$SNRSER$\mathrm{SF}=8$,
ML MC$\mathrm{SF}=8$, ML Approx. (43)$\mathrm{SF}=8$, FFT MC$\mathrm{SF}=8$,
FFT Approx. (43)$\mathrm{SF}=9$, ML MC$\mathrm{SF}=9$, ML Approx.
(43)$\mathrm{SF}=9$, FFT MC$\mathrm{SF}=9$, FFT Approx. (43) (a)

$-20$$-19$$-18$$-17$$-16$$-15$$-14$$-13$$-12$$-11$$-10$$-9$$-8$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$SNRSER$\mathrm{SF}=8$,
ML MC$\mathrm{SF}=8$, ML Approx. (43)$\mathrm{SF}=8$, FFT MC$\mathrm{SF}=8$,
FFT Approx. (43)$\mathrm{SF}=9$, ML MC$\mathrm{SF}=9$, ML Approx.
(43)$\mathrm{SF}=9$, FFT MC$\mathrm{SF}=9$, FFT Approx. (43) (b)

Figure 9: LB SER performance in the AWGN channel at SF=8,9: (a) $N=2$, (b)
$N=4$.

In Fig. 8, we compare the SER using ML decoder in the AWGN channel at
$\text{SF}\in\left\\{7,8,9\right\\}$ for normal LoRa and LB with the different
number of quantization phases, where, as a benchmark, the curves for normal
LoRa are calculated by computing integrals numerically:

$\int_{0}^{\infty}\left[1-\left[1-\exp\left[-\frac{l^{2}}{2\sigma^{2}}\right]^{M-1}\right]\right]\,f_{\mathcal{L}_{\left(a,a\right)}^{\text{ML}}}\left(l\right)\mathrm{d}l,$
(70)

and the curves for LB are from Monte Carlo simulation. It can be observed that
the SER for LB using ML decoder has a strong similarity to normal LoRa for
various values of the SFs and number of loads, even though the LB waveforms
are non-orthogonal. We verified the results for small lower values of
$\mathrm{SF}$ and still find the error performance similar to each other. The
results suggest that existing works discussing the error performance of
orthogonal LoRa can be also applied to LB using ML when considering detectors.

Fig. 9 presents comparisons of derived LB SER approximation in (43) with the
theoretical SER calculated by numerically computing the integral in (33) at
$\text{SF}=8,9$ in the AWGN channel. The figures also include the results from
the Monte Carlo simulation. As shown in the figure, the derived expression in
(43) exhibits a tight approximation to the theoretical results. Furthermore,
comparing the SER using ML and FFT decoders, the SER gap between them differs
hugely depending on the number of quantization phases. Improved SER using FFT
decoder is evident when increasing $N$. While the gap is around $1$ dB at SER
of $10^{-3}$ for $\text{SF}=9$ when using only 4 quantization phases, there
are nearly no differences between the performance of two decoders with $16$
quantization phases, indicating that the error performance with $N\geq 4$ is
similar to orthogonal LoRa for both ML and FFT detectors in AWGN channel.

###  V-C Error Performance in Fading Channels With Fixed Transmit Power and
Limited Average Power

$6$$8$$10$$12$$14$$16$$18$$20$$22$$24$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$Tag
moving closer to TxSNRSER$d_{2}=\,\quad d_{1}$, MC$d_{2}=\,\quad d_{1}$,
Approx. $d_{2}=\,\quad d_{1}$, MC$d_{2}=\,\quad d_{1}$, Approx.
$d_{2}=\,\,4\,\,d_{1}$, MC$d_{2}=\,\,4\,\,d_{1}$, Approx.
$d_{2}=\,\,4\,\,d_{1}$, MC$d_{2}=\,\,4\,\,d_{1}$, Approx. $d_{2}=16\,d_{1}$,
MC$d_{2}=16\,d_{1}$, Approx. $d_{2}=16\,d_{1}$, MC$d_{2}=16\,d_{1}$, Approx.
Figure 10: LB SER performance in double Nakagami-m fading channels using ML
decoder for $\textrm{SF}=7$ and $N=2$ (dashed), $N=4$ (solid).
$6$$8$$10$$12$$14$$16$$18$$20$$22$$24$$10^{-3}$$10^{-2}$$10^{-1}$SNRSERSF=7,
$N=2$, MCSF=7, $N=2$, Num (66)SF=7, $N=4$, MCSF=7, $N=4$, Num (66)SF=8, $N=2$,
MCSF=8, $N=2$, Num (66)SF=8, $N=4$, MCSF=8, $N=4$, Num (66) Figure 11: SER
performance for LB with water-filling technique using FFT decoder at
$\text{SF}\in\left\\{7,8\right\\}$ and $N\in\left\\{2,4\right\\}$.

In Fig. 10, we plot the SER performance of LB over double Nakagami-m fading
channels at $\text{SF}=7$ and $N=2$$,4$ where we fix the total distance of the
reflected path between the Tx and Rx through the tag, i.e., $d=d_{1}+d_{2}$.
In addition, we consider a path loss model where the energy is inversely
proportional to the squared distance, i.e.,
$\mathbb{V}\left[\left|h_{1}\right|\right]\varpropto 1/d_{1}^{2}$ and
$\mathbb{V}\left[\left|h_{2}\right|\right]\varpropto 1/d_{2}^{2}$. In this
simulation, we consider the FFT decoder. It is shown that a better SER
performance is achieved when the tag is placed closer to the Tx, which is
similar to the experimental results in [15]. Also, the comparison of the
derived SER approximation (62) to the curves calculated by numerical
integrating (47) and the Monte Carlo simulation presents an accurate
approximation.

On the other hand, for the model with limited average power, we show the error
performance of LB with a water-filling power allocation scheme in Fig. 11
using FFT decoder at $\text{SF}\in\left\\{7,8\right\\}$ and
$N\in\left\\{2,4\right\\}$. The theoretical curves from numerical integrating
(66) are consistent with the Monte Carlo simulations.

##  VI Conclusion

In this paper, we provide the first mathematical expression of LB signals with
a finite number of loads. Based on the expression, we derived the closed-form
expressions for the power spectrum of LB, showing the staircase-shaped
spectrum. To satisfy the wireless transmission regulations, measures such as
increasing the number of phases or decreasing the transmit power may be
needed. An analytical approximation of SER for LB in both AWGN and double
Nakagami-m fading channels using two different decoders is derived. The
results suggest that the SER performance of LB using ML decoder with a small
number of phases is similar to LoRa modulation, while the SER performance
using FFT decoder is worse. On the other hand, SER performance using FFT
decoder will improve with more quantization phases. In the double fading
scenarios, a longer communication range can be achieved by setting the
backscatter tag closer to Tx.

## Appendix A
Proof that maximum likelihood decoder is optimum for LoRa backscatter

The Tx sends the real part of the passband signal with symbol $a$ to the Rx:

$R\left(t\right)=\operatorname{Re}\left[r_{a}\left(t\right)e^{j2\pi
f_{\textrm{c}}t}\right],$ (71)

where
$r_{a}\left(t\right)=\left|h\right|\sqrt{E_{\textrm{s}}}x_{a}\left(t\right)e^{j\theta_{h}}+\omega\left(t\right)$
is the corresponding baseband signal, $\theta_{h}$ is the phase of complex
channel gain $h$, and $f_{\textrm{c}}$ is the carrier frequency. We assume
that the channel amplitude $\left|h\right|$ and the delay $\tau$ are known to
Rx. The symbol decision is based on choosing the maximum among a set of
posteriori probabilities[55]

$\displaystyle\hat{a}$ $\displaystyle=\mathop{\arg\max}\limits_{0\leq a\leq
M-1}P_{a}\,p\left(R\left(t\right)\big{|}x_{a}\left(t\right),\left|h\right|,\tau\right)$
(72) $\displaystyle=\mathop{\arg\max}\limits_{0\leq a\leq
M-1}\int_{0}^{2\pi}p\left(R\left(t\right)\big{|}x_{a}\left(t\right),\left|h\right|,\tau,\theta_{h}\right)p\left(\theta_{h}\right)\mathrm{d}\theta_{h},$

where $P_{a}$ is the probability of transmitting symbol $a$,
$p\left(\theta_{h}\right)$ is the PDF of $\theta_{h}$. We consider
equiprobable transmit symbols and uniformly distributed phases as is typical
for Rayleigh, Rician, and Nakagami fading channels. The probability
$p\left(R\left(t\right)\big{|}x_{a}\left(t\right),\left|h\right|,\tau,\theta_{h}\right)$
is a joint Gaussian PDF

$\displaystyle
p\left(R\left(t\right)\big{|}x_{a}\left(t\right),\left|h\right|,\tau,\theta_{h}\right)=$
$\displaystyle\quad\quad\quad\quad
K\exp\left(\operatorname{Re}\left\\{\frac{\left|h\right|}{N_{0}}e^{-j\theta_{h}}y_{a}\left(\tau\right)\right\\}-\frac{\left|h\right|^{2}E_{a}}{N_{0}}\right),$
(73)

where $K$ is an integral constant and

$\displaystyle\begin{aligned}
y_{a}\left(\tau\right)=\int_{0}^{T_{s}}r_{a}\left(t+\tau\right)x_{a}^{*}\left(t\right)\mathrm{d}t\end{aligned}$
(74)

is the complex correlation of the received signal and the signal waveform of
symbol $a$ and

$\displaystyle\begin{aligned}
E_{a}=\frac{1}{2}\int_{0}^{T_{s}}\left|x_{a}\left(t\right)\right|^{2}\mathrm{d}t\end{aligned}$
(75)

is the energy of the waveform $x_{a}\left(t\right)$. Thus, we obtain

$\displaystyle
p\left(R\left(t\right)\big{|}x_{a}\left(t\right),\left|h\right|,\tau\right)=$
(76)
$\displaystyle\exp\left(-\frac{\left|h\right|^{2}E_{a}}{N_{0}}\right)\frac{K}{2\pi}\int_{0}^{2\pi}\exp\left(\operatorname{Re}\left\\{\frac{\left|h\right|}{N_{0}}e^{-j\theta_{h}}y_{a}\left(\tau\right)\right\\}\right)\mathrm{d}\theta_{h}.$

Substitute (76) into (72) and remove the terms that are independent of $a$, we
obtain

$\displaystyle\hat{a}$ $\displaystyle=\mathop{\arg\max}\limits_{0\leq a\leq
M-1}\int_{0}^{2\pi}\exp\left(\operatorname{Re}\left\\{e^{-j\theta_{h}}y_{a}\left(\tau\right)\right\\}\right)\mathrm{d}\theta_{h}$
$\displaystyle=\mathop{\arg\max}\limits_{0\leq a\leq
M-1}\int_{0}^{2\pi}\exp\left(\left|y_{a}\left(\tau\right)\right|\cos\left[\theta_{h}-\arg\left(y_{a}\left(\tau\right)\right)\right]\right)\mathrm{d}\theta_{h}$
$\displaystyle=\mathop{\arg\max}\limits_{0\leq a\leq
M-1}I_{0}\left(\left|y_{a}\left(\tau\right)\right|\right),$ (77)

where the term $E_{a}$ is independent of $a$ because the waveforms of
different symbols have the same energy. Since $I_{0}$ is a monotonically
increasing function, the optimum decision is made by choosing the symbol that
has the maximum cross-correlation with the received signal.

## Appendix B
How to obtain $\left\\{t_{m}\right\\}_{m=1}^{\psi}$

$t_{m}$ are times when the value of
$\widetilde{Q}_{N}\left[f\left(t\right)\right]$ make discrete changes. The set
of $t_{m}$ can be obtained by solving

$f\left(t\right)=\frac{i\pi}{2^{N-1}},\quad i\in\mathbb{Z},\quad
t\in\left[0,M/B\right],$ (78)

where the two parts of $f\left(t\right)$ are quadratic functions so that the
roots can be easily obtained. It is feasible to directly solve
$f_{1}\left(t\right)={i\pi}/{2^{N-1}}$ and
$f_{2}\left(t\right)={i\pi}/{2^{N-1}}$ individually in their respective
regions. However, it is worth noting that we have the equality

$f_{2}\left(t\right)=f_{1}\left(t-\frac{M}{B}\right).$ (79)

We see that $f_{2}\left(t\right)$ is a shifted version of
$f_{1}\left(t\right)$, which can be leveraged to simplify the solution. More
precisely, the roots of $f_{2}\left(t\right)={i\pi}/{2^{N-1}}$ in the region
$\frac{M-a}{B}\leq t\leq\frac{M}{B}$ are corresponding to the roots of
$f_{1}\left(t\right)={i\pi}/{2^{N-1}}$ in region $\frac{-a}{B}\leq t\leq 0$.
Hence, all the roots can be obtained by calculating the roots of
$f_{1}\left(t\right)={i\pi}/{2^{N-1}}$ for
$t\in\left[-\frac{a}{B},\frac{M-a}{B}\right]$ using root formula for quadratic
equation

$\displaystyle
z_{i}^{\pm}=\frac{\left(M-2a\right)\pm\sqrt{\left(M-2a\right)^{2}+iM2^{3-N}}}{2B},$
$\displaystyle\lceil-\frac{\left(M-2a\right)^{2}}{2^{3-n}M}\rceil\leq
i\leq\lfloor\frac{a\left(M-a\right)}{2^{1-n}M}\rfloor.$ (80)

Next, add $M/B$ to $z_{i}^{\pm}$ if they are negative. Also, $t=M/B$ should be
added to the roots set because it overlaps with root $t=0$ before being
shifted. Finally, rearrange all roots from smallest to largest, namely,
$t_{1}\leq t_{2}\leq\cdots\leq t_{\psi}$.

## Appendix C
Proof that $\mathcal{L}^{D}_{\left(a,i\right)}$ follow Rician distribution
with shape parameter
$\kappa^{D}_{\left(a,i\right)}=\frac{\left|h\sqrt{E_{\textrm{s}}}\xi^{D}_{\left(a,i\right)}\right|^{2}}{2\sigma^{2}}$

The Rician distribution is the probability distribution of the magnitude of a
circularly-symmetric bivariate normal random variable. A r.v.
$\mathcal{L}=\sqrt{\mathcal{X}^{2}+\mathcal{Y}^{2}}$ follows the Rician
distribution with shape parameter $\kappa={\left|v\right|^{2}}/{2\sigma^{2}}$
if $\mathcal{X}\sim\mathcal{N}\left(v\sin{\theta},\sigma^{2}\right)$ and
$\mathcal{Y}\sim\mathcal{N}\left(v\cos{\theta},\sigma^{2}\right)$, where
$\mathcal{X}$ and $\mathcal{Y}$ are independently distributed normal r.v.s and
$\theta=v/\left|v\right|$. $\mathcal{L}^{D}_{\left(a,i\right)}$ is the
amplitude of a complex r.v.

$\displaystyle\mathcal{L}^{D}_{\left(a,i\right)}$
$\displaystyle=\left|h\sqrt{E_{\textrm{s}}}\xi^{D}_{\left(a,i\right)}+\mathcal{W}^{D}_{\left(a,i\right)}\right|$
(81)
$\displaystyle=\left|\left|h\sqrt{E_{\text{s}}}\xi_{\left(a,i\right)}^{D}\right|e^{j\beta}e^{-j\beta}+\mathcal{W}^{D}_{\left(a,i\right)}e^{-j\beta}\right|,$
(82)

where $e^{j\beta}$ is the angle of
$h\sqrt{E_{\textrm{s}}}\xi^{D}_{\left(a,i\right)}$, and the second term is the
rotated version of complex Gaussian r.v.. Thus, we define

$\displaystyle\mathcal{X}$
$\displaystyle=\left|h\sqrt{E_{\textrm{s}}}\xi^{D}_{\left(a,i\right)}\right|+\operatorname{Re}\left[\mathcal{W}^{D}_{\left(a,i\right)}e^{-j\beta}\right],$
(83) $\displaystyle\mathcal{Y}$
$\displaystyle=\operatorname{Im}\left[\mathcal{W}^{D}_{\left(a,i\right)}e^{-j\beta}\right],$
(84)

where
$\mathcal{L}^{D}_{\left(a,i\right)}=\left|\mathcal{X}+j\mathcal{Y}\right|$.
Since
$\mathcal{X}\sim\mathcal{N}\left(\left|h\sqrt{E_{\textrm{s}}}\xi^{D}_{\left(a,i\right)}\right|,\sigma^{2}\right)$
and $\mathcal{Y}\sim\mathcal{N}\left(0,\sigma^{2}\right)$,
$\mathcal{L}^{D}_{\left(a,i\right)}$ follows Rician distribution with shape
parameter
$\kappa^{D}_{\left(a,i\right)}=\frac{\left|h\sqrt{E_{\textrm{s}}}\xi^{D}_{\left(a,i\right)}\right|^{2}}{2\sigma^{2}}$.

## References

  * [1] F. Boccardi, R. W. Heath, A. Lozano, T. L. Marzetta, and P. Popovski, “Five disruptive technology directions for 5G,” _IEEE Communications Magazine_ , vol. 52, no. 2, pp. 74–80, Feb. 2014.
  * [2] S. Dang, O. Amin, B. Shihada, and M.-S. Alouini, “What should 6G be?” _Nature Electronics_ , vol. 3, no. 1, pp. 20–29, 2020.
  * [3] L. S. Vailshery, “Iot connected devices worldwide 2019-2030,” Jun 2022. [Online]. Available: https://www.statista.com/statistics/1183457/iot-connected-devices-worldwide/
  * [4] M. Chiani and A. Elzanaty, “On the LoRa modulation for IoT: Waveform properties and spectral analysis,” _IEEE Internet of Things Journal_ , vol. 6, no. 5, pp. 8463–8470, Oct. 2019.
  * [5] U. Raza, P. Kulkarni, and M. Sooriyabandara, “Low power wide area networks: An overview,” _IEEE Communications Surveys & Tutorials_, vol. 19, no. 2, pp. 855–873, Secondquarter 2017.
  * [6] D. Pianini, A. Elzanaty, A. Giorgetti, and M. Chiani, “Emerging distributed programming paradigm for cyber-physical systems over LoRAWANs,” in _2018 IEEE Globecom Workshops (GC Wkshps)_ , Dec. 2018, pp. 1–6.
  * [7] S. Sobhi, A. Elzanaty, M. Y. Selim, A. M. Ghuniem, and M. F. Abdelkader, “Mobility of LoRaWAN gateways for efficient environmental monitoring in pristine sites,” _Sensors_ , vol. 23, no. 3, 2023. [Online]. Available: https://www.mdpi.com/1424-8220/23/3/1698
  * [8] F. Rezaei, C. Tellambura, and S. Herath, “Large-scale wireless-powered networks with backscatter communications—a comprehensive survey,” _IEEE Open Journal of the Communications Society_ , vol. 1, pp. 1100–1130, July 2020.
  * [9] C. Song, Y. Ding, A. Eid, J. G. D. Hester, X. He, R. Bahr, A. Georgiadis, G. Goussetis, and M. M. Tentzeris, “Advances in wirelessly powered backscatter communications: From antenna/RF circuitry design to printed flexible electronics,” _Proceedings of the IEEE_ , vol. 110, no. 1, pp. 171–192, Jan. 2022.
  * [10] N. Van Huynh, D. T. Hoang, X. Lu, D. Niyato, P. Wang, and D. I. Kim, “Ambient backscatter communications: A contemporary survey,” _IEEE Communications Surveys & Tutorials_, vol. 20, no. 4, pp. 2889–2922, May 2018.
  * [11] G. Wang, F. Gao, R. Fan, and C. Tellambura, “Ambient backscatter communication systems: Detection and performance analysis,” _IEEE Transactions on Communications_ , vol. 64, no. 11, pp. 4836–4846, Nov. 2016.
  * [12] V. Liu, A. Parks, V. Talla, S. Gollakota, D. Wetherall, and J. R. Smith, “Ambient backscatter: Wireless communication out of thin air,” _ACM SIGCOMM Computer Communication Review_ , vol. 43, no. 4, pp. 39–50, Oct. 2013\.
  * [13] Y. H. Al-Badarneh, A. Elzanaty, and M.-S. Alouini, “On the performance of spectrum-sharing backscatter communication systems,” _IEEE Internet of Things Journal_ , vol. 9, no. 3, pp. 1951–1961, Feb. 2022.
  * [14] V. Iyer, V. Talla, B. Kellogg, S. Gollakota, and J. Smith, “Inter-technology backscatter: Towards internet connectivity for implanted devices,” in _Proceedings of the 2016 ACM SIGCOMM Conference_ , ser. SIGCOMM ’16. New York, NY, USA: Association for Computing Machinery, Aug. 2016, p. 356–369.
  * [15] V. Talla, M. Hessar, B. Kellogg, A. Najafi, J. R. Smith, and S. Gollakota, “LoRa backscatter: Enabling the vision of ubiquitous connectivity,” _Proceedings of the ACM on Interactive, Mobile, Wearable and Ubiquitous Technologies_ , vol. 1, no. 3, pp. 1–24, Sept. 2017.
  * [16] X. Tang, G. Xie, and Y. Cui, “Self-sustainable long-range backscattering communication using RF energy harvesting,” _IEEE Internet of Things Journal_ , vol. 8, no. 17, pp. 13 737–13 749, Sept. 2021.
  * [17] M. Katanbaf, A. Weinand, and V. Talla, “Simplifying backscatter deployment: Full-duplex LoRa backscatter,” in _18th USENIX Symposium on Networked Systems Design and Implementation (NSDI 21)_ , 2021, pp. 955–972.
  * [18] H. Li, X. Tong, Q. Li, and X. Tian, “XORLoRa: LoRa backscatter communication with commodity devices,” in _2020 IEEE 6th International Conference on Computer and Communications (ICCC)_ , Dec. 2020, pp. 706–711.
  * [19] Y. Peng, L. Shangguan, Y. Hu, Y. Qian, X. Lin, X. Chen, D. Fang, and K. Jamieson, “Plora: A passive long-range data network from ambient LoRa transmissions,” in _Proceedings of the 2018 Conference of the ACM Special Interest Group on Data Communication_ , 2018, pp. 147–160.
  * [20] G. Huang, P. Yang, H. Zhou, Y. Yan, X. He, and X. Li, “Freeback: Blind and distributed rate adaptation in LoRa-based backscatter networks,” in _2021 IEEE Wireless Communications and Networking Conference (WCNC)_ , Mar.-Apr. 2021, pp. 1–6.
  * [21] X. Guo, L. Shangguan, Y. He, J. Zhang, H. Jiang, A. A. Siddiqi, and Y. Liu, “Efficient ambient LoRa backscatter with on-off keying modulation,” _IEEE/ACM Transactions on Networking_ , pp. 1–14, Nov. 2021.
  * [22] M. Lazaro, A. Lazaro, and R. Villarino, “Feasibility of backscatter communication using LoRAWAN signals for deep implanted devices and wearable applications,” _Sensors_ , vol. 20, no. 21, p. 6342, Nov. 2020.
  * [23] N. Fasarakis-Hilliard, P. N. Alevizos, and A. Bletsas, “Coherent detection and channel coding for bistatic scatter radio sensor networking,” _IEEE Transactions on Communications_ , vol. 63, no. 5, pp. 1798–1810, May 2015.
  * [24] J. Qian, F. Gao, G. Wang, S. Jin, and H. Zhu, “Noncoherent detections for ambient backscatter system,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 3, pp. 1412–1422, march 2017.
  * [25] J. K. Devineni and H. S. Dhillon, “Ambient backscatter systems: Exact average bit error rate under fading channels,” _IEEE Transactions on Green Communications and Networking_ , vol. 3, no. 1, pp. 11–25, March 2019.
  * [26] J. Qian, Y. Zhu, C. He, F. Gao, and S. Jin, “Achievable rate and capacity analysis for ambient backscatter communications,” _IEEE Transactions on Communications_ , vol. 67, no. 9, pp. 6299–6310, Sept. 2019.
  * [27] G. Ferré and A. Giremus, “LoRa physical layer principle and performance analysis,” in _2018 25th IEEE International Conference on Electronics, Circuits and Systems (ICECS)_ , Dec. 2018, pp. 65–68.
  * [28] J. G. Proakis and M. Salehi, _Digital Communications_ , 5th ed. New York, NY, USA: McGraw-Hill, 2007.
  * [29] C. Ferreira Dias, E. Rodrigues de Lima, and G. Fraidenraich, “Bit error rate closed-form expressions for LoRa systems under Nakagami and Rice fading channels,” _Sensors_ , vol. 19, no. 20, p. 4412, Oct. 2019.
  * [30] V. Savaux and G. Ferré, “Simple asymptotic BER expressions for LoRa system over Rice and Rayleigh channels,” in _2021 Wireless Telecommunications Symposium (WTS)_ , Apr. 2021, pp. 1–4.
  * [31] T. Elshabrawy and J. Robert, “Closed-form approximation of LoRa modulation BER performance,” _IEEE Communications Letters_ , vol. 22, no. 9, pp. 1778–1781, Sept. 2018.
  * [32] J. Courjault, B. Vrigneau, O. Berder, and M. R. Bhatnagar, “A computable form for LoRa performance estimation: Application to Ricean and Nakagami fading,” _IEEE Access_ , vol. 9, pp. 81 601–81 611, Apr. 2021.
  * [33] G. Baruffa, L. Rugini, L. Germani, and F. Frescura, “Error probability performance of chirp modulation in uncoded and coded lora systems,” _Digital Signal Processing_ , vol. 106, p. 102828, Nov. 2020.
  * [34] S. An, H. Wang, Y. Sun, Z. Lu, and Q. Yu, “Time domain multiplexed lora modulation waveform design for iot communication,” _IEEE Communications Letters_ , vol. 26, no. 4, pp. 838–842, April 2022.
  * [35] M. J. Faber, K. M. van der Zwaag, W. G. V. dos Santos, H. R. d. O. Rocha, M. E. V. Segatto, and J. A. L. Silva, “A theoretical and experimental evaluation on the performance of LoRa technology,” _IEEE Sensors Journal_ , vol. 20, no. 16, pp. 9480–9489, Aug. 2020.
  * [36] O. Afisiadis, M. Cotting, A. Burg, and A. Balatsoukas-Stimming, “On the error rate of the LoRa modulation with interference,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 2, pp. 1292–1304, Feb. 2020.
  * [37] T. Elshabrawy and J. Robert, “Analysis of BER and coverage performance of LoRa modulation under same spreading factor interference,” in _2018 IEEE 29th Annual International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC)_ , Sept. 2018, pp. 1–6.
  * [38] D. Croce, M. Gucciardo, S. Mangione, G. Santaromita, and I. Tinnirello, “Impact of LoRa imperfect orthogonality: Analysis of link-level performance,” _IEEE Communications Letters_ , vol. 22, no. 4, pp. 796–799, Apr. 2018.
  * [39] F. Benkhelifa, Y. Bouazizi, and J. A. McCann, “How orthogonal is LoRa modulation?” _IEEE Internet of Things Journal_ , vol. 9, no. 20, pp. 19 928–19 944, May 2022.
  * [40] A. S. Ali, L. Bariah, S. Muhaidat, P. Sofotasios, and M. Al-Qutayri, “Performance analysis of LoRa-backscatter communication system in awgn channels,” in _2021 4th International Conference on Advanced Communication Technologies and Networking (CommNet)_ , Dec. 2021, pp. 1–5.
  * [41] “Electromagnetic compatibility and radio spectrum matters (ERM); shortrange devices (SRD); radio equipment to be used in the 25 mhz to1000 mhz frequency range with power levels ranging up to 500 mw;part 1: Technical characteristics and test methods,” Eur. Telecommun.Stand. Inst., Sophia Antipolis, France, Standard, 2012.
  * [42] J. W. Harris and H. Stöcker, _Handbook of mathematics and computational science_. Springer Science & Business Media, 1998.
  * [43] T. Elshabrawy, P. Edward, M. Ashour, and J. Robert, “On the different mathematical realizations for the digital synthesis of lora-based modulation,” in _European Wireless 2019; 25th European Wireless Conference_ , May 2019, pp. 1–6.
  * [44] E. A. Lee and D. G. Messerschmitt, _Digital communication_. Springer Science & Business Media, 2012.
  * [45] S. Benedetto and E. Biglieri, _Principles of digital transmission: with wireless applications_. Springer Science & Business Media, 1999.
  * [46] M. Simon and M. Alouini, “A unified approach to the performance analysis of digital communication over generalized fading channels,” _Proceedings of the IEEE_ , vol. 86, no. 9, pp. 1860–1877, 1998.
  * [47] M. Abramowitz and I. A. Stegun, _Handbook of mathematical functions with formulas, graphs, and mathematical tables_. US Government printing office, 1964, vol. 55.
  * [48] H. Gudbjartsson and S. Patz, “The Rician distribution of noisy MRI data,” _Magnetic resonance in medicine_ , vol. 34, no. 6, pp. 910–914, 1995.
  * [49] A. Subhash, S. Kalyani, Y. H. Al-Badarneh, and M.-S. Alouini, “On the asymptotic performance analysis of the k-th best link selection over non-identical non-central chi-square fading channels,” _IEEE Transactions on Communications_ , vol. 70, no. 11, pp. 7191–7206, Nov. 2022.
  * [50] G. K. Karagiannidis, N. C. Sagias, and P. T. Mathiopoulos, “N$\ast$Nakagami: A novel stochastic model for cascaded fading channels,” _IEEE Transactions on Communications_ , vol. 55, no. 8, pp. 1453–1458, Aug. 2007.
  * [51] I. S. Gradshteyn and I. M. Ryzhik, _Table of integrals, series, and products_. Academic press, 2014.
  * [52] P. Rabinowitz and G. Weiss, “Tables of abscissas and weights for numerical evaluation of integrals of the form$\int_{0}^{\infty}e^{-x}x^{n}f\left(x\right)\mathrm{d}x$,” _Mathematical Tables and Other Aids to Computation_ , pp. 285–294, 1959.
  * [53] Z.-H. Yang and S.-Z. Zheng, “The monotonicity and convexity for the ratios of modified Bessel functions of the second kind and applications,” _Proceedings of the American Mathematical Society_ , vol. 145, no. 7, pp. 2943–2958, Jan. 2017.
  * [54] A. Goldsmith, _Wireless communications_. Cambridge university press, 2005.
  * [55] M. K. Simon and M.-S. Alouini, _Digital communication over fading channels_. John Wiley & Sons, 2005, vol. 95.

  *[LB]: LoRa backscatter
  *[LPWAN]: low power wide area networks
  *[SER]: symbol error rate
  *[ML]: maximum likelihood
  *[FFT]: fast Fourier transform
  *[AWGN]: additive white Gaussian noise
  *[ETSI]: European Telecommunications Standards Institute
  *[IoT]: Internet of Things
  *[5G]: 5th-generation cellular networks
  *[eMBB]: enhanced mobile broadband
  *[URLLC]: ultra reliable low latency communications
  *[mMTC]: massive machine type communications
  *[LoRaWAN]: low power long range wide area network
  *[NB-IoT]: narrowband Internet of things
  *[BC]: backscatter communication
  *[Tx]: transmitter
  *[Rx]: receiver
  *[CSS]: chirp spread spectrum
  *[Txs]: transmitter
  *[PDF]: probability density function
  *[RFEH]: RF energy harvesting
  *[OOK]: ON-OFF Keying
  *[FSK]: frequency shift keying
  *[CSI]: channel state information
  *[BER]: bit error rate
  *[BPSK]: binary phase shift keying
  *[SNR]: signal-to-noise ratio
  *[CDF]: cumulative distribution function
  *[DFT]: discrete Fourier transform
  *[r.v.]: random variable
  *[r.v.s]: random variable
  *[SEP]: symbol error probability
  *[EVT]: extreme value theorem
  *[LOS]: Line-of-Sight