1. Introduction {#sec1-sensors-18-01781}
===============

Wireless sensor networks (WSNs) are subject to the constraint of energy storage at each mobile node \[[@B1-sensors-18-01781],[@B2-sensors-18-01781]\]. Saving on energy consumption or extending the battery life for sensor nodes has become an important research issue in wireless sensor networks. Recently, energy harvesting (EH) has attracted enormous attention from researchers as a promising cost-effective technique to maximize energy efficiency of a wireless network \[[@B3-sensors-18-01781],[@B4-sensors-18-01781],[@B5-sensors-18-01781]\], especially in wireless sensor networks \[[@B6-sensors-18-01781],[@B7-sensors-18-01781],[@B8-sensors-18-01781],[@B9-sensors-18-01781]\]. Various energy harvesting sources have been studied, such as natural sources (solar \[[@B10-sensors-18-01781]\], wind \[[@B11-sensors-18-01781]\], thermal \[[@B12-sensors-18-01781]\], etc.), strongly coupled magnetic resonances \[[@B13-sensors-18-01781]\], etc. Among them, radio frequency (RF) energy from the ambient transmitters is the most popular source for energy harvesting since it can be received more effectively from RF signals \[[@B4-sensors-18-01781]\]. Because RF signals can convey both energy and information simultaneously, an RF-based energy harvesting technique, called simultaneous wireless information and power transfer (SWIPT), is becoming a more and more promising research topic for WSNs \[[@B14-sensors-18-01781],[@B15-sensors-18-01781]\]. The idea of SWIPT was first introduced in 2008 in the seminal paper of Varshney \[[@B16-sensors-18-01781]\], in which authors proposed an ideal receiver design that is able to simultaneously observe the information and extract power from the same received signal. Zhou et al. \[[@B17-sensors-18-01781]\] connected Varshney's idea to practice by proposing two realizable receiver architectures design: time switching (TS) and power splitting (PS). The performance analysis of these two EH protocols was conducted rigorously by Nasir et al. \[[@B18-sensors-18-01781]\]. The authors derived analytical expressions for the outage probability and the ergodic capacity of one-way amplify-and-forward (AF) relay networks over a Rayleigh fading channel.

The applications of SWIPT to wireless sensor networks have been investigated in recent works. In \[[@B7-sensors-18-01781]\], the authors proposed a distributed iteration algorithm to solve an energy-efficient cooperative transmission problem for SWIPT in clustered WSNs. In \[[@B9-sensors-18-01781]\], the authors employed wireless energy harvesting techniques and realistic energy converters in dense and often randomly deployed WSNs, and quantified the potential energy gains that can be achieved in the network. Peng et al. \[[@B19-sensors-18-01781]\] considered a wireless energy harvesting two-way relay network (TWRN) using power splitting protocol, where the effects of practical hardware impairments were taken into consideration. Mouapi and Hakem \[[@B20-sensors-18-01781]\] proposed a new approach to defining the specifications of a stand-alone wireless node based on an RF EH system and implemented a hardware circuit to illustrate their energy optimization method. The performance of wireless powered sensor networks for Internet of Things (IoT) application was studied in \[[@B4-sensors-18-01781]\], where the authors proposed the optimal power allocation to maximize the system throughput and also derived the closed-form of that solution. However, in those papers, the channel gains are assumed either to be constant or to be a Rayleigh distributed random variable.

Wireless communication systems in the real radio environments are not only affected by the short-term fading (multipath), but are also subject to the shadowing effects. The vital issue in studying the performance of energy-harvesting-based wireless networks in these conditions is the outage probability analysis. An important statistical characteristic that can describe the behavior of the wireless channels is the probability density function (PDF) of the signal-to-noise-ratio (SNR) at the receiver output. To derive this PDF for different radio propagation environments is sometimes a difficult mathematical task, especially for complicated channel models such as Nakagami-m or Rician channels. In fact, not only the papers we mentioned about SWIPT for WSN, but most of the other results on outage performance up to now \[[@B21-sensors-18-01781],[@B22-sensors-18-01781]\] also focus on Rayleigh fading channels, where we can exploit the fact that the square of channel gain magnitude is exponentially distributed.

As a general comment, it is noted that very few publications about energy harvesting for Rician fading channel exist in the open technical literature. In fact, together with Rayleigh channel, Rician channel should be also considered as one of the small-scale fading models for WSNs, especially for the cases of relatively short range power transfer distance and with existence of a strong line-of-sight (LOS) path \[[@B23-sensors-18-01781]\]. Recently, Zhao et al. \[[@B24-sensors-18-01781]\] have derived the capacity expressions for wireless powered communication systems over a Rician fading channel. However, there is no relay in this study. The source directly harvests the energy from the power beacon. In \[[@B25-sensors-18-01781]\], the authors provided the throughput analysis of relay networks with two energy harvesting protocols (continuous and discrete) over a Rician fading channel, but were only limited in the case of perfect hardware. Furthermore, the paper only provided the integral form of the throughput for the continuous case, so it is not a computationally friendly result. Mishra et al. \[[@B26-sensors-18-01781]\] did provide impressive results on joint optimization of power allocation, power splitting for EH, and relay placement for SWIPT over a Rician channel. Nevertheless, they only use an approximation expression of outage probability to formulate the optimization problem, and did not consider the hardware impairment at the nodes. In practice, the transceiver hardware is imperfect due to phase noise, I/Q (In-phase/Quadrature) imbalance and amplifier nonlinearities \[[@B27-sensors-18-01781]\]. The modeling of hardware impairment in system performance analysis has been presented in many works, for instance, in \[[@B28-sensors-18-01781]\], where the authors analyzed the performance of dual-hop relaying systems in hardware impairment condition, in terms of the capacity, throughput and symbol error rate (SER). We also proposed and evaluated an energy harvesting-based spectrum access model in cognitive radio network with hardware impairment \[[@B29-sensors-18-01781]\]. Again, these works about hardware impairment considered the Rayleigh channels only.

In our current work, we tackle the problem of investigating the performance analysis of SWIPT for WSN over a Rician fading channel, which takes into account the hardware impairment at source and relay nodes. Specifically, we consider an AF two-way half-duplex energy harvesting relay network model suffering from hardware impairments at all nodes over the Rician fading channels. The exact analytical expressions of the achievable throughput, outage probability, and the exact-form expression for the PDF of SNR at each destination node of a half-duplex AF bidirectional wireless sensor networks over a Rician fading channel are derived rigorously. The main contributions of this paper can be described in more details as follows:The exact form expression of outage probability and achievable throughput at each destination node with imperfect hardware and in Rician fading environment are derived mathematically.We derive the exact-form cumulative distribution function (CDF) of the SNR at each destination node, and use this result to derive the integral exact-form of the SER at destination nodes.We also conduct the asymptotic analysis and provide the approximation of all performance factors mentioned above at high $P/N_{0}$ regime.The analytical results are all confirmed by Monte Carlo simulations. Using the simulation results, the effect of various system parameters on the system performance is carefully studied.

The rest of this paper is organized as follows. [Section 2](#sec2-sensors-18-01781){ref-type="sec"} describes the system model and the EH protocol that is used in this paper. [Section 3](#sec3-sensors-18-01781){ref-type="sec"} provides the detailed performance analysis of the system, including exact analysis and asymptotic analysis. The numerical results to validate the analysis are presented in [Section 4](#sec4-sensors-18-01781){ref-type="sec"}. Finally, conclusions are drawn in [Section 5](#sec5-sensors-18-01781){ref-type="sec"}.

2. System Model {#sec2-sensors-18-01781}
===============

2.1. System Model Description {#sec2dot1-sensors-18-01781}
-----------------------------

We consider a bidirectional network with two sources and one relay as illustrated in [Figure 1](#sensors-18-01781-f001){ref-type="fig"}. The transmission model follows the principles of analog network coding \[[@B30-sensors-18-01781],[@B31-sensors-18-01781]\]. This concept is the extension of linear network coding to multihop wireless networks. Traditionally, two sources transmit their data to the relay node in two different time slots to avoid interference. However, with analog coding, the data can be transmitted simultaneously to the relay. The relay node forwards this noisy sum of two signals during the next time slot. Because each node already knows one of the signals by virtue of having forwarded it earlier, it can cancel its own part in the received signals and decode the message sent by the other source.

In our model, every terminal has only one antenna and operates in a half-duplex mode. Let $S_{1}$ and $S_{2}$ denote the first node and the second node that are going to exchange their messages, respectively, and *R* denote the relay node. We denote the channel gain between node $S_{i}$ and the relay *R* as $g_{i}$, for $i = 1,2$. Here, both channels are assumed to be Rician fading channels and be reciprocal. In this research, we take into account the hardware impairment at both source nodes $S_{1}$, $S_{2}$, and relay node *R*. Regarding the communication between two nodes, we assume that the direct link between two nodes is very weak, and the communication relies solely with the help of relay. However, the relay has only the energy to serve their own purpose, so it needs to harvest the energy from the two nodes of interest before forwarding the information messages. As in many previous publications on two-way half-duplex channels, we assume that the two sources as well as the relay know the channel gains.

2.2. Energy Harvesting and Information Transfer Protocols {#sec2dot2-sensors-18-01781}
---------------------------------------------------------

The energy harvesting protocol for this system is illustrated in [Figure 2](#sensors-18-01781-f002){ref-type="fig"}. In this protocol, the transmission is divided into blocks of length *T*. Each transmission block consists of three time slots. In the first time slot, which has the duration of $\alpha T$, the relay harvests energy from the nodes $S_{1}$ and $S_{2}$. The two remaining time slots are used for information transferring. First, $S_{1}$ and $S_{2}$ simultaneously transmit to *R* at the second time slot with the transmitted power P. Then, at the third time slot, *R* amplifies the signal that it received during the second time slot and broadcasts it to $S_{1}$ and $S_{2}$. Both the second and third time slots have the same length of $\left( 1 - \alpha \right)T/2$.

Let *P* denote the average power transmitted by each node $S_{1}$ and $S_{2}$ during the first time slot. After the first time slot, the amount of the harvested energy at node *R* can be computed as $$E_{r} = \left( {P\left| g_{1} \right|^{2} + P\left| g_{2} \right|^{2}} \right)\eta\alpha T,$$ where $\eta$ is the energy harvesting efficiency.

Let $s_{1}$ and $s_{2}$ denote the information-bearing symbols transmitted from $S_{1}$ and $S_{2}$, respectively. Again, assume that both $s_{1}$ and $s_{2}$ have the same average power *P*, i.e. $E\left\lbrack \left| s_{i} \right|^{2} \right\rbrack = P$. Thus, the received signal at *R* at the second time slot can be written as \[[@B28-sensors-18-01781]\] $$y_{r} = \sum\limits_{i = 1}^{2}{g_{i}\left( s_{i} + \eta_{i} \right) + n_{r}},i = 1,2,$$ where:$\eta_{i}$ denotes the hardware distortion noise at $S_{i}$ with zero mean and variance $\kappa^{2}P$. Here, $\kappa$ is sufficient to characterize the aggregate level of impairments of the channel \[1\].$n_{r}$ is the additive white Gaussian noise (AWGN) at *R* with zero mean and variance $\sigma_{r}^{2}$.

At the second time slot, *R* amplifies the signal $y_{r}$ with an amplifying coefficient $\beta$ and then retransmits it to $S_{1}$ and $S_{2}$. The received signal $y_{i}$ at node $S_{i}$ at the third time slot is $$y_{i} = g_{i}\left( \beta y_{r} + \eta_{r} \right) + n_{i},i = 1,2,$$ where $\eta_{r}$ denotes the distortion noise with zero mean and variance $\kappa_{r}^{2}P_{r}$, $n_{i}$ is the additive white Gaussian noise (AWGN) at $S_{i}$ with zero mean and variance $\sigma_{i}^{2}$, for $i = 1,2$.

According to the law of energy conservation, the total energy that node *R* uses to transmit the amplified signal to $S_{1}$ and $S_{2}$ must be equal to the energy that *R* received from $S_{1}$ and $S_{2}$ at the first time slot, $E_{r}$. Now, because the transmission duration of the relay is $\frac{\left( 1 - \alpha \right)T}{2}$, the average power of the signal transmitted by *R* at the third time slot can be computed as $$P_{r} = \frac{E_{r}}{\left( 1 - \alpha \right)T/2} = \psi\left\lbrack {P\left| g_{1} \right|^{2} + P\left| g_{2} \right|^{2}} \right\rbrack,$$ where $\psi = \frac{2\eta\alpha}{1 - \alpha}$.

In order to ensure that the transmission power at *R* is $P_{r}$, the amplifying coefficient $\beta$ can be chosen as $$\beta = \frac{r}{y_{r}} = \sqrt{\frac{P_{r}}{X_{1}P + X_{1}P\kappa^{2} + X_{2}P + X_{2}P\kappa^{2} + \sigma_{r}^{2}}},$$ where $X_{i} = \left| g_{i} \right|^{2}$ denotes the square of channel gain amplitude between $S_{i}$ and *R*.

Note that the relay node *R* does not need to estimate the hardware noise parameters as well as the individual channel gains $g_{i}$ because the denominator of Equation ([5](#FD5-sensors-18-01781){ref-type="disp-formula"}) is the received symbol power during the 2nd time slot. Now, because of the similarity between the roles of $S_{1}$ and $S_{2}$, we can focus only on node $S_{1}$. In fact, the received signal $y_{1}$ can be rewritten as

y

1

=

β

g

1

2

s

1

\+

β

g

1

2

η

1

\+

β

g

1

g

2

s

2

\+

β

g

1

g

2

η

2

\+

β

g

1

n

r

\+

g

1

η

r

\+

n

1

.

This signal contains both messages $s_{1}$ and $s_{2}$, while only $s_{2}$ is the desired signal at $s_{1}$. Since node $S_{1}$ perfectly knows its own transmitted symbol $s_{1}$, it can eliminate the corresponding the self-interference term $\beta\left| g_{1} \right|^{2}s_{1}$ from $y_{1}$.

From this fact, Equation ([6](#FD6-sensors-18-01781){ref-type="disp-formula"}) can be rewritten as

y

1

=

β

g

1

g

2

s

2

︸

s

i

g

n

a

l

\+

β

g

1

2

η

1

\+

β

g

1

g

2

η

2

\+

β

g

1

n

r

\+

g

1

η

r

\+

n

1

︸

n

o

i

s

e

.

Therefore, the end-to-end signal-to-noise-ratio at $S_{1}$ for detection of the symbol $s_{1}$ is given by

γ

1

=

E

s

i

g

n

a

l

2

E

n

o

i

s

e

2

=

β

2

g

1

2

g

2

2

P

β

2

g

1

4

P

κ

2

\+

β

2

g

1

2

g

2

2

P

κ

2

\+

β

2

g

1

2

σ

r

2

\+

g

1

2

P

r

κ

r

2

\+

σ

1

2

.

After doing some algebra and using the fact that $\frac{\sigma_{1}^{2}\sigma_{r}^{2}}{P_{r}} \approx 0$, combining with Equation ([6](#FD6-sensors-18-01781){ref-type="disp-formula"}), we obtain:$$\gamma_{1} = \frac{X_{1}X_{2}}{X_{1}X_{2}a + X_{1}^{2}a + X_{1}b + c},$$ where $a = \kappa^{2} + \kappa_{r}^{2}\left( 1 + \kappa^{2} \right),b = \frac{\left( 1 + \kappa_{r}^{2} \right)}{P/\sigma_{r}^{2}},c = \frac{\left( 1 + \kappa^{2} \right)}{\psi P/\sigma_{1}^{2}}$ (assume that $\sigma_{1}^{2} = \sigma_{2}^{2} = \sigma_{r}^{2} = N_{0}$).

3. System Performance {#sec3-sensors-18-01781}
=====================

3.1. Outage Probability {#sec3dot1-sensors-18-01781}
-----------------------

Since the random variable (RV) $X_{i}$ is Rician distributed for $i = 1,2$, the probability density function (PDF) of $X_{i}$ can be found as $$f_{X_{i}}\left( x \right) = \frac{\left( K + 1 \right)e^{- K}}{\lambda_{i}}e^{- \frac{(K + 1)x}{\lambda_{i}}}I_{0}\left( {2\sqrt{\frac{K\left( K + 1 \right)x}{\lambda_{i}}}} \right),$$ where $\lambda_{i}$ is the mean value of the RV $X_{i},i = 1,2$, *K* is the Rician K-factor defined as the ratio of the power of the line-of-sight (LOS) component to the scattered components and $I_{0}\left( \cdot \right)$ is the zero-th order modified Bessel function of the first kind.

Using the equality $I_{0}\left( x \right) = \sum\limits_{l = 0}^{\infty}\frac{x^{2l}}{2^{2l}\left( l! \right)^{2}}$ \[[@B32-sensors-18-01781]\], Equation ([10](#FD10-sensors-18-01781){ref-type="disp-formula"}) can be rewritten as $$f_{X_{i}}\left( x \right) = \zeta_{i}e^{- K}\sum\limits_{l = 0}^{\infty}{\frac{\left( \zeta_{i}K \right)^{l}}{\left( l! \right)^{2}}x^{l}e^{- \zeta_{i}x}},$$ where $\zeta_{i} = \frac{K + 1}{\lambda_{i}}$.

The cumulative density function (CDF) of the RV $X_{i}$ ($i = 1,2$) can be derived as in \[[@B33-sensors-18-01781]\]:$$F_{X_{i}}\left( x \right) = \int\limits_{0}^{x}{f_{X_{i}}\left( t \right)dt = 1 - e^{- K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{m = 0}^{l}\frac{K^{l}\zeta_{i}^{m}}{l!m!}}}x^{m}e^{- \zeta_{i}x}.$$

The outage probability of this system is defined as the probability that the end-to-end SNR falls below a desired threshold---let us say, $\gamma_{th}$. In this case, it can be derived as $$P_{out\_ 1} = F_{\gamma_{1}}\left( \gamma_{th} \right) = \Pr\left( {\gamma_{1} < \gamma_{th}} \right) = \Pr\left( {\frac{X_{1}X_{2}}{X_{1}X_{2}a + X_{1}^{2}a + X_{1}b + c} < \gamma_{th}} \right),$$ where $\gamma_{th} = 2^{2R} - 1$, *R* is the source transmission rate:$$\begin{matrix}
{\Pr\left( {\frac{X_{1}X_{2}}{X_{1}X_{2}a + X_{1}^{2}a + X_{1}b + c} < \gamma_{th}} \right) = \Pr\left\{ {X_{2}\left( X_{1} - \gamma_{th}X_{1}a \right) < \gamma_{th}X_{1}^{2}a + \gamma_{th}X_{1}b + \gamma_{th}c} \right\},} \\
{= \left\{ \begin{matrix}
{\Pr\left\{ {X_{2} < \frac{\gamma_{th}X_{1}^{2}a + \gamma_{th}X_{1}b + \gamma_{th}c}{X_{1} - \gamma_{th}X_{1}a}} \right\},{if}\; 1 - \gamma_{th}a \geq 0,} \\
{1,{if}\; 1 - \gamma_{th}a < 0.} \\
\end{matrix} \right.} \\
\end{matrix}$$

In this subsection, we assume that $1 - \gamma_{th}a$ is positive, which is reasonable in practice. From Equation ([14](#FD14-sensors-18-01781){ref-type="disp-formula"}), we can derive the exact-form expression for the outage probability of the system. This is stated by the following theorem.

*The exact outage probability for node $S_{1}$ of the proposed half-duplex bidirectional wireless sensor network using time-switching energy harvesting strategy over a Rician fading channel can be expressed as* $$\begin{matrix}
{{P_{out\_ 1} = 1 - 2\zeta_{1}e^{- 2K}e^{- b_{1}\zeta_{2}{\widetilde{\gamma}}_{th}}}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{n = 0}^{m}{\sum\limits_{p = 0}^{n}\frac{{K^{l + k}\zeta_{1}^{k}\left( \zeta_{2}{\widetilde{\gamma}}_{th} \right)^{\frac{3m - 2n + p + k + 1}{2}}}a^{m - n}b^{p}c^{\frac{m - p + k + 1}{2}}}{l!p!\left( n - p \right)!\left( m - n \right)!\left( k! \right)^{2}\left( \zeta_{1} + a\zeta_{2}{\widetilde{\gamma}}_{th} \right)^{\frac{m - 2n + p + k + 1}{2}}}}}}}}} \\
{\times K_{m - 2n + p + k + 1}\left( {2\sqrt{\left( \zeta_{1} + a\zeta_{2}{\widetilde{\gamma}}_{th} \right)\zeta_{2}{\widetilde{\gamma}}_{th}c}} \right),} \\
\end{matrix}$$ *where $a = \kappa^{2} + \kappa_{r}^{2}\left( 1 + \kappa^{2} \right),b = \frac{\left( 1 + \kappa_{r}^{2} \right)}{P/\sigma_{r}^{2}},c = \frac{\left( 1 + \kappa^{2} \right)}{\psi P/\sigma_{i}^{2}},\zeta_{i} = \frac{K + 1}{\lambda_{i}}$, and ${\widetilde{\gamma}}_{th} = \frac{\gamma_{th}}{1 - a\gamma_{th}}$.*

See Appendix A. ☐

*The outage probability at node $S_{2}$ can also be obtained by exchanging the roles of $\zeta_{1}$ and $\zeta_{2}$ in Equation* ([15](#FD15-sensors-18-01781){ref-type="disp-formula"})*.*

3.2. Achievable Throughput {#sec3dot2-sensors-18-01781}
--------------------------

With the outage probability obtained from Theorem 1, the achievable throughput of the considered system can be given by $$\tau_{i} = \left( 1 - P_{out\_ i} \right).\frac{R}{2}.\left( 1 - \alpha \right),i \in \left\{ 1,2 \right\}.$$

3.3. SER Analysis {#sec3dot3-sensors-18-01781}
-----------------

In this section, we derive the formula for symbol error rate (SER) at the destination. In fact, the theoretical SER at the receiving node can be computed in terms of the end-to-end SNR at that node \[[@B34-sensors-18-01781]\]. By taking into account the randomness of the SNR in our model, the SER at node $S_{i}$ can be derived by taking the expectation with respect to the random variable SNR: $${SER}_{i} = E\left\lbrack {\omega Q\left( \sqrt{2\theta\gamma_{i}} \right)} \right\rbrack,i \in {1,2},$$ where $Q\left( t \right) = \frac{1}{\sqrt{2\pi}}\int\limits_{t}^{\infty}{e^{- x^{2}/2}dx}$ is the Gaussian Q-function; ω and θ are constants that are specific for each modulation type---in particular, $\omega = 1,\theta = 1$ for Binary Phase Shift keying (BPSK) and $\omega = 1,\theta = 2$ for Quadrature Phase Shift Keying (QPSK). Hence, before obtaining the SER formula, the CDF of $\gamma_{i}$ is needed. Let $F_{\gamma_{i}}\left( x \right)$ be the CDF of $\gamma_{i}$. By exchanging the order of integration and changing variable, we come up with the following formula: $${SER}_{i} = \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{0}^{\infty}{\frac{e^{- \theta x}}{\sqrt{x}}F_{\gamma_{i}}\left( x \right)dx}.$$

Now, we can claim the following theorem on SER.

*The exact symbol error rate for node $S_{1}$ of the proposed half-duplex bidirectional wireless sensor network using time-switching energy harvesting strategy over a Rician fading channel can be expressed as:* $$\begin{array}{l}
{{SER_{1}} = \frac{\omega}{2} - \zeta_{1}e^{- 2K}\omega\sqrt{\theta}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{n = 0}^{m}{\sum\limits_{p = 0}^{n}\frac{{K^{l + k}\zeta_{1}^{k}}a^{m - n}b^{p}c^{\frac{m - p + k + 1}{2}}}{l!p!\left( n - p \right)!\left( m - n \right)!\left( k! \right)^{2}}}}}}} \times} \\
{\times \int\limits_{0}^{\infty}{\frac{e^{- \theta x - b\zeta_{2}{\widetilde{\gamma}}_{th}{(x)}}}{\sqrt{\pi x}}K_{m - 2n + p + k + 1}\left( {2\sqrt{\left\lbrack \zeta_{1} + a\zeta_{2}{\widetilde{\gamma}}_{th}\left( x \right) \right\rbrack\zeta_{2}{\widetilde{\gamma}}_{th}\left( x \right)c}} \right)\frac{\left\lbrack \zeta_{2}{\widetilde{\gamma}}_{th}\left( x \right) \right\rbrack^{\frac{3m - 2n + p + k + 1}{2}}}{\left\lbrack \zeta_{1} + a\zeta_{2}{\widetilde{\gamma}}_{th}\left( x \right) \right\rbrack^{\frac{m - 2n + p + k + 1}{2}}}dx},} \\
\end{array}$$ *where ${\widetilde{\gamma}}_{th}\left( x \right) = \frac{x}{1 - ax}$.*

See Appendix B. ☐

3.4. Asymptotic Analysis {#sec3dot4-sensors-18-01781}
------------------------

The asymptotic analysis is significant to provide further insights into the impact of hardware impairments on the network performance. It can also be used to verify the correctness of the exact analysis. In this section, the asymptotic outage probability and SER at the high SNR regime are going to be derived rigorously.

### 3.4.1. Outage Probability {#sec3dot4dot1-sensors-18-01781}

As the $P/N_{0}$ approaches infinity, it is obvious to see that the SNR in Equation ([9](#FD9-sensors-18-01781){ref-type="disp-formula"}) is asymptotically equal to $$\gamma_{1}^{\infty} = \frac{X_{2}}{aX_{2} + aX_{1}} = \frac{X_{2}}{a\left( X_{1} + X_{2} \right)},$$and, similarly, $$\gamma_{2}^{\infty} = \frac{X_{1}}{a\left( X_{1} + X_{2} \right)},$$ where $a = \kappa^{2} + \kappa_{r}^{2}\left( 1 + \kappa^{2} \right)$.

By using the similar approach as in the exact analysis of outage probability, we can get the results as stated in the following theorem.

*The asymptotic outage probability and achievable throughput for node $S_{1}$ of the proposed half-duplex bidirectional wireless sensor network using time-switching energy harvesting strategy over a Rician fading channel can be provided as* $$P_{out\_ 1}^{\infty} = 1 - \zeta_{1}e^{- 2K}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\frac{K^{l + k}{\zeta_{1}^{k}\zeta_{2}^{m}}}{l!m!\left( k! \right)^{2}}\left( {a{\widetilde{\gamma}}_{th}} \right)^{m}}\frac{\Gamma\left( m + k + 1 \right)}{\left( \zeta_{2}a{\widetilde{\gamma}}_{th} + \zeta_{1} \right)^{m + k + 1}}}}}$$ *and* $$\tau_{1} \triangleq \left( 1 - P_{out\_ 1}^{\infty} \right)\frac{R}{2} = \frac{\zeta_{1}Re^{- 2K}}{2}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\frac{K^{l + k}{\zeta_{1}^{k}\zeta_{2}^{m}}}{l!m!\left( k! \right)^{2}}\left( {a{\widetilde{\gamma}}_{th}} \right)^{m}}\frac{\Gamma\left( m + k + 1 \right)}{\left( \zeta_{2}a{\widetilde{\gamma}}_{th} + \zeta_{1} \right)^{m + k + 1}}}}},$$ *where $a = \kappa^{2} + \kappa_{r}^{2}\left( 1 + \kappa^{2} \right),b = \frac{\left( 1 + \kappa_{r}^{2} \right)}{P/\sigma_{r}^{2}},c = \frac{\left( 1 + \kappa^{2} \right)}{\psi P/\sigma_{i}^{2}},\zeta_{i} = \frac{K + 1}{\lambda_{i}},{\widetilde{\gamma}}_{th} = \frac{\gamma_{th}}{1 - a\gamma_{th}}$, and $\Gamma\left( \cdot \right)$ is the complete Gamma function, which is defined by $\Gamma\left( z \right) \triangleq \int_{0}^{\infty}{x^{z - 1}e^{- x}dx}$.*

See Appendix C. ☐

### 3.4.2. SER Analysis {#sec3dot4dot2-sensors-18-01781}

As $P/N_{0}$ approaches infinity, the symbol-error-rate at node $S_{1}$ becomes $${SER}_{1}^{\infty} = \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{0}^{\infty}{\frac{e^{- \theta x}}{\sqrt{x}}F_{\gamma_{1}}^{\infty}\left( x \right)dx},$$where $F_{\gamma_{1}}^{\infty}\left( x \right)$ is the CDF of $\gamma_{1}$ as $P/N_{0}$ goes to infinity.

The result of our asymptotic SER analysis is presented in Theorem 4.

*The asymptotic symbol error rate for node $S_{1}$ of the proposed half-duplex bidirectional wireless sensor network using a time-switching energy harvesting strategy over a Rician fading channel can be provided as the following:* *If $\xi_{1} = \xi_{2}$:* $$\begin{array}{l}
{{SER}_{1}^{\infty} = \frac{\omega}{2}} \\
{- \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{1}^{k}\zeta_{2}^{m}a^{m + v}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}\zeta_{1}^{m + 1}\theta^{m + v + \frac{1}{2}}}}}}}\gamma\left( {m + v + \frac{1}{2},\frac{\theta}{a}} \right).} \\
\end{array}$$*If $\xi_{1} \neq \xi_{2}$:* $$\begin{array}{l}
{SER_{1}^{\infty} = \frac{\omega}{2} - \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\sum\limits_{p = 0}^{\infty}{\left( m + k + 1 \right)_{p}\left\lbrack \left( \xi_{1} - \xi_{2} \right) \right\rbrack^{p}}\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}}}}}} \\
{\times \frac{K^{l + k}\zeta_{2}^{m}a^{m + v + p}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}p!\zeta_{1}^{m + p + 1}\theta^{m + v + p + \frac{1}{2}}}\gamma\left( {m + v + p + \frac{1}{2},\frac{\theta}{a}} \right),} \\
\end{array}$$ *where $a = \kappa^{2} + \kappa_{r}^{2}\left( 1 + \kappa^{2} \right),b = \frac{\left( 1 + \kappa_{r}^{2} \right)}{P/\sigma_{r}^{2}},c = \frac{\left( 1 + \kappa^{2} \right)}{\psi P/\sigma_{i}^{2}},\zeta_{1} = \frac{K + 1}{\lambda_{1}},{\widetilde{\gamma}}_{th} = \frac{\gamma_{th}}{1 - a\gamma_{th}}$, $\Gamma\left( \cdot \right)$ is the complete Gamma function, $\gamma\left( \cdot , \cdot \right)$ is the incomplete Gamma function, and $\left( \cdot \right)_{p}$ is the Pochhammer symbol.*

See Appendix D. ☐

3.5. Optimal Time-Switching Factor {#sec3dot5-sensors-18-01781}
----------------------------------

It is not difficult to learn that there is always a trade-off between the amount of energy used for transmission and the duration of the transmission in the considered protocol. Specifically, if more time is allocated for energy harvesting, a higher available transmission power can be obtained, which may lead to a higher throughput. However, at the same time, less time resources are left for signal transmission, which may lead to the decrease of the transmission rate. Hence, there exists an optimal time-switching factor $\alpha^{*}$ that provides the best throughput performance.

Given the throughput expression obtained in Equation ([16](#FD16-sensors-18-01781){ref-type="disp-formula"}), the optimal time switching can be obtained by solving Equation $\frac{d\tau_{i}\left( \alpha \right)}{d\alpha} = 0$. However, due to the complicated infinite series and Bessel function in each throughput expression, this optimization problem can hardly admit a closed-form solution.

Here, we apply an iterative algorithm to solve this problem numerically. In particular, the Golden section search algorithm, which has been used in many global optimization problems in communications (for example, in \[[@B35-sensors-18-01781]\]), is chosen for this work. For a detailed algorithm as well as its related theory, please refer to \[[@B36-sensors-18-01781]\].

4. Numerical Results and Discussion {#sec4-sensors-18-01781}
===================================

For the purpose of validation, the correctness of the derived outage probability and SER expressions as well as investigation of the effect of various parameters on the system performance, a set of Monte Carlo simulations are conducted and presented in this section. For each simulation, we first provide the graphs of the outage probability and throughput obtained by the analytical formulas. Secondly, we plot the same outage probability and throughput curves that result from Monte Carlo simulation. To do this, we generate $10^{5}$ random samples of each channel gain, which are Rician distributed. Using these random samples, the SNR at destination node $S_{1}$ is calculated and compared with the threshold value *γ*. The outage probability occurs if this SNR falls below the threshold. By taking the number of cases that $SNR < \gamma$ divided by the number of samples, we can estimate the outage probability and then the throughput of system. The analytical curve and the simulation one should match together to verify the correctness of our analysis.

The hardware impairment parameters are chosen as $\kappa = \kappa_{r} = 0.1$. The ideal hardware impairment situation ($\kappa = \kappa_{r} = 0$) is also considered as a benchmark performance for simulation. The channel gains are considered as Rician fading with $\lambda_{1} = \lambda_{2} = 0.5$ and with the Rician K-factors equal to 3 for both channels. The transmit power are set to the same value $P_{1} = P_{2} = P$ for both two sources, so that the ratio $P/N_{0}$ varies in the range from 0 to 50 dB. The energy harvesting efficiency is set to be 0.7. The source transmission rate is chosen as 1.5 bps/Hz. From the Shannon's theorem on capacity of the channel, we can calculate the $SNR$ threshold as $\gamma = 2^{2R} - 1$. All simulation parameters are listed in [Table 1](#sensors-18-01781-t001){ref-type="table"}.

4.1. Effects of Various Parameters on the System Performance {#sec4dot1-sensors-18-01781}
------------------------------------------------------------

[Figure 3](#sensors-18-01781-f003){ref-type="fig"} and [Figure 4](#sensors-18-01781-f004){ref-type="fig"} show the effect of $P/N_{0}$ on the outage probability and throughput of the proposed system, respectively. For this simulation, the utilized parameter settings are: $\kappa = \kappa_{r} = 0.1$ or $0.2,\alpha = 0.5$ and $\eta = 0.7$. We choose $\alpha = 0.5$ to consider the case that the duration of energy harvesting and the duration of transmission are balanced. The case $\kappa = 0$ (no hardware impairment) is also introduced for comparison. The first observation is that the outage probability and throughput obtained from mathematical analysis match with the corresponding Monte Carlo simulations. Regarding the effect of *κ*, the outage probability decreases and the throughput increases as *κ* varies from 0 to 0.2. When $P/N_{0}$ increases, the outage probability and throughput approach the corresponding asymptotic values obtained from analysis. Furthermore, the lower the value of *κ*, the faster the outage probability and throughput converge to their asymptotic values.

The effect of hardware impairment level on the outage probability and the achievable throughput at each node is presented more thoroughly in [Figure 5](#sensors-18-01781-f005){ref-type="fig"} and [Figure 6](#sensors-18-01781-f006){ref-type="fig"}. Here, $P/N_{0}$ is set at 20 dB and the transmission rate is fixed at 1.5 bps. Three values of *α* are chosen: 0.2, 0.5, and 0.8, corresponding to three cases: the energy harvesting duration is dominant, there is a balance between energy harvesting and information transmission, and the information transmission duration is dominant.

Again, it is observed that the exact-form expressions of outage probability and throughput obtained by the analysis coincide with the ones that are obtained by Monte Carlo simulations. From the numerical results, it is evident that the achievable throughput decreases and the outage probability increases significantly at each destination node when the impairment level *κ* increases. In addition, the outage probability tends to reduce at higher time-switching factor. This can be explained because the larger value of *α* means more power is used for data transmission. However, this doesn't mean that the throughput is better for larger *α*. In [Figure 6](#sensors-18-01781-f006){ref-type="fig"}, the throughput performance is improved when *α* increases from 0.2 to 0.5, but then degraded when *α* increases from 0.5 to 0.8.

[Figure 7](#sensors-18-01781-f007){ref-type="fig"} and [Figure 8](#sensors-18-01781-f008){ref-type="fig"} illustrate more clearly the effect of time-switching factor on the outage and throughput performance. In this simulation, the parameters are chosen as $\eta = 0.7$ and $P/N_{0} = 20$ dB (this value is chosen because it is in the middle range of $P/N_{0}$). The transmission rate varies among three values: 0.5 bps/Hz, 1 bps/Hz and 1.5 bps/Hz, while the time-switching factor varies in the range $\left( 0,1 \right)$. The results confirm what we mentioned just above. There should be a unique value of *α* that maximizes the throughput. This is because, when we increase *α* initially, there is more power used for transmission, so the outage probability is reduced and the throughput increases correspondingly. However, when *α* keeps increasing, the duration of transmission is also reduced, hence, less data is transmitted during a given time interval. As a result, the throughput performance becomes worse.

4.2. Effect of Various Parameters on SER {#sec4dot2-sensors-18-01781}
----------------------------------------

The purpose of the following simulations is to confirm the correctness of the SER formulas provided in the analysis. First, [Figure 9](#sensors-18-01781-f009){ref-type="fig"} presents the effect of the hardware impairment level on the SER performance. In this simulation, the time-switching factor is chosen as $\alpha = 0.5$ with the same reason as in [Section 4.1](#sec4dot1-sensors-18-01781){ref-type="sec"}, the transmission rate is fixed at 1 bps/Hz, and the ratio $P/N_{0}$ varies in the range from 0 dB to 40 dB. From the results, it is showed that SER decreases to the asymptotic value when the ratio $P_{s}/N_{0}$ increases. The Monte Carlo simulation curves overlap with the corresponding analysis curves. This confirms the validity of our analysis. When the hardware impairment level goes higher, the SER also has a larger value, as expected. Furthermore, the SER performance of the QPSK scheme is better than the one of the BPSK scheme in the same simulation condition. This can be explained because the QPSK modulation scheme can transmit two bits in one symbol while the BPSK scheme can only transmit one bit per symbol. Hence, if we have the same constraint on the transmission rate for both methods, the required SNR for maintaining good communication would be smaller for the QPSK method. As a result, the outage probability and the SER for QPSK modulation would be smaller than the ones of BPSK.

In a similar way, the influence of the time-switching factor on the symbol-error-rate at the destination node is illustrated in [Figure 10](#sensors-18-01781-f010){ref-type="fig"}. The simulation parameters are $\kappa = 0.1,R = 1$ bps/Hz, and $P/N_{0}$ varies from 0 dB to 40 dB. Again, the simulation curves match perfectly with the corresponding analysis curves. The SER tends to approach its asymptotic value, and QPSK modulation still provides the better SER performance than BPSK for both values of *α*. Note that, for this simulation, the asymptotic values do not depend on the time-switching factor *α*. However, *α* surely has an effect on the immediate value of SER. In fact, the SER performance should be better with the value of *α* that is in the middle of its range. For example, in [Figure 10](#sensors-18-01781-f010){ref-type="fig"}, the SER value for $\alpha = 0.5$ is less than the one with $\alpha = 0.2$. The explanation is the same as the case of [Figure 8](#sensors-18-01781-f008){ref-type="fig"}.

4.3. Optimal Time-Switching Factor {#sec4dot3-sensors-18-01781}
----------------------------------

As mentioned in [Section 3](#sec3-sensors-18-01781){ref-type="sec"}, the optimal time-switching factor to maximize the achievable throughput of the considered system can be found numerically by using an iterative algorithm such as the golden section search method \[[@B36-sensors-18-01781]\]. [Figure 11](#sensors-18-01781-f011){ref-type="fig"} plots the optimal value $\alpha^{*}$ for various values of the ratio $P/N_{0}$ at different hardware impairment levels.

It can be observed that the optimal time-switching factor decreases as the ratio $P/N_{0}$ increases. This is because, for large $P/N_{0}$, the outage probability tends to reduce, so it is not necessary to use a large amount of energy to transmit data. Reversely, we need to spend more time resources to increase the throughput of the system.

On the other hand, we can learn from this simulation that the optimal *α* does not change much for different hardware impairment levels. Especially, for small *κ*, the value of $\alpha^{*}$ is nearly the same as the one for a perfect hardware case.

5. Conclusions {#sec5-sensors-18-01781}
==============

Recent development in wireless sensor networks have led to an exponential growth of the energy demand for operating the networks, which raises a question about how to efficiently use the available energy in the wireless environment. This paper rigorously analyzes the performance of a half-duplex AF bidirectional sensor network in which the relay node is equipped with time-switching-based energy harvesting protocol. The channel considered in this paper is a Rician fading channel. We also take into account the hardware impairment at source and relay nodes. We derive both exact and asymptotic forms of the outage probability, achievable throughput, as well as the symbol-error-rate at each destination node. The analysis results are validated by Monte Carlo simulation. From the results of this work, we can gain an insightful understanding of the effect of various parameters on the system performance. Furthermore, the optimal time-switching factor (i.e., the best energy harvesting strategy) is also founded by numerical algorithms.

Tan N. Nguyen created the main idea and developed the performance evaluation by analysis and simulation for the basic case that $\lambda_{1} = \lambda_{2}$. Tran Hoang Quang Minh set up the simulation configurations, performed the numerical experiments and wrote the paper. Phuong T. Tran extended the analysis for the general case. Phuong T. Tran and Miroslav Voznak were working as the advisors of Tan N. Nguyen to discuss and advise the main idea and performance analysis together.

This research received funding from the grant No. SP2018/59 conducted by VSB--Technical University of Ostrava, Czech Republic.

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

The following abbreviations are used in this manuscript:

AF

amplify-and-forward

AWGN

additive white Gaussian noise

BPSK

binary phase shift keying

CDF

cumulative distribution function

EH

energy harvesting

I/Q

In-phase/Quadrature

LOS

line-of-sight

PDF

probability density function

PS

power splitting

QPSK

quadrature phase shift keying

RF

radio frequency

RV

random variable

SER

symbol-error-rate

SNR

signal-to-noise ratio

SWIPT

simultaneous wireless information and power transfer

TS

time switching

TWRN

two-way relay network

WSN

wireless sensor network

In this section, the detailed proof of Theorem 1 is provided. Without loss of generality, we can consider the node $S_{1}$ only. From Equation ([14](#FD14-sensors-18-01781){ref-type="disp-formula"}), we can fix $X_{1}$ at some value $x_{1}$ and compute the probability in Equation ([14](#FD14-sensors-18-01781){ref-type="disp-formula"}) in terms of $x_{1}$. Then, we take the average of the obtained function with respect to the distribution of $X_{1}$ to get the outage probability. Using this approach, Equation ([14](#FD14-sensors-18-01781){ref-type="disp-formula"}) can be rewritten as $$P_{out\_ 1} = \int\limits_{0}^{\infty}{F_{X_{2}}\left( \frac{\gamma_{th}x_{1}^{2}a + \gamma_{th}x_{1}b + \gamma_{th}c}{x_{1} - \gamma_{th}x_{1}a} \right){f_{X_{1}}\left( x_{1} \right)}dx_{1}}.$$

Let us denote ${\widetilde{\gamma}}_{th} \triangleq \frac{\gamma_{th}}{1 - a\gamma_{th}}$. Then, we have $\frac{\gamma_{th}x_{1}^{2}a + \gamma_{th}x_{1}b + \gamma_{th}c}{x_{1} - \gamma_{th}x_{1}a} = {\widetilde{\gamma}}_{th}\left( {x_{1}a + \frac{c}{x_{1}} + b} \right)$. Now, Equation ([A1](#FD27-sensors-18-01781){ref-type="disp-formula"}) can be rewritten as $$\begin{array}{l}
{P_{out\_ 1} =} \\
{= 1 - \int\limits_{0}^{\infty}{{\sum\limits_{l = 0}^{\infty}{\sum\limits_{m = 0}^{l}\frac{e^{- K}K^{l}\left( \zeta_{2}{\widetilde{\gamma}}_{th} \right)^{m}}{l!m!}\left( {x_{1}a + \frac{c}{x_{1}} + b_{1}} \right)^{m}}}e^{- \zeta_{2}{\widetilde{\gamma}}_{th}{({ax_{1} + \frac{c}{x_{1}} + b})}}}\zeta_{1}e^{- K}\sum\limits_{k = 0}^{\infty}\frac{\left( \zeta_{1}Kx_{1} \right)^{k}e^{- \zeta_{1}x_{1}}}{\left( k! \right)^{2}}dx_{1}} \\
{= 1 - \zeta_{1}e^{- 2K}\int\limits_{0}^{\infty}{{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\frac{K^{l + k}\zeta_{2}^{m}{\widetilde{\gamma}}_{th}^{m}}{l!m!\left( k! \right)^{2}}\left( {x_{1}a + \frac{c}{x_{1}} + b} \right)}}}}^{m}e^{- {\zeta_{2}{\widetilde{\gamma}}_{th}}{({x_{1}a + \frac{c}{x_{1}} + b})}}\left( \zeta_{1}x_{1} \right)^{k}e^{- {\zeta_{1}x_{1}}}dx_{1}}.} \\
\end{array}$$

By applying the formula $\left( x + y \right)^{m} = \sum\limits_{n = 0}^{m}\begin{pmatrix}
m \\
n \\
\end{pmatrix}x^{m - n}y^{n}$ and letting $\xi_{12} \triangleq \zeta_{1}e^{- 2K}e^{- b\zeta_{2}{\widetilde{\gamma}}_{th}}$ , we obtain: $$\begin{matrix}
P_{out\_ 1} & {= 1 - \xi_{12}\int\limits_{0}^{\infty}{{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\frac{K^{l + k}\zeta_{2}^{m}{\widetilde{\gamma}}_{th}^{m}}{l!m!\left( k! \right)^{2}}\sum\limits_{n = 0}^{m}{\left( \begin{array}{l}
m \\
n \\
\end{array} \right)\left( {x_{1}a} \right)^{m - n}\left( {\frac{c}{x_{1}} + b} \right)^{n}}}}e^{- {\zeta_{2}{\widetilde{\gamma}}_{th}}{({x_{1}a + \frac{c}{x_{1}}})}}\left( \zeta_{1}x_{1} \right)^{k}}}e^{- {\zeta_{1}x_{1}}}dx_{1}}} \\
 & {= 1 - \xi_{12}\int\limits_{0}^{\infty}{{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{n = 0}^{m}{\frac{{K^{l + k}{\zeta_{1}^{k}\zeta_{2}^{m}}{\widetilde{\gamma}}_{th}^{m}}a^{m - n}}{l!n!\left( m - n \right)!\left( k! \right)^{2}}{x_{1}}^{m - n + k}\sum\limits_{p = 0}^{n}{\left( \begin{array}{l}
n \\
p \\
\end{array} \right)\left( \frac{c}{x_{1}} \right)^{n - p}b^{p}}}}}e^{- {\zeta_{2}{\widetilde{\gamma}}_{th}}{({x_{1}a + \frac{c}{x_{1}}})}}}}e^{- {\zeta_{1}x_{1}}}dx_{1}}} \\
 & {= 1 - \xi_{12}\int\limits_{0}^{\infty}{{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{n = 0}^{m}{\sum\limits_{p = 0}^{n}{\frac{{K^{l + k}{\zeta_{1}^{k}\zeta_{2}^{m}}{\widetilde{\gamma}}_{th}^{m}}a^{m - n}b^{p}c^{n - p}}{l!p!\left( n - p \right)!\left( m - n \right)!\left( k! \right)^{2}}{x_{1}}^{m - 2n + k + p}}}}}e^{- {\zeta_{2}{\widetilde{\gamma}}_{th}}{({x_{1}a + \frac{c}{x_{1}}})}}}}e^{- {\zeta_{1}x_{1}}}dx_{1}}} \\
 & {= 1 - \xi_{12}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{n = 0}^{m}{\sum\limits_{p = 0}^{n}\frac{{K^{l + k}{\zeta_{1}^{k}\zeta_{2}^{m}}{\widetilde{\gamma}}_{th}^{m}}a^{m - n}b^{p}c^{n - p}}{l!p!\left( n - p \right)!\left( m - n \right)!\left( k! \right)^{2}}}}}}}\int\limits_{0}^{\infty}{{x_{1}}^{m - 2n + k + p}e^{- {(\zeta_{1} + a\zeta_{2}{\widetilde{\gamma}}_{th})}x_{1} - \zeta_{2}{\widetilde{\gamma}}_{th}\frac{c}{x_{1}}}dx_{1}}.} \\
\end{matrix}$$

Now, by letting $\zeta_{12} \triangleq \zeta_{1} + a\zeta_{2}{\widetilde{\gamma}}_{th}$ and using Equation (3.471.9) in \[[@B32-sensors-18-01781]\], we obtain the outage probability: $$\begin{array}{ll}
P_{out\_ 1} & {= 1 - \xi_{12}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{n = 0}^{m}{\sum\limits_{p = 0}^{n}\frac{{K^{l + k}{\zeta_{1}^{k}\zeta_{2}^{m}}{\widetilde{\gamma}}_{th}^{m}}a^{m - n}b^{p}c^{n - p}}{l!p!\left( n - p \right)!\left( m - n \right)!\left( k! \right)^{2}}}}}}}2\left( \frac{\zeta_{2}{\widetilde{\gamma}}_{th}c}{\zeta_{12}} \right)^{\frac{m - 2n + p + k + 1}{2}}K_{m - 2n + p + k + 1}\left( {2\sqrt{\zeta_{12}\zeta_{2}{\widetilde{\gamma}}_{th}c}} \right)} \\
 & {= 1 - 2\xi_{12}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{n = 0}^{m}{\sum\limits_{p = 0}^{n}\frac{{K^{l + k}\zeta_{1}^{k}\left( \zeta_{2}{\widetilde{\gamma}}_{th} \right)^{\frac{3m - 2n + p + k + 1}{2}}}a^{m - n}b^{p}c^{\frac{m - p + k + 1}{2}}}{l!p!\left( n - p \right)!\left( m - n \right)!\left( k! \right)^{2}\zeta_{12}^{\frac{m - 2n + p + k + 1}{2}}}}}}}}K_{m - 2n + p + k + 1}\left( {2\sqrt{\zeta_{12}\zeta_{2}{\widetilde{\gamma}}_{th}c}} \right),} \\
\end{array}$$ where $K_{v}\left( \cdot \right)$ is the $\nu^{th}$-order modified Bessel function of the second kind.

Again, without loss of generality, we can only consider the case $i = 1$ (the SER at node 1). Obviously, the CDF of the SNR at node $S_{1}$, denoted by $F_{\gamma_{1}}\left( x \right)$, is the outage probability at $S_{1}$ when the threshold SNR is equal to *x*. Thus, it can be written as $$F_{\gamma_{1}}\left( x \right) = \Pr\left( {\frac{X_{1}X_{2}}{X_{1}X_{2}a + X_{1}^{2}a + X_{1}b + c} < x} \right) = \left\{ \begin{array}{l}
{\Pr\left( {X_{2} < \frac{xX_{1}^{2}a + xX_{1}b + xc}{X_{1} - xX_{1}a}} \right),{if}\; x \leqslant \frac{1}{a},} \\
{1,{if}\; x > \frac{1}{a}.} \\
\end{array} \right.$$

Then, we can substitute Equation ([A5](#FD31-sensors-18-01781){ref-type="disp-formula"}) into Equation ([19](#FD19-sensors-18-01781){ref-type="disp-formula"}) to get $$\begin{array}{r}
{{SER}_{1} = \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{0}^{\infty}{\frac{e^{- \theta x}}{\sqrt{x}}F_{\gamma_{1}}\left( x \right)dx} = \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{0}^{1/a}{\frac{e^{- \theta x}}{\sqrt{x}}F_{\gamma_{1}}\left( x \right)dx} + \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{1/a}^{\infty}{\frac{e^{- \theta x}}{\sqrt{x}}dx}.} \\
\end{array}$$

By using the results from the previous version, the CDF $F_{\gamma_{1}}\left( x \right)$ for $x \leq 1/a$ is given by $$\begin{matrix}
{{F_{\gamma_{1}}\left( x \right) = 1 - 2\zeta_{1}e^{- 2K}e^{- b\zeta_{2}{\widetilde{\gamma}}_{th}}}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{n = 0}^{m}{\sum\limits_{p = 0}^{n}\frac{{K^{l + k}\zeta_{1}^{k}\left( \zeta_{2}{\widetilde{\gamma}}_{th} \right)^{\frac{3m - 2n + p + k + 1}{2}}}a^{m - n}b^{p}c^{\frac{m - p + k + 1}{2}}}{l!p!\left( n - p \right)!\left( m - n \right)!\left( k! \right)^{2}\left( \zeta_{1} + a\zeta_{2}{\widetilde{\gamma}}_{th} \right)^{\frac{m - 2n + p + k + 1}{2}}}}}}}}} \\
{\times K_{m - 2n + p + k + 1}\left( {2\sqrt{\left( \zeta_{1} + a\zeta_{2}{\widetilde{\gamma}}_{th} \right)\zeta_{2}{\widetilde{\gamma}}_{th}c}} \right).} \\
\end{matrix}$$

Let us denote $\Xi\left( x \right) \triangleq 1 - F_{\gamma_{1}}\left( x \right)$. Then, Equation ([A6](#FD32-sensors-18-01781){ref-type="disp-formula"}) can be rewritten as $$\begin{array}{ll}
{SER}_{1} & {= \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{0}^{1/a}{\frac{e^{- \theta x}}{\sqrt{x}}\left\{ {1 - \Xi\left( x \right)} \right\} dx} + \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{1/a_{1}}^{\infty}{\frac{e^{- \theta x}}{\sqrt{x}}dx} = \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{0}^{\infty}{\frac{e^{- \theta x}}{\sqrt{x}}dx} - \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{0}^{1/a}{\frac{e^{- \theta x}}{\sqrt{x}}\Xi\left( x \right)dx}} \\
 & {= \frac{\omega}{2} - \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{0}^{1/a}{\frac{e^{- \theta x}}{\sqrt{x}}\Xi\left( x \right)dx}.} \\
\end{array}$$

Then, Equation ([19](#FD19-sensors-18-01781){ref-type="disp-formula"}) can be obtained by substituting $\Xi\left( x \right) \triangleq 1 - F_{\gamma_{1}}\left( x \right)$ and Equation ([A7](#FD33-sensors-18-01781){ref-type="disp-formula"}) into Equation ([18](#FD18-sensors-18-01781){ref-type="disp-formula"}).

The PDF of the random variable $X_{i}$ is given by $f_{X_{i}}\left( x \right) = \zeta_{i}e^{- K}\sum\limits_{l = 0}^{\infty}{\frac{\left( \zeta_{i}K \right)^{l}}{\left( l! \right)^{2}}x^{l}e^{- \zeta_{i}x}}$. Hence, the asymptotic outage probability at node $S_{1}$ can be written as $$\begin{array}{ll}
P_{out\_ 1}^{\infty} & {= \Pr\left( \gamma_{1}^{\infty} < \gamma_{th} \right) = \Pr\left( {\frac{X_{2}}{X_{1} + X_{2}} < {a\gamma_{th}}} \right) = \Pr\left( {X_{2} < \frac{a\gamma_{th}X_{1}}{1 - a\gamma_{th}}} \right)} \\
 & {= \int\limits_{0}^{\infty}{F_{X_{2}}\left( \frac{a\gamma_{th}x_{1}}{1 - a\gamma_{th}} \right)}f_{X_{1}}\left( x_{1} \right)dx_{1}} \\
 & {= 1 - \int\limits_{0}^{\infty}{e^{- K}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\frac{K^{l}\zeta_{2}^{m}}{l!m!}\left( {a{\widetilde{\gamma}}_{th}x_{1}} \right)^{m}}}}e^{- \zeta_{2}a{\widetilde{\gamma}}_{th}x_{1}}\zeta_{1}e^{- K}\sum\limits_{k = 0}^{\infty}\frac{\left( \zeta_{1}K \right)^{k}}{\left( k! \right)^{2}}x_{1}^{k}e^{- \zeta_{1}x_{1}}dx_{1}}} \\
 & {= 1 - \zeta_{1}e^{- 2K}\int\limits_{0}^{\infty}{{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\frac{K^{l + k}{\zeta_{1}^{k}\zeta_{2}^{m}}}{l!m!\left( k! \right)^{2}}\left( {a{\widetilde{\gamma}}_{th}x_{1}} \right)^{m}}}}}e^{- \zeta_{2}a{\widetilde{\gamma}}_{th}x_{1}}x_{1}^{k}e^{- \zeta_{1}x_{1}}dx_{1}}} \\
 & {= 1 - \zeta_{1}e^{- 2K}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\frac{K^{l + k}{\zeta_{1}^{k}\zeta_{2}^{m}}}{l!m!\left( k! \right)^{2}}\left( {a{\widetilde{\gamma}}_{th}} \right)^{m}}\int\limits_{0}^{\infty}{x_{1}^{m + k}.e^{- {(\zeta_{2}a{\widetilde{\gamma}}_{th} + \zeta_{1})}x_{1}}dx_{1}}}}}.} \\
\end{array}$$

Using the Formula (3.381.4) in \[[@B32-sensors-18-01781]\], we obtain the asymptotic formula: $$\begin{aligned}
P_{out\_ 1}^{\infty} & {= 1 - \zeta_{1}e^{- 2K}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\frac{K^{l + k}{\zeta_{1}^{k}\zeta_{2}^{m}}}{l!m!\left( k! \right)^{2}}\left( {a{\widetilde{\gamma}}_{th}} \right)^{m}}\frac{\Gamma\left( m + k + 1 \right)}{\left( \zeta_{2}a{\widetilde{\gamma}}_{th} + \zeta_{1} \right)^{m + k + 1}}}}},} \\
\end{aligned}$$ where $\Gamma\left( \cdot \right)$ is the complete gamma function.

The asymptotic throughput can be obtained by substituting Equation ([A9](#FD35-sensors-18-01781){ref-type="disp-formula"}) into the definitive formula of the throughput.

The CDF of $\gamma_{1}$ as high SNR regime can be written as $$F_{\gamma_{1}^{\infty}}\left( x \right) = \Pr\left( \gamma_{1}^{\infty} < x \right) = \Pr\left( {\frac{X_{2}}{aX_{2} + aX_{1}} < x} \right) = \left\{ \begin{array}{l}
{\Pr\left( {X_{2} < \frac{axX_{1}}{1 - ax}} \right),{if}\; x \leqslant \frac{1}{a},} \\
{1,{if}\; x > \frac{1}{a}.} \\
\end{array} \right.$$

Now, we can rewrite the formula of SER by substituting Equation ([A10](#FD36-sensors-18-01781){ref-type="disp-formula"}) into Equation ([24](#FD24-sensors-18-01781){ref-type="disp-formula"}) $$SER_{1}^{\infty} = \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{0}^{\infty}{\frac{e^{- \theta x}}{\sqrt{x}}F_{\gamma_{1}}^{\infty}\left( x \right)dx} = \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{0}^{1/a}{\frac{e^{- \theta x}}{\sqrt{x}}F_{\gamma_{1}}^{\infty}\left( x \right)dx} + \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{1/a}^{\infty}{\frac{e^{- \theta x}}{\sqrt{x}}dx}.$$

For $x \leq \frac{1}{a}$, we can apply the result of Theorem 3 with ${\widetilde{\gamma}}_{th}\left( x \right) = \frac{x}{1 - ax}$ to get $$\begin{array}{ll}
{F_{\gamma_{1}}^{\infty}\left( x \right)} & {= \Pr\left( \gamma_{1}^{\infty} < x \right) = 1 - \zeta_{1}e^{- 2K}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\frac{K^{l + k}{\zeta_{1}^{k}\zeta_{2}^{m}}}{l!m!\left( k! \right)^{2}}\left( {a{\widetilde{\gamma}}_{th}\left( x \right)} \right)^{m}}\frac{\Gamma\left( m + k + 1 \right)}{\left\lbrack \zeta_{2}a{\widetilde{\gamma}}_{th}\left( x \right) + \zeta_{1} \right\rbrack^{m + k + 1}}}}}} \\
 & {= 1 - \zeta_{1}e^{- 2K}{\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}\frac{K^{l + k}{\zeta_{1}^{k}\zeta_{2}^{m}a^{m}}}{l!m!\left( k! \right)^{2}}.\frac{x^{m}\left( 1 - ax \right)^{k + 1}\Gamma\left( m + k + 1 \right)}{\left\lbrack \left( \zeta_{2} - \zeta_{1} \right)ax + \zeta_{1} \right\rbrack^{m + k + 1}}}}}} \\
 & {= 1 - \zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{1}^{k}\zeta_{2}^{m}a^{m + v}\left( - 1 \right)^{v}x^{m + v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}\left\lbrack \left( \zeta_{2} - \zeta_{1} \right)ax + \zeta_{1} \right\rbrack^{m + k + 1}}}}}}} \\
 & {= 1 - \zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{2}^{m}a^{m + v}\left( - 1 \right)^{v}x^{m + v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}\zeta_{1}^{m + 1}\left\lbrack {1 + \frac{\left( \zeta_{2} - \zeta_{1} \right)ax}{\zeta_{1}}} \right\rbrack^{m + k + 1}}}}}}.} \\
\end{array}$$

The aymptotic SER is now obtained by substituting Equation ([A12](#FD38-sensors-18-01781){ref-type="disp-formula"}) into Equation ([A11](#FD37-sensors-18-01781){ref-type="disp-formula"}): $$\begin{array}{l}
{SER_{1}^{\infty} = \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\int\limits_{0}^{\infty}{\frac{e^{- \theta x}}{\sqrt{x}}dx}} \\
{- \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{1}^{k}\zeta_{2}^{m}a^{v - k - 1}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}\left( \zeta_{2} - \zeta_{1} \right)^{m + k + 1}}}}}}\int\limits_{0}^{\infty}\frac{x^{m + v - \frac{1}{2}}e^{- \theta x}dx}{\left\lbrack {x + \frac{\zeta_{1}}{\left( \zeta_{2} - \zeta_{1} \right)a}} \right\rbrack^{m + k + 1}}} \\
{= \frac{\omega}{2} - \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{1}^{k}\zeta_{2}^{m}a^{v - k - 1}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}\left( \zeta_{2} - \zeta_{1} \right)^{m + k + 1}}}}}}\int\limits_{0}^{1/a}\frac{x^{m + v - \frac{1}{2}}e^{- \theta x}dx}{\left\lbrack {x + \frac{\zeta_{1}}{\left( \zeta_{2} - \zeta_{1} \right)a}} \right\rbrack^{m + k + 1}}.} \\
\end{array}$$

We consider three cases:

**Case 1:** ξ 2 =ξ 1 .

Equation ([A13](#FD39-sensors-18-01781){ref-type="disp-formula"}) can be rewritten as $$\begin{array}{l}
{{SER}_{1}^{\infty} = \frac{\omega}{2} - \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{1}^{k}\zeta_{2}^{m}a^{m + v}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}\zeta_{1}^{m + 1}}}}}}\int\limits_{0}^{1/a}{x^{m + v - \frac{1}{2}}e^{- \theta x}dx}} \\
{= \frac{\omega}{2} - \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{1}^{k}\zeta_{2}^{m}a^{m + v}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}\zeta_{1}^{m + 1}}}}}}\theta^{- m - v - \frac{1}{2}}\gamma\left( {m + v + \frac{1}{2},\frac{\theta}{a}} \right),} \\
\end{array}$$ where $\gamma\left( \alpha,z \right) \triangleq \int_{0}^{z}{e^{- t}t^{\alpha - 1}dt}$ is the incomplete gamma function as defined in Equation (8.350.1) in \[[@B32-sensors-18-01781]\] and the last equality comes from the Formula (3.381.1) in \[[@B32-sensors-18-01781]\].

**Case 2:** ξ 2 \<ξ 1 .

For all $x \in \left\lbrack 0,\frac{1}{a} \right\rbrack$, we have $x < \frac{1}{a}.\frac{\xi_{1}}{\xi_{1} - \xi_{2}}$. Hence, we can apply the Taylor's series expansion of $\frac{1}{\left\lbrack {1 - \frac{\left( \xi_{1} - \xi_{2} \right)ax}{\xi_{1}}} \right\rbrack^{m + k + 1}}$ as in \[[@B37-sensors-18-01781]\]: $$\frac{1}{\left\lbrack {1 - \frac{\left( \xi_{1} - \xi_{2} \right)ax}{\xi_{1}}} \right\rbrack^{m + k + 1}} = \sum\limits_{p = 0}^{\infty}{\frac{\left( m + k + 1 \right)_{p}}{p!}\left\lbrack \frac{\left( \xi_{1} - \xi_{2} \right)ax}{\xi_{1}} \right\rbrack^{p}},$$ where $\left( \cdot \right)_{p}$ is the Pochhammer symbol, which is defined as $$\left( \alpha \right)_{p} \triangleq \left\{ \begin{array}{l}
{1,{if}\mspace{720mu} p = 0,} \\
{\alpha\left( \alpha + 1 \right)...\left( \alpha + p - 1 \right),{if}\mspace{720mu} p \neq 0.} \\
\end{array} \right.$$

Now, we can obtain the following equation by substituting Equation ([A15](#FD41-sensors-18-01781){ref-type="disp-formula"}) into Equation ([A13](#FD39-sensors-18-01781){ref-type="disp-formula"}) and applying (3.381.1) in \[[@B32-sensors-18-01781]\]: $$\begin{array}{l}
{SER_{1}^{\infty} = \frac{\omega}{2}} \\
{- \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{1}^{k}\zeta_{2}^{m}a^{v - k - 1}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}\left( \zeta_{2} - \zeta_{1} \right)^{m + k + 1}}}}}}\int\limits_{0}^{1/a}\frac{x^{m + v - \frac{1}{2}}e^{- \theta x}dx}{\left\lbrack {x + \frac{\zeta_{1}}{\left( \zeta_{2} - \zeta_{1} \right)a}} \right\rbrack^{m + k + 1}}} \\
{= \frac{\omega}{2} - \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{2}^{m}a^{m + v}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}\zeta_{1}^{m + 1}}}}}}\int\limits_{0}^{1/a}\frac{x^{m + v - \frac{1}{2}}e^{- \theta x}dx}{\left\lbrack {1 + \frac{\left( \zeta_{2} - \zeta_{1} \right)ax}{\zeta_{1}}} \right\rbrack^{m + k + 1}}} \\
{= \frac{\omega}{2} - \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{2}^{m}a^{m + v}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}\zeta_{1}^{m + 1}}}}}}} \\
{\times \int\limits_{0}^{1/a}{x^{m + v - \frac{1}{2}}e^{- \theta x}\sum\limits_{p = 0}^{\infty}{\frac{\left( m + k + 1 \right)_{p}}{p!}\left\lbrack \frac{\left( \xi_{1} - \xi_{2} \right)ax}{\xi_{1}} \right\rbrack^{p}}dx}} \\
{= \frac{\omega}{2} - \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\sum\limits_{p = 0}^{\infty}{\left( m + k + 1 \right)_{p}\left\lbrack \left( \xi_{1} - \xi_{2} \right) \right\rbrack^{p}}\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{2}^{m}a^{m + v + p}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}p!\zeta_{1}^{m + p + 1}}}}}}} \\
{\times \int\limits_{0}^{1/a}{x^{m + v + p - \frac{1}{2}}e^{- \theta x}dx}} \\
{= \frac{\omega}{2} - \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\sum\limits_{p = 0}^{\infty}{\left( m + k + 1 \right)_{p}\left\lbrack \left( \xi_{1} - \xi_{2} \right) \right\rbrack^{p}}\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}\frac{K^{l + k}\zeta_{2}^{m}a^{m + v + p}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}p!\zeta_{1}^{m + p + 1}\theta^{m + v + p + \frac{1}{2}}}}}}}} \\
{\times \gamma\left( {m + v + p + \frac{1}{2},\frac{\theta}{a}} \right).} \\
\end{array}$$

**Case 3:** ξ 2 \>ξ 1 .

Using the same approach as in Case 2, with the Taylor's series expansion $\frac{1}{\left\lbrack {1 + \frac{\left( \xi_{2} - \xi_{1} \right)ax}{\xi_{1}}} \right\rbrack^{m + k + 1}} = \sum\limits_{p = 0}^{\infty}{\left( - 1 \right)^{p}\frac{\left( m + k + 1 \right)_{p}}{p!}\left\lbrack \frac{\left( \xi_{2} - \xi_{1} \right)ax}{\xi_{1}} \right\rbrack^{p}}$, we get the following result: $$\begin{array}{l}
{SER_{1}^{\infty} = \frac{\omega}{2} - \frac{\omega\sqrt{\theta}}{2\sqrt{\pi}}\zeta_{1}e^{- 2K}\sum\limits_{l = 0}^{\infty}{\sum\limits_{k = 0}^{\infty}{\sum\limits_{m = 0}^{l}{\sum\limits_{v = 0}^{k + 1}{\sum\limits_{p = 0}^{\infty}{\left( m + k + 1 \right)_{p}\left\lbrack \left( \xi_{1} - \xi_{2} \right) \right\rbrack^{p}}\begin{pmatrix}
{k + 1} \\
v \\
\end{pmatrix}}}}}} \\
{\times \frac{K^{l + k}\zeta_{2}^{m}a^{m + v + p}\left( - 1 \right)^{v}\Gamma\left( m + k + 1 \right)}{l!m!\left( k! \right)^{2}p!\zeta_{1}^{m + p + 1}\theta^{m + v + p + \frac{1}{2}}}\gamma\left( {m + v + p + \frac{1}{2},\frac{\theta}{a}} \right).} \\
\end{array}$$

![System model.](sensors-18-01781-g001){#sensors-18-01781-f001}

![The energy harvesting and information processing in the system model.](sensors-18-01781-g002){#sensors-18-01781-f002}

![Outage probability versus $P/N_{0}$.](sensors-18-01781-g003){#sensors-18-01781-f003}

![Achievable throughput versus $P/N_{0}$](sensors-18-01781-g004){#sensors-18-01781-f004}

![Outage probability versus *κ*.](sensors-18-01781-g005){#sensors-18-01781-f005}

![Achievable throughput versus *κ*.](sensors-18-01781-g006){#sensors-18-01781-f006}

![Outage probability versus *α*.](sensors-18-01781-g007){#sensors-18-01781-f007}

![Achievable throughput versus *α*.](sensors-18-01781-g008){#sensors-18-01781-f008}

![Effect of *κ* and modulation scheme on SER.](sensors-18-01781-g009){#sensors-18-01781-f009}

![Effect of *α* and modulation scheme on SER.](sensors-18-01781-g010){#sensors-18-01781-f010}

![Effect of *α* and modulation scheme on SER.](sensors-18-01781-g011){#sensors-18-01781-f011}

sensors-18-01781-t001_Table 1

###### 

Simulation parameters.

  Symbol                  Parameter Names                Values
  ----------------------- ------------------------------ ---------------
  $\eta$                  Energy harvesting efficiency   0.7
  $\lambda_{1}$           Mean of $|g_{1}|^{2}$          0.5
  $\lambda_{2}$           Mean of $|g_{2}|^{2}$          0.5
  *K*                     Rician K-factor                3
  $P/N_{0}$               Source-power-to-noise ratio    0--50 dB
  $\kappa = \kappa_{r}$   Hardware impairment levels     0, 0.15, 0.25
  *R*                     Source transmission rate       1.5 bps/Hz
