PEP Analysis of AF Relay NOMA Systems Employing Order Statistics of Cascaded Channels

Abstract

The precise error performance analysis is challenging for non-orthogonal multiple access (NOMA) systems due to nonlinear successive interference cancellation (SIC) processing among NOMA users. In this paper, the pairwise error probability (PEP) performance of different users is investigated for relay NOMA simultaneous wireless information and power transfer (SWIPT) systems. By employing the order statistics theory, we obtain the ordered probability density function of the cascaded channel through Source-to-Relay-to-User links. Then we derive the analytical closed-form PEP expressions for NOMA users. To obtain the approximate closed-form PEP, we explore the finite series representation of the power of the modified Bessel function to replace the integrand terms. Monte Carlo simulation results show that the approximate analytical PEP of each user is basically in agreement with the simulated PEP. Furthermore, on the basis of the closed-form PEP, the influence of relevant system parameters on the error performance is examined via numerical simulations, which manifests that the choice of power allocation coefficients should be balanced between the users’ channel conditions and the demanded quality of service.

1. Introduction

Compared to the current fourth generation (4G) wireless networks, the fifth generation (5G) networks are expected to support massive connectivity of devices so as to satisfy the growth of internet of things. The channel multiple access schemes determine the number of users to share the wireless resources, which mainly include orthogonal multiple access (OMA) and non-orthogonal multiple access (NOMA) [1]. It is difficult for OMA to meet the spectral efficiency requirements of 5G services. Power-domain NOMA and code-domain NOMA are two categories of NOMA schemes. For the former, multiple users are allocated with different power levels at the same time and channel frequency. Specifically, higher power is assigned to the users with worse channel gains; and lower power is allocated to the ones with better channel gains. In terms of compatibility with existing OMA schemes, power-domain NOMA is more promising.

Cooperative communication adopts relay strategies to improve system reliability, amplify-and-forward (AF) and decode-and-forward (DF) are two common types of processing protocols by the relay terminals [2]. A prior study [3] derived the outage probability and diversity multiplexing trade-off curve for an uplink cooperative system with DF protocol. Cooperative NOMA transmission has been proposed to enhance the overall performance of communication systems, as well as reducing the transmit power [4]. The status of recent research in power-domain NOMA-based cooperative networks are surveyed in [5]. In order to prolong the lifetime of energy constrained relay systems, NOMA can be integrated with simultaneous wireless information and power transfer (SWIPT) [6,7].

For AF relay systems, data detection at the destination only requires the information of the cascaded channel through the source-to-relay-to-destination link. Therefore, the estimation algorithms of the cascaded channel are proposed in [8,9]. Additionally, the statistical properties of the cascaded relay fading channels are investigated in [10], where the cascaded channel is modeled as a double Gaussian channel, i.e., a product of two complex Gaussian channels.

There have been numerous studies that focus on the performance analysis of cooperative NOMA systems [11,12,13]. One prior study [11] derived the exact and asymptotic outage probabilities for NOMA AF systems. The authors in [12] studied NOMA-based downlink AF relaying network under Nakagami-m fading and derived the closed-form expressions of the outage probability, with imperfect channel state information taken into account. In [13], a two-stage relay selection scheme was proposed and a closed-form expression of the outage probability was derived. Almost all the previous literature concentrated on evaluating the outage probability. There are few reports about the error performance analysis for cooperative NOMA systems. In [14], the error rate was analyzed in non-cooperative NOMA systems over Nakagami-m fading channels.

To the best of our knowledge, there is no closed-form solution for the pairwise error probability (PEP) expressions of NOMA users based on order statistics of the cascaded channel for relay NOMA communication systems in the literature. In this paper, we mainly analyze the PEP performance of NOMA systems with multiple AF relay–user pairs. Our specific contributions of this paper are summarized as follows:

• We derive the approximate closed-form PEP expressions of the first and the kth users over Rayleigh fading channels, respectively. In particular, the ordered probability density function (PDF) of the cascaded channel is utilized to solve the integrals. Simulation results reveal the consistency between the derived closed-form PEP expressions and their corresponding Monte Carlo simulations.
• We attempt to substitute $K_1^p\left(x\right)$, namely, the pth power of modified Bessel function $K_1\left(x\right)$, with an approximately equivalent finite series representation since the presence of the power of the modified Bessel function of the second kind, $K_1^p\left(x\right)$, makes the integrals about the PEP of NOMA users intractable. In existing literature, there is no appropriate finite series representation of $K_1^p\left(x\right)$. Therefore, we derive an effective and simplified series representation of $K_1^p\left(x\right)$ based on the series representation of $K_1\left(x\right)$ in [15]. The major difficulty lies in that a smaller x value will result in a significant truncation error between the approximated ${\tilde{K}}_1^p\left(x\right)$ and the actual $K_1^p\left(x\right)$ once the choice of the finite order is improper. To some extent, the effectiveness of the derived closed-form PEP expression is dependent on the accurate representation of finite series of $K_1^p\left(x\right)$.
• The impacts of some related parameters, such as power allocation coefficients of NOMA users and power splitting factor related to energy harvesting, can be evaluated through the simulations by using the approximate closed-form PEP expressions.

The rest of the paper is organized as follows. First, the system model is given in Section 2 for the considered cooperative NOMA system. The analytical PEP expression of the first user is derived in Section 3. In Section 4, we present the approximate finite series representation of $K_1^p\left(x\right)$. The closed-form PEP of the kth successive interference cancellation (SIC) useris derived in Section 5. The numerical and simulation results are provided for analytical confirmation in Section 6. Finally, the paper is concluded in Section 7.

2. System Model

We consider a downlink relay NOMA system with K users assisted by K relays as depicted in Figure 1, where all are equipped with a single antenna. NOMA technique enables various users to manipulate individual information under the conditions of identical time and frequency resources. The base station (BS) broadcasts the superimposed information to K users with the help of respective relays, which are equipped with the function of simultaneous wireless information and power transfer (SWIPT). We adopt a half-duplex amplify-and-forward (AF) relay with no direct link due to severe obstruction or shadow. Considering energy harvesting, the power splitter (PS) scheme at the relay is adopted, which can split the receiver signal power into two parts comprising of energy harvesting and information processing [16].

At the BS with known power allocation coefficient ${\alpha }_k$ of the kth user, the superimposition signal of NOMA users is transmitted, which is denoted as

$$s={\displaystyle \sum _{k=1}^K}\sqrt{{\alpha }_k{P}_t}{x}_k,$$

(1)

where $P_t$ is the total power at the BS. $x_k$ is the data symbol of the kth user. Moreover, ${\alpha }_1>{\alpha }_2>\cdots >{\alpha }_K$ and ${\sum }_{k=1}^K{\alpha }_k=1$.

The users accomplish signal reception through two time slots. In the first time slot, BS broadcasts the signal to the relay, thus the received signal at the relay is

$$y_R^k=h_R^{k}s+n_r,$$

(2)

where $h_R^k\sim \mathbb{CN}(0,{\lambda }_r)$ represents the channel coefficients from the source to relay. $n_r\sim \mathbb{CN}(0,{\sigma }_{n_r}^2)$ indicates the additive white Gaussian noise at the relay.

After power splitting at the relay, the signals about energy harvesting and information forwarding are written as, respectively

$$y_R^{E,k}=\sqrt{\rho }(h_R^{k}s+n_r),$$

(3)

$$y_R^{I,k}=\sqrt{1-\rho }(h_R^{k}s+n_r)+n_c,$$

(4)

where $\rho \in (0,1)$ is a PS factor, and $n_c\sim \mathbb{CN}(0,{\sigma }_{n_c}^2)$ is the conversion noise from RF to baseband.

In the second time slot, with the relay operating the AF protocol and forwarding the information signal, the received signal of the kth user can be written as

$$y_U^k=G_A{h}_U^k{y}_R^{I,k}+n_k=G_A^{\prime }{h}_R^k{h}_U^{k}s+G_A^{\prime }{h}_U^k{n}_r+G_A{h}_U^k{n}_c+n_k,$$

(5)

where $G_A^{\prime }=G_A\sqrt{1-\rho }$, and $h_U^k\sim \mathbb{CN}(0,{\lambda }_k)$ represents the channel coefficients from the relay to users. Moreover, $n_k\sim \mathbb{CN}(0,{\sigma }_{n_c}^2)$ is the additive white Gaussian noise at the destination. Assume ${\lambda }_r={\lambda }_k=\lambda$ and ${\sigma }_{n_r}^2={\sigma }_{n_c}^2={\sigma }_{n_k}^2={\sigma }_n^2$. It can be observed that the knowledge of individual channels $h_R^k$ and $h_U^k$ is not required for data detection. As a result, the cascade channel $h_R^k{h}_U^k$ is considered in the following analysis. The amplification factor $G_A$ is denoted as

$${\mathrm{G}}_A=\sqrt{{\displaystyle \frac{P_r}{(1-\rho ){\left|h_R^k\right|}^2{P}_t+(2-\rho ){\sigma }_n^2}}},$$

(6)

where $P_r=\eta \rho P_t{\left|h_R^k\right|}^2$ is the harvested power at the relay in the second time slot. $\eta \sim (0,1]$ denotes energy conversion efficiency.

Considering the fixed-gain amplification factor, $G_A$ can be rewritten as

$$G_A=\sqrt{{\displaystyle \frac{\eta \rho \lambda P_t}{(1-\rho )\lambda P_t+(2-\rho ){\sigma }_n^2}}}.$$

(7)

For each receiver, successive interference cancellation (SIC) is utilized to remove the interference and implement the accurate signal detection. Specifically, for the kth user, the manipulation of SIC is to successively detect the signals $x_1,\cdots ,x_{k-1}$ and subtract them from the received signal. The signals of the remaining users are treated as interference.

3. PEP Derivation of the First User

The SIC process is not required for the first user. Plugging (1) into (5), the received signal of the first user can be formulated as

$$y_U^1=G_A^{\prime }\left(\sqrt{{\alpha }_1{P}_t}{x}_1+{\displaystyle \sum _{k=2}^K}\sqrt{{\alpha }_k{P}_t}{x}_k\right)h_R^1{h}_U^1+G_A^{\prime }{h}_U^1{n}_r+G_A{h}_U^1{n}_c+n_1,$$

(8)

where ${\sum }_{k=2}^K\sqrt{{\alpha }_k{P}_t}{x}_k$ is viewed as the interference component during the process of detecting the first user. The conditional pairwise error probability (PEP) of the first user can be expressed as

$$\mathrm{Pr}(x_1\to {\check{x}}_1|h_R^1{h}_U^1)=\mathrm{Pr}({\left|y_U^1-G_A^{\prime }\sqrt{{\alpha }_1{P}_t}{h}_R^1{h}_U^1{\check{x}}_1\right|}^2\le {\left|y_U^1-G_A^{\prime }\sqrt{{\alpha }_1{P}_t}{h}_R^1{h}_U^1{x}_1\right|}^2),$$

(9)

where $x_1\to {\check{x}}_1$ denotes the event in which $x_1$ is erroneously detected into the symbol ${\check{x}}_1$.

Inserting (8) into (9), (9) can be rewritten as

$$\begin{array}{ll}\mathrm{Pr}& (x_1\to {\check{x}}_1|h_R^1{h}_U^1)\hfill \\ & =\mathrm{Pr}\left(2\mathfrak{Re}\left\{h_R^1{h}_U^1{\check{\Delta }}_1{n}_{to}^{*}\right\}\le -{\left|h_R^1{h}_U^1\right|}^2{G}_A^{\prime }\left[\sqrt{{\alpha }_1{P}_t}{\left|{\check{\Delta }}_1\right|}^2+2\mathfrak{Re}\left\{{\check{\Delta }}_1{\displaystyle \sum _{k=2}^K}\sqrt{{\alpha }_k{P}_t}{x}_k^{*}\right\}\right]\right),\hfill \end{array}$$

(10)

where $n_{to}^{*}=G_A^{\prime }{h_U^1}^{*}{n}_r^{*}+G_A{h_U^1}^{*}{n}_c^{*}+n_1^{*}$, ${\check{\Delta }}_1=x_1-{\check{x}}_1,x_1\ne {\check{x}}_1$. Moreover, $2\mathfrak{Re}\left\{h_R^1{h}_U^1{\check{\Delta }}_1{n}_{to}^{*}\right\}\sim \mathcal{N}\left(0,2{\sigma }_n^2|h_R^1{h}_U^1{|}^2|{\check{\Delta }}_1{|}^2(1+|G_A^{\prime }{|}^2|h_U^1{|}^2+|G_A{|}^2|h_U^1{|}^2)\right)$, which further turns Equation (10) into

$$\mathrm{Pr}(x_1\to {\check{x}}_1|{\omega }_1)=Q\left({\displaystyle \frac{{\omega }_1{\xi }_1}{{\mu }_1}}\right),$$

(11)

where ${\omega }_1=\left|h_R^1{h}_U^1\right|$ is the amplitude of the cascaded channel. $Q(·)$ denotes the Q-function. ${\xi }_1$ and ${\mu }_1$ are given by

$${\xi }_1=G_A^{\prime }\left[\sqrt{{\alpha }_1{P}_t}{\left|{\check{\Delta }}_1\right|}^2+2\mathfrak{Re}\left\{{\check{\Delta }}_1{\displaystyle \sum _{k=2}^K}\sqrt{{\alpha }_k{P}_t}{x}_k^{*}\right\}\right],$$

(12)

$${\mu }_1=\sqrt{2}{\sigma }_n\left|{\check{\Delta }}_1\right|\sqrt{1+(2-\rho )|G_A{|}^2{\left|h_U^1\right|}^2}.$$

(13)

Without taking the instantaneous channel into account, ${\mu }_1$ can be further written as

$${\mu }_1=\sqrt{2}{\sigma }_n\left|{\check{\Delta }}_1\right|\sqrt{1+(2-\rho )|G_A{|}^2\lambda }.$$

(14)

Generally, the ordered PDF of the cascaded channel ${\omega }_k=\left|h_R^k{h}_U^k\right|$ regarding the kth user is represented as [17]

$$f_k\left({\omega }_k\right)=A_{k}f\left({\omega }_k\right){\left[F\left({\omega }_k\right)\right]}^{k-1}{[1-F\left({\omega }_k\right)]}^{K-k},$$

(15)

where $f\left({\omega }_k\right)$ and $F\left({\omega }_k\right)$ represent the PDF and cumulative distribution function (CDF) of ${\omega }_k$, and

$$A_k={\displaystyle \frac{K!}{(k-1)!(K-k)!}}.$$

(16)

Assume that all the separate channels undergo Rayleigh fading, which means ${\omega }_k$ is the amplitude of the product of two complex Gaussian random variables [18], so the PDF and CDF of ${\omega }_k$ are expressed as

$$f\left({\omega }_k\right)={\displaystyle \frac{4{\omega }_k}{{\lambda }_r{\lambda }_k}}{K}_0\left({\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}\right),$$

(17)

$$F\left({\omega }_k\right)=1-{\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}{K}_1\left({\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}\right),$$

(18)

where $K_0(·)$ and $K_1(·)$ are the modified Bessel function of the second kind with zero order and first order, respectively.

Consequently, with ${\lambda }_r={\lambda }_k=\lambda$ assumed, the ordered PDF of the first user is given by

$$f_1\left({\omega }_1\right)=A_1{\displaystyle \frac{4{\omega }_1}{{\lambda }^2}}{K}_0\left({\displaystyle \frac{2\omega }{\lambda }}\right){\left({\displaystyle \frac{2{\omega }_1}{\lambda }}{K}_1\left({\displaystyle \frac{2{\omega }_1}{\lambda }}\right)\right)}^{K-1}.$$

(19)

It follows that the exact PEP of the first user is solved as

$$\begin{array}{cc}\hfill \mathrm{Pr}(x_1\to {\check{x}}_1)=& {\int }_0^{\infty }{f}_1\left({\omega }_1\right)Q\left({\displaystyle \frac{{\omega }_1{\xi }_1}{{\mu }_1}}\right)d{\omega }_1\hfill \\ \hfill \stackrel{x=2{\omega }_1/\lambda }{=}& K{\int }_0^{\infty }{x}^K{K}_0\left(x\right)K_1{\left(x\right)}^{K-1}Q\left(ax\right)dx,\hfill \end{array}$$

(20)

where $A_1=K$, $a=\lambda {\xi }_1/\left(2{\mu }_1\right)$.

Proposition 1.

The approximate closed-form PEP of the first user can be represented as

$$\mathrm{Pr}(x_1\to {\check{x}}_1)\approx {\displaystyle \frac{1}{2}}-{\displaystyle \frac{a}{\sqrt{2\pi }}}exp\left({\displaystyle \frac{K^2}{4a^2}}\right){\displaystyle \sum _{q=0}^Q}{C}_q{a}^{-(q+1)}\Gamma (q+1)D_{-(q+1)}\left({\displaystyle \frac{K}{a}}\right),$$

(21)

where $\Gamma (·)$ is the gamma function, $D_v\left(z\right)$ denotes a parabolic Cylinder function, and

$$C_q=\left\{\begin{array}{c}1,q=0\hfill \\ {\displaystyle \frac{1}{q}}{\displaystyle \sum _{z=1}^q}(Kz-q+z){\displaystyle \sum _{\ell =z}^{Q^{\prime }}}\Lambda (1,\ell ,z)C_{q-z},q\ge 1.\hfill \end{array}\right.$$

(22)

Proof of Proposition 1.

Equation (20) can be further formulated as

$$\begin{array}{cc}\hfill \mathrm{Pr}(x_1\to {\check{x}}_1)=& K{\int }_0^{\infty }{x}^K{K}_0\left(x\right)K_1{\left(x\right)}^{K-1}Q\left(ax\right)dx\hfill \\ \hfill =& {\displaystyle \frac{-1}{2}}{\int }_0^{\infty }\mathrm{erfc}\left({\displaystyle \frac{a}{\sqrt{2}}}x\right)\mathrm{d}\left(x^K{K}_1^K\left(x\right)\right)\hfill \\ \hfill =& -{\displaystyle \frac{1}{2}}\left[\mathrm{erfc}\left({\displaystyle \frac{a}{\sqrt{2}}}x\right)x^K{K}_1^K\left(x\right){|}_0^{\infty }-{\int }_0^{\infty }{x}^K{K}_1^K\left(x\right)\mathrm{d}\left(\mathrm{erfc}\left({\displaystyle \frac{a}{\sqrt{2}}}x\right)\right)\right]\hfill \\ \hfill =& {\displaystyle \frac{1}{2}}-{\displaystyle \frac{a}{\sqrt{2\pi }}}{\int }_0^{\infty }{x}^K{K}_1^K\left(x\right)exp\left(-{\displaystyle \frac{a^2}{2}}{x}^2\right)\mathrm{d}x.\hfill \end{array}$$

(23)

In (23), $Q\left(ax\right)$ is replaced by ${\displaystyle \frac{1}{2}}\mathrm{erfc}\left(ax/\sqrt{2}\right)$. It can be easily derived that ${\left(x^K{K}_1^K\left(x\right)\right)}^{\prime }=-K{x}^K{K}_0\left(x\right)K_1{\left(x\right)}^{K-1}$. Moreover, the first term on the right side of the third equal sign can be computed as

$$\left\{\begin{array}{c}\underset{x\to \infty }{lim}\mathrm{erfc}\left({\displaystyle \frac{a}{\sqrt{2}}}x\right)=0\hfill \\ \underset{x\to \infty }{lim}{x}^K{K}_1^K\left(x\right)=\underset{x\to \infty }{lim}{x}^K{\left(\sqrt{{\displaystyle \frac{\pi }{2x}}}{e}^{-x}\right)}^K=0,\hfill \end{array}\right.$$

(24)

$$\left\{\begin{array}{c}\underset{x\to 0}{lim}\mathrm{erfc}\left({\displaystyle \frac{a}{\sqrt{2}}}x\right)=1\hfill \\ \underset{x\to 0}{lim}{x}^K{K}_1^K\left(x\right)=\underset{x\to 0}{lim}{x}^K{x}^{-K}=1.\hfill \end{array}\right.$$

(25)

In the above limit derivations, we exploit the asymptotic expressions of $K_1\left(x\right)$, namely $K_1\left(x\right)\approx x^{-1}$ when $x\to 0$, while $K_1\left(x\right)\approx \sqrt{{\displaystyle \frac{\pi }{2x}}}{e}^{-x}$ when $x\to \infty$ [19].

Therefore, based on Equation (34) in Section 4, $K_1{\left(x\right)}^K$ can be further expressed as

$$K_1^K\left(x\right)\approx {\displaystyle \frac{e^{-Kx}}{x^K}}{\displaystyle \sum _{q=0}^Q}{C}_q{x}^q,$$

(26)

where $C_q$ is given by (22) in the next section. It is worth noting that $C_q$ can be stored as a vector with constants in advance.

Subsequently, plugging (26) into (23) yields

$$\begin{array}{cc}\hfill \mathrm{Pr}(x_1\to {\check{x}}_1)\approx & {\displaystyle \frac{1}{2}}-{\displaystyle \frac{a}{\sqrt{2\pi }}}{\displaystyle \sum _{q=0}^Q}{C}_q{\int }_0^{\infty }{x}^{q}exp\left(-Kx-{\displaystyle \frac{a^2}{2}}{x}^2\right)\mathrm{d}x\hfill \\ \hfill =& {\displaystyle \frac{1}{2}}-{\displaystyle \frac{a}{\sqrt{2\pi }}}exp\left({\displaystyle \frac{K^2}{4a^2}}\right){\displaystyle \sum _{q=0}^Q}{C}_q{a}^{-(q+1)}\Gamma (q+1)D_{-(q+1)}\left({\displaystyle \frac{K}{a}}\right),\hfill \end{array}$$

(27)

where the last equation is obtained according to ([20], Equation (3.462.1)). □

4. Approximate Finite Series Representation of ${\mathit{K}}_{\mathbf{1}}^{\mathit{p}}\left(\mathit{x}\right)$

In order to facilitate the integral in (23), it is essential to explore the equivalent series representation of $K_1^p\left(x\right)$. According to [15] (Equation (24)), the series representation of $K_1\left(x\right)$ can be given by

$$K_1\left(x\right)={\displaystyle \frac{e^{-x}}{x}}{\displaystyle \sum _{q=0}^{\infty }}{\displaystyle \sum _{\ell =q}^{\infty }}\Lambda (1,\ell ,q)x^q,$$

(28)

with the coefficients being

$$\Lambda (1,\ell ,q)={\displaystyle \frac{{(-1)}^{q+1}\Gamma (\ell -{\displaystyle \frac{1}{2}})L(\ell ,q)}{2^{2-q}\Gamma (\ell +{\displaystyle \frac{3}{2}})\ell !}},$$

(29)

where $L(\ell ,q)$ denotes the Lah numbers ([15], Equation (8)), namely, $L(\ell ,q)=\left(\genfrac{}{}{0pt}{}{\ell -1}{q-1}\right){\displaystyle \frac{\ell !}{q!}}$, for $\ell ,q>0$. The other values can be taken as: $L(0,0)=1$; for $\ell >0$, $L(\ell ,0)=0$ and $L(\ell ,1)=\ell !$. Let $a_q={\sum }_{\ell =q}^{\infty }\Lambda (1,\ell ,q)$, (28) can be simplified by

$$K_1\left(x\right)={\displaystyle \frac{e^{-x}}{x}}{\displaystyle \sum _{q=0}^{\infty }}{a}_q{x}^q.$$

(30)

Hence, it is straightforward to obtain the pth power of $K_1\left(x\right)$ as

$$K_1^p\left(x\right)={\displaystyle \frac{e^{-px}}{x^p}}{\left({\displaystyle \sum _{q=0}^{\infty }}{a}_q{x}^q\right)}^p,p\ge 1.$$

(31)

Utilizing the mathematical operation regarding power series raised to powers [21], the power of $K_1\left(x\right)$ can be further converted into

$$K_1^p\left(x\right)={\displaystyle \frac{e^{-px}}{x^p}}{\displaystyle \sum _{m=0}^{\infty }}{c}_m{x}^m,$$

(32)

where

$$c_m=\left\{\begin{array}{c}{\left(a_0\right)}^p,m=0\hfill \\ {\displaystyle \frac{1}{m{a}_0}}{\displaystyle \sum _{z=1}^m}(pz-m+z)a_z{c}_{m-z},m\ge 1\hfill \end{array}\right.$$

(33)

with $a_0={\displaystyle \sum _{\ell =0}^{\infty }}\Lambda (1,\ell ,0)$, $a_z={\displaystyle \sum _{\ell =z}^{\infty }}\Lambda (1,\ell ,z)$.

From the perspective of practical implementation, it is normal to consider the approximate representation in terms of finite series terms of the Bessel function $K_1\left(x\right)$ [15]. However, there is no existing work about the finite series representation of $K_1^p\left(x\right)$. In the following, we will examine whether $K_1^p\left(x\right)$ can be truncated by only using a small series order with high precision instead of infinite series in (32).

It is apparent that there exist two infinity series in the expression of $K_1^p\left(x\right)$ in Equation (32). The outer series can be approximated by the finite order Q, namely,

$$K_1^p\left(x\right)={\displaystyle \frac{e^{-px}}{x^p}}{\displaystyle \sum _{m=0}^Q}{c}_m{x}^m+ϵ={\tilde{K}}_1^p\left(x\right)+ϵ,$$

(34)

where $ϵ$ indicates the truncation error, which is negligible with a small Q, as illustrated in

Table 1. The Truncation Error of $K_1^p\left(x\right)$ with different choices of $Q^{\prime }$ and Q.

${\mathit{K}}_1^{\mathit{p}}\left(\mathit{x}\right)$
Truncation Errror $\mathit{ϵ}$	${\mathit{K}}_{\mathbf{1}}\left(\mathit{x}\right)$	${\mathit{K}}_{\mathbf{1}}^{\mathbf{2}}\left(\mathit{x}\right)$	${\mathit{K}}_{\mathbf{1}}^{\mathbf{3}}\left(\mathit{x}\right)$
$(Q^{\prime },Q)=(5,5)$	$9.4640\times {10}^{-4}$	0.0015	0.0022
$(Q^{\prime },Q)=(8,8)$	$4.8373\times {10}^{-4}$	$8.9121\times {10}^{-4}$	0.0013
$(Q^{\prime },Q)=(10,1)$	0.0154	$7.6091\times {10}^{-4}$	0.0011
$(Q^{\prime },Q)=(10,2)$	0.0173	$7.0637\times {10}^{-4}$	$9.8693\times {10}^{-4}$
$(Q^{\prime },Q)=(12,1)$	0.0158	$6.3690\times {10}^{-4}$	$8.2680\times {10}^{-4}$
$(Q^{\prime },Q)=(12,2)$	0.0204	$6.1105\times {10}^{-4}$	$7.8361\times {10}^{-4}$

. Moreover, the smaller Q is, the less coefficients are required to store. The series $a_z$ within the coefficient $c_m$ is truncated by the parameter $Q^{\prime }$, namely,

$$a_q\approx {\displaystyle \sum _{\ell =q}^{Q^{\prime }}}\Lambda (1,\ell ,q),q=0,\cdots ,Q^{\prime },$$

(35)

where $Q^{\prime }$ and Q are positive integers, $Q^{\prime }\ge Q$, and $a_0\approx {\sum }_{\ell =0}^{Q^{\prime }}\Lambda (1,\ell ,0)=1$ [15].

Remark 1.

It is worth noting that two identical finite orders are chosen to replace two infinity series about $K_1\left(x\right)$ in [15] (Equation (24)). Unlike that, we attempt to exploit two different finite orders to truncate the infinite series of $K_1^p\left(x\right)$, which can improve the accuracy when $p>1$. Apparently, the former can be regarded as our special case.

Remark 2.

It is likely for Q and $Q^{\prime }$ to take different values for different p in order to minimize the truncation error. We need to search for proper values about Q and $Q^{\prime }$ for each p. Once these parameters are fixed, the coefficients $c_m$ can be stored in advance as shown in

Table 2. Coefficients $c_m$ in (34) with different $Q^{\prime }$ and Q.

${\mathit{c}}_{\mathit{m}}$	${\mathit{c}}_0$	${\mathit{c}}_1$	${\mathit{c}}_2$	${\mathit{c}}_3$	${\mathit{c}}_4$
$K_1\left(x\right):(Q^{\prime },Q)=(4,4)$	1	0.8889	−0.3429	0.1016	−0.0106
$K_1^2\left(x\right):(Q^{\prime },Q)=(10,2)$	1	1.9048	−0.5023
$K_1^2\left(x\right):(Q^{\prime },Q)=(12,2)$	1	1.92	−0.6471
$K_1^3\left(x\right):(Q^{\prime },Q)=(10,2)$	1	2.8571	0.6070
$K_1^3\left(x\right):(Q^{\prime },Q)=(12,2)$	1	2.88	0.4117

.

The accuracy of series approximation of $K_1^p\left(x\right)$ exerts a direct impact on the analytical PEP derivation. We define the normalized truncation errors as $ϵ={\sum }_x^{}|{K_1\left(x\right)}^p-{\tilde{K}}_1^p\left(x\right)|/{\sum }_x^{}\left|{K_1\left(x\right)}^p\right|$ in the simulation, where the true ${K_1\left(x\right)}^p$ is generated by MATLAB command "besselk", and ${\tilde{K}}_1^p\left(x\right)={\displaystyle \frac{e^{-px}}{x^p}}{\sum }_{m=0}^Q{c}_m{x}^m$. In Table 1, take discrete values on $x=0.01:0.1:2$ for example, where $0.1$ is the adopted sample interval. From Table 1, we observe that the case with $Q^{\prime }=12$ and $Q=2$ has the least truncation error for $p=2,3$. In addition, $K_1^p\left(x\right)$ with $Q^{\prime }=Q$ is better than $Q^{\prime }>Q$ for $p=1$ whereas $K_1^p\left(x\right)$ exhibits the smaller truncation error with $Q^{\prime }>Q$ than the case $Q^{\prime }=Q$ for $p>1$.

For specific $K_1^p\left(x\right)$, $c_m$ can be generated in advance given $Q^{\prime }$ and Q, as shown in Table 2. It is worth noting that the coefficients $c_0$ and $c_1$ in the cases $Q^{\prime }=10,12$ and $Q=1$ are identical to the corresponding $c_0$ and $c_1$ with the case $Q=2$, hence, the former cases are not shown in this table. Moreover, Table 2 also illustrates $c_m$ in (34) with $p=1$ and $Q^{\prime }=Q$ has completely the same coefficients as $a_{1,k,q}$ listed in [15] (Table 1). In any case, it is noticed that $c_0$=1 remains invariant.

5. PEP Derivation of the $\mathit{k}$th SIC User

From the second user to the last user, the SIC is performed by each user. During each SIC process, the output of the kth user is

$$y_U^k=G_A^{\prime }\left(\sqrt{{\alpha }_k{P}_t}{x}_k+{\displaystyle \sum _{j=k+1}^K}\sqrt{{\alpha }_j{P}_t}{x}_j+{\displaystyle \sum _{u=1}^{k-1}}\sqrt{{\alpha }_u{P}_t}{\widehat{\Delta }}_u\right)h_R^k{h}_U^k+G_A^{\prime }{h}_U^k{n}_r+G_A{h}_U^k{n}_c+n_k,$$

(36)

where ${\widehat{\Delta }}_u=x_u-{\widehat{x}}_u$ is the interference cancellation error of the previously detected signal ${\widehat{x}}_u$.

After the same manipulation as above, the conditional PEP of the kth user is given by

$$\mathrm{Pr}(x_k\to {\check{x}}_k|{\omega }_k)=Q\left({\displaystyle \frac{{\omega }_k{\xi }_k}{{\mu }_k}}\right),$$

(37)

where ${\omega }_k=\left|h_R^k{h}_U^k\right|$, ${\xi }_k$ and ${\mu }_k$ are represented as

$${\xi }_k=G_A^{\prime }\left\{\sqrt{{\alpha }_{k}P}{\left|{\check{\Delta }}_k\right|}^2+2\mathfrak{Re}\left\{{\check{\Delta }}_k\left[{\displaystyle \sum _{u=1}^{k-1}}\sqrt{{\alpha }_u{P}_t}{\widehat{\Delta }}_u^{*}+{\displaystyle \sum _{j=k+1}^K}\sqrt{{\alpha }_j{P}_t}{x}_j^{*}\right]\right\}\right\},$$

(38)

$${\mu }_k=\sqrt{2}{\sigma }_n\left|{\check{\Delta }}_k\right|\sqrt{1+(2-\rho )|G_A{|}^2\lambda },$$

(39)

where ${\check{\Delta }}_k=x_k-{\check{x}}_k$. It is worth highlighting that since all previous user symbols have been detected, ${\xi }_K$ for the last user is given by

$${\xi }_K=G_A^{\prime }\left\{\sqrt{{\alpha }_K{P}_t}{\left|{\check{\Delta }}_K\right|}^2+2\mathfrak{Re}\left\{{\check{\Delta }}_K\left[{\displaystyle \sum _{u=1}^{K-1}}\sqrt{{\alpha }_u{P}_t}{\widehat{\Delta }}_u^{*}\right]\right\}\right\},$$

(40)

where ${\widehat{\Delta }}_u$ is zero if the perfect SIC is taken into account.

From (19), the ordered PDF of the kth user is

$$\begin{array}{rr}{f}_k\left({\omega }_k\right)=A_k{\displaystyle \frac{4{\omega }_k}{{\lambda }_r{\lambda }_k}}{K}_0\left({\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}\right)& {\left(1-{\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}{K}_1\left({\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}\right)\right)}^{k-1}\\ & \times {\left({\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}{K}_1\left({\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}\right)\right)}^{K-k}.\end{array}$$

(41)

Taking advantage of binomial expansion, we have

$${\left(1-{\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}{K}_1\left({\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}\right)\right)}^{k-1}={\displaystyle \sum _{\ell =0}^{k-1}}\left(\genfrac{}{}{0pt}{}{k-1}{\ell }\right){\left(-{\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}{K}_1\left({\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}\right)\right)}^{\ell }.$$

(42)

Then (41) can be rewritten as

$$f_k\left({\omega }_k\right)=A_k{\displaystyle \sum _{\ell =0}^{k-1}}\left(\genfrac{}{}{0pt}{}{k-1}{\ell }\right){(-1)}^{\ell }{\displaystyle \frac{4{\omega }_k}{{\lambda }_r{\lambda }_k}}{\left({\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}\right)}^{K^{\prime }-1}{K}_0\left({\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}\right)K_1^{K^{\prime }-1}\left({\displaystyle \frac{2{\omega }_k}{\sqrt{{\lambda }_r{\lambda }_k}}}\right),$$

(43)

where $K^{\prime }=K-k+\ell +1$.

With ${\lambda }_r={\lambda }_k=\lambda$ assumed, it turns out from (43)

$$f_k\left({\omega }_k\right)=A_k{\displaystyle \sum _{\ell =0}^{k-1}}\left(\genfrac{}{}{0pt}{}{k-1}{\ell }\right){(-1)}^{\ell }{\displaystyle \frac{4{\omega }_k}{{\lambda }^2}}{\left({\displaystyle \frac{2{\omega }_k}{\lambda }}\right)}^{K^{\prime }-1}{K}_0\left({\displaystyle \frac{2{\omega }_k}{\lambda }}\right)K_1^{K^{\prime }-1}\left({\displaystyle \frac{2{\omega }_k}{\lambda }}\right).$$

(44)

Therefore, the exact PEP of the kth user can be given by

$$\begin{array}{cc}\hfill \mathrm{Pr}(x_k\to {\check{x}}_k)& ={\int }_0^{\infty }{f}_k\left({\omega }_k\right)Q\left({\displaystyle \frac{{\omega }_k{\xi }_k}{{\mu }_k}}\right)d{\omega }_k\hfill \\ \hfill \stackrel{x=2{\omega }_k/\lambda }{=}& A_k{\displaystyle \sum _{\ell =0}^{k-1}}\left(\genfrac{}{}{0pt}{}{k-1}{\ell }\right){(-1)}^{\ell }{\int }_0^{\infty }{x}^{K^{\prime }}{K}_0\left(x\right)K_1^{K^{\prime }-1}\left(x\right)Q\left(bx\right)dx,\hfill \end{array}$$

(45)

where $b=\lambda {\xi }_k/\left(2{\mu }_k\right)$.

Proposition 2.

The approximate closed-form PEP of the kth user can be derived as

$$\mathrm{Pr}(x_k\to {\check{x}}_k)\approx {\displaystyle \frac{1}{2}}-A_k{\displaystyle \frac{b}{\sqrt{2\pi }}}{\displaystyle \sum _{\ell =0}^{k-1}}{T}_{\ell }exp\left({\displaystyle \frac{{K^{\prime }}^2}{4b^2}}\right){\displaystyle \sum _{q=0}^Q}{M}_q{b}^{-(q+1)}\Gamma (q+1)D_{-(q+1)}\left({\displaystyle \frac{K^{\prime }}{b}}\right),$$

(46)

where $T_{\ell }={\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{k-1}{\ell }\right){(-1)}^{\ell }}{K^{\prime }}}$, and

$$M_q=\left\{\begin{array}{c}1,q=0\hfill \\ {\displaystyle \frac{1}{q}}{\displaystyle \sum _{z=1}^q}\left(K^{\prime }z-q+z\right){\displaystyle \sum _{i=z}^{Q^{\prime }}}\Lambda (1,i,z)M_{q-z},otherwise.\hfill \end{array}\right.$$

(47)

Proof of Proposition 2.

The integral term in Equation (45) can be derived as

$$\begin{array}{ll}{I}_k& ={\int }_0^{\infty }{x}^{K-k+\ell +1}{K}_0\left(x\right)K_1^{K^{\prime }-1}\left(x\right)Q\left(bx\right)dx\\ & ={\displaystyle \frac{1}{2}}{\int }_0^{\infty }{x}^{K^{\prime }}{K}_0\left(x\right)K_1^{K^{\prime }-1}\left(x\right)\mathrm{erfc}\left({\displaystyle \frac{bx}{\sqrt{2}}}\right)dx\hfill \\ & ={\displaystyle \frac{1}{2K^{\prime }}}-{\displaystyle \frac{b}{K^{\prime }\sqrt{2\pi }}}{\int }_0^{\infty }{x}^{K^{\prime }}{K}_1^{K^{\prime }}\left(x\right)exp\left(-{\displaystyle \frac{b^2}{2}}{x}^2\right)\mathrm{d}x.\hfill \end{array}$$

(48)

Therefore, the exact PEP expression can be obtained as

$$\begin{array}{cc}\hfill \mathrm{Pr}(x_k\to {\check{x}}_k)& =A_k{\displaystyle \sum _{\ell =0}^{k-1}}\left(\genfrac{}{}{0pt}{}{k-1}{\ell }\right){(-1)}^{\ell }\hfill \\ & \times \left({\displaystyle \frac{1}{2K^{\prime }}}-{\displaystyle \frac{b}{K^{\prime }\sqrt{2\pi }}}{\int }_0^{\infty }{x}^{K^{\prime }}{K}_1^{K^{\prime }}\left(x\right)exp\left(-{\displaystyle \frac{b^2}{2}}{x}^2\right)\mathrm{d}x\right).\hfill \end{array}$$

(49)

The above closed-form expression has a similar manipulation process to (23). It should be pointed out that the approximate series representation of $K_1^{K^{\prime }}\left(x\right)$ is given by

$$K_1^{K^{\prime }}\left(x\right)\approx {\displaystyle \frac{e^{-K^{\prime }x}}{x^{K^{\prime }}}}{\displaystyle \sum _{q=0}^Q}{M}_q{x}^q$$

(50)

with the coefficients $M_q$ manifested in (47). By plugging (50) into (48), we have

$$\begin{array}{cc}\hfill I_k& \approx {\displaystyle \frac{1}{2K^{\prime }}}-{\displaystyle \frac{b}{K^{\prime }\sqrt{2\pi }}}{\displaystyle \sum _{q=0}^Q}{M}_q{\int }_0^{\infty }{x}^{q}exp\left(-K^{\prime }x-{\displaystyle \frac{b^2}{2}}{x}^2\right)\mathrm{d}x\hfill \\ & ={\displaystyle \frac{1}{2K^{\prime }}}-{\displaystyle \frac{b}{K^{\prime }\sqrt{2\pi }}}exp\left({\displaystyle \frac{{K^{\prime }}^2}{4b^2}}\right){\displaystyle \sum _{q=0}^Q}{M}_q{b}^{-(q+1)}\Gamma (q+1)D_{-(q+1)}\left({\displaystyle \frac{K^{\prime }}{b}}\right).\hfill \end{array}$$

(51)

Consequently, inserting (51) into (45) yields

$$\mathrm{Pr}(x_k\to {\check{x}}_k)={\displaystyle \frac{N_k}{2}}-A_k{\displaystyle \frac{b}{\sqrt{2\pi }}}{\displaystyle \sum _{\ell =0}^{k-1}}{T}_{\ell }exp\left({\displaystyle \frac{{K^{\prime }}^2}{4b^2}}\right){\displaystyle \sum _{q=0}^Q}{M}_q{b}^{-(q+1)}\Gamma (q+1)D_{-(q+1)}\left({\displaystyle \frac{K^{\prime }}{b}}\right),$$

(52)

where $N_k=A_k{\displaystyle \sum _{\ell =0}^{k-1}}{\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{k-1}{\ell }\right){(-1)}^{\ell }}{K^{\prime }}}$. Moreover, leveraging on [20] (Equation 0.160), the summation terms in $N_k$ can be rewritten as

$$\begin{array}{cc}\hfill {\displaystyle \sum _{\ell =0}^{k-1}}{\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{k-1}{\ell }\right){(-1)}^{\ell }}{K^{\prime }}}& ={\displaystyle \sum _{\ell =0}^{k-1}}\left(\genfrac{}{}{0pt}{}{k-1}{\ell }\right){(-1)}^{\ell }{\displaystyle \frac{\Gamma (\ell +K-k+1)}{\Gamma (\ell +K-k+2)}}\hfill \\ & =B(k,K-k+1).\hfill \end{array}$$

(53)

Combined with the key property of Beta function $B(a,b)={\displaystyle \frac{(a-1)!(b-1)!}{(a+b-1)!}}$ [22], it follows that

$$N_k=A_{k}B(k,K-k+1)=A_k{\displaystyle \frac{(k-1)!(K-k)!}{K!}}=1.$$

(54)

 □

The proof of (46) is completed.

6. Simulation Results

This section validates the derived PEP expressions in the preceding sections. Consider a two-user NOMA SWIPT system over Rayleigh flat fading channels. BPSK modulation is adopted. The channel variance is assumed to be $\lambda =1$. In addition, other parameters are fixed as $({\alpha }_1,{\alpha }_2)=(0.72,0.28)$, $\rho =0.6$, and $\eta =0.8$, unless otherwise specified. The simulated curves are conducted through Monte Carlo simulations.

It should be pointed out that the aforementioned parabolic Cylinder function $D_v\left(z\right)$ ([20], Equation (9.240)) can be converted into

$$D_v\left(z\right)=2^{{\displaystyle \frac{1}{4}}+{\displaystyle \frac{v}{2}}}{z}^{-{\displaystyle \frac{1}{2}}}{W}_{{\displaystyle \frac{1}{4}}+{\displaystyle \frac{v}{2}},-{\displaystyle \frac{1}{4}}}\left({\displaystyle \frac{z^2}{2}}\right),$$

(55)

where W(·) is the Whittaker function, which is available in MATLAB toolbox.

In Figure 2, we evaluate the derived closed-form PEP (21) for the first user compared with the simulated PEP. It is worth pointing out that the PEP of the first user is associated with $x_2$ when ${\check{\Delta }}_1$ is fixed. The solid lines represent the simulated PEP, while the markers denote the analytical PEP. The average analytical PEP is obtained by taking the average of both cases $x_2=\pm 1$. As seen in Figure 2a,b, the accuracy of the closed-formed PEP is related to the value $(Q,Q^{\prime })$. When $(Q,Q^{\prime })=(4,4)$, the analytical PEP with $x_2=1$ is evidently inconsistent with the simulated PEP at SNR = 30 dB, which is due to the approximation error of finite series representation of $K_1{\left(x\right)}^K$. However, the analytical PEP almost agrees well with the simulated PEP when $(Q,Q^{\prime })=(2,8)$.

Figure 3 verifies the accuracy of the derived analytical closed-form PEP expression for the second user. Whichever modulation is adopted, the closed-form expression (46) in Proposition 2 is consistent with the simulated PEP acquired by random Monte Carlo simulation over Rayleigh fading channels. To preclude the influence of the approximate error of finite series representation of $K_1{\left(x\right)}^{K^{\prime }}$, (49) is provided for more reliable reference. Specifically, for (49), we use the existing MATLAB “Besselk” function to implement $K_1{\left(x\right)}^{K^{\prime }}$ as a benchmark. As anticipated, the simulated PEP and the analytical in Equation (49) completely overlap. Meanwhile, Equation (46) Proposition 2 is totally in agreement with the other two curves when $(Q,Q^{\prime })=(4,4)$.

The impact of several parameter values on PEP is straightforwardly presented. The following three figures are simulated at SNR = 20 dB, and the average analytical PEP is evaluated for the first user. First, we observe the performance of PEP with varying PS factors $\rho$ for two users as seen in Figure 4. Fortunately, $\rho$ = 0.5 or 0.6 simultaneously corresponds to the optimal performance of two users.

In order to evaluate the effect of power allocation coefficients, we alter ${\alpha }_1$ of the first user. For the second user, its power allocation coefficient will make a corresponding change due to ${\alpha }_1+{\alpha }_2=1$. As depicted in Figure 5, the first user with a weaker channel condition performs gradually better as ${\alpha }_1$ slowly increases. On the other hand, the second user with stronger channel quality becomes worse due to its reduced ${\alpha }_2$. The choice of power allocation coefficient depends on the required quality of service for different users.

The effect of energy conversion efficiency on the error performance is explored as shown in Figure 6. For any user, only a slight change in PEP happens when $\eta$ varies from 0.5 to 0.9.

7. Conclusions

An approximate closed-form PEP analysis is derived for cooperative NOMA SWIPT systems over Rayleigh fading channels. Based on the ordered statistics of the composite channel across the Source-Relay-User link, the PDF and CDF of $|h_R^k{h}_U^k|$ are related to the modified Bessel functions of the second kind, which renders the PEP analysis very challenging. In order to obtain the closed-form PEP, the finite series representation of $K_1{\left(x\right)}^k$ is introduced. We verify the effectiveness of the approximately analytical PEP for NOMA users, which is expressed by finite series representation. Numerical results substantiate that our method agrees well with Monte Carlo simulations. It is worth pointing out that for more than three users, it needs some time to search the suitable pairs of $(Q,Q\prime )$ to match the finite series representation of $K_1^p\left(x\right)$. Specifically, assuming K = k = 3, the accuracy of PEP calculations for the third user depends on the joint precision of finite series approximations of $K_1\left(x\right)$, $K_1^2\left(x\right)$, and $K_1^3\left(x\right)$, which can be seen from Equation (49). Our future work will explore an effective approach to determine the joint pairs of $(Q,Q\prime )$ for different $K_1^p\left(x\right)$ to adapt more users.