1

Inference from Randomized Transmissions by
Many Backscatter Sensors
Guangxu Zhu, Seung-Woo Ko and Kaibin Huang

arXiv:1707.08535v2 [cs.IT] 10 Jan 2018

Abstract
Attaining the vision of Smart Cities requires the deployment of an enormous number of sensors
for monitoring various conditions of the environment ranging from air quality to traffic. Backscatter
sensors have emerged to be a promising solution for two reasons. First, transmissions by backscattering
allow sensors to be powered wirelessly by radio-frequency (RF) waves, overcoming the difficulty in
battery recharging for billions of sensors. Second, the simple backscatter hardware leads to low-cost
sensors suitable for large-scale deployment. On the other hand, backscatter sensors with limited signalprocessing capabilities are unable to support conventional algorithms for multiple access and channel
training. Thus, the key challenge in designing backscatter sensor networks is to enable readers to
accurately detect sensing values given simple ALOHA random access, primitive transmission schemes,
and no knowledge of channel states and statistics. We tackle this challenge by proposing the novel
framework of backscatter sensing (BackSense) featuring random encoding at backscatter sensors and
statistical inference at readers. Specifically, assuming the widely used on/off keying for backscatter
transmissions, the practical random-encoding scheme causes the on/off transmission of a sensor to be
randomized and follow a distribution parameterized by the sensing values. Facilitated by the scheme,
statistical inference algorithms are designed to enable a reader to infer sensing values from randomized
transmissions by multiple backscatter sensors. The specific design procedure involves the construction
of Bayesian networks, namely deriving conditional distributions for relating unknown parameters and
variables (including sensing values, noise power, sensing measurements, number of active sensors) to
signals observed by the reader. Then based on the Bayesian networks and the well-known expectationmaximization (EM) principle, inference algorithms are derived to recover sensing values. Simulation of
the BackSense system demonstrates high accuracy in reader inference despite the mentioned limitations
of backscatter sensors, which grows with increasing numbers of received symbols and reader antennas.

I. I NTRODUCTION
Realizing the visions of Internet-of-Things (IoT) and Smart Cities will require the deployment
of billions of wireless IoT sensors in our society for automating a wide range of applications such
G. Zhu, S.-W. Ko and K. Huang are with the Department of Electrical and Electronic Engineering, The University of Hong
Kong, Pok Fu Lam, Hong Kong (e-mail: gxzhu@eee.hku.hk, swko@eee.hku.hk, haungkb@eee.hku.hk).

2

UAV mounted reader
Continuous sinusoidal wave (CW)
Backscatter transmission

Power beacon (PB)

Fig. 1.

Cluster of
sensors

A BackSense system comprising a cluster of sensors deployed at a remote or hazardous location, a PB wirelessly

powering the sensors, and a UAV mounted reader for collecting sensing data.

as health care, smart homes, pollution and traffic monitoring, and industrial control. Backscatter
transmission has emerged to be a promising solution for tackling some key challenges on
deploying large-scale sensor networks. The main feature of a backscatter radio is to transmit by
backscattering and modulating an incident radio-frequency (RF) wave [1]. This allows backscatter
sensors to be powered wirelessly by a reader, overcoming the difficulty of battery recharging
for a massive number of sensors. Furthermore, requiring no oscillators and RF components,
backscatter sensors can be manufactured to have a small form factor and low cost. Thereby,
they are suitable for large-scale deployment. Last, with recent advancements in multi-antenna
beamforming and low-power electronics, the ranges for backscatter links have been increased
from several meters in the classic RFID applications to tens of meters in state-of-the-art systems
[2], [3]. This makes it possible to collect measurement data from backscatter sensors using
versatile mobile readers mounted on vehicles and unmanned aerial vehicles (UAVs) (see Fig. 1).
Motivated by the potential of backscatter sensor networks, this work presents a novel design
framework, called randomized backscatter sensing (BackSense), for sensing-data uploading and
detection based on random access and machine learning. To be specific, to resolve transmission
collisions, statistical-inference algorithms are designed for readers to directly infer sensing values
from collided signals by exploiting their spatial correlation and partial information of their
distributions. The algorithms build on a proposed scheme of randomized on/off transmissions
targeting backscatter sensors.

3

A. Backscatter Communications and Networking
Backscatter communications are expected to play a key role in next-generation sensor systems
e.g., IoT and Smart Cities. These applications motivate growing research on developing advanced
techniques for backscatter communications that are more complex than those for the traditional
RFID applications [1], [4]–[8]. The modulation and space-time coding schemes were designed
in [1] for enhancing the rate and reliability of backscatter communication links. Traditionally,
a reader transmits a single sinusoidal wave for powering a backscatter tag but such a design is
sub-optimal in terms of energy harvesting efficiency. Instead, a waveform design superimposing
multiple sinusoidal waves was proposed in [4] for enhancing the efficiency, via optimization
based on the harvester-response function. For backscatter tags with low-complexity, non-coherent
detection schemes are preferred due to the simple transmitter architecture. Specific schemes
based on on/off keying or frequency-shift keying were proposed in [5] and [6], respectively. The
channel-utilization efficiency for backscatter communications tends to be low as the primitive
backscatter architecture is incapable of supporting sophisticated transmission designs such as
high-order modulation or spatial multiplexing. An attempt was made in [7] on improving the
efficiency by designing full-duplex backscatter communication networks. In this design, fullduplexing is realized by double (coherent and non-coherent) modulating backscattered signals
and using time-hopping spread spectrum to mitigate interference. Last, to extend the ranges as
well as simultaneoully power a large number of sensors, a novel architecture for backscatter
communication networks was recently proposed in [8] where tags are powered by distributed
power beacons (PBs) instead of readers. Then the tradeoff between network coverage and PB
density was studied therein using a stochastic-geometric network model.
A key design issue for backscatter sensor networks is multiple access by sensors. The natural
design approach is to adopt classic multiple-access schemes including random access [9], [10],
time division multiple access (TDMA) [11], [12], space-division multiple access (SDMA) [13],
[14], and code division multiple access (CDMA) [15], [16]. As mentioned, for large-scale sensors
networks, these schemes have the drawbacks of excessive network overhead due to protocols
and algorithms for orthogonalizing transmissions as well as long latency and low spectrumutilization efficiency. Moreover, the required signal processing (e.g., FFT, channel estimation
and feedback) are impractical for backscatter sensors with a primitive hardware architecture.
On the other hand, ALOHA-type random access does not have the above drawbacks but suffers
from performance degradation caused by transmission collisions. A novel approach for resolving
collisions is proposed in [16]. Specifically, assuming sparse random access by backscatter tags,

4

collided signals are treated as a sparse code and recovered using the techniques of compressive
sensing and rateless coding. The design is effective only in the scenario of dense but mostly
inactive sensors and thus cannot support the deployment of dense sensor networks with many
active sensors.
B. Machine Learning for Wireless Sensor Networks
Machine learning provides a rich set of techniques for learning and prediction of data. Recent
breakthroughs in the field especially the area of artificial intelligence motivated researchers to
apply relevant techniques for bringing intelligence to communication systems and networks (e.g.,
UAV assisted networking [17] and molecular communication [18]). In particular, a thrust of the
research focuses on applying machine-learning techniques to streamline various operations of
sensors networks such as routing [19], abnormality detection [20], positioning and network
topology discovery [21], and sensing-data compression and feature extraction [22].
The existing research most relevant to the current work is the application of statisticallearning techniques to detection of sensing values from signals transmitted by wireless sensors
based on random access [23]–[27]. In [23], [24], the sensors are designed such that their noisy
measurements received by the reader provide sufficient statistics of the measured environmental
parameter, thereby allowing it to be inferred using the classic maximum-likelihood algorithm.
This requires each sensor to estimate and inverse its channel to the reader and furthermore
map the noisy measurement to a pre-determined set of orthogonal waveforms. Such operations,
however, are too complex for a backscatter sensor that usually supports only primitive signal
processing and modulates signals using the simple on/off keying. In a series of other research,
advanced machine-learning techniques were applied to enhance the performance of the carriersensing multiple access (CSMA) scheme [25]–[27]. In [25], cognitive transmitting sensors rely
on statistical inference to infer the channel status (idle or busy) at the receiver based on local
measurements. Such information allows the sensors to intelligently switch between active and
sleep modes to avoid collision and reduce energy consumption. A different kind of cognitive
sensor is designed in [26] based on reinforcement learning to have the capability of adapting its
duty cycle according to the traffic load and channel state, thereby reducing its energy consumption
as well as enhancing network throughput. Another technique, namely neural network, is applied
in [27] to learn the optimal scheduling policy that minimizes a schedule cycle such that all
sensors can be scheduled within each cycle without any collision. The existing designs of
cognitive transmitters are impractical for backscatter sensors again due to their low complexity

5

and limited signal-processing capability. The more practical approach, especially for largescale sensor networks, should be one reducing sensor complexity and compensating for it by
centralized machine learning at readers or in the cloud.
C. Contributions and Organization
Based on this approach, we design techniques for efficient implementation of backscatter
sensor networks. Specifically, adopting the network architecture proposed in [8], we consider
a system illustrated in Fig. 1 where a cluster of backscatter sensors are powered wirelessly by
a PB for transmission to a mobile reader (also called an information collector) mounted on
e.g., a UAV. The BackSense framework for designing such a system satisfies several practical
constraints. First, the transmission of each sensor is based on on/off keying. Second, multiple
access by sensors is based on the simple ALOHA-type random access (without carrier sensing
and scheduling). Last, both sensors and reader have neither channel state information (CSI)
nor parameters of channel distribution. The only knowledge the reader has about the channels is
their distribution type. Under these constraints, the BackSense framework is designed for efficient
sensing-data uploading and accurate detection. To the best of the authors’ knowledge, this current
work represents the first attempt on applying machine learning to the design of backscatter-sensor
systems. The resultant BackSense framework comprises two key components: one is randomized
transmissions by sensors and the other statistical inference at the reader, described as follows.
Randomized transmission for a backscatter sensor: Under the constraint of on/off keying,
a novel random-encoding scheme for sensor transmission is proposed for embedding the sensing
value into the distribution of backscattered signals. Mathematically, the sensor state (backscatter
or not) is distributed as a Bernoulli random variable (r.v.) whose distribution is parameterized
by the sensing value. In other words, the sensing value governs the transmission probability.
The randomized transmission is realized by the proposed encoder design that one-to-one maps
the sensing value using a sigmoid function (one with a “S” shape) to yield a variable with a
normalized range and then comparing the result with a uniform r.v. in the same range, generating
a sequence of binary bits. The bits switch the connection of the sensor antenna with either of
two load impedances and as the result, randomly turn backscatter transmission on/off, thereby
modulating the bits by on/off keying.
Statistical inference at the reader: Statistical inference is capable of estimating parameters
of a signal distribution based on signal observation. Using the theory and given randomized
transmission by sensors, the reader is designed to infer their transmission probabilities (or equiva-

6

lently their sensing values) from observations over time and antennas. To this end, algorithms for
reader inference are designed based on the expectation-maximization (EM) framework. A simpler
version of the algorithms are designed assuming no measurement noise and the assumption is
relaxed subsequently. To design the algorithms, a Bayesian network is constructed that relates
the signals received by the reader to the sensing values via a number of fixed or random variables
including the mapped sensing values, number of active sensors, variances of channel noise and
channel gains. Their relations are specified by their conditional distributions as derived. Building
on the Bayesian network, the specific EM-based algorithms targeting the BackSense system are
developed by deriving two iterative steps, called the Expectation and Maximization steps, for
both the criteria of maximum likelihood (ML) and maximum a posteriori (MAP).
The inference algorithms are extended to the case with measurement noise. As a result,
the Bayesian network in the preceding case should be modified to include a set of new r.v.,
modelling measurement noise, as many as the number of sensors. This dramatically increases
the dimensionality of the latent-space, the space of latent variables (r.v. unobserved by the reader).
As the result, the corresponding EM-based algorithms become computationally demanding and
thus impractical. To overcome the difficulty, we exploit an approximate implementation of the
EM-framework based on the variational-inference method [28]. Its basic idea for complexity
reduction is to restrict the approximate posterior distribution of the latent variables to take a
factorized form. This converts the required a larger number of nested integrals over the latentvariable space in the EM framework to parallel integrals, enabling computationally-efficient
implementation. Based on the method, we develop the EM-based algorithms for inference at the
reader for the ML and MAP criteria.
The remainder of the paper is organized as follows. Section II introduces the system model
and problem formulation. Section III presents the randomized transmission scheme at tags. The
EM-based algorithms for reader inference are designed in Sections IV and V for the cases
without and with measurement noise, respectively. Simulation results are provided in Section
VI, followed by concluding remarks in Section VII.
II. S YSTEM M ODEL AND P ROBLEM F ORMULATION
A. System Model
We consider a sensing system (see Fig. 2) comprising a PB, a reader equipped with M
antennas, and N single-antenna backscatter sensors. Provisioned with reliable power supply, the
PB is equipped with an antenna array and able to power the sensors by energy beamforming.

7

The reader uses the antennas for receiving simultaneous signals transmitted by the sensors. In the
following, the models of backscatter sensor, wireless channels and sensing values are described.
Provisioned with a backscatter antenna, each sensor reflects a fraction of continuous wave
(CW) back to the reader and harvests the energy of the remaining fraction for powering the sensor
circuit. In the process, the sensor modulates the reflected CW as follows. First of all, note that
the reflection coefficient S, the ratio between incident and reflected CWs, depends on the level of
mismatch between the antenna and load impedances. Mathematically, S =

ZL −ZA ∗
,
ZL +ZA ∗

where ZL and

ZA are load and antenna impedances, respectively [1]. The sensor modulates the backscattered
CW by adapting the reflection coefficient that determines the phase and magnitude of the wave
[1]. The coefficient variation can be implemented by switching over a set of load impedances
(see Fig. 2). We consider the modulation scheme based on on/off keying that is widely used
for backscatter transmission for its simplicity. The scheme requires switching between two load
impedances Z1 and Z2 with Z1 6= ZA∗ and Z2 = ZA∗ . In other words, for a sensor, switching to
Z1 turn on backscatter transmission (S =

Z1 −ZA ∗
Z1 +ZA ∗

= ρ̄) and to Z2 turns transmission off (S = 0).

The channel model is described as follows. Two kinds of channels are cascaded. One is from
the PB to the sensors for energy transfer (ET) and the other is from the sensors to the reader
for information transfer (IT). First, given energy beamforming and assuming that the size of
the sensor cluster is much smaller than its distance to the PB, ET channels can be modelled
as line-of-sight (LOS) channels with identical path losses. Second, with longer distances and
isotropic wave propagation, IT channels are assumed to be characterized by rich scattering and
hence modelled as Rayleigh fading. Specifically, time is slotted into symbol durations and we
consider block fading where channels coefficients are independent and identically distributed
(i.i.d.) over different time slots. The IT channel in slot i is denoted as a M × N matrix H(i) ,
(i)

whose (`, n)th element, denoted by h`,n , represents the coefficient from sensor n to the reader’s
(i)

(i)

antenna `. It is assumed that all channel coefficients {h`,n } are i.i.d. CN (0, σh2 ) r.v.. Let Sn be
the reflection coefficient of sensor n in slot i. The vector collecting the reflection coefficients for
(i)

(i)

(i)

all N sensors in slot i is then denoted by s(i) = [S1 , S2 , · · · , SN ]T . Given s(i) , the received
(i)

(i)

(i)

signal at the reader in slot i, denoted by y(i) = [y1 , y2 , · · · , yM ]T , is as follows1
p
i = 1, · · · , L,
y(i) = GPt H(i) s(i) + w(i) ,
(i)

(i)

(1)

(i)

where w(i) = [w1 , w2 , · · · , wM ]T is the additive white Gaussian noise (AWGN) with the entries

following i.i.d. CN (0, σw2 ) distributions, L denotes the length of the observation period in slot,
1

The interference from PB to the reader can be easily canceled as it is a CW and thus is not considered.

8

G captures the energy beamforming gain, and Pt is the transmit power of the PB. Without loss
of generality, G and Pt are set as one to simplify notation.
Assumption 1 (CSI Free). No CSI of individual backscatter channels is available at both of the
reader and the sensors. The reader only has the knowledge of distribution type of H(i) and w(i)
(i.e., complex Gaussians), but not the distribution parameters σh2 and σw2 .
The sensing values measured by different sensors are spatially correlated (e.g., temperature,
humidity, and wind strength). The model of their joint distribution is described as follows. Let the
sensing values to be measured by N sensors be denoted by x = [x1 , x2 , · · · , xN ]T . It is assumed

that the prior knowledge of the distribution of x is available at reader by estimation using
historical measurement data. For tractability, x is assumed to follow the multivariate Gaussian

distribution x ∼ N (µx , Σx ) with mean µx and covariance matrix Σx . The probability density
function (PDF) of x is then given by



1
T −1
− (x − µx ) Σx (x − µx ) .
p(x) =
1 exp
N
2
(2π) 2 |Σx | 2
1

(2)

Assumption 2 (Time-Scale Differentiation). The environmental features usually vary much
slower than wireless channels especially given reader mobility. Thereby, the sensing values x
are assumed to remain unchanged during the whole observation window of L time slots, while
wireless channels are i.i.d. over slots.
In addition, the sensor measurement of a sensing value x, denoted as x̃, is corrupted by
measurement noise: x̃ = x + ∆ where the r.v. ∆ represents the noise and assumed to follow the
N (0, δ 2 ) distribution. Moreover, measurement noises at different sensors and slots are i.i.d..
B. Problem Formulation
The conventional deterministic transmission design that first quantizes the sensing value and
then transmits the output bits is unsuitable for a backscatter sensor without an analog-to-digital
converter and under the constraint of on/off keying. To overcome the limitations, the proposed
BackSense design framework relies on a randomized design for sensor transmission and matching
statistical-inference algorithms for the reader to recover the transmitted sensing values. Two
corresponding design problems are formulated as follows.
1) Randomized Transmission Problem for Backscatter Sensors: Consider an arbitrary symbol
duration. Let Sn with Sn ∈ {0, ρ̄} denote the symbol transmitted by sensor n. The symbols

9

from N sensors are grouped as a vector s = [S1 , S2 , · · · , SN ]T . The design objective is to
randomize sensor transmissions such that the sensing value vector x is encoded into a sequence

of realizations of the random vector {s(i) }, whose distribution function, p(s(i) ), is parameterized

by the sensing values in x. Designing the mentioned random-encoding scheme is challenging.

First, it should be implementable under the constraint of on/off keying and using the simple
backscatter hardware architecture. Next, the scheme should yield a distribution function p(s(i) )
that allows tractable design of statistical-inference algorithms for the reader.
2) Statistical Inference Problem for Reader: The algorithms for statistical inference at the
reader are designed targeting randomized transmissions by backscatter sensors from solving the
preceding design problem. Furthermore, under the constraints of ALOHA-type random access
and no CSI, the algorithms should be capable of scaling up the inference accuracy with growing
correlation between sensing values and the number of received observations. To this end, we
design the inference algorithms based on the ML or MAP criterion. Let y(i) denote the symbol
vector received by the multi-antenna reader in the slot i. The set of data received over L slots is
represented by D. Moreover, let Φ = {x, σh2 , σw2 , δ 2 } denote the set of all the unknown parameters
including the desired sensing value x. Then the ML problem formulation is given as follows:
max log p(D|Φ) ⇐⇒ max

(ML)

Φ

Φ

L
X
i=1

log p(y(i) |Φ).

(3)

If the prior distribution of Φ is available, the MAP algorithms can be obtained using the following
MAP-to-ML conversion:
(MAP)

max log p(x|D) ⇐⇒ max log
Φ

Φ

⇐⇒ max
Φ

L
X
i=1

p(D|Φ)p(Φ)
p(D)
log p(y(i) |Φ) + log p(Φ).

(4)

The challenge for designing the inference algorithms based on the ML or MAP criteria arises
from the existence of latent variables. In statistics, a latent variable is defined as a variable not
directly observed but has a direct or indirect impact on the observed variables. Denote z as the
vector that collects all the involved latent variables, which capture the randomized effect of the
measurement noise, random transmitted symbols at sensors, and channel noise and coefficients as
specified in the sequel. Note that, via the latent variables, the design of random-encoding scheme
in the preceding problem formulation affects the design of inference algorithms. The difficulty
of solving (3) or (4) lies in the computation of the likelihood function therein which requires
R
marginalization over the latent variables: p(y|x) = p(y, z|x)dz. Essentially, the marginalization

10
Power beacon
Continuous sinusoidal wave (CW)
Backscatter signal
Sensing Measurement
values noise

x̃1
RN
generator

x̂1

Z2

Z1

xN

(i)
y1

Antenna 1
···

···
N

x̃2
···

(i)

< ✓1 ?

(i)
s1

2

Inferred
sensing values

Reader

x̃1
Sensor 2

x2
···

1

Fig. 2.

1

1

0
(i)
✓1

1

x1

Sensor 1

Sensor N

x̃N

Statistical
inference
(ML or MAP)

x̂2
···

Random encoding scheme

(i)

yM

Antenna M

x̂N

The design of BackSense system.

leads to undesired integrals inside the logarithm operation before the likelihood function, making
the direct optimization on (3) or (4) intractable. Furthermore, one can observe that given a large
set of latent variables (corresponding to a high-dimension space for z), such marginalization
involving many nested integrals can be computationally demanding. We tackle the challenge by
applying theories of EM and variational inference in the algorithmic design.
III. BACK S ENSE D ESIGN : R ANDOMIZED S ENSOR T RANSMISSION
In this section, the design of random-encoding scheme is presented that solves the design
problem formulated in Section II-B1. Its implementation on the BackSense architecture is illustrated in Fig. 2. Then the distributions of the randomized sensor signals are analyzed. The
results facilitate the design of reader inference algorithms in the following two sections.
A. random-encoding scheme
As shown in Fig. 2, the random-encoding scheme consists of the following two key steps.
Step 1) Range normalization. Distributed as a continuous variable in real domain, a measured
sensing value has an infinite range. The goal of the current step is to map it to a new variable
with a unit range so as to facilitate random encoding in Step 2. Consider the sensor measurement
x̃ of a sensing value x. Let θ be a variable within the range [0, 1]. Then x̃ can be one-to-one
mapped to θ using the well-known sigmoid function (a function characterized by a S-shaped
curve), denoted as F:
θ = F(x̃) =

1
,
1 + exp − x̃−µ
σ

(5)

11

where µ and σ are mean and variance of x, respectively. This gives the name of θ as the mapped
sensing value. It is assumed that the parameters µ and σ are known by the sensor and reader
via estimation using historical sensing data. Consider the sensing values and their measurements
for all N sensors, represented by the vectors x and x̃, respectively, with x̃ = x + ∆. The vector
∆ groups i.i.d. measurement-noise r.v. The element-wise mapping using (5) is represented by
F(x̃), where the mapping of the n-th element is specified by the corresponding parameters µn

and σn . The result is denoted as θ, namely θ = F(x̃). Since the mapping F is one-to-one, the
mapped sensing values θ are used in place of x̃ in the subsequent analysis and design.

Step 2) Encoding by comparison with a uniform distributed r.v. Consider an arbitrary sensor.
As shown in Fig. 2, the mapped sensing value θ is compared with a r.v., denoted as χ, uniformly
distributed in the same range of [0, 1]. The r.v. is generated independently for different symbol
durations. Then, the comparison yields a random sequence of binary bits for turning on/off
backscattering of the sensor (see Fig. 2). The r.v. χ can be generated locally using a simple
circuit (see e.g., [29], [30]). Let S denote an arbitrary symbol transmitted by the sensor using
on/off keying. Then S = ρ̄ if χ < θ or otherwise S = 0. Thereby, the distribution of S is
parameterized by the mapped sensing value θ as follows:

 ρ̄, w.p θ
S=
 0, w.p 1 − θ.

(6)

Remark 1 (How to decode the sensing value?). Given random encoding, the sensing value can
be decoded at the reader by inferring the statistics of many observed symbols transmitted by
sensors as illustrated by the following toy example case where all N sensors with identical
(i)

and fixed channel gains: hn = h for all 1 ≤ n ≤ N , 1 ≤ i ≤ L and have uniform mapped
sensing value θ to deliver. As a result, the reader receives a set of L symbols transmitted by
P
(i) (i)
(i)
(i)
sensors, { N
n=1 hn Sn + w }, where {w } represent the set of i.i.d. samples of the channel-

noise process. If L is sufficiently large and applying the law of large numbers, then the mapped
sensing value can be recovered as
L

N

1 X X (i) (i)
θ = lim
hn Sn + w(i) .
L→∞ ρ̄hN L
i=1 n=1
The measured sensing value is then obtained as x̃ = F −1 (θ). The above illustration shows how

the sensing value can be decoded in a simplified scenario and the performance of sensing-data
uploading improves with the growth of observation duration. However, the target scenario is
more complex involving practical factors including transmissions over multi-antenna channels,

12

unknown and fluctuating channel states, spatial variation of sensing values, and finite numbers
of observed symbols. This is the reason why the powerful statistical inference is needed for
recovering sensing values (see Sections IV and V).
B. Distributions of Randomized Transmissions
To facilitate statistical inference at the reader, we derive in this subsection the distributions of
some r.v. arising from randomized transmissions using the proposed random-encoding scheme.
1) Distribution of mapped sensing values: First, consider the case without measurement noise.
For the current case, the mapping between θ and x is one-to-one and thus θ can be considered
as an equivalence of x. Using the mapping in (5), the prior distribution of θ, denoted as p(θ),
can be derived based on that of x, i.e., p(x) given in (2), as presented below.
Lemma 1. For the case of heterogeneous sensing values without measurement noise, the distribution of θ is specified by
p(θ) =

N
Y

σn
θ − θn2
n=1 n

!


1
T −1
− r(θ) Σx r(θ) ,
N
1 exp
2
(2π) 2 |Σx | 2


1

where r(θ) is an element-wise function of θ given by




T
1
1
r(θ) = −σ1 log
.
− 1 , · · · , −σN log
−1
θ1
θN

(7)

(8)

Proof: See Appendix A.
For the simplified case of uniform sensing values (all elements of x are equal), the expression
given in (2) no longer holds as the covariance matrix Σx becomes singular. Alternatively, note
that x is specified by the scalar uniform value x only and whose distribution can be modelled by a
univariate Gaussian x ∼ N (µ0 , σ02 ), where µ0 and σ02 denote the mean and variance respectively.


2
0)
Thus the PDF of x is then given by p(x) = √ 1 2 exp − (x−µ
.
2
2σ
2πσ0

0

Accordingly, following the similar procedure as for deriving Lemma 1, the distribution of



x−µ0
θ = 1/ 1 + exp − σ0
can be derived as follows:
2 !
log 1θ − 1
1
√ exp −
p(θ) =
.
(9)
2
(θ − θ2 ) 2π
Next, consider the case with measurement noise. For this case, the mapped sensing-value

vector θ is a function of both the sensing values x and measurement noise ∆. The randomized
effect of the measurement noise on the mapping between θ and x is captured by the conditional

13

distribution of θ (a latent variable) given x (a parameter) which is required in the subsequent
reader inference. The result is derived using a similar method as Lemma 1 and given as follows.
Lemma 2. For the case of heterogeneous sensing values with measurement noise, the conditional
distribution function of θ given x is given as:



 PN
!
1
2
N
−
1
+
µ
−
x
)
(−σ
log
Y
n
n
n
n=1
θn
σn
1
.
p(θ|x) =
exp −
N
2
2
θ
−
θ
2δ
2 δN
n
(2π)
n
n=1

(10)

2) Distribution of the number of active sensors: A latent variable that affects the distribution
of received signal at the reader is the number of active (on) sensors, denoted as T (i) for the
i-th slot (see Fig. 3 and 4). Given on/off keying, T (i) is the sum over a sequence of Bernoulli
P
(i)
(i)
r.v. that represent the i-th symbols transmitted by all sensors: T (i) = N
n=1 Sn where {Sn }
follow the distribution in (6). Consider the case where the mapped sensing values differ over

sensors, denoted as θ1 , θ2 , · · · , θN . The r.v. is a Poisson-Binomial r.v. that describes the numbers
of successes in N independent Bernoulli trails with different individual success probabilities,

denoted by PB(N, θ1 , θ2 , · · · , θN ). The corresponding probability mass function (PMF) is given
as [31]:

N
N
1 X −`m Y
c
[1 + (c` − 1)θn ],
m = 0, 1, · · · , N,
(11)
N + 1 `=0
n=1

√
2πj
where c = exp N +1 with j = −1. For the special case of uniform mapped sensing values

p(T (i) = m|θ) =

denoted as θ, the distribution of T (i) reduces to the Binomial distribution:
 
N m
(i)
p(T = m|θ) =
θ (1 − θ)N −m ,
m = 0, 1, · · · , N.
m

(12)

IV. BACK S ENSE D ESIGN : R EADER I NFERENCE WITH I DEAL M EASUREMENTS
This section addresses the design problem of reader inference formulated in Section II-B2.
To this end, statistical inference algorithms are designed based on the EM framework, which
enable the reader to recover the sensing values form the received signals. Ideal measurements
are assumed and the effect of measurement noise on the design is addressed in the next section.
We start with a brief introduction to the basic principle of the EM framework. Then a Bayesian
network is constructed for specifying the statistical dependencies between the observations and
unknowns (including the latent variables and model parameters). Based on the Bayesian network,
EM algorithms targeting BackSense are designed for both the criteria of ML and MAP.

14

A. Principle of the EM Framework
Due to the existence of latent variables, directly solving the potentially non-convex ML and
MAP problems formulated in Section II-B2 is intractable. Alternatively, the EM framework
tempts to find local-optimal solutions by an iterative procedure, involving iterations between
two main steps, i.e., the Expectation step (E-step) and the Maximization step (M-step) (see e.g.,
[32]). In the E-step, the expectation of the complete-data log likelihood is used as a surrogate of
the required incomplete-data log likelihood as defined later. Then in the M-step, the surrogate is
maximized instead of the actual incomplete-data log likelihood. The use of the surrogate makes
finding the solution more tractable without compromizing its local-optimality. For exposition,
the details are discussed in the sequel based on the ML criterion while the same principle also
applies to the MAP criterion.
Assume that we have L i.i.d. observations. We denote all of the observed variables as Y =
[y(1) , y(2) , · · · , y(L) ], all of the latent variables as Z = [z(1) , z(2) , · · · , z(L) ], and all of the
parameters as Φ. By introducing an auxiliary distribution q(Z) defined over the latent variables
Z, the incomplete-data log likelihood can be decomposed by
log p(Y|Φ) = L(q, Φ) + KL(q||p),
(13)
n
o
R
where L(q, Φ) = q(Z) log p(Y,Z|Φ)
dZ is the expectation of the complete-data log likelihood
q(Z)
o
n
R
dZ is known as Kullback-Leibler (KL) divergence.
and the KL(q||p) = − q(Z) log p(Z|Y,Φ)
q(Z)
Note that the KL divergence measures the “similarity” between the two distributions of q and
p. Hence the measure is a non-negative function attaining the minimum when q and p are
identically distributed. Therefore, it follows from (13) that L(q, Φ) lower bounds log p(Y|Φ).
One can see that the bound reaches equality, namely KL(q||p) = 0, by letting q(Z) = p(Z|Y, Φ).
This essentially gives the E-step. The resultant tight lower bound can then serve as a surrogate
of log p(Y|Φ) that is optimized to provide updated parameters Φ, constituting the M-step. As a
result, the maximization of (13) can be solved by repeating the following EM cycle, where the
subscript (t) is used to denote the iteration index:
E−step :

M−step :

q (t+1) (Z) ← arg max L(q, Φ(t) )
q

(14)

=⇒ q (t+1) (Z) = p(Z|Y, Φ(t) ).

(15)

Φ(t+1) ← arg max L(q (t+1) , Φ)
Φ
Z
(t+1)
=⇒ Φ
← arg max p(Z|Y, Φ(t) ) log p(Y, Z|Φ(t) )dZ.

(16)

Φ

(17)

15

The above EM cycle starts from some initial values for the parameters Φ and keep repeating until
convergence. The convergence of the EM framework to at least a local optimum is guaranteed
since both the E-step (14) and M-step (16) lead to a non-decreasing objective L(q, Φ) at each
EM cycle [32].

B. Construction of Bayesian Network
Designing the BackSense inference algorithms begins with constructing a Bayesian network
that specifies the set of parameters and latent variables affecting the received signal as well as
their interdependence characterized by their conditional distributions. The set of parameters are
unknown to the readers and thus the target of inference. The parameter set Φ comprises the
mapped sensing values θ, channel variance σh2 and channel-noise variance σw2 : Φ = {θ, σh2 , σw2 }.
For the current case without measurement noise, there exists only one type of latent variables,

namely the set of randomized symbols transmitted by sensors, that affect the observed variable,
namely the received symbols over multiple antennas.
Consider the i-th slot. To complete the construction of the Bayesian network, the conditional
distribution of the received symbol vector y(i) given the transmitted symbol vector s(i) is derived
as follows. Note that the elements of the additive noise and the channel matrix in (1) are
complex Gaussians. It follows from (1) that given s(i) , y(i) is a sum of complex Gaussian vectors
that also follows the complex Gaussian distribution. Specifically, conditioned on s(i) , y(i) ∼
P
(i)
2
CN (0, σh2 N
n=1 Sn + σw ). One can see that given on/off keying for sensor transmissions, the
P
(i)
(i)
where T (i) is the number of active tags and ρ̄ is the reflection
sum σh2 N
n=1 Sn is equal to ρ̄T
coefficient for an active sensor. Without loss of generality, we set ρ̄ = 1 in the subsequent analysis
for ease of notation. Mathematically,
p(y

(i)

(i)

| s ) = p(y

(i)

(i)

|T )=

1

ky(i) k22
exp −
2(σh2 T (i) + σw2 )




.
(18)
M
[2π(σh2 T (i) + σw2 )]
This suggests that the latent variable s(i) can be replaced with T (i) without changing the distribution of y(i) . The distribution of T (i) is derived earlier as given in (11) and (12). Note that
the variable replacement reduces the dimensions of latent space from N to just 1, reducing the
complexity of inference. Based on the above discussion, the Bayesian network corresponding to
the i-th received symbol is illustrated in Fig. 3.
C. Reader ML Inference with Uniform Sensing Values
In this subsection, we consider the simple case of uniform sensing values θ and the ML
criterion for designing the reader-inference algorithm. The design is subsequently extended to

16

Unknown
parameter set

Sensing
values

Mapped
measurement

x

✓

2
w

Latent variable
(# of active tags)

Channel
variance

2
h

y(i)

T (i)

Variable set

Fig. 3.

Noise
variance

Observed variable
(Received signals)

Bayesian network for BackSense inference without measurement noise.

the case of heterogeneous sensing values as well as the MAP criterion. Based on the EM
framework introduced in Section IV-A, the algorithm for inferring θ is designed for BackSense
by driving the corresponding E-step and M-step shortly. The resultant procedure is summarized
in Algorithm 1. Given the inferred sensing value θ̂ and the mapping in (5), the inferred sensing
value is then obtained as

x̂ = σ0 log

1
θ̂


− 1 + µ0 .

(19)

a) Derivation of E-step
Consider the Bayesian network in Fig. 3. Given the parameter set Φuni = {θ, σh2 , σw2 }, the required
posterior distribution of the latent variable T (i) can be computed via Bayes’ Law as
fm
,
q(T (i) = m) = p(T (i) = m | y(i) , Φuni ) = P
n fn

(20)

where fm is defined as
fm =

1
M

[2π(σh2 m + σw2 )]


exp −

ky(i) k22
2(σh2 m + σw2 )

 
N m
θ (1 − θ)N −m , m = 0, 1, · · · , N. (21)
m
(i)

This completes the derivation of the E-step as indicated in (15). For ease of notation, let qT (m)
denote q(T (i) = m) in the rest of the paper.
b) Derivation of M-step
According to the M-step in (17), the parameter set Φuni should be optimized by solving
Φ∗uni

= arg max
Φuni

L X
N
X
i=1 m=0

(i)

qT (m) log p(y(i) , T (i) = m | Φuni ),

(22)

17

Algorithm 1 EM algorithm for ML reader inference
Initialization:
Initialize the model parameters θ, σh2 , σw2 properly.
Iteration:
(i)

1) E-step: Update the posterior distribution of the latent variable, {qT (m)}, according to (20).
(i)

2) M-step: Update the parameter θ and σh2 , σw2 by substituting the latest values of {qT (m)} into
(25) and (29), respectively.
Until Convergence.

where the joint distribution over the observed and latent variables is given by
p(y(i) , T (i) = m | Φuni ) = fm ,

m = 0, 1, · · · , N,

(23)

where fm is defined in (21). Recall that the parameter set Φuni = {θ, σh2 , σw2 }. A close observation
of (22) and (23) reveals that θ and {σw2 , σh2 } are decoupled in the objective function in (22) and
thereby they can be optimized separately.

Firstly, we derive the updating formula for the parameter θ by considering those terms in (22)
that are related to θ only, the corresponding optimization problem is given as
∗

θ = arg max
θ∈(0,1)

L X
N
X
i=1 m=0

(i)

qT (m)[m log θ + (N − m) log(1 − θ)],

(24)

Proposition 1.A. The value of θ that solves the optimization problem in (24) is given by
L

N

1 X X (i)
θ =
q (m)m.
LN i=1 m=1 T
∗

(25)

Proof: See Appendix B.
The above result yields the formula for updating the parameter θ in the M-step.
Next, to derive the updating formula for the parameters σw2 and σh2 , group the relevant terms
in the objective function in (22). The corresponding optimization problem is


L X
N
X
ky(i) k22
(i)
2
2
.
max
qT (m) −M log(σh m + σw ) −
2 ,σ 2
2(σh2 m + σw2 )
σh
w
i=1 m=0

(26)

The problem is non-convex and thus it is difficult to find its solution in closed form. To overcome
the difficulty, we develop a sub-optimal approach that relaxes the problem to a quasiconvex one
and solving it provides closed-form updates for the target parameters σh2 and σw2 . For this purpose,
one can observe that the objective function depends on σh2 and σw2 via a weighted sum, defined

18

as Σm = σh2 m + σw2 . By substituting Σm , the optimization problem in (26) is converted to the
following quasiconvex problem:
max
{Σm }

L X
N
X

(i)
qT (m)

i=1 m=0



ky(i) k22
.
−M log(Σm ) −
2Σm

(27)

Proposition 1.B. The optimal values of {Σm } that solve the problem in (27) are given by
PL (i)
q (m)ky(i) k22
∗
,
m = 0, 1, · · · , N.
(28)
Σm = i=1 PT L (i)
2M i=1 qT (m)
Proof: See Appendix C.
Note that in general, it is infeasible to find the values of σh2 and σw2 that satisfy (N + 1)
equations: Σm = σh2 m + σw2 with m = 0, 1, · · · , N . Thus, we resort to finding approximate
values, denoted as (σw2∗ , σh2∗ ), by linear regression:
(σh2∗ , σw2∗ ) = arg min2

2 ,σ )
(σw
h

N
X
m=0

σw2 + mσh2 − Σ∗m

2

PN

N
∗
N
m=0 (m − 2 )(Σm − Σ̄)
, Σ̄ − σh2∗
PN
N 2
2
m=0 (m − 2 )

=

!
(29)

where Σ̄ denotes the mean of the optimal values of {Σ∗m } derived in (28). The above result

gives the formula for updating the parameters σh2 and σw2 for the M-step. Combining this result
and that in Proposition 1.A completes the derivation of the M-step.
D. Extension: Reader MAP Inference
The EM algorithm based on the ML criterion as shown in Algorithm 1 is extended to the
MAP criterion as follows. Comparing (3) and (4), one can see that the objective function for the
MAP problem differs from the ML counterpart only by the extra term log p(Φ), corresponding
to the “log-prior” of the unknown parameters. Thus, deriving the MAP algorithm requires only
modifying the M-step of the ML counterpart that updates the estimation of the parameters.
Specifically, the M-step for the MAP algorithm can be written as
Φ∗uni

= arg max

Recall that Φuni =

Φuni

L X
N
X
i=1 m=0

{θ, σh2 , σw2 }.

(i)

qT (m) log p(y(i) , T (i) = m | Φuni ) + log p(Φuni ).

(30)

The prior distribution of θ is derived as shown in (9). However,

for those of σh2 and σw2 , in practice, such information is difficult to obtain from historical inference
results or estimate due to lack of training signals. A common technique to tackle the difficulty

19

is to replace the required prior distributions with some constants, known as non-informative
priors (see e.g., [33]). To some extent, the technique reduces the performance gain of MAP
interference over ML inference. Using the technique and prior distribution of θ in (9), the
optimization problem for updating θ in the M-step for the MAP inference is modified from the
ML counterpart in (24) as
max
θ∈(0,1)

L X
N
X
i=1 m=0

(i)
qT (m)[m log θ


+ (N − m) log(1 − θ)] − log θ(1 − θ) −

log

2
−1
.
2

1
θ

(31)

Proposition 2. The optimal value of θ that solves the optimization problem in (31) satisfies


1
(LN −2)θ∗
e
=a ∗ −1 ,
(32)
θ
P P

(i)
L
N
q
where a is a constant and given by a = exp
(m)m
−
1
.
i=1
m=0 T
The proof is straightforward and omitted here for brevity. Though θ∗ has no closed form, it can
be computed by a simple numerical search. On the other hand, due to the use of non-informative
priors, the procedure for updating the other parameters, i.e., σh2 and σw2 , remains the same as the
ML counterpart.
In summary, the iterative algorithm for MAP inference for the current case can be modified
from Algorithm 1 by replacing (25) in Step 2) with solving the equation in (32).
Remark 2. (ML versus MAP) A comparison between (22) and (30) reveals that, for MAP
inference, the reader’s knowledge of parameters’ prior distributions servers as an additional
regularization term in the objective function to avoid overfitting to the potential outliers in the
data set. This makes the performance more robust even when there are errorneous data points
in the data set or its size is not large enough.
E. Extension: Heterogeneous Sensing Values
In the preceding sub-sections, we designed inference algorithms for the simple case of uniform
sensing values. They are extended in this sub-section to the general case of heterogeneous sensing
values. The corresponding parameter θ is a vector comprising N elements instead of being a
scalar in the preceding case. Their joint distribution is derived earlier as shown in (7). Given their
relation shown in Fig. 3, the change in θ causes the number of active tags, T (i) , to follow the
Poisson-Binomial distribution in (11) instead of the Binomial distribution in the case of uniform
sensing values. The above changes in parametric distributions require the EM algorithms for
reader inference to be modified as follows.

20

a) Derivation of E-step
According to (15), given initial values of the unknown parameters Φcor = {θ, σh2 , σw2 }, the E-step
calculates the posterior distribution of the latent variable by Bayes’ law as follows.
gm
,
q(T (i) = m) = p(T (i) = m | y(i) , Φcor ) = P
n gn

(33)

where gm with m = 0, 1, · · · , N is defined as


N
N
ky(i) k22
1
1 X −`m Y
exp −
gm =
[1 + (c` − 1)θn ].
c
2
M
2
2
2
2(σh m + σw ) N + 1 `=0
[2π(σh m + σw )]
n=1

(34)

b) Derivation of M-step
The formulas for updating the parameters σw2 and σh2 remain the same as in (29). Only that for
updating θ needs to be modified as follows. Similar to (24), by considering only those terms
related to θ in the M-step optimization, the update of θ requires to solve the following two
optimization problems based on the ML and MAP criteria, respectively:
(ML)
(MAP)

max g̃(θ),
θ
"

(35)

#
1
max g̃(θ) −
log(θn − θn2 ) − r(θ)T Σ−1
x r(θ) ,
θ
2
n=1
N
X

(36)

where g̃(θ) is defined as
g̃(θ) =

L X
N
X
i=1 m=0

(i)
qT (m) log

( N
X
`=0

c−`m

N
Y

)
[1 + (c` − 1)θn ] .

(37)

n=1

It can be observed from the above two optimization problems that the mapped sensing values
{θn } are coupled with each other, making the problem non-convex and the direct solutions
intractable. In the area of statistical inference, a common approach for tackling this challenge,
known as the generalized EM (GEM) framework (see e.g., [34]), is to update the parameters
in a tractable way to increase the M-step objective in (35) or (36) (e.g., using gradient ascent)
instead of exactly maximizing it. Using this technique, each complete EM cycle of the GEM
algorithm is guaranteed to increase the value of the log likelihood objective until converging to
a local optimum [34]. Specifically, the iterative equation for updating θ can be derived as shown
below:
θ (t+1) = θ (t) + α∇L(θ (t) ),

(38)

21

where α is the step size and ∇L(θ) denotes the derivative of the objective function in either
(35) or (36), representing the gradient direction, and whose the n-th element is given by:
PN −`m `
Q
L X
N
`
X
(c − 1) N
`=0 c
n0 6=n [1 + (c − 1)θn0 ]
(i)
[∇L(θ)]n (ML) =
qT (m)
.
(39)
PN −`m QN
`
`=0 c
n=1 [1 + (c − 1)θn ]
i=1 m=0
h
i
 −1

σn
1 − 2θn
−
Σx r(θ) n .
(40)
[∇L(θ)]n (MAP) = ∇θ L̃(θ) (ML) −
2
θn (1 − θn ) (θn − θn )
n
In summary, the iterative algorithm for reader inference in the current case can be modified

from Algorithm 1 by simply replacing (20) in the E-step with (33) and (25) in M-step with (38).
V. BACK S ENSE D ESIGN : R EADER I NFERENCE WITH N OISY M EASUREMENTS
In the preceding section, the mapped sensing value vector θ is a parameter in the Bayesian
network (see Fig. 3) since it is one-to-one mapped from the sensing-value vector x and remains
fixed throughout the observation duration. In this section, in the presence of measurement
noise, θ is a random vector conditioned on x and varies from symbol-to-symbol, making it
a latent variable in the Bayesian network. This causes the dimensionality of the latent space
increasing linearly with the number of sensors. To overcome the resultant prohibitive complexity
of statistical inference for a large number of sensors, we design practical approximate EM
algorithms using the method of variational inference. The details are presented targeting ML
inference. The extension to the MAP inference is straightforward and omitted for brevity.
A. Construction of Bayesian Network
With the addition of measurement noise, the corresponding Bayesian network is modified from
the noise-free counterpart (see Fig. 3) as shown in Fig. 4 and described as follows. First of all,
the mapped sensing vector θ (i) is changed from a parameter to a latent variable. Moreover, an
extra parameter δ 2 , namely the measurement-noise variance, is included that affects the latent
variable θ (i) together with the sensing-value vector x. The corresponding conditional distribution
of θ (i) for given δ 2 and x is derived earlier as shown in (10). On the other hand, the conditional
distributions of the other latent variable T (i) and the observed variable y(i) are identical as their
measurement-noise-free counterparts given in (11) and (18), respectively.
B. Reader Inference Based on Approximate EM Implementation
1) Principle of the Variational Inference: The method of variation inference is one lowcomplexity approximation of the EM framework. The key idea is to apply the mean field theory

22

Unknown
Parameter set

Variable set

Fig. 4.

Sensing
values

Measurement
noise variance

x

2

Noise
variance

Channel
variance

2
w

2
h

✓ (i)

T (i)

y(i)

Latent variable A
(Mapped measurement)

Latent variable B
(# of active tags)

Observed variable
(Received signals)

Bayesian network for BacksSense with the measurement noise.

[35] that constrains the variational distribution q(Z) in the EM framework [see (14): E-step
and (16): M-step] to belong to a family of distribution functions charaterized by the factorized
Q
form q(Z) = k qk (zk ), where zTk is the k-th row of Z and qk is the corresponding factor
function [28]. When a Bayesian network can be decoupled into sub-networks or the complete
data likelihood can be factorized into product form, the approximation can reduce integrals over
a high-dimensional latent space to multiple parallel integrals with low dimensionality and hence
Q
low complexity. By substituting the approximate form of q(Z) = k qk (zk ), the two steps of
EM framework in (15) and (17) can be modified and written as follows [28]:
!
Z
Y
Y

Variational E−step : q (t+1) (Z) = b
exp
log p Y, Z|Φ(t)
qk (zk )dzk , (41)
`

Variational M−step : Φ(t+1) ← arg max
Φ

k6=`

Z

q (t+1) (Z) log p(Y, Z|Φ(t) )dZ,

(42)

where b is a normalization constant whose value can be computed by inspection. One can
observe from (41) and (42) that the high-dimensional integrals therein can be reduced into low
dimensional ones if the term p Y, Z|Φ(t) can be decomposed as the product of factors that are
functions of non-overlapping subsets of the rows of Z.
2) Inference Algorithm Derivation: Based on the Bayesian network in Fig. 4 and applying
the mean field approximation, the variational distribution in the EM framework is constrained
Q
(i)
to take the fully factorized form q(z(i) ) = q(T (i) ) n q(θn ). Then according to (41) we have
n h
io
q ∗ (T (i) ) = b0 exp Eθ log p(y(i) , T (i) , θ (i) | Φnoi ) ,
(43)
n
h
io
q ∗ (θn(i) ) = bn exp ET,θ\θn log p(y(i) , T (i) , θ (i) | Φnoi ) ,
∀n,
(44)
where {bn } are normalization constants. The explicit expressions of the approximate marginal

distributions of individual latent variables shown in (43) and (44) can be derived based on

23

the following procedure. Considering (43), apply the conditional distributions derived earlier
as shown in (10), (11), (18) and together with the chain rule, followed by separating the terms
depending on T in (43) from those on θn . Similar procedure can be applied to (44). As the result,
we can obtain the said marginal distributions as shown below, giving the variational E-step.
Proposition 3.A. The variational E-step involves calculation of the approximate marginal distributions of latent variables given as:
ga (m)
q ∗ (T (i) = m) = P
,
k ga (k)

q ∗ (θn(i) = zn ) = R 1
0

gb (zn )
gb (zn )dzn

, ∀n,

(45)

where ga (m) and gb (zn ) are given as


F1 (m)
ky(i) k22
ga (m) = 2
exp −
,
σw + mσh2
2(σw2 + mσh2 )
(



2 )
1
F2 (zn )
1
exp − 2 −σn log
gb (zn ) =
− 1 + µ n − xn
,
zn (1 − zn )
2δ
zn
wherein we define two functions for ease of notation:
n  X N
o
YN
F1 (m) = exp Eθ log
c−`m
[1 + (c` − 1)θk ]
.
`=0
k=1
n
 XN
o
YN
F2 (zn ) = exp ET,θ\θn log
c−`T
[1 + (c` − 1)θk ]
`=0

k=1

(46)
.

(47)

θn =zn

Next, we derive the variational M-step. To this end, substitute the marginal distributions of
latent variables in (45) into (42), and solve the resultant optimization problems. Thereby, we
have the following set of equations for parametric updating, constituting the variational M-step.
Proposition 3.B. For the variational M-step, the parameters {x, δ 2 , σh2 , σw2 } are updated as:
L Z

1X
∗
n = 1, 2, · · · , N.
q(θn(i) )σn log 1/θn(i) − 1 dθn(i)
xn = µ n −
L i=1
L
N Z
1 XX
2∗
δ =
q(θn(i) )(−σn log(1/θn(i) − 1) + µn − xn )2 dθn(i) .
LN i=1 n=1
PL
q(T (i) = m)ky(i) k22
∗
Σm = i=1 PL
,
m = 0, 1, · · · , N
2M i=1 q(T (i) = m)
!
PN
N
∗
(m
−
)(Σ
−
Σ̄)
N
m
m=0
2
(σh2∗ , σw2∗ ) =
, Σ̄ − σh2∗
PN
N 2
2
m=0 (m − 2 )
Finally, the main steps of the approximate EM algorithm is summarized in Algorithm 2.

24

Algorithm 2 Summary of the approximate EM algorithm using variational inference
Initialization:
(i)

Initialize the model parameters x, δ 2 , σh2 , σw2 , and the marginal posteriors q(T (i) ) and {q(θn )};
Iteration:

1) E-step: Compute the approximate marginal posteriors of latent variables iteratively by cycling
the variational distributions presented in (45) in Proposition 3.A for rounds until convergence;
2) M-step: Update the parameters according to the formulas in Proposition 3.B;
Until Convergence.
Remark 3. (Complexity Reduction by Variational Inference) It can be observed from Proposition 3.B that only one-dimension integrals are required in the variational M-step. However, high
dimension integrals are still needed in the variational E-step (see Proposition 3.A) for the reason
that given the parameters, the joint distribution function of observation and latent variables cannot
be decomposed into the fully factorized form. Despite this, the method of variation inference
still leads to substantial complexity reduction with respect to the exact EM inference where both
the E-step and M-step involve N -dimension integrals over the full latent space.

VI. S IMULATION R ESULTS
For simulation, the performance metric is the inference error, defined by

kx̂−xk
.
kxk

The simulation

parameters are set as follows unless specified otherwise. The average transmit signal-and-noise
ratio (SNR), defined as Pt /σw2 , is set to be 10 dB. The variances of the wireless channel and
the additive Gaussian noise are σh2 = σw2 = 1. For the prior distribution of x, the mean is
µn = 25, ∀n, and the covariance matrix Σx can be decomposed as Σx = D1/2 Cx D1/2 , where

D is a diagonal matrix with the diagonal elements collecting the variances of the elements of
x, namely [D]n,n = σn2 , and we set σn2 = 1, ∀n. Moreover, the correlation matrix Cx is assumed
to have the linear exponent autoregressive structure [36]:

1
ρ
ρ2 · · ·

 ρ
1
ρ
···

Cx =  .
.
.
..
..
 ..

ρN −1 ρN −2 ρN −3 · · ·

ρN −1



ρN −2
..
.




,



1

in which the parameter ρ ∈ [0, 1) controls the level of correlation between the adjacent sensor

tags, the correlation between tags further away is assumed to decrease with the growth of their
spatial separation following a geometric sequence 1, ρ, ρ2 , · · · .

25

A. Effect of Observation Duration
Consider the scenario without the measurement noise with N = 4 and ρ = 0.5. Increasing
the observation duration, namely L slots (or symbol durations), generates a larger data set and
thereby improves the inference accuracy. The gain is evaluated in Fig 5(a). It is observed that the
inference errors for both the cases of ML and MAP inference are monotonic decreasing functions
of L. MAP inference outperforms ML especially in the regime with a relatively short observation
duration (small L) and the difference diminishes as L increases. In this data-deficient regime,
ML tends to overfit the potential outliners arising from noise corruption. On the other hand,
overfitting is mitigated by MAP due to a regularization term that attempts to align inference
results with the prior distribution (see Remark 2). The performance convergence between ML
and MAP in the data-sufficient regime (large L) is due to the increasingly dominant contribution
from the likelihood function to the posterior objective as indicated in (4), which reduces that
of the regularization imposed by the prior distribution. Furthermore, the inference accuracy is
also observed to increase with the number of reader’s receive antenna M . The multi-antenna
gain arises from the diversity gain resulting from the increased observation dimension. The same
observations also hold for the scenario with measurement noise.
B. Effect of Spatial Correlation of Sensing Values
The impact of the correlation coefficient ρ on the inference accuracy is evaluated in Fig.
5(b). The measurement noise is assumed negligible, and we set M = 4 and L = 100. One can
observe from the curves that spatial correlation between sensing values helps reader inference by
improving its accuracy. In particular, MAP inference sees a much more significant performance
improvement compared with the ML counterpart as ρ increases. The reason is that the correlation
information is explicitly encoded in the prior distribution and exploited by the MAP inference
but not ML. In other words, the accessibility of the prior distribution of the environment variables
plays a critical role in determining how much gain can be exploited from the spatial correlation
between them. The above observation hold for different N settings. However, increasing N leads
to a decline in decoding accuracy as the number of unknown to be inferred by the reader grows.
C. Effect of Average Transmit SNR
Fig. 5(c) displays the curves of inference error versus the average transmit SNR for different
values of M and L. The figure only shows the curves for MAP inference without measurement
noise. Simulation results for other cases (e.g., ML inference and the presence of measurement

26

0.2

0.06

(1)

Inference Error

ML inference with M = 4
MAP inference with M = 4
ML inference with M = 8
MAP inference with M = 8

(2)
Inference Error

(1)
(2)
(3)
(4)

0.1

(1)
(2)
0.01

(3)
0.01

(4)

(1) ML inference with N = 8
(2) MAP inference with N = 8
(3) ML inference with N = 4
(4) MAP inference with N = 4

(3)
(4)
0.002
20

40

60

80

100

120

140

160

180

0.001

200

0

0.1

0.2

0.3

Size of data set

0.6

0.7

0.8

0.9

1

(b) Effect of Spatial Correlation of Sensing Values

3 × 10-2

10-1

Noise-limited region

Data-limited region

(1)

10-2

Inference Error

Inference Error

(2)
(3)
(4)

10-3

(1)
(2)
(3)
(4)
0

MAP
MAP
MAP
MAP

2

4

inference
inference
inference
inference
6

8

with
with
with
with
10

M=4,
M=8,
M=4,
M=8,
12

10-2

(1), (2), (3), (4)

10-3

L=100
L=100
L=300
L=300

14

16

18

20

22

24

26

28

30

Average Transmit SNR (dB)

Fig. 5.

0.5

Correlation

(a) Effect of Observation Duration

3 × 10-4

0.4

10-4
0.01

(1) ML inference by approximate EM implementation
(2) ML inference by exact EM implementation
(3) MAP inference by approximate EM implementation
(4) MAP inference by exact EM implementation

0.1

1

10

Variance of measurement noise

(c) Effect of Average Transmit SNR
(d) Effect of Measurement Noise
The effects of system parameters on the accuracy of reader’s statistical inference.

noise) are omitted as they all show similar trends. One can observe that as the SNR increases, the
inference error sees fast decrease in the low-to-moderate SNR regime. However when the SNR
is large (exceeds about 20 dB), the accuracy cannot be further improved by simply increasing
the transmission power. The reason is that in the high SNR regime the bottleneck of accurate
inference is no longer the noise but the data-set size. Increasing the size (or equivalently the
observation duration) can lower the error-saturation level. Furthermore, it is also noted that
increasing the number of receive antenna M can considerably alleviate the distortion from noise
when SNR is low due to the said diversity gain, but still suffer from the same error floor in the
high SNR regime.
D. Tightness of Variational-Inference Approximation
Consider the scenario with presence of measurement noise, resulting in the dimensionality
of the latent space linearly growing with the number of sensors. The performance of the exact

27

EM implementation is compared with its approximate implementation by variational inference
in Fig. 5(d) over varying the variance of measurement noise. It is observed that, for both the
cases of ML and MAP inference, the approximation leads to only a marginal performance
loss compared with exact EM but leads to reduced computation complexity. This confirms
the effectiveness of variational inference. Moreover, one can also observe that the inference
performance is satisfactory even in the regime of strong measurement noise (e.g., the error is
less than 0.1 when δ 2 = 10). This shows the robustness of reader inference against the corruption
of measurement noise.
VII. C ONCLUDING R EMARKS
The work presents a novel design framework for backscatter sensor systems. The design
features the application of machine learning to tackle hardware and signal processing limitations
of backscatter sensors as well as reducing overhead (due to e.g., scheduling protocol and
channel estimation and feedback). This work represents an initial investigation of a promising
approach of leveraging powerful techniques from machine learning to implement intelligent
sensor networks. There are many interesting directions warranting further research. For instance,
for large-scale networks with mobile readers (e.g., for Smart Cities), machine learning can be
applied to jointly design reader inference and path planning. As another example, inference from
large-scale sensing data can be centralized in the cloud using advanced machine learning such
as deep learning. Alternatively, a hierarchical framework can be developed for data inference
combining centralized and distributed inference. All such interesting directions call for novel
designs integrating machine learning and wireless communications in the same vein as the current
work but in potentially more complex ways.
A PPENDIX
A. Proof of Lemma 1
Starting from the cumulative distribution function (CDF) of θ, denoted by Fθ (u), by definition
and the mapping function shown in (5) we have


1
1

 < u1 , · · · ,

 < uN 
Fθ (u) = Prob 
x1 −µ1
xN −µN
1 + exp − σ1
1 + exp − σN






1
1
= Prob x1 < −σ1 log
− 1 + µ1 , · · · , xN < −σN log
− 1 + µN
u1
uN






1
1
= Fx −σ1 log
− 1 + µ1 , · · · , −σN log
− 1 + µN ,
(48)
u1
uN

28

in which Fx (·) represent the CDF of x. Then the PDF of θ can be evaluated by taking derivative
of (48) with respect to u = [u1 , u2 , · · · , uN ]T , namely






∂
1
1
p(θ) =
Fx −σ1 log
− 1 + µ1 , · · · , −σN log
− 1 + µN
∂u
u1
uN
u=θ






N
Y σn
1
1
− 1 + µ1 , · · · , −σN log
− 1 + µN ,
=
f −σ1 log
2 x
θ
−
θ
θ
θ
n
1
N
n
n=1

(49)

where fx (·) = p(x) denotes the PDF of x as given by (2). Then the desired result can be obtained
by substituting (2) into (49), completing the proof.
B. Proof of Proposition 1.A
Let L̃(θ) denote the objective function shown in (24), then the first and second derivatives of

L̃(θ) with respect to θ can be expressed as follows


L X
N
X
m N −m
(i)
0
qT (m)
L̃ (θ) =
−
.
θ
1−θ
i=1 m=0


L X
N
X
N −m
m
(i)
00
L̃ (θ) = −
qT (m) 2 +
.
2
θ
(1
−
θ)
i=1 m=0

(50)

(51)

(i)

Note that qT (m) is non-negative as it is a probability measure, hence it is easy to see that
L̃00 (θ) < 0, ∀θ ∈ (0, 1), suggesting L̃(θ) is a concave function over its domain θ ∈ (0, 1). As

a result, the maximum of L̃(θ) can be achieved by setting its first derivative to zero, leading to
the desired result in (25). However, one should take note that the zero of L̃0 (θ) may not lie in

the range of (0, 1), thereby requiring further validation of the solution. Fortunately, the result
shown in (25) is guaranteed to meet the constraint θ ∈ (0, 1). To see this, let’s rewrite (25) as

P  P
(i)
q
(m)m
. Note that the term inside the bracket can be identified as
θ∗ = L1 Li=1 N1 N
m=0 T
the expectation of T (i) with respect to its posterior distribution q (i) normalized by its maximum
value N which should fall in the range (0, 1) for each i. Hence a further average over i still
keeps the value within (0, 1), which completes the proof.
C. Proof of Proposition 1.B
Denoting L̃(Σm ) as the objective function shown in (27). Then taking the first and second

derivatives of L̃(Σm ) with respect to Σm gives
∇Σm L̃(Σm ) =

L
X
i=1

(i)
qT (m)



M
ky(i) k2
−
−
Σm
2Σ2m


.

(52)

29

∇2Σm L̃(Σm ) = −

L
X

(i)

qT (m)

i=1



(i)

M Σm − ky k
Σ3m


2
=

Σm −

(i)
(i) k2
2
i=1 qT (m)ky
PL
(i)
M i=1 qT (m)

PL

Σ3m /(M

(i)
i=1 qT (m))

PL

.

(53)

Although ∇2Σm L̃(Σm ) is not non-positive for all Σm ∈ R+ , in other words, L̃(Σm ) is not a

concave function on its domain Σm ∈ R, we can still prove that the maximum of L̃(Σm ) is
attained at the zero of its first derivative due to its quasi-concavity as shown below. Particularly,
note that there is only one root for ∇Σm L̃(Σm ) = 0, namely Σ∗m derived in (28). Furthermore,

as can be seen from (53) that the second derivative at that critical point is strictly less than zero,
i.e., ∇2Σm L̃(Σ∗m ) < 0. The above observations plus the fact that L̃(Σm ) is a continuous function
and differentiable everywhere on its domain Σm > 0 lead to the desired result.
R EFERENCES
[1] C. Boyer and S. Roy, “Backscatter communication and RFID: Coding, energy, and MIMO analysis,” IEEE Trans. Commun.,
vol. 62, no. 3, pp. 770–785, Mar. 2014.
[2] G. Yang, C. K. Ho, and Y. L. Guan, “Multi-antenna wireless energy transfer for backscatter communication systems,”
IEEE J. Sel. Areas Commun., vol. 33, no. 12, pp. 2974–2987, Dec. 2015.
[3] J. Kimionis, A. Bletsas, and J. N. Sahalos, “Increased range bistatic scatter radio,” IEEE Trans. Commun., vol. 62, no. 3,
pp. 1091–1104, Mar. 2014.
[4] B. Clerckx, Z. B. Zawawi, and K. Huang, “Wirelessly powered backscatter communications: Waveform design and snrenergy tradeoff,” to appear in IEEE Commun. Lett.
[5] J. Qian, F. Gao, G. Wang, S. Jin, and H. Zhu, “Noncoherent detections for ambient backscatter system,” IEEE Trans.
Wireless Commun., vol. 16, no. 3, pp. 1412–1422, Mar. 2017.
[6] P. N. Alevizos, A. Bletsas, and G. N. Karystinos, “Noncoherent short packet detection and decoding for scatter radio sensor
networking,” IEEE Trans. Commun., vol. 65, no. 5, pp. 2128–2140, May 2017.
[7] W. Liu, K. Huang, X. Zhou, and S. Durrani, “Full-duplex backscatter interference networks based on time-hopping spread
spectrum,” IEEE Trans. Wireless Commun., vol. 16, no. 7, pp. 4361–4377, Jul. 2017.
[8] K. Han and K. Huang, “Wirelessly powered backscatter communication networks: Modeling, coverage, and capacity,”
IEEE Trans. Wireless Commun., vol. 16, no. 4, pp. 2548–2561, Apr. 2017.
[9] B. Zhen, M. Kobayashi, and M. Shimizu, “Framed ALOHA for multiple RFID objects identification,” IEICE Trans.
Commun., vol. E88-B, no. 3, pp. 991–999, Mar. 2005.
[10] M. T. Isik and O. B. Akan, “PADRE: Modulated backscattering-based passive data retrieval in wireless sensor networks,”
in Proc. IEEE Wireless Commun. Networking Conf. (WCNC), Budapest, Hungary, Apr. 2009, pp. 1–6.
[11] Y. Zheng, M. Li, Y. Zheng, and M. Li, “Read bulk data from computational RFIDs,” IEEE/ACM Trans. Netw., vol. 24,
no. 5, pp. 3098–3108, Oct. 2016.
[12] D. T. Hoang, D. Niyato, P. Wang, D. I. Kim, and L. B. Le, “Optimal data scheduling and admission control for backscatter
sensor networks,” IEEE Trans. Commun., vol. 65, no. 5, pp. 2062–2077, May 2017.
[13] C. Angerer, R. Langwieser, and M. Rupp, “RFID reader receivers for physical layer collision recovery,” IEEE Trans.
Commun., vol. 58, no. 12, pp. 3526–3537, Dec. 2010.
[14] S. Kim, S. Kwack, S. Choi, and B. G. Lee, “Enhanced collision arbitration protocol utilizing multiple antennas in RFID
systems,” in Proc. 17th Asia Pacific Conf. on Commun., Sabah, Malaysia, Oct. 2011, pp. 925–929.

30

[15] T. Demeechai and S. Siwamogsatham, “Using CDMA to enhance the MAC performance of ISO/IEC 18000-6 Type C,”
IEEE Commun. Lett., vol. 15, no. 10, pp. 1129–1131, Oct. 2011.
[16] J. Wang, H. Hassanieh, D. Katabi, and P. Indyk, “Efficient and reliable low-power backscatter networks,” in Proc. ACM
SIGCOMM, Helsinki, Finland, Aug. 2012, pp. 61–72.
[17] J. Chen, U. Yatnalli, and D. Gesbert, “Learning radio maps for UAV-aided wireless networks: A segmented regression
approach,” in Proc. IEEE Int. Conf. Commun. (ICC), Paris, France, May 2017.
[18] N. Farsad and A. Goldsmith, “Detection algorithms for communication systems using deep learning,” arXiv preprint
arXiv:1705.08044, 2017.
[19] J. Barbancho, C. León, F. J. Molina, and A. Barbancho, “A new QoS routing algorithm based on self-organizing maps for
wireless sensor networks,” Telecommun. Syst.,, vol. 36, no. 1-3, pp. 73–83, Mov. 2007.
[20] B. Krishnamachari and S. Iyengar, “Distributed Bayesian algorithms for fault-tolerant event region detection in wireless
sensor networks,” IEEE Trans. Comput., vol. 53, no. 3, pp. 241–250, Mar. 2004.
[21] M. R. Morelande, B. Moran, and M. Brazil, “Bayesian node localisation in wireless sensor networks,” in Proc. IEEE Int.
Conf. Acousti., Speech Signal Process., Mar. 2008, pp. 2545–2548.
[22] J. Kho, A. Rogers, and N. R. Jennings, “Decentralized control of adaptive sampling in wireless sensor networks,” ACM
Trans. Sen. Netw., vol. 5, no. 3, pp. 19:1–19:35, May 2009.
[23] G. Mergen and L. Tong, “Type based estimation over multiaccess channels,” IEEE Trans. on Signal Process., vol. 54,
no. 2, pp. 613–626, Feb. 2006.
[24] B. Chen, L. Tong, and P. K. Varshney, “Channel-aware distributed detection in wireless sensor networks,” IEEE Signal
Processing Mag., vol. 23, no. 4, pp. 16–26, Jul. 2006.
[25] I. Glaropoulos, M. Laganà, V. Fodor, and C. Petrioli, “Energy efficient COGnitive-MAC for sensor networks under WLAN
co-existence,” IEEE Trans. Wireless Commun., vol. 14, no. 7, pp. 4075–4089, Jul. 2015.
[26] Z. Liu and I. Elhanany, “RL-MAC: a reinforcement learning based mac protocol for wireless sensor networks,” International
Journal of Sensor Networks, vol. 1, no. 3-4, pp. 117–124, 2006.
[27] Y.-J. Shen and M.-S. Wang, “Broadcast scheduling in wireless sensor networks using fuzzy Hopfield neural network,”
Expert systems with applications, vol. 34, no. 2, pp. 900–907, 2008.
[28] M. J. Beal, Variational algorithms for approximate Bayesian inference.

University of London United Kingdom, 2003.

[29] W. Che, H. Deng, W. Tan, and J. Wang, “A random number generator for application in RFID tags,” in Networked RFID
Systems and Lightweight Cryptography.

Springer Berlin Heidelberg, 2008, pp. 279–287.

[30] G. K. Balachandran and R. E. Barnett, “A 440-nA true random number generator for passive RFID tags,” IEEE Trans.
Circuits Syst. I, Reg. Papers, vol. 55, no. 11, pp. 3723–3732, Dec. 2008.
[31] M. Fernández and S. Williams, “Closed-form expression for the poisson-binomial probability density function,” IEEE
Trans. Aerosp. Electron. Syst., vol. 46, no. 2, pp. 803–817, Apr. 2010.
[32] M. R. Gupta and Y. Chen, “Theory and use of the EM algorithm,” Foundations and Trends R in Signal Processing, vol. 4,
no. 3, pp. 223–296, 2011.
[33] K. P. Murphy, Machine learning: a probabilistic perspective.

MIT press, 2012.

[34] R. M. Neal and G. E. Hinton, “A view of the EM algorithm that justifies incremental, sparse, and other variants,” in
Learning in graphical models.

Springer, 1998, pp. 355–368.

[35] L. P. Kadanoff, “More is the same; phase transitions and mean field theories,” Journal of Statistical Physics, vol. 137, no.
5-6, pp. 777–797, 2009.
[36] S. L. Simpson, L. J. Edwards, K. E. Muller, P. K. Sen, and M. A. Styner, “A linear exponent AR (1) family of correlation
structures,” Statistics in medicine, vol. 29, no. 17, pp. 1825–1838, 2010.

