1

Optimal Spectrum Sensing Policy in RF-Powered
Cognitive Radio Networks

arXiv:1803.09193v1 [] 25 Mar 2018

Hae Sol Lee, Muhammad Ejaz Ahmed, and Dong In Kim

Abstract—An orthogonal frequency division multiple access
(OFDMA)-based primary user (PU) network is considered, which
provides different spectral access/energy harvesting opportunities
in RF-powered cognitive radio networks (CRNs). In this scenario,
we propose an optimal spectrum sensing policy for opportunistic
spectrum access/energy harvesting under both the PU collision
and energy causality constraints. PU subchannels can have different traffic patterns and exhibit distinct idle/busy frequencies, due
to which the spectral access/energy harvesting opportunities are
application specific. Secondary user (SU) collects traffic pattern
information through observation of the PU subchannels and
classifies the idle/busy period statistics for each subchannel. Based
on the statistics, we invoke stochastic models for evaluating SU
capacity by which the energy detection threshold for spectrum
sensing can be adjusted with higher sensing accuracy. To this end,
we employ the Markov decision process (MDP) model obtained
by quantizing the amount of SU battery and the duty cycle
model obtained by the ratio of average harvested energy and
energy consumption rates. We demonstrate the effectiveness of
the proposed stochastic models through comparison with the
optimal one obtained from an exhaustive method.
Index Terms—OFDMA, RF-powered cognitive radio networks
(CRNs), traffic patterns, spectral access, energy harvesting,
spectrum sensing.

I. I NTRODUCTION
Recently, there have been many proposals on designing
efficient circuits and devices for radio frequency (RF) energy
harvesting suitable for low-power wireless applications [1] [4]. The maximum power available for RF energy harvesting
at a free space distance of 40 meters is known to be 7µW and
1µW for 2.4GHz and 900MHz frequency, respectively [5].
With the RF energy harvesting capability, wireless, especially
mobile device can operate perpetually without periodic energy
replenishment. In cognitive radio networks (CRNs), secondary
users (SUs) equipped with RF energy harvesting device can
opportunistically not only access primary user (PU) channels
but also harvest RF energy carried on PU channels through
spectrum sensing. Hence, selecting PU channel for harvesting
or transmitting through accurate spectrum sensing is a crucial
component for SUs to achieve an optimal performance.
Most works on RF energy harvesting in CRNs rely on a
predefined assumption on PU channel idle time distribution [6]
- [15] and focus on optimizing SU spectral access based on the
battery level. Exploiting multiple channels will offer different
idle channel distributions and PU’s signal strengths. Therefore,
it will give SU more chances to choose between transmitting data and harvesting energy, which in turn improves the
transmission and energy harvesting efficiency, as demonstrated
in [16] and [17]. However, in practice, the channel idle time
distribution depends on specific traffic patterns carried over the

PU channel [18]. Therefore, it is of paramount importance for
SU to be aware of PU traffic patterns so that it can adapt its
harvesting/transmission strategies accordingly. The challenge,
however, lies in efficient classification of PU traffic patterns
based on the applications’ fingerprints (traffic features).
Existing solutions for the traffic patterns identification fall
into the following three categories: 1) port-based, 2) signaturebased, and 3) deep packet inspection. But these approaches
do not perform well under the dynamic nature of the traffic
patterns. In this paper, we propose Dirichlet process mixture
model (DPMM) to efficiently classify various applications
(traffic patterns). The DP is a family of Bayesian nonparametric (BNP) models which are mostly used for density
estimation, clustering, model selection/averaging, etc. The DPs
are nonparametric which means that the number of hidden
traffic applications is not known in advance. Due to the
nonparametric nature of these models, they do not require
the number of clusters (applications) to be known a priori.
Moreover, such models can adapt dynamically over time as the
number of traffic patterns grows. The proposed DPMM traffic
classification is unsupervised, and based only on observations
without any control overhead.
Based on the classification of PU traffic patterns, we propose
an optimal spectrum sensing policy for opportunistic spectral
access/energy harvesting in RF-powered CRNs. Towards this,
SU collects the information on distinct traffic patterns through
observation of the PU subchannels and identifies the PU traffic
patterns by classifying the received data packets into distinct
features. Following the DPMM approach developed in [19],
[20], we can classify the PU subchannels with idle/busy period
statistics which will be used for optimal spectral access/energy
harvesting. For this, we need to obtain the appropriate energy
detection thresholds associated with distinct traffic patterns, so
as to establish the optimal spectrum sensing policy with higher
sensing accuracy.
Suppose the energy detection threshold is low, then SU will
likely identify the PU subchannel as busy due to noise/cochannel interference. Hence, the probability of the PU subchannel being recognized as idle is low, even if the PU
subchannel in fact is idle, resulting in less transmission opportunity for SU. On the other hand, if the threshold is high,
the probability of the PU subchannel being recognized as
busy is low, causing SU to transmit aggressively. This will
result in collision of PU and SU transmissions, reducing their
transmission efficiency, and also incur the energy depletion
of SU device due to its frequent transmissions. Therefore,
we invoke two stochastic models for evaluating SU capacity,
and then derive an optimal energy detection threshold for

2

spectrum sensing to maximize the SU capacity. To this end, we
employ the Markov decision process (MDP) model obtained
by quantizing the amount of SU battery and the duty cycle
model obtained by considering the average energy harvesting
and energy consuming rates. We confirm that the SU capacity
obtained from two stochastic models can be close to the
optimal capacity obtained from an exhaustive method.
The rest of the paper is organized as follows. Section II
describes the system model along with the key assumptions.
In Section III, we give an overview of traffic classification
based on the DPMM approach. In Section IV, we formulate an
optimization problem using the duty cycle model for deriving
an optimal energy detection threshold, while in Section V we
propose a stochastic model based on the MDP for SU with
some statistic information about PU subchannels. Finally, the
performance obtained from analysis is examined in Section
VI, and concluding remarks are given in Section VII.
II. S YSTEM M ODEL
We consider a RF-powered CRN as shown in Fig. 1, where
SU is equipped with RF energy harvesting capability and
performs opportunistic transmission or energy harvesting by
accessing the corresponding PU subchannel. We assume a
PU network which employs orthogonal frequency division
multiple access (OFDMA) with a total of Nc subchannels and
synchronous time-slot based communication across PUs. Here,
the frequency band is divided into several non-overlapping
narrow frequency subbands assigned to different PUs.

Fig. 1. Illustration of cognitive radio network with energy harvesting.

We assume the subchannel of each PU shows independent
idle/busy time statistics, varying with K traffic sources. If
SU can identify K traffic sources on PU subchannels, it
will increase opportunities for SU transmission and energy
harvesting. To identify K traffic sources and classify their
traffic patterns, we consider the three features, such as the
packet length, packet interarrival time, and variance in packet
length. The three features are observed by SUs by inspecting
the packet header from PU traffic. It is assumed that SUs
collaborate with each other in sharing these features via a
common control channel.

Traffic classification is important in the problem under
consideration. Each application follows a unique pattern, and
recognizing this pattern is important for the following two
reasons: 1) we can predict how frequently the energy arrives
(packet interarrival). 2) we can predict for how long the
subchannel is occupied. For each traffic application those
values are different, so estimating those values would save us
energy by avoiding undue spectrum sensing because in that
case we do not know the pattern of the traffic application. If
traffic classification is removed, we would waste undue energy
by blindly sensing channels. In our approach we can use this
energy to transmit data instead.
We assume the packet arrival rate on subchannel c follows
a random process with mean λc . After clustering subchannels,
SU identifies the harvest subchannel, denoted by ch , for energy
harvesting such that λch is maximum. Similarly, the transmit
subchannel, denoted by ct , is identified for transmission such
that λct is minimum. From the subchannel with maximum
λch , SU can harvest more energy as ch is the subchannel with
frequent energy/packet arrivals. Similarly, for the subchannel
with lowest packet arrivals, i.e., ct , a possible PU collision
is less likely, compared to the other subchannels. We define
the current state (idle/busy) of the selected subchannels for
harvesting and transmission, i.e., Sch ∈ {0(idle), 1(busy)},
Sct ∈ {0(idle), 1(busy)}. Both subchannels coexist unless
all Nc subchannels have the same traffic pattern. If they all
have the same traffic pattern and have the same subchannel
states, the two subchannels are randomly selected. Given the
subchannels ch and ct , we define Pr(Sch = 0) = pci h and
Pr(Sct = 0) = pci t as the probabilities of the corresponding
subchannels for harvesting and transmission being idle.
Note that the energy harvesting can be performed over
multiple PU subchannels, and the channel gain affects largely
the performance of harvesting and transmission, both of which
need to be addressed further. Moreover, the energy harvesting using cooperation among multiple PU subchannels can
increase the rate capacity of SU, especially under the energy
causality in the RF-powered CRN considered herein. Due to
the complexity of determining the optimal energy detection
threshold through the traffic classification, thereby the optimal
sensing policy, if consider such multi-channel cooperation,
in this paper we demonstrate an improvement in the rate
capacity by choosing the best subchannel either for energy
harvesting or for transmission. Such multi-channel cooperation
will be an interesting issue for further extending the framework
developed in this paper to improve the rate capacity of SU with
self-powering.
A. SU Battery Model
The idle and busy probabilities on each PU subchannel for
harvesting and transmission can be estimated through the PU
traffic patterns identification. Here, we assume that SU battery
is charged by the energy harvesting which stores energy into
a rechargeable battery of finite capacity Bmax ∈ R+ for a
non-negative real number R+ . As shown in Fig. 2 below, SU
selects active mode or sleep mode for which the action can be
denoted as at ∈ {0(sleep), 1(active)} in the slot t.

3

subchannel. The PU signal and noise are assumed to be modeled as independent circularly symmetric complex Gaussian
(CSCG) random processes with mean zero and variances σp2
2
and σw
, respectively. Then, from [24], the probabilities of
false alarm Pf (ǫ) and detection Pd (ǫ) for the transmission
subchannel ct are evaluated as
Fig. 2. The frame structure of SU with energy harvesting. at represents the
mode selection of SU in a slot t.

In Fig. 2, we denote a time slot by Tslot where each Tslot
is divided into sensing and transmission time durations. We
represent the sensing time duration with Ts , the transmission
time duration with Tt = Tslot − Ts , and the residual battery
level in a slot t with Bt . For simplicity, we assume that SU
always has data to transmit.
SU selects the mode (active or sleep) according to the
battery level observed in each time slot. Then, by comparing
the residual battery level Bt at current slot t with sufficient
energy for transmission, we define the action as
(
1, Bt ≥ es + et
(1)
at =
0, Bt < es + et
where es is the energy consumption for sensing and et the
energy consumption for transmission.
In active mode, SU performs spectrum sensing with the
energy es = Ts Ps consumed over Ts with sensing power Ps .
We represent the sensing outcome as observations, i.e., ot ∈
{0(idle), 1(busy)}. If SU observation is busy (ot = 1) after
sensing the subchannel ct , SU does not transmit data in the
transmission phase (Tt ) of a slot. On the contrary, if ot = 0,
SU transmits data that consumes et = Tt (Pt /η + Pnc ) in the
transmission phase where Pt is the SU transmit power, Pnc
the non-ideal circuit power, and η (0 < η ≤ 1) the efficiency
of power amplifier [21]. Note that if ot = 1 in active mode,
SU can still harvest energy from the subchannel ch . But we
do not account for this as clustering can secure ct with higher
pci t , which can be ignored here for tractable analysis. On the
other hand, in sleep mode, SU turns off its transceiver (except
an energy harvesting circuit) until next slot begins.
To find the amount of residual energy in battery in each
time slot, the harvested energy Eth and energy consumption
Etc in the slot t are defined as
Eth = g(1 − at )ϕPp Tt Sch

(2)

Etc

(3)

= at [es + (1 − ot )et ]

where Pp is the PU transmit power and g is the subchannel
gain between PU and SU. Here we assume that SU receives
the subchannel gain information via a common control subchannel. Therefore, the residual energy in battery in the next
slot t + 1 can be updated as
Bt+1 = min(Bt − Etc + Eth , Bmax ).

(4)

B. Spectrum Sensing and Duty Cycle Model
The sensing probability varies with specific energy detection threshold and subchannel condition associated with PU

Pf (ǫ) = Pr(ot = 1|Sct = 0, at = 1)
p 

ǫ
−1
Ns
=Q
2
σw
Pd (ǫ|g) = Pr(ot = 1|Sct = 1, at = 1, g)



p
ǫ
= Q  gσ2
− 1  Ns 
2
( σ2p + 1)σw

(5)

(6)

w

where ǫ ∈ R+ is a detection threshold for the energy detector,
Ns denotes the number of samples, and Q(x) Q-function.
In order to analyze a stochastic performance of SU, we need
to find how often the current battery level has enough energy
to transmit, i.e., the probability of SU being active Pr(at =
1). Note that we take into account the available energy to
formulate the SU’s action policy in (7) below. For this we
consider the duty-cycling behavior between active mode and
sleep mode to formulate Pr(at = 1). Thus, we define Pr(at =
1) = Tactive /(Tactive + Tsleep ), where Tactive and Tsleep are
the average times spent in active and sleep modes, respectively.
Assuming that the energy harvested in sleep mode must equal
the energy consumed during active mode, the probability of
SU being active can be expressed as
ρh
Pa (ǫ|g) = P r(at = 1) =
(7)
ρh + ρc (ǫ|g)
ρh = g ξpcoh
ρc (ǫ|g) = es +

n

 o

1 − Pf (ǫ) pci t + 1 − Pd (ǫ|g) pcot et

(8)
(9)

where ξ = ϕPp Tt and pco is the probability of the subchannel
c ∈ {ch , ct } being busy with pco = 1 − pci . In the above, we
have assumed that during sleep mode, SU consumes no energy
but harvests the energy with rate ρh . On the other hand, during
active mode, SU consumes the energy at rate ρc . Hence, the
probability of SU being active based on the duty cycle model
above represents the ratio of the energy harvesting rate to the
sum of the energy harvesting and energy consuming rates.
III. T RAFFIC C LASSIFICATION
Sensing multiple PU subchannels may be required to
identify idle/busy subchannels among them, which carry K
traffic applications (patterns), without any prior information.
This may not lead to a performance gain of SU under
the energy causality in the RF-powered CRN because of
undue sensing time and high energy consumption costs. To
overcome this, we attempt to classify traffic patterns with
observed features and select the best subchannel either for
energy harvesting or for transmission accordingly. In this
paper, we consider the following three traffic features with
regard to application fingerprinting: packet length (pl ), packet

4

interarrival time (pt ), and variance in packet length (∆), which
are sufficient for acquiring the channel state information, as
detailed in [22] and [23]. Then, the packet length vector
Plen is represented as Plen = {pl1 , pl2 , · · · , plN }, where
N is the number of packet length samples. Similarly, Pint
and ∆ are represented as Pint = {pt1 , pt2 , · · · , ptN } and
∆ = {[var(pl )]W1 = 0, [var(pl )]W2 , · · · , [var(pl )]WN },
respectively. Here, [var(pl )]Wn is the temporal variance of
packet length in a window Wn of size n, spanning over
[1, . . . , n](st)th observations.
With regard to some application backgrounds, we provide
more detailed description of each feature as follows:
1) Packet Length: Packet lengths for different traffic payloads are likely to be different. For example, the UDP
packet size is longer, the gaming packet size may vary
depending on game dynamics, and the VoIP packet has
smaller packet lengths to minimize jitter. So we use
the packet length as one feature point for identifying
different traffic applications.
2) Packets Interarrival Time: Packet interarrival times for
different applications also vary depending upon the
requirements of applications. For example, in VoIP, the
inter-arrival time is small to avoid annoying effects
caused by jitter.
3) Variance in Packet Length: The packet lengths may
change in every connection. For example, by investigating real wireless traces in [26], [27], for gaming data,
we have observed that packet lengths vary significantly
during the communication.
Let xn be the feature vector of the nth training feature point,
given by


T

xn = pln , ptn , [var(pl )]Wn , X = x1 , · · · , xN .

The matrix of observations X is observed by SU by examining
the packet header of PU.

‫ܩ‬௢

ߠ௞

ߙ

ߨ௞

λ

Z
‫ݖ‬௡

λ

ܰ

In practice, we do not know how many traffic applications
are active in the PU network, and therefore, we would like to
learn it from the data observed. The Dirichlet processes can
be employed to have a mixture model (DPMM) with infinite
components (applications) which can be viewed as taking the
limit of the finite mixture model for K → ∞. The theoretical
model of data generation for the DPMM is
∼
∼

p(θk )
G

G

∼

DP (α, Go )

A. DPMM Representation
The model defined above is a theoretical one. In order
to realize DPMM models, Stick-breaking process, Chinese
restaurant process, and Polya-urn process are used. Here we
represent the DPMM in (10) using the Stick-breaking process.
Consider two infinite collections of independent random varii.i.d.
i.i.d.
ables Vk ∼ Beta(1, α) and θk∗ ∼ Go for k = {1, 2, . . . }.
The Stick-breaking process of G is given by
πk = Vk

k−1
Y

(1 − Vj )

(11)

j=1

G=

∞
X

πk δθk∗

(12)

k=1

where V = {V1 , V2 , . . .} with V0 = 0. The mixing weights
{πk } are given by breaking a stick of unit length into infinitely
small segments. In the DPMM, the vector π represents the
infinite vector of mixing weights and {θ1∗ , θ2∗ , . . . } are the
atoms which correspond to mixing components. Since Zn is
the cluster assignment random variable to the feature point xn ,
the data for the DPMM is generated as
i.i.d.

௡

Fig. 3. Graphical model of an DPMM. Observations are represented by the
shaded node. Nodes are random variables, edges are dependencies, and plates
are replications.

xn
θk

P∞
∗
where G =
k=1 πk δθk ∼ DP (α, Go ), Go is the base
∗
distribution and δθk is used as a short notation for δ(θ = θk∗ )
which is a delta function that takes 1 if θ = θk∗ and 0
otherwise. θk are the cluster parameters sampled from G where
k ∈ {1, 2, . . . }. The generative distribution p(θk ) is governed
by the cluster parameters θk and is used to generate the
observations xn . Then, the multimodal
probability distribution
P
∗
can be defined as p(xn ) = ∞
k=1 πk p(·|δθk ) which is called
the mixture distribution with mixing weights πk and mixing
components p(·|δθk∗ ). The graphical model of the DPMM is
shown in Fig. 3. α is the scalar hyperparameter of the DPMM
and affects the number of clusters obtained. Zn is the cluster
assignment variable such that the feature point xn belongs
to the kth cluster. Larger the value of α, the more clusters,
while smaller the value of α, the fewer clusters. Note that
the value of α indicates the strength of belief in Go . A large
value means that most of the samples will be distinct and have
values concentrated on Go .

(10)

1) Draw Vk |α ∼ Beta(1, α), k ∈ {1, 2, . . . }.
i.i.d.
2) Draw θk∗ ∼ G0 , k ∈ {1, 2, . . . }.
3) For the feature point xn , do:
(a) Draw Zn |{v1 , v2 , . . . } ∼ Multinomial(π).
(b) Draw xn |zn ∼ p(xn |θz∗n ).
Here we restrict ourselves to the DPMM for which the
observable data are drawn from Normal distribution and where
the base distribution for the DP is the corresponding conjugate
distribution.
B. Inference for DPMM
Since the Dirichlet processes are nonparametric, we cannot
use the EM algorithm to estimate the random variables {Zn }
(which store the cluster assignments) for our DPMM model
in (10) due to the fact that EM is generally used for inference
in a mixture model, but here G is nonparametric, making EM

5

difficult. Hence, in order to estimate these assignment variables
in the paradigm of Bayesian nonparametrics, there exist two
candidate approaches for inferences: First, a sampling-based
approach uses Markov chain Monte Carlo (MCMC) to cluster
the traffic patterns X. The MCMC based sampling (also known
as Gibbs sampling) approach is more accurate in classifying
feature points. The second approach is based on variational
methods which convert inference problems into optimization
problems [20]. The main idea that governs variational inference is that it formulates the computation of marginal or
conditional probability in terms of an optimization problem
that depends on a number of free parameters, i.e., variational
parameters. We discuss both approaches in the following
subsections with their detailed comparison in Table I, which
are confirmed with the simulation results.
1) Collapsed Gibbs Sampling for Traffic Classification: In
our model, the data follow multivariate normal distribution,
Nk (~µk , Σk ), where the parameters are 3-dimensional mean
vector and covariance matrix. The conjugate distributions for
mean vector µ
~ k and covariance matrix Σk are given by
µk ∼ N (~µ0 , Σk /κ0 ) and Σk ∼ Inverse-Wishartν0 (Λ−1
~
0 ),
respectively. The Dirichlet hyperparameters, here symmetric
α/K, encodes our beliefs about how uniform/skewed the class
mixture weights are. The parameters to the Normal times
Inverse-Wishart prior, Λ−1
0 , ν0 , κ0 imply our prior knowledge
regarding the shape and position of the mixture densities.
For instance, ~µ0 specifies where we believe the mean of the
mixture densities are expected to be, where κ0 is the number
of pseudo-observations we are willing to ascribe to our belief.
The hyper-parameters Λ−1
0 , ν0 behave similarly for the mixture
density covariance.
Collapsed Gibbs sampler requires to select the base distribution Go which is a conjugate prior of the generative
distribution p(xn |θz∗n ), in order to solve analytically and be
able to sample directly from p(Zn |Z−n , X). The posterior
distribution under our model is
!
N
Y
P (Z, Θ, π, α) ∝ P (X|Z, Θ)P (Θ|G0 )
P (zn |π)
n=1

· P (π|α)P (α).

(13)

By integrating-out certain parameters, the posterior distribution is given by [10]
P (zi = k|Z−n , X, Θ, π, α) ∝ P (xn |zn , Θ)P (zn |Z−n , α).
(14)
For the first term in the above equation,
we
use
multivariate
i
h
(κn +1)
, since we
µn , Λκnn (ν
Student-t distribution, i.e., tνn −2 ~
n −2)
chose Inverse-Wishart as conjugate prior for Σk and Normal
distribution for µ
~ k , where the second term is called Chinese
restaurant process and is given by

mk

, if k ≤ K+ ,
n
−
1+α
(15)
P (zn = k|Z−n ) =
α

, if k > K+
n−1+α

where Z−n = Z/zn , K + is the number
PNof classes containing
at least one data point, and mk =
n=1 I(zi = k) is the
number of data points in class k.

The steps involved in collapsed Gibbs sampling are enumerated below as:
1) Initialize the cluster assignments {zn } randomly.
2) Repeat until convergence:
(a) Randomly select xn .
(b) Fix all other zn for every n 6= n: Z−n .
(c) Sample zn ∼ p(Zn |Z−n , X) from (15).
(d) If zn > K, then update K = K + 1.
An overall framework for PU traffic pattern classification
and application-specific optimization is shown in Fig. 4.
PU traffic pattern classification
PU network
traffic

Preprocessing
Feature
extraction
DPMM data
generation

Gibbs
sampling

inference
Variational
inference

Traffic classification
Idle/busy distribution
estimation

PU application-specific optimal
harvesting/transmission strategy

PU application-specific optimization
Fig. 4. A flow diagram for classification of PU traffic patterns.

2) Variational Bayes Inference for Traffic Classification:
The main idea that governs variational inference is that it
formulates the computation of marginal or conditional probability in terms of an optimization problem that depends on
the number of free parameters, i.e., variational parameters.
In other words, we choose a family of distributions over the
latent variables with its own set of variational parameters ν,
i.e., q(W1:M |ν). We are interested in finding the setting of
the parameters that makes our approximation q closest to the
posterior distribution. Consequently, we can use q with the
fitted parameters in place of the posterior.
Here we assume N observations, X = {xn }N
n=1 and M
latent variables, i.e., W = {Wm }M
.
The
fixed
parameters
m=1
H could be the parametrization of the distribution over the

6

observations or the latent variables. With the given notations,
the posterior distribution under Bayesian paradigm is
p(W|X, H) = R

p(W, X|H)
.
p(W, X|H)
W

(16)

The posterior density in (16) is in an intractable form (often
involving integrals) which cannot easily be solved analytically.
Thus we rely on an approximate inference method, i.e., variational Bayes inference. Our goal is to find an approximation
of the posterior distribution p(W|X, H) as well as the model
evidence p(X|H). We introduce a distribution q(W) defined
over the latent variables and observe that for any choice of
q(W), the following decomposition holds
log p(X|H) = L(q) + KL(q||p)

(17)

(

(18)

where
)
p(X, W|H)
L(q) = q(W) log
dW
q(W)
)
(
Z
p(W|X, H)
dW.
KL(q||p) = − q(W) log
q(W)
Z

(Eq [log p(Zn |V)] + Eq [log p(xn |Zn )])

n=1

− Eq [log q(V, θ , Z)].
∗

(20)

To maximize the bound, we must find a family of variational
distributions that approximate the distributions of infinitedimensional random measure G, where G is expressed in
terms of V = {V1 , V2 , . . . } and θ ∗ = {θ1∗ , θ2∗ , . . . }. The
factorized family of variational distributions for mean-field
inference can be expressed as
q(v, θ ∗ , z) =

K−1
Y
k=1

1 We

q(vk |ζk )

K
Y

k=1

q(θk∗ |εk )

N
Y

(22)

Since multiple parameter options exist for each latent variable
under the variational distribution, we need to optimize the
bounds in (20) based on ν above.
Optimizing by employing the coordinate ascent algorithm
from [9], we optimize the bounds in (20) with respect to the
variational parameters ν in (22). From the coordinate ascent
algorithm, we acquire the following statistics: the number of
traffic patterns (K) and their corresponding parameters θk∗ .
These statistics are utilized by SU for optimal harvesting and
transmission strategy in the sequel.
C. Comparison of Inference Algorithms
Table I presents the detailed comparison of the two algorithms. The collapsed Gibbs sampling is more accurate
than the variational Bayes inference, but the former has some
limitations such as slow convergence [19] and difficult to
diagnose convergence. We confirm the comparison results, in
terms of accuracy and time complexity, through simulations.
The variational inference is biased (underfitting), whereas the
Gibbs sampling’s bias diminishes as the number of runs for the
Markov chain increases. For non-conjugate prior distributions,
the latter is preferred which is much faster than the former.
The latter is deterministic, which means that we always obtain
the same optimal value, given the same starting value and
an objective function without huge local optima problems.
The latter is quicker since it approximates the posterior using
optimization with free variables.
TABLE I
C OMPARISON OF INFERENCE ALGORITHMS .

log p(X|H) ≥ Eq [log p(V|H)] + Eq [log p(θ ∗ |G0 )]
N
X

ν = {ζ1 , . . . , ζK−1 , ε1 , . . . , εK , ϑ1 , . . . , ϑN }.

(19)

From (19), we see that KL(q||p) is the Kullback-Leibler
divergence between q(W) and the posterior distribution
p(W|X, H), where the KL divergence satisfies KL(q||p) ≥ 0
with equality if and only if q(W) = p(W|X, H), i.e., q(·) is
the true posterior. It follows from (17) that L(q) ≤ log p(X|θ),
in other words, L(q) is a lower bound on p(X|θ). Therefore,
we can maximize the lower bound L(q) by optimizing with
respect to q(W), which is equivalent to minimizing the KL
divergence. We consider a restricted family of distributions
q(W) and then seek the member from this family for which
the KL divergence is minimized.
The latent variables1 for DPMM are the stick lengths, atoms,
and cluster assignment variables: W = {V, θ ∗ , Z} and the
hyperparameters are the scaling parameter α and the parameter
for the conjugate base distribution G0 , H = {α, G0 }. Thus,
the marginal distribution of the data (evidence) lower bound
is evaluated as

+

where q(vk |ζk ) are beta distributions parameterized with ζk ,
q(θk∗ |εk ) are exponential family distributions with natural
parameter εk , and q(zn |ϑn ) are multinomial distributions with
parameter σn . The latent variables in Fig. 3 are governed by
the same distribution, but, following the fully factorized variational variables in the mean-field variational approximation,
given there is an independent distribution for each variable.
Thus, the variational parameters are defined by

q(zn |ϑn ) (21)

n=1

use the variational inference to pick a family of distributions over the
latent variables W with its own variational parameters ν, and then set ν to
render q(·|ν) close to the posterior of interest.

Item

Gibbs sampling

Variational inference

Speed
Biasness
Computational requirement
For non-conjugate dist.
Deterministic
Accuracy
Convergence diagnoses
Converge to true posterior
Inference approach
Approximates

Slower
Not biased
Higher
Not preferred
No
Higher
Difficult
Yes
MCMC sampling
Integrals

Faster
Biased
Lower
Preferred
Yes
Lower
Easy
No
Optimization
Data distribution

IV. O PTIMAL E NERGY D ETECTION T HRESHOLD
E STIMATION
After selecting the subchannels for harvesting/transmission,
SU uses ch to harvest energy and ct to transmit information.
Specifically, the goal is to adjust the detection threshold ǫ of

7

the SU energy detector under the energy causality and PU
collision constraints. For instance, increasing the detection
threshold results in frequent SU transmissions, as a result,
the increased probability of accessing the occupied spectrum,
which may result in collision with PU transmission. Also, it
incurs excessive energy usage which is not good for energyconstrained SU. On the other hand, lowering the detection
threshold ǫ reduces unnecessary sensing and transmission
actions and consequently saves energy for future transmission. However, it decreases the probability of accessing the
unoccupied spectrum, causing the achievable rate loss to SU.
Thus, it is of vital importance to finely tune the SU energy
detection threshold ǫ for optimal sensing subject to the design
constraints, i.e., PU protection and energy causality.
A. Problem Formulation
The achievable rate capacity of SU is defined as C =
W log2 (1 + SN R) for the signal-to-noise ratio (SN R) of the
SU link when the PU transmit subchannel ct of bandwidth W
is idle. Then the SU rate capacity is expressed as
Tt
R (ǫ|g) =
C Pr(at = 1, ot = 0, Sct = 0|g)
Tslot


Tt
(23)
C 1 − Pf (ǫ) Pa (ǫ|g)pci t .
=
Tslot
Here we have used the relation Pr(at = 1, ot = 0, Sct =
0|g) = Pr(ot = 0|at = 1, Sct = 0) Pr(at = 1|Sct =
0, g) Pr(Sct = 0) based on Bays’ rule. The rate capacity
above converges to a specific value due to the Q-function
characteristic of Pf and Pa as the energy detection threshold
ǫ increases.
In this case, however, the PU performance is degraded due
to the collision of PU and SU transmissions as the latter
becomes more aggressive. Therefore, we should put some
constraint on the collision probability, which is evaluated as
Pc (ǫ|g) = Pr(at = 1, ot = 0|Sct = 1, g)


= 1 − Pd (ǫ|g) Pa (ǫ|g).

(24)

Now we can formulate an optimization problem to find an
appropriate value of ǫ, which leads to an optimal spectrum
sensing policy for maximizing the SU rate capacity as
max R (ǫ|g) s.t. Pc (ǫ|g) ≤ P c
ǫ

(25)

where P c is the target probability of collision with which PU
can be sufficiently protected.
B. Distinctions of RF-Powered CRNs from General CRNs
As defined in (7), the probability of action for SU in
RF-powered CRNs is a function of its battery level, given
the energy causality (i.e., self-powering) is applied to SU
for joint opportunistic energy harvesting and transmission,
unlike general CRNs. Therefore, as formulated in (23) - (25),
the optimal sensing policy in the RF-powered CRNs should
take into account the energy state, unlike that in the general
CRNs. For this, we have acquired the channel state information
through the traffic classification developed in this paper, which
is a crucial factor for determining the optimal sensing policy
in the RF-powered CRNs.

C. Optimal Energy Detection Threshold
To find an optimal energy detection threshold
ǫ in (25), we

define an objective function as O(ǫ|g) = 1 − Pf (ǫ) Pa (ǫ|g)
which is affected only by ǫ in (23), and then reformulate the
optimization problem above as
max O (ǫ|g) s.t. Pc (ǫ|g) ≤ P c .
ǫ

(26)

In addition, we define the constraint function as Φ(ǫ, P c ) =
Pc (ǫ|g) − P c to obtain a proper threshold range of [0, ǫc ],
where ǫc is the solution of the following equation:
Φ(ǫc , P c ) = 0.

(27)

We will find the following Propositions 1 and 2 useful in
obtaining the solution ǫc of (27).
Proposition 1. The collision probability Pc (ǫ|g) in (24) and
the object function O(ǫ|g) in (26) can be classified into two
types of function f (ǫ|g) as follows:
1) f (ǫ1 |g) < f (ǫ2 |g) for ǫ1 < ǫ2 ≤ ǫm ,
f (ǫ3 |g) > f (ǫ4 |g) for ǫm ≤ ǫ3 < ǫ4 ,
limǫ→∞ f (ǫ|g) = γ1 where ǫm = arg maxǫ f (ǫ|g).
2) f (ǫ1 |g) < f (ǫ2 |g) for ǫ1 < ǫ2 ,
limǫ→∞ f (ǫ|g) = γ1
where
ρh
γ1 =
.
(28)
ρh + e s + e t
Proof. See Appendix A.
Proposition 2. The constraint function in (27) yields a unique
solution ǫc for γ1 > P c . If γ1 ≤ P c , we have two or no
solution.
According to Proposition 2, we have the constraint range
[0, ǫc ] if γ1 > P c and otherwise, the constraint range is [0, ∞).
Following the IEEE 802.22 WRAN, if P c = 0.1, we consider
γ1 ≤ P c as an extreme case corresponding to an extremely
low energy harvesting rate, namely ρh ≪ (es + et ) in (28).
Thus, we may choose the constraint range to be [0, ǫc ] as the
only possible option.
To find ǫc , we resort to the secant method as a root-finding
algorithm where the constraint function can be approximated
by a secant line through two points of the function. Starting
with the two initial iterates ǫ0 and ǫ1 , the next iterate ǫ2 is obtained by computing the value at which the secant line passing
through the two points (ǫ0 , Φ(ǫ0 , P c )) and (ǫ1 , Φ(ǫ1 , P c )) as
Φ(ǫ1 , P c ) − Φ(ǫ0 , P c )
(ǫ1 − ǫ0 ) + Φ(ǫ1 , P c ) = 0,
ǫ2 − ǫ1

(29)

which yields the solution
ǫ2 = ǫ1 − Φ(ǫ1 , P c )

ǫ1 − ǫ0
.
Φ(ǫ1 , P c ) − Φ(ǫ0 , P c )

(30)

Hence, we can derive the recurrence relation as
ǫk =

ǫk−2 Φ(ǫk−1 , P c ) − ǫk−1 Φ(ǫk−2 , P c )
.
Φ(ǫk−1 , P c ) − Φ(ǫk−2 , P c )

(31)

Since the value of P c is small, a solution can be obtained
quickly by setting ǫ0 and ǫ1 close to zero. Then, we iterate

8

until |ǫk − ǫk−1 | becomes very small, which is described in
Algorithm 1. With the constraint rage [0, ǫc ] fixed, we should
be able to find an optimal energy detection threshold which
maximizes the objective function O(ǫ).
Algorithm 1 Optimization algorithm
1: Initialize ǫ0 and ǫ1 .
2: for k = 1, 2, 3, ... do
3:
Update ǫk using the recurrence relation (31)
4:
if |ǫk − ǫk−1 | is sufficiently small then
5:
ǫc = ǫk return ǫc
6:
end if
7: end
8: if ∇O(ǫc ) > 0 then
9:
ǫc is the optimal solution.
10: else
11:
Initialize ǫ0 = ǫc .
12:
for k = 1, 2, 3, ... do
13:
Update ǫk using the recurrence relation (32)
14:
if |ǫk − ǫk−1 | is sufficiently small then
15:
ǫ∗ = ǫk return ǫ∗
16:
end if
17:
end
18: end if

tion threshold play crucial role in obtaining the optimal sensing
policy for SU.
PU traffic pattern classification
Bayesian nonparametric
subchannel clustering

PU network
traffic

ܿ௧ and ܿ௛

‫ܤ‬௧ ൒ ݁௦ ൅ ݁௧

Yes

No

Spectrum sensing ܿ௧ with ߳
obtained by Algorithm 1

Update using (4)

Harvest energy
from ܿ௛

‫݋‬௧ ൌൌ Ͳ

We use the gradient descent method to find the maximum
value in the constraint range with the recurrence relation as

Yes

ǫk = ǫk−1 + β ∇O(ǫk−1 )

(32)

for the step size β. We use a fixed value for β to avoid the
complication of calculation and find an optimal one. In the first
step of the recurrence, we assume ǫ0 = ǫc . From Proposition 1,
if ∇O(ǫc ) > 0, ǫc is an optimal point as the objective function
is an increasing function in [0, ǫc ]. If ∇O(ǫc ) < 0, we continue
the recurrence process to find the maximum point given the
objective function is a type-1 function with unique maximum
point. Algorithm 1 for finding the optimal ǫ∗ is stated above.
Fig. 5 illustrates a whole process for subchannel clustering,
spectral access and energy harvesting. The three traffic features
are used to classify the PU traffic patterns through the BNP
subchannel clustering. Based on the obtained idle and busy period statistics from the output of MCMC (Gibbs sampling), we
find the corresponding subchannels ct and ch with maximum
energy harvesting and transmission probabilities, respectively.
If the residual energy in SU battery Bt is less than the
required energy for transmission, ch is selected to harvest
energy and otherwise, ct to transmit data. Then, SU senses
the selected subchannel using the optimal energy detection
threshold obtained from the gradient descent method above.
After sensing, if the sensing result is ot = 0, SU transmits
data and otherwise, turns off the transmission. The residual
energy in battery is then updated using (4). To enable this, SU
has to acquire information of pci h and pci t from the subchannel
clustering, which in turn influences the energy harvesting
and consuming rates ρh and ρc , respectively. Therefore, the
statistics information obtained through the accurate clustering
process and resulting sensing parameter of the energy detec-

Transmit
PU application-specific optimization
Fig. 5. A flow diagram for spectral access and energy harvesting.

V. M ARKOV D ECISION P ROCESS F ORMULATION
Unlike the duty-cycling behavior between active mode and
sleep mode, we introduce the MDP model to accurately predict
the actions of SU which vary with the evolution of the residual
energy in battery. Since idle/busy time distributions (to estimate harvesting and transmission opportunities) are obtained
from clustering traffic patterns, we need to take decision in
a real time. Even if we know the distribution parameters, we
still need to automate the decision process because the SU harvesting and transmission is a real-time process. The parameters
obtained would help us to achieve early convergence. For this
purpose, the MDP which is Markov-chain based approach,
is incorporated given idle/busy time distributions for traffic
applications. Therefore, the MDP renders a good solution for
online/real-time decision making.
A. State Space of Battery
SU decides whether to harvest energy or transmit based
on the battery level. The event of at = 1 means that the
amount of the residual energy in battery is sufficient enough
to transmit. To find the probability of having this event, we use

9

the MDP formulation assuming the state is the discretization of
the battery capacity, and evaluate the steady-state probabilities
of the battery level with sufficient energy to transmit. We
discretize thejcurrentk residual energy in battery bt in Nb levels
denotes the maximum amount of energy
where Nb = Bmax
eq
quanta that can be stored in battery. Here, one energy quantum
t)
corresponds to eq = (esn+e
where nτ represents the number
τ
of states that enter harvesting mode. In general, if Nb is
sufficiently high, the discrete model can be considered as a
good approximation of the continuous one. Then, (4) can be
rewritten in terms of energy quanta as
bt+1 = min(bt − ect + eht , Nb )
(33)
j hk
l cm
E
E
where eht = eqt and ect = eqt . Here, the floor is used to
have a conservative harvesting performance, while the ceiling
to assure a required energy consumption. Thus, the worst case
of the battery level is assumed.

B. Transition Probability Matrix

1

•• − 1

•• + 1

••

!

−1

Fig. 6. The battery state transition with nτ energy harvesting states for Nb state Markov chain model.

We define the transition probability matrix U for the Nb state Markov chain model in Fig. 6 as
h

U(ǫ|g) = Uh (ǫ|g) | Ua (ǫ|g)



u0,0
u1,0
..
.







 unτ −1,0

Uh (ǫ|g) = 
 un ,0
τ


 un +1,0
τ


..

.
uNb −1,0

0
u1,1
..
.
..
.
..
.
..
.
..
.
...

...
...
..
.
0
0
0
0
0
..
.

0
0
0
uNb −1−nτ ,Nb −1
0
0
0
0

...

0
uNb −1,Nb −1











 . (36)









Here, Uh (ǫ|g) denotes the Nb × nτ matrix whose current
battery level induces SU to enter harvesting mode, while
Ua (ǫ|g) the Nb × (Nb − nτ ) matrix whose current battery
level induces SU to enter active mode. The components of
these matrices are defined as
ui,i = 1[0≤E h <eq ] pcot + pci t , (i < nτ )
t
ui,j = 1[(i−j)eq ≤E h <(i−j+1)eq ] pcot ,

(37)

t

In the Markov chain model with Nb states as shown in
Fig. 6, the harvesting state i ∈ {0, 1, ..., nτ − 1} changes
to state j (j ≥ i) through energy harvesting as the current
battery level is insufficient for transmission. The active state
i ∈ {nτ , nτ + 1, ..., Nb − 1} with sufficient energy to transmit
will change to (i − nτ )-state or fail to transfer and return to
i-state again.

0

Ua (ǫ|g) =

u0,n
0

0
u
1,n
τ +1



0
0


0
0


0
0


0
0

 unτ −,nτ
0


0
unτ +1,nτ +1



0
0
0
0

...
...
..
.
..

.

0
unτ −1,nτ −1

..

.

unτ ,nτ −1

..

.

unτ +1,nτ −1
..
.

..

.
. . . uNb −1,nτ −1

ui,i = Pf (ǫ) pci t + Pd (ǫ|g) pcot , (n ≤ i ≤ Nb − 1).
(40)
Here, ui,i (i < nτ ) is the probability that energy harvesting
is successful when ch is busy but there is insufficient energy
to reach a higher battery level, or ch is idle. ui,j (0 < j <
nτ , j < i < Nb ) is the one that energy harvesting is successful
when ch is busy and the state changes from j to i. ui,i+n
(i < nτ ) is the one that ct is idle with no false alarm and
successful transmission, or ct is busy with missed detection
and collision. ui,i (n ≤ i ≤ Nb − 1) is the one that ct is idle
with false alarm, or ct is busy with no missed detection.
C. Steady-State Probability and Optimal Threshold Algorithm
We define the steady-state probability vector of the Nb -state
Markov chain as Π = [π0 , π1 , ..., πNb −1 ], where Π is the left
eigenvector of U(ǫ) corresponding to the unit eigenvalue as
ΠU(ǫ|g) = Π.

iT
0
0

ui,i+nτ

(0 < j < nτ , j < i < Nb )
(38)
 ct

 ct 
= 1 − Pf (ǫ) pi + 1 − Pd (ǫ|g) po , (i < n)
(39)

(34)

To derive the steady-state probability vector Π, we need to
make the necessary assumption below.
Assumption 1. The maximum energy quanta Nb should be
sufficient enough to satisfy the following conditions:



















(41)

(35)

eht + nτ − 1 ≤ Nb − 1
 h
Et
≤ N b − nτ
eq
Eth
< N b − nβ + 1
eq


gϕPp Tt
Nb >
+ 1 nτ − 1.
es + et

(42a)
(42b)
(42c)
(42d)

10

Assumption 1 implies that the maximum state Nb − 1 must be
greater than the sum of the harvesting energy quanta eht and
the maximum number of harvesting state nτ − 1. This means
that the maximum state number should always be greater than
the maximum allowable state due to harvesting. Thus, (33)
can be rewritten as
bt+1 = bt − ect + eht .

VI. R ESULTS
The simulation results for the proposed duty cycle and
MDP based stochastic models are also presented. To show
the effectiveness of the proposed scheme by comparing how
close to actual capacity, we define
β=

(43)

With Assumption 1, we define the number of energy quanta
charged through energy harvesting as

nκ = nκ ∈ N ∩ {0} nκ eq ≤ Eth < (nκ + 1)eq


gnτ ϕPp Tt
.
(44)
=
es + et
Proposition 3. Then, the steady-state probability vector Π can
be evaluated as

τ

(0 ≤ i < nτ )
 nα α+nτ τ ,
κ
πi = nκ κ+nτ τ ,
(45)
(nτ ≤ i < nτ + nκ )


0,
(nτ + nκ ≤ i < Nb )

where

κ = pcoh




τ = 1 − Pf (ǫ) pci t + (1 − Pd (ǫ|g) pcot .

(46)
(47)

Proof. See Appendix B.

Using the steady-state probability vector Π obtained by
Proposition 3, the probability of SU entering active mode
PaM (ǫ|g) can be derived as

NX
nκ
b −1
nκ κ
nτ κ
M
πi =
Pa (ǫ|g) =
= nκ 
.
(48)
nκ κ + nτ τ
nτ κ + τ
i=n

Nt
Tt
1 X
C at (1 − ot )(1 − Sct )
Nt t=1 Tslot

(52)

for the actual capacity obtained from simulation based on the
Monte-Carlo method where Nt is the number of simulation
iterations using the exhaustive search for energy detection
threshold. We use real wireless traces available online [26]
[27] for 3G network. We utilize three sources (UDP, VoIP,
Game) in our data set. Unless otherwise stated, the values
of the parameters used here are listed in Table II, which are
mainly drawn from [25].
TABLE II
S IMULATION PARAMETERS
Symbol
W
Tslot
Ts
Tt
B0
Ps
Pt
Pnc
η
ϕ
2
σp2 /σw
Pc
Bmax

Definition
Bandwith
Slot duration
Sensing duration
Transmission duration
Initial energy
Sensing power
Transmit power
Non-ideal circuit power
Efficiency of power amplifier
Energy harvesting efficiency
SNR of PU signal at SU transmitter
Target probability of collision
Maximum capacity of battery

Value
1MHz
10ms
2ms
98ms
0J
110 mW
50 mW
115.9 mW
-5.65 dB
0.2
-10 dB
0.1
1 mJ

τ

If nτ is sufficiently large, we can approximate
in (25), which yields
PaM (ǫ|g) =

gϕPp Tt κ
.
gϕPp Tt κ + (es + et )τ

nκ
nτ

∼
=

ϕPp Tt
es +et

(49)

To find the optimal energy detection threshold, we define
the MDP objective function using PaM (ǫ|g) obtained by the
MDP as OM (ǫ) = (1 − Pf (ǫ))PaM (ǫ|g), and then express the
optimization problem again as
max OM (ǫ|g) s.t. Pc (ǫ|g) ≤ P c .
ǫ

(50)

Proposition 4. The MDP objective function OM (ǫ|g) can be
classified into two types of function f M (ǫ|g) as follows:
1) f M (ǫ1 |g) < f M (ǫ2 |g) for ǫ1 < ǫ2 ≤ ǫm ,
f M (ǫ3 |g) > f M (ǫ4 |g) for ǫm ≤ ǫ3 < ǫ4 ,
limǫ→∞ f (ǫ|g) = γ2 where ǫm = arg maxǫ f M (ǫ|g).
2) f M (ǫ1 |g) < f M (ǫ2 |g) for ǫ1 < ǫ2 ,
limǫ→∞ f M (ǫ|g) = γ2
where
gϕPp Tt κ
.
(51)
γ2 =
gϕPp Tt κ + es + et
Proof. See Appendix C.
By Proposition 4 we can optimize in the same way as
Algorithm 1 in Section IV.

A. Traffic-Awareness via Clustering
Before evaluating the performance of the proposed scheme,
we confirm the performance of traffic patterns clustering that
results from the proposed MCMC based sampling and the
variational inference. In Fig. 7, we compare the accuracy of
the MCMC based sampling, variational inference, and Kmeans (as baseline) clustering algorithms when data points
(observations) are generated from 3 different (UDP, VoIP,
Game) traffic patterns being mixed. The K-means algorithm
is one for grouping a given data into K clusters by minimizing
the dispersion of the distance between each clusters. Unlike
the other approaches, the K-means algorithm cannot estimate
the number of clusters, and it should be performed only
by assuming a fixed number of traffic sources. We see that
the Bayesian approaches offer higher accuracy than the Kmeans algorithm. This is because the Bayesian approaches are
to approximate a prior probability and a likelihood function
derived from a statistical model for the observed data, whereas
the K-means considers only the differences in observed traffic
values. In the variational inference, we observe some errors
compared to the MCMC method because we approximate the
latent variables assumed by the mean-field theory. In Fig. 8,
we compare the two Bayesian approaches in terms of their
elapsed times. We notice that the elapsed times increase as the

11

MCMC (Gibbs Sampling)
Variational Inference
K−Means (Parametric with K=3)

Accurately identficating K (%)

100

90

80

70

60

50
90

150

300
450
Data points (mixed with different traffics)

600

750

Fig. 7. Clustering accuracy of the proposed MCMC and variational inference,
and K-means algorithms.
Fig. 9. SU capacity versus SU transmit power Pt and energy detection
c
c
threshold ǫ. (top: VoIP with pi h = 0.2 bottom: Game with pi h = 0.5)

100
MCMC (Gibbs Sampling)
Variational Inference

Elapsed time (sec)

80

0.1

60

Optimal Point (Simulation)
0.12

Optimal Point (Theory: Duty Cycle)

0.4

Optimal Point (Theory: MDP)

40

0.14
0.35

20

0
90

150

300
450
Data points (mixed with different traffics)

600

750

SU transmit power Pt (W)

0.16

Fig. 8. Elapsed time of the MCMC and variational inference algorithms.

0.3
0.18
0.25
0.2

0.2
0.22

0.24

0.15

0.26

0.1

0.28
0.05
0.95

number of data points increases, and it is confirmed that the
variational inference shows less elapsed time than the MCMC
method. Hence, if we can derive a set of equations used to
iteratively update the parameters well, the former converges
faster than the latter requiring a large amount of sampling
work.

1

1.05

Energy detection threshold ε

1.1

1.15

1.2

(a)
0.1
Optimal Point (Simulation)
0.12

Optimal Point (Theory: Duty Cycle)

0.4

Optimal Point (Theory: MDP)
0.14
0.35

B. SU Optimal Sensing Threshold
Fig. 9 shows the SU achievable rate capacity based on the
duty cycle model with varying energy detection threshold ǫ and
SU transmit power Pt . Note that the probability of the harvest
subchannel ch being idle was measured as pci h = (0.2, 0.5) for
the (VoIP, Game) traffic applications, respectively. We observe
an interesting trade-off in choice of the variables ǫ and Pt .
For small ǫ, the false alarm probability increases, resulting
in low transmission probability. To the contrary, for large ǫ,
it decreases and SU more likely transmits data if its residual
energy is enough for transmission. We notice that the capacity
converges to a specific value proportional to γ1 in (28), which
is the ratio of the energy harvesting rate to the sum of the
energy harvesting and consuming rates. Hence, an optimal ǫ
balances the sensing accuracy trade-off. Likewise, there is the

SU transmit power Pt (W)

0.16
0.3
0.18
0.25
0.2

0.2
0.22

0.24

0.15

0.26

0.1

0.28
0.05
0.95

1

1.05

Energy detection threshold ε

1.1

1.15

1.2

(b)
Fig. 10. An optimal point of SU transmit power Pt and energy detection
c
threshold ǫ for (VoIP, Game) traffic with ((a) VoIP pi h = 0.2 (b) Game
ch
pi = 0.5)

12

0.9

0.8

SU capacity (bits/s/Hz)

0.7

Opimal
MCMC + Duty Cycle
Variational Inference + Duty Cycle
MCMC + MDP
Variational Inference + MDP

0.6

0.5

0.4

0.3

0.2

0.1

0
0.9

0.95

1

1.05
1.1
Energy detection threshold ε

1.15

Fig. 11. SU rate capacity versus energy detection threshold ǫ when Pt =0.2W.

0.9

0.85

SU capacity (bits/s/Hz)

energy causality trade-off in the SU transmit power. For large
Pt , the probability of SU being active decreases, while for
small Pt , it increases but the SNR decreases. We see the VoIP
traffic offering better performance than the Game traffic, as
the former yields higher harvesting rate with continuous and
short intervals between voice packets. However, the latter with
pci h = 0.5 shows the packet intervals changing dynamically,
resulting in low harvesting rate.
Fig. 10 shows an optimal point of the actual capacity in (52),
the duty cycle and MDP models for (VoIP, Game) traffic. The
black line shows an optimal ǫ obtained from (26) according
to Pt , where the optimal ǫ decreases slightly as Pt increases.
This is because a tight energy causality due to the increased
Pt requires less transmission opportunity. We can see that the
optimal point of VoIP traffic has a larger value of Pt than that
of Game traffic. In VoIP traffic, the increased Pt results in
low transmission opportunity but increases the SNR of SU,
and a small value of pci h guarantees the energy causality. It
means that VoIP traffic subchannel idle/busy statistics, which
show higher harvesting opportunity than Game traffic, offset
the tight energy causality due to the increased Pt . The optimal
point of the duty cycle and MDP models offer almost the same
performance as the actual capacity from simulation. In Figs.
9 and 10, it is evident that the subchannels carrying distinct
traffic patterns exhibit different harvesting rates, and hence the
appropriate values of ǫ should be determined considering both
the energy causality and PU collision constraints.

Opimal
MCMC + Duty Cycle
Variational Inference + Duty Cycle
MCMC + MDP
Variational Inference + MDP

0.8

0.75

0.7

C. SU Performance by Clustering Algorithms

0.65
0.1

0.12

0.14

0.16

0.18
0.2
0.22
SU transmit power P (W)

0.24

0.26

0.28

t

We evaluate the SU achievable rate capacity by the proposed
clustering algorithms with respective threshold optimization.
We assume 3 different (UDP, VoIP, Game) traffic sources
with 10 subchannels, respectively. In this setting, SU selects
the subchannels ct and ch from Nc = 30 subchannels with
450 data points generated by using MCMC and variational
inference, respectively. Then, the optimal energy detection
threshold for the rate capacity is determined based on the duty
cycle and MDP models. Fig. 11 shows the SU achievable
rate capacity for varying energy detection threshold ǫ. The
black line represents actual capacity using ǫ obtained from
(26) with accurate clustering information. We see that the
MCMC is closer to the optimal line than variational inference,
while the duty cycle and MDP models offer almost the same
performance in optimizing ǫ. This clearly shows the higher
sensing accuracy in selecting ct and ch subchannels of the
MCMC than variational inference.
Fig. 12 shows the SU achievable rate capacity obtained by
using the optimal value ǫ∗ for varying SU transmit power Pt .
We see that the variational inference reaches maximum when
Pt = 0.24W , but the optimal actual capacity does slightly later
when Pt = 0.26W , like the MCMC. It means that accurate
clustering information leads to higher energy harvesting rate,
which allows SU to increase the residual energy in battery.
Hence, SU can increase the maximum achievable rate capacity
with higher transmit power Pt .

Fig. 12. SU rate capacity versus SU transmit power Pt with optimal ǫ∗ .

VII. C ONCLUSION
We have proposed an optimal spectrum sensing policy for
maximizing the SU capacity in OFDMA based CRNs which
is powered by energy harvesting. SU collected traffic pattern
information through observation of PU subchannels and classified the idle/busy period statistics for each PU subchannel
using the MCMC and variational inference algorithms. Based
on these statistics, we developed the stochastic SU capacity
models which are the duty cycle based one defined by the
times spent in active and sleep mode, and the MDP model
based on the evolution of the residual energy in battery. The
energy detection threshold was optimized to maximize the SU
capacity while satisfying the energy causality and PU collision constraints according to traffic patterns. We have shown
the performance trade-off of the BNP subchannel clustering
algorithms by comparing the accuracy and elapsed time of
algorithms. It was shown that SU can optimize the stochastic
model by selecting the threshold referring to the idle/busy
period statistics of PU subchannels. It was also shown that the
proposed duty cycle and MDP model achieve similar capacity
to that of the actual capacity from simulation based on the
Monte-Carlo method.

13

ACKNOWLEDGMENT
This work was supported by the National Research Foundation of Korea (NRF) Grant funded by the Korean Government
under Grant 2014R1A5A1011478.
A PPENDIX A
P ROOF OF PROPOSITION 1
We define A(ǫ|g) = [1 − Pd (ǫ|g)] and B(ǫ|g) = Pa (ǫ|g)
for the proof of the collision probability Pc (ǫ|g) in (24).
df (ǫ|g)
dA(ǫ|g)
dB(ǫ|g)
=
B(ǫ|g) + A(ǫ|g)
.
(53)
dǫ
dǫ
dǫ
The increasing function A(ǫ|g) and the decreasing function
B(ǫ|g) satisfy the following properties:
lim A(ǫ|g) < lim B(ǫ|g)

(54)

lim A(ǫ|g) > lim B(ǫ|g)

(55)

ǫ→0

ǫ→0

ǫ→∞

ǫ→∞

A(ǫe |g) = B(ǫe |g)
(56)
dB(ǫ|g)
dA(ǫ|g)
>−
.
(57)
dǫ
dǫ
The same results can be obtained even if A(ǫ) = [1 − Pf (ǫ)]
for the proof of the objective function O(ǫ|g) in (26).
Proposition 5. If 0 < ǫ < ǫe ,

df (ǫ|g)
dǫ

> 0.

Proof. We have A(ǫ|g) < B(ǫ|g) from (54) and (56), and
using (57), we have the following inequalities:


dB(ǫ|g)
dA(ǫ|g)
(58a)
B(ǫ|g) > B(ǫ|g) −
dǫ
dǫ


dB(ǫ|g)
dA(ǫ|g)
(58b)
B(ǫ|g) > A(ǫ|g) −
dǫ
dǫ
df (ǫ|g)
> 0.
(58c)
dǫ
Proposition 6. If ǫe < ǫ and
where ∆ ∈ R+ .

df (ǫ|g)
dǫ

< 0,

df (ǫ+∆|g)
dǫ

< 0

Proof.
dA(ǫ + ∆|g)
dA(ǫ|g)
<
dǫ
dǫ
dA(ǫ+∆|g)
dǫ

A(ǫ|g)
dA(ǫ+∆|g)
dǫ

A(ǫ + ∆|g)

<
<

dA(ǫ|g)
dǫ

A(ǫ|g)
dA(ǫ|g)
dǫ

A(ǫ|g)

.

(59a)
(59b)
(59c)

Inequality (59a) is satisfied by ǫi < ǫe where ǫi is the inflection
point of A(ǫ|g), and (59c) is satisfied by the increasing
function property A(ǫ + ∆|g) > A(ǫ|g)
dA(ǫ|g)
dǫ

A(ǫ|g)
dB(ǫ+∆|g)
dǫ

B(ǫ|g)
dA(ǫ+∆|g)
dǫ

A(ǫ + ∆|g)

Inequality (60a) is satisfied by ǫe < ǫ, and (60b) is satisfied
by the decreasing function property B(ǫ + ∆|g) < B(ǫ|g).
Finally, from (59c), the inequality (60c) holds.
(ǫ|g)
= 0 has one
From Propositions 5 and 6, the equation df dǫ
solution or no solution, so that the function f (ǫ|g) converges
to γ1 = (limǫ→0 A(ǫ|g))(limǫ→0 B(ǫ|g)).

A PPENDIX B
P ROOF OF PROPOSITION 3
From (41), we have
Π (I − U(ǫ|g)) = 0
T

(61a)

T

(61b)

T

(61c)

(I − U(ǫ|g)) Π = 0
G(ǫ|g)Π = 0.
T

We define G(ǫ|g) = (I − U(ǫ|g)) . If Assumption 1 holds,
then G(ǫ|g) is given by
h
i
G(ǫ|g) = G1 (ǫ|g) G2 (ǫ|g)
(62)












G1 (ǫ|g) = 












α
..

.

0

0
..
.

0
−α

0
..

0

0
..
.

..

.

0
..
.
0
..

α
0
0

.

−β






 0


G2 (ǫ|g) = 

 0

 β






.

−α
0























(63)






−β 

0 


.

0 





β

(64)

The steady-state vector ΠT is the null vector of U(ǫ|g), and
hence the vector in (45) is a kind of the null vector.
A PPENDIX C
P ROOF OF PROPOSITION 4

<

− dB(ǫ+∆|g)
dǫ
B(ǫ|g)

(60a)

<

− dB(ǫ+∆|g)
dǫ
B(ǫ + ∆|g)

We define A(ǫ) = (1 − Pf (ǫ)) and B(ǫ|g) = Pa (ǫ|g) for
the proof of the MDP objective function OM (ǫ|g) case.

(60b)

<

− dB(ǫ+∆|g)
dǫ
.
B(ǫ + ∆|g)

(60c)

dA(ǫ)
dB(ǫ|g)
df M (ǫ|g)
=
B(ǫ|g) + A(ǫ)
.
(65)
dǫ
dǫ
dǫ
The increasing function A(ǫ) and the decreasing function
B(ǫ|g) satisfy the following properties (54) - (57). From

14

M

Proposition 5, we can derive that if 0 < ǫ < ǫe , df dǫ(ǫ|g) > 0.
Also, from Proposition 6, we can derive that if ǫe < ǫ
M
M
and df dǫ(ǫ|g) < 0, df (ǫ+∆|g)
< 0 where ∆ ∈ R+ .
dǫ
df M (ǫ|g)
= 0 has one solution or
Thus, the equation
dǫ
no solution, and the function f M (ǫ|g) converges to γ2 =
(limǫ→0 A(ǫ))(limǫ→0 B(ǫ|g)).
R EFERENCES
[1] T. Chen, Y. Yang, H. Zhang, H. Kim, and K. Horneman, “Network energy
saving technologies for green wireless access networks,” IEEE Wireless
Commun., vol. 18, no. 5, pp. 30-38, Oct. 2011.
[2] C. Han et al., “Green radio: Radio techniques to enable energy-efficient
wireless networks,” IEEE Commun. Mag., vol. 49, no. 6, pp. 46-54, Jun.
2011.
[3] X. Lu, P. Wang, D. Niyato, D. I. Kim, and Z. Han, “Wireless networks
with RF energy harvesting: A contemporary survey,” IEEE Commun.
Surveys & Tutorials, vol. 17, no. 2, pp. 757-789, Second Quarter 2015.
[4] X. Lu, P. Wang, D. Niyato, D. I. Kim, and Z. Han, “Wireless charging
technologies: Fundamentals, standards, and network applications,” IEEE
Commun. Surveys & Tutorials, vol. 18, no. 2, pp. 1413-1452, Second
Quarter 2016.
[5] A. M. Zungeru, L. Ang, S. Prabaharan, and K. P. Seng, “Radio frequency
energy harvesting and management for wireless sensor networks,” Green
Mobile Devices and Networks: Energy Optimization and Scavenging
Techniques, ch. 13, pp. 341-368, CRC Press 2012.
[6] S. Lee, R. Zhang, and K. Huang, “Opportunistic wireless energy harvesting in cognitive radio networks,” IEEE Trans. Wireless Commun., vol.
12, pp. 4788-4799, Sep. 2013.
[7] S. Park, H. Kim, and D. Hong, “Cognitive radio networks with energy
harvesting,” IEEE Trans. Wireless Commun., vol. 12, pp. 1386-1397, Mar.
2013.
[8] Q. Zhao, L. Tong, A. Swami, and Y. Chen, “Decentralized cognitive
MAC for opportunistic spectrum access in ad hoc networks: A POMDP
framework,” IEEE J. Select. Area. Commun., vol. 25, no. 3, pp. 589-600,
2007.
[9] D. Blei and M. Jordan, “Variational inference for Dirichlet process
mixtures,” Bayesian Analysis, vol. 1, no. 1, pp. 121-144, Aug. 2006.
[10] F. Wood and M. J. Black, “A nonparametric bayesian alternative to spike
sorting,” Journal of Neuroscience Methods, vol. 173, no. 1, pp. 1-12, Jun.
2008.
[11] T. L. Griffiths and Z. Ghahramani, “Infinite latent feature models and the
Indian buffet process,” Gatsby Computational Neuroscience Unit, Tech.
Rep. 2005-001, 2005.
[12] S. Park, J. Heo, B. Kim, W. Chung, H. Wang, and D. Hong, “Optimal
mode selection for cognitive radio sensor networks with RF energy
harvesting,” IEEE Proc. PIMRC 2012, pp. 2155-2159, 2012.
[13] S. Park, H. Kim, and D. Hong, “Cognitive radio networks with energy
harvesting,” IEEE Trans. Wireless Commun., vol. 12, no. 3, pp. 13861397,
Mar. 2013.
[14] C. M. Bishop, “Mixtures of gaussians,” in Pattern recognition and
machine learning 1st ed., Springer-Verlag New York, ch. 9, sec. 2, pp.
430-435.
[15] A. Sultan, “Sensing and transmit energy optimization for an energy
harvesting cognitive radio,” IEEE Wireless Commun. Letters, vol. 1, no.
5, pp. 500-503, Oct. 2012.
[16] D. T. Hoang, D. Niyato, P. Wang, D. I. Kim, “Opportunistic channel
access and RF energy harvesting in cognitive radio network,” IEEE
Journal on Selected Areas in Communications - Cognitive Radio Series,
vol. 32, pp. 2039-2052, Nov. 2014.
[17] D. T. Hoang, D. Niyato, P. Wang, D. I. Kim, “Performance optimization
for cooperative multiuser cognitive radio networks with RF energy
harvesting capability,” IEEE Trans. Wireless Commun., vol. 14, pp. 36143629, July 2015.
[18] M. E. Ahmed, J. B. Song, Z. Han, and D. Y. Suh, “Sensing-transmission
edifice using Bayesian nonparametric traffic clustering in cognitive radio
networks,” IEEE Trans. Mobile Computing, vol. 13, pp. 2141-2155, Sep.
2014.
[19] M. E. Ahmed, D. I. Kim, J. Y. Kim, and Y. A. Shin, “Energy-arrivalaware detection threshold in wireless-powered cognitive radio networks,”
IEEE Trans. Vehic. Technol., vol. 66, pp. 9201-9213, Oct. 2017.
[20] M. E. Ahmed, D. I. Kim, and K. W. Choi, “Traffic-aware optimal
spectral Access in wireless powered cognitive radio networks,” IEEE
Trans. Mobile Computing, vol. 17, pp. 734-745, Mar. 2018.

[21] J. Xu, and R. Zhang, “Throughput optimal policies for energy harvesting
wireless transmitters with non-ideal circuit power,” IEEE J. Select. Area.
Commun., vol. 32, pp. 322-332, Feb. 2014.
[22] M. E. Ahmed, J. B. Song, N. T. Nguyen, and Z. Han, “Nonparametric
Bayesian identification of primary users’ payloads in cognitive radio
networks,” IEEE Proc. ICC 2012, pp. 1586-1591, June 2012.
[23] M. E. Ahmed, J. B. Song, Z. Han, and D. Y. Suh, “Sensing-transmission
edifice using Bayesian nonparametric traffic clustering in cognitive radio
networks,” IEEE Trans. Mobile Computing, vol. 13, no. 9, pp. 2141-2155,
Sept. 2014.
[24] W. S. Chung, S. S. Park, S. M. Lim, and D. S. Hong, “Spectrum sensing
optimization for energy-harvesting cognitive radio systems,” IEEE Trans.
Wireless Commun., vol. 13, pp. 2601-2613, May 2014.
[25] Y. Pei, Y. C. Liang, K. C. Teh, and K. H. Li, “Energy-efficient design of
sequential channel sensing in cognitive radio networks: Optimal sensing
strategy, power allocation, and sensing order,” IEEE J. Select. Area.
Commun., vol. 29, pp. 1648-1659, Aug. 2011.
[26] [Online]. Available FTP: http://crawdad.cs.dartmouth.edu/meta.php?name=snu/wowviaw
[27] [Online]. Available FTP: http://crawdad.cs.dartmouth.edu/meta.php?name=kaist/wibro

