# Analysis of HARQ-IR over Time-Correlated Rayleigh Fading Channels

Zheng Shi, Haichuan Ding, Shaodan Ma, and Kam-Weng Tam Manuscript received
January 14, 2015; revised May 29, 2015; accepted July 19, 2015. The associate
editor coordinating the review of this paper and approving it for publication
was M. Elkashlan.Zheng Shi, Shaodan Ma and Kam-Weng Tam are with the
Department of Electrical and Computer Engineering, University of Macau, Macao
<EMAIL_ADDRESS><EMAIL_ADDRESS>kentam@umac.mo).Haichuan
Ding was with University of Macau, and is now with the Department of
Electrical and Computer Engineering, University of Florida, U.S.A. (email:
dhcbit@gmail.com).This work was supported by the Research Committee of
University of Macau under grants: MYRG078(Y1-L2)-FST12-MSD and
MYRG101(Y1-L3)-FST13-MSD.

###### Abstract

In this paper, performance of hybrid automatic repeat request with incremental
redundancy (HARQ-IR) over Rayleigh fading channels is investigated. Different
from prior analysis, time correlation in the channels is considered. Under
time-correlated fading channels, the mutual information in multiple HARQ
transmissions is correlated, making the analysis challenging. By using
polynomial fitting technique, probability distribution function of the
accumulated mutual information is derived. Three meaningful performance
metrics including outage probability, average number of transmissions and long
term average throughput (LTAT) are then derived in closed-forms. Moreover,
diversity order of HARQ-IR is also investigated. It is proved that full
diversity can be achieved by HARQ-IR, i.e., the diversity order is equal to
the number of transmissions, even under time-correlated fading channels. These
analytical results are verified by simulations and enable the evaluation of
the impact of various system parameters on the performance. Particularly, the
results unveil the negative impact of time correlation on the outage and
throughput performance. The results also show that although more transmissions
would improve the outage performance, they may not be beneficial to the LTAT
when time correlation is high. Optimal rate design to maximize the LTAT is
finally discussed and significant LTAT improvement is demonstrated.

###### Index Terms:

Hybrid automatic repeat request, incremental redundancy, time correlation,
Rayleigh fading channels, outage probability.

##  I Introduction

Over the last decade, wireless data traffic has experienced an explosive
growth. Contrary to this boom in data traffic, transmission reliability and
throughput of wireless channels are limited by unideal propagation environment
[1]. As a combination of forward error control and automatic repeat request
(ARQ), hybrid automatic repeat request (HARQ) has been proved as an effective
technique to improve transmission reliability and boost system throughput. It
thus has been adopted in various wireless standards, such as High Speed Packet
Access (HSPA) and Long Term Evolution (LTE) [2, 3]. Generally, there exist
three kinds of HARQ schemes, i.e., Type-I HARQ, HARQ with chase combining
(HARQ-CC) and HARQ with incremental redundancy (HARQ-IR). In Type-I HARQ, the
erroneously received packets are discarded and each retransmitted packet is
decoded independently. In HARQ-CC and HARQ-IR, the previously failed packets
are stored and combined with the packets received in subsequent
retransmissions for decoding. Specifically, the same packet is retransmitted
in each transmission attempt in HARQ-CC scheme, while redundant information is
incrementally transmitted in each HARQ round in the case of HARQ-IR. By
exploiting additional coding gain, HARQ-IR is able to achieve a higher link
throughput than Type-I HARQ and HARQ-CC [2]. The focus of this paper is thus
turned to HARQ-IR scheme.

Performance of HARQ-IR under various wireless systems has been investigated in
the literature [4, 5, 6, 7, 8, 9, 10]. To name a few, delay and throughput of
a HARQ-IR enabled multicast system are investigated and scaling laws with
respect to the number of users are discovered based on an upper bound of
outage probability in [4]. Aiming at throughput maximization, rate allocation
and adaptation for HARQ-IR are discussed in [5]. For quasi-static fading
environments (i.e., channel coefficients are constant during multiple HARQ
rounds), power efficiency of HARQ-IR is concerned and optimal power allocation
is obtained to minimize outage-limited average transmission power in [6].
Considering the same quasi-static environments as [6], average rate of HARQ-IR
enabled spectrum sharing networks is analyzed in [7]. By expressing outage
probability through k-fold convolution, throughput of network coded HARQ-IR
with arbitrary number of users is derived in [8]. In [9], an optimal rate
adaptation policy is proposed for cooperative HARQ-IR with a multi-bit
feedback channel. Dynamic programming is then employed to find the optimal
rate for maximizing the throughput of the outage-constrained transmission.
Moreover, based on dominant term approximation and upper bounds of outage
probability, energy-delay-tradeoff (EDT) is analyzed for both one-way and two-
way relaying systems with HARQ-IR in [10]. Noticing the lack of exact
analytical result on outage probability of HARQ-IR, an analytical approach is
proposed to derive the outage probability in a closed form through the
generalized Fox’s H function in [11]. Unfortunately, all of the prior studies
are conducted for either quasi-static fading channels or fast fading channels
(i.e., channel coefficients in multiple HARQ rounds are independent and
identical distributed). The results are not applicable to time-correlated
fading channels, which usually occur when the transceiver has low-to-medium
mobility [12, 13]111In a dense scattering environment, the time-correlation
between two channel amplitudes with time spacing of $\tau$ is quantified by
$\rho={J_{0}}^{2}\left({2\pi{f_{c}}\tau v{c^{-1}}}\right)$ where
$J_{0}(\cdot)$ denotes the zero-th order Bessel function of the first kind,
$f_{c}$ represents the carrier frequency, $c$ is the speed of light and $v$
refers to the moving speed of a mobile terminal [14]. Taking a 3GPP LTE system
as an example, the successive transmissions are not carried out in adjacent
time slots, and the time spacing between two successive HARQ transmissions is
$\tau=8$ms [15]. When the LTE system is operated at a carrier frequency of
$f_{c}=2.6$GHz [16], the time correlation coefficient $\rho$ between the
channel amplitudes in two HARQ transmissions is 0.83 and 0.45 for moving
speeds of 5 km/h and 10 km/h, respectively..

Considering the wide occurrence of time correlation in fading channels in
practice, it is necessary and meaningful to analyze the performance of HARQ-IR
operating over time-correlated fading channels. However the analysis is very
challenging because of the difficulty in handling a product of multiple
correlated random variables (RVs). Notice that the analysis is essentially
different from [12, 13] where HARQ-CC is analyzed and the sum of multiple
correlated RVs is concerned. In this paper, we consider HARQ-IR operating over
time-correlated Rayleigh fading channels. The accumulated mutual information
after multiple HARQ rounds is first expressed as a logarithm function of a
product of multiple shifted correlated signal-to-noise ratios. By using
polynomial fitting technique, probability distribution function (PDF) of the
accumulated mutual information is derived as a product of Gamma distribution
and a correction polynomial. Outage probability, average number of
transmissions and long term average throughput (LTAT) are then obtained in
closed-forms. Moreover, diversity order of HARQ-IR is also analyzed. It is
proved that full diversity can be achieved, i.e., the diversity order is equal
to the number of transmissions, even under time-correlated fading channels.
The impact of channel time correlation on the performance is also investigated
and optimal rate design to maximize the throughput is finally discussed. The
results reveal that time correlation of the channel causes negative effect on
the outage and throughput performance. More HARQ rounds may not be beneficial
to the throughput under highly correlated channels, although they do improve
the outage performance.

The remainder of this paper is organized as follows. In Section II, a point-
to-point HARQ-IR enabled system operating over time-correlated Rayleigh fading
channels is introduced. In Section III, outage probability, average number of
transmission and LTAT are derived in closed-forms by using polynomial fitting
technique. Section IV analyzes the diversity order of HARQ-IR over time-
correlated fading channels. The analytical results are verified through Monte-
Carlo simulations, and the impact of time correlation on the performance of
HARQ-IR and optimal rate design are then discussed in Section V. Section VI
finally concludes this paper.

##  II System Model

Consider a point-to-point system with one source and one destination, as shown
in Fig. 1. To enhance the transmission reliability, HARQ-IR protocol is
adopted here. Notice that most of the prior research on HARQ-IR is carried out
over quasi-static or fast fading channels [10, 11, 17]. Different from the
prior analysis, time-correlated fading channels are considered in this paper.
Specifically, the HARQ-IR protocol and channel model are introduced in the
following.

Figure 1: System model.

###  II-A HARQ-IR Protocol

Following the HARQ-IR protocol, prior to transmit a message with $b$ bits, the
source first encodes the message into $M$ packets, each with $L$ symbols. The
$M$ packets are denoted as $B_{1},B_{2},\cdots,B_{M}$, as shown in Fig. 1.
Thus the maximum allowable number of transmissions for this message is limited
to $M$. Then the source transmits the $M$ packets one by one in multiple HARQ
rounds till the destination succeeds to decode the message. If the destination
succeeds/fails to decode the message, a positive/negative acknowledgement
(ACK/NACK) message will be fed back to the source. At the source, after
receiving an NACK message from the destination, the subsequent packet will be
delivered in the next HARQ round until the maximum allowable number of
transmissions is reached or an ACK message is received. Once either of these
two events happens, the source will initiate the transmission of a new message
following the same procedure.

At the destination, the received packets are corrupted by fading channels and
additive white Gaussian noises. The corrupted packets associated with
$B_{1},B_{2},\cdots,B_{M}$ are denoted as $\\{C_{1},C_{2},...,C_{M}\\}$. The
channel decoder attempts to recover the message based on all the previously
received packets. More specifically, after $k$ HARQ rounds, the received
packets from $C_{1}$ to $C_{k}$ are utilized for decoding the message at the
destination. If the message can be successfully recovered, an ACK message will
be fed back to the source. Otherwise, a feedback of failure notification will
be sent to the source, namely, NACK message.

###  II-B Channel Model

Denote the modulated signal from the $l$th packet $B_{l}$ as ${\bf{x}}_{l}$.
The signal received at the destination in the $l$th HARQ round is written as

${\bf{y}}_{l}=h_{l}{\bf{x}}_{l}+{\bf{n}}_{l}$ (1)

where ${\bf{n}}_{l}$ represents complex additive white Gaussian noise with
zero mean and variance $\mathfrak{N}_{l}$, i.e.,
${\bf{n}}_{l}\sim\mathcal{CN}({\bf{0}},\mathfrak{N}_{l}{{\bf{I}}})$, $h_{l}$
denotes the block Rayleigh fading channel coefficient in the $l$th HARQ round,
i.e., the magnitude of $h_{l}$ obeys a Rayleigh distribution, such that
${\left|{{h_{l}}}\right|}\sim Rayleigh({2{\sigma_{l}}^{2}})$ and the
expectation of the squared channel magnitude is
${\mathrm{E}}\\{{{\left|{{h_{l}}}\right|}^{2}}\\}=2\sigma_{l}^{2}$. The PDF of
$|h_{l}|$ is given by

${f_{\left|{{h_{l}}}\right|}}\left(x\right)=\frac{x}{{{\sigma_{l}}^{2}}}\exp\left({-\frac{{{x^{2}}}}{{2{\sigma_{l}}^{2}}}}\right),x\in[0,+\infty).$
(2)

In this paper, time correlation of the channels is considered. Herein, a
widely used correlated channel model [18] is adopted, that is, the channel
coefficients $\left|{\bf{h}}\right|=\\{|h_{1}|,|h_{2}|,\cdots,|h_{K}|\\}$ are
modeled as a multivariate Rayleigh distribution with generalized correlation.
The joint PDF of $\left|{\bf{h}}\right|$ follows as 222In fact, the joint PDF
of channel amplitudes (II-B) is consistent with the conditional PDF of channel
power gains given in [12, Eq. 5]. However, our channel model is different from
[19, Eq. 11] where fast fading channels with imperfect channel state
information are considered.

$\displaystyle{f_{\left|{\bf{h}}\right|}}\left({\left|{{h_{1}}}\right|={x_{1}},\cdots,\left|{{h_{K}}}\right|={x_{K}}}\right)=\prod\limits_{k=1}^{K}{\frac{1}{{{\sigma_{k}}^{2}\left({1-{\lambda_{k}}^{2}}\right)}}}$
$\displaystyle\times\int\nolimits_{t=0}^{\infty}{{e^{-\left({1+\sum\limits_{k=1}^{K}{\frac{{{\lambda_{k}}^{2}}}{{1-{\lambda_{k}}^{2}}}}}\right)t}}\prod\limits_{k=1}^{K}{{x_{k}}{e^{-\frac{{{x_{k}}^{2}}}{{2{\sigma_{k}}^{2}\left({1-{\lambda_{k}}^{2}}\right)}}}}}}$
$\displaystyle\times{}_{0}{F_{1}}\left({;1;\frac{{{x_{k}}^{2}{\lambda_{k}}^{2}t}}{{2{\sigma_{k}}^{2}{{\left({1-{\lambda_{k}}^{2}}\right)}^{2}}}}}\right)dt,\,|{\lambda_{k}}|<1,$
(3)

where ${}_{0}{F_{1}}\left({\cdot}\right)$ denotes the confluent hypergeometric
limit function and $\lambda_{k}$ indicates the correlation degree of the
channels. Under this model, the time correlation coefficient between channels
associated with the $k$th and the $l$th HARQ rounds is given as

$\displaystyle{\rho_{{{{{k}}}},{{{{l}}}}}}$
$\displaystyle=\frac{{{\rm{E}}\left({{{\left|{{h_{k}}}\right|}^{2}}{{\left|{{h_{l}}}\right|}^{2}}}\right)-{\rm{E}}\left({{{\left|{{h_{k}}}\right|}^{2}}}\right){\rm{E}}\left({{{\left|{{h_{l}}}\right|}^{2}}}\right)}}{{\sqrt{{\rm{Var}}\left({{{\left|{{h_{k}}}\right|}^{2}}}\right){\rm{Var}}\left({{{\left|{{h_{l}}}\right|}^{2}}}\right)}}}$
$\displaystyle={\lambda_{k}}^{2}{\lambda_{l}}^{2}<1,\quad k,l\in[1,K],$ (4)

where the notation ${\rm{Var}}(x)$ denotes the variance of $x$. Notice that
this channel model is applicable to the time-correlated channels with
correlation coefficients of ${\rho_{k,l}}<1$. For fully correlated Rayleigh
fading channels (i.e., quasi-static Rayleigh fading channels), the channel
remains static in all HARQ rounds, i.e., $|h_{1}|=|h_{2}|=\cdots=|h_{M}|\sim
Rayleigh(2{\sigma}^{2})$ and $\rho_{k,l}=1$. The analysis of HARQ-IR over
fully correlated fading channels has been conducted in [6, 7] and our analysis
focuses on the time-correlated channels with ${\rho_{k,l}}<1$.

From (1), the received signal-to-noise ratio (SNR) at the destination in the
$l$th HARQ round is written as

$\gamma_{l}=\frac{{{{\left|{h_{l}}\right|}^{2}}P_{l}}}{{\mathfrak{N}_{l}}}$
(5)

where $P_{l}$ is the transmitted signal power in the $l$th HARQ round. The
accumulated mutual information at the destination after $K$ HARQ rounds is
then given by

$I_{K}^{IR}=\sum\limits_{l=1}^{K}{{I_{l}}}$ (6)

where $I_{l}$ represents the mutual information acquired from the $l$th HARQ
round and is given as

${I_{l}}={\log_{2}}\left({1+\gamma_{l}}\right).$ (7)

##  III Performance Analysis

###  III-A Performance Metrics

To investigate the performance of HARQ-IR over time-correlated Rayleigh fading
channels, three widely adopted metrics including outage probability, average
number of transmissions and long term average throughput (LTAT) are discussed
here.

####  III-A1 Outage Probability

In each HARQ round, the destination combines the current received packet with
all the previously received packets for joint decoding. When the accumulated
mutual information at the destination is less than the transmission rate
$\mathcal{R}$, an outage (i.e., the failure of the decoding) would occur. The
outage probability after $K$ HARQ rounds $P_{out}^{IR}\left(K\right)$ is then
given as

$P_{out}^{IR}\left(K\right)=\Pr\left({I_{K}^{IR}<\mathcal{R}}\right).$ (8)

This outage probability can well approximate the error probability when a
capacity achieving coding is adopted, and is of great importance in the
analysis of HARQ schemes [11].

####  III-A2 Average Number of Transmissions

HARQ scheme is a combination of forward error control and automatic repeat
request. To enhance the transmission reliability, each message may be
retransmitted through multiple HARQ rounds. When the channel condition is
good, few retransmissions are sufficient for successful decoding, while more
retransmissions are needed over a poor channel. From statistical point of
view, it is meaningful to know the average transmission time for each message,
which can been well characterized by the average number of transmissions.
Given the maximum allowable number of transmissions $M$, the average number of
transmissions $\bar{\mathcal{N}}$ is expressed as [20]

$\bar{\mathcal{N}}=1+\sum\limits_{K=1}^{M-1}{P_{out}^{IR}\left(K\right)}.$ (9)

####  III-A3 LTAT

As an effective metric to characterize the system throughput of HARQ schemes,
the LTAT given the transmission rate $\mathcal{R}$ and the maximum number of
transmissions $M$ is defined as [21]

$\bar{\mathcal{R}}=\frac{{\mathcal{R}\left({1-P_{out}^{IR}\left(M\right)}\right)}}{{\bar{\mathcal{N}}}}.$
(10)

Clearly, the average number of transmissions and LTAT only depend on the
outage probability, when the transmission rate and the maximum number of
transmissions are given. Meanwhile, the outage probability is equivalent to
the cumulative distribution function (CDF)
${{F_{I_{K}^{IR}}}\left(\mathcal{R}\right)}$ of the accumulated mutual
information $I_{K}^{IR}$, i.e.,
$P_{out}^{IR}\left(K\right)={{F_{I_{K}^{IR}}}\left(\mathcal{R}\right)}$. As
further proved in [22], the essential parameter to characterize the
performance of HARQ schemes is the CDF of the accumulated mutual information
in each round. It will then be particularly investigated in the next
subsection.

###  III-B Analysis of Outage Probability

By substituting (5) and (7) into (8), the outage probability becomes

$\displaystyle P_{out}^{IR}\left(K\right)$
$\displaystyle=\Pr\left({{{\log}_{2}}\left({\prod\limits_{l=1}^{K}{\left({1+{\gamma_{l}}}\right)}}\right)<{\cal
R}}\right)$ $\displaystyle=\int_{0}^{\cal
R}{{f_{I_{K}^{IR}}}\left(x\right)dx}$ (11)

where ${{f_{I_{K}^{IR}}}\left(x\right)}$ stands for the PDF of $I_{K}^{IR}$.
Since the considered channel is Rayleigh fading, it is easy to get that the
received SNR ${\gamma_{l}}$ follows exponential distribution with a PDF of

${f_{{\gamma_{l}}}}\left(x\right)=\frac{1}{{{2\sigma_{l}^{\prime}}{{}^{2}}}}\exp\left({-\frac{{{x}}}{{2{\sigma_{l}^{\prime}}{{}^{2}}}}}\right),x\in[0,+\infty)$
(12)

where
${\sigma_{l}^{\prime}}={\left({{{{P_{l}}}}/{{{\mathfrak{N}_{l}}}}}\right)^{\frac{1}{2}}}{\sigma_{l}}$.
Due to the time correlation in the channels, the SNRs ${\gamma_{l}}$ are
correlated. By making simple substitutions of variables on (II-B), the joint
distribution of
${\boldsymbol{\gamma}_{1:K}}=\left\\{{{\gamma_{1}},{\gamma_{2}},\cdots,{\gamma_{K}}}\right\\}$
can be readily derived as

$\displaystyle{{f_{\boldsymbol{\gamma}_{1:K}}}\left({{\gamma_{1}}={x_{1}},\cdots,{\gamma_{K}}={x_{K}}}\right)=\prod\limits_{k=1}^{K}{\frac{1}{{2{\sigma_{k}^{\prime}}{{}^{2}}\left({1-{\lambda_{k}}^{2}}\right)}}}}$
$\displaystyle\times{\int\nolimits_{t=0}^{\infty}{{{\rm{e}}^{-\left({1+\sum\limits_{k=1}^{K}{\frac{{{\lambda_{k}}^{2}}}{{1-{\lambda_{k}}^{2}}}}}\right)t}}\prod\limits_{k=1}^{K}{{e^{{-\frac{{{x_{k}}}}{{2{\sigma_{k}^{\prime}}{{}^{2}}\left({1-{\lambda_{k}}^{2}}\right)}}}}}}}}$
$\displaystyle\times{}_{0}{F_{1}}\left({;1;\frac{{{x_{k}}{\lambda_{k}}^{2}t}}{{2{\sigma_{k}^{\prime}}{{}^{2}}{{\left({1-{\lambda_{k}}^{2}}\right)}^{2}}}}}\right)dt.$
(13)

It is clear from (III-B) that the distribution of a product of time-correlated
shifted-exponential RVs is necessary to derive the outage probability. As
reported in the literature, Mellin transform can be exploited to derive the
distribution of the product of independent RVs [23, 24, 25, 26, 27, 28, 11].
Unfortunately, it is inapplicable for the case with correlated RVs due to the
involvement of multiple integral. In fact, the presence of time correlation
makes the derivation of the exact distribution of ${I_{K}^{IR}}$ intractable.
To proceed with the analysis, we resort to find a good approximation of the
distribution of ${I_{K}^{IR}}$ based on polynomial fitting technique which
will be introduced in the following.

As shown in [21], the accumulated mutual information in HARQ-IR systems over
independent fading channels can be well approximated as a Gamma RV by using
Laguerre series. Inspired by this result, the PDF of the accumulated mutual
information $I_{K}^{IR}$ over correlated Rayleigh fading channels can be
written as the product of a Gamma PDF $\varphi(x)$ and a correction term
$\psi\left(x\right)$ as 333Notice that when the channels are independent
fading, the correction term can be approximated as $\psi\left(x\right)\approx
1$ and the PDF is reduced as ${f_{I_{K}^{IR}}}(x)\approx\varphi(x)$ [21].

${f_{I_{K}^{IR}}}(x)=\varphi(x)\psi\left(x\right).$ (14)

In (14), the Gamma PDF $\varphi(x)$ serves as a basis function and is given by

$\varphi(x)=\frac{{{x^{\zeta-1}}{e^{-\frac{x}{\theta}}}}}{{{\theta^{\zeta}}\Gamma\left(\zeta\right)}},x\geq
0$ (15)

where $\Gamma(\cdot)$ denotes Gamma function, and the parameters $\zeta$ and
$\theta$ are determined by matching the first two moments of
$f_{I_{K}^{IR}}(x)$ with that of $\varphi(x)$, thus leading to [29]

$\zeta=\frac{{{{\mathcal{M}}^{2}}\left(1\right)}}{{{\mathcal{M}}\left(2\right)-{{\mathcal{M}}^{2}}\left(1\right)}}$
(16)
$\theta=\frac{{{\mathcal{M}}\left(2\right)-{{\mathcal{M}}^{2}}\left(1\right)}}{{{\mathcal{M}}\left(1\right)}}$
(17)

where ${\mathcal{M}}(i)$ denotes the $i$th moment with respect to
$f_{I_{K}^{IR}}(x)$. As shown in Appendix A, the $i$th moment can be derived
using Gaussian Quadrature as

$\displaystyle{\mathcal{M}}\left(i\right)\approx\frac{{i!}}{{1+\sum\limits_{k=1}^{K}{\frac{{{\lambda_{k}}^{2}}}{{1-{\lambda_{k}}^{2}}}}}}\sum\limits_{\sum\limits_{l=1}^{K}{i_{l}}=i,{i_{l}}\geq
0}{\frac{1}{{{i_{1}}!{i_{2}}!\cdots{i_{K}}!}}}$
$\displaystyle\times\sum\limits_{{q_{k}}\in\left[{1,{N_{Q}}}\right],k\in\left[{1,K}\right]}{\prod\limits_{k=1}^{K}{{\varrho_{{q_{k}}}}{{\log}_{2}}^{{i_{k}}}\left({1+{w_{k}}{\xi_{{q_{k}}}}}\right)}}$
$\displaystyle\times\Psi_{2}^{\left(K\right)}\left({1;1,1,\cdots,1;{\varpi_{1}}{\xi_{{q_{1}}}},{\varpi_{2}}{\xi_{{q_{2}}}},\cdots,{\varpi_{K}}{\xi_{{q_{K}}}}}\right).$
(18)

where
${\varpi_{k}}=\frac{{{\lambda_{k}}^{2}}}{{1-{\lambda_{k}}^{2}}}{\left({1+\sum\nolimits_{l=1}^{K}{\frac{{{\lambda_{l}}^{2}}}{{1-{\lambda_{l}}^{2}}}}}\right)^{-1}}$,
${w_{k}}=2{\sigma_{k}^{\prime}}^{2}\left({1-{\lambda_{k}}^{2}}\right)$,
$N_{Q}$ is the quadrature order, the weights $\varrho_{p}$ and the abscissas
${{\xi_{p}}}$ are tabulated in [30, Table 25.9], and
$\Psi_{2}^{\left(K\right)}(\cdot)$ is defined as the confluent form of
Lauricella hypergeometric function [31, Definition A.20], [32]. Specifically,
${\mathcal{M}}(0)=1$.

On the other hand, the correction term $\psi\left(x\right)$ in (14) is used to
compensate the difference between ${f_{I_{K}^{IR}}}(x)$ and the basis function
$\varphi(x)$. Apparently, it is very difficult to derive the exact expression
for $\psi\left(x\right)$. However, $\psi\left(x\right)$ can be generally
approximated as a polynomial $\hat{\psi}_{N}(x)\in\mathbb{P}_{N}$ with degree
$N$ by means of polynomial fitting technique [33], where the fitting error is
characterized as $e\left(x\right)=\psi\left(x\right)-\hat{\psi}_{N}(x)$. The
remaining problem here is then to find the optimal polynomial
$\hat{\psi}_{N}(x)$ which can minimize the mean square error (MSE), i.e.,
${\rm{E}}\\{{e\left(x\right)^{2}}\\}{\rm{=}}\int_{0}^{\infty}{\varphi\left(x\right){{\left({\psi\left(x\right)-\hat{\psi}_{N}(x)}\right)}^{2}}dx}$.

It is widely known that any polynomial can be written as a unique linear
combination of orthogonal polynomials. Denote the monic orthogonal polynomials
$\boldsymbol{\mathcal{P}}(x)=[\mathcal{P}_{0}(x),\mathcal{P}_{1}(x),\cdots,\mathcal{P}_{N}(x)]^{\rm{T}}$
with respect to a measure $d\mu(x)$ be the basis in the space of polynomials
of degree less than or equal to $N$, where
${{\mathcal{P}}_{n}}\left(x\right)=\sum\nolimits_{k=0}^{n}{{C_{n,k}}{x^{k}}}\in{{\mathbb{P}}_{n}}$
and $C_{n,n}=1$. As pointed out in [34, Theorem 1.27],
$\boldsymbol{\mathcal{P}}(x)$ is uniquely determined given the measure
$d\mu(x)$. Moreover, the monic orthogonal polynomials obey the orthogonality

$\left\langle{{\mathcal{P}_{n}}\left(x\right),{\mathcal{P}_{k}}\left(x\right)}\right\rangle={\delta_{n,k}}{{\cal
D}_{n}}=\left\\{{\begin{array}[]{*{20}{c}}{{{\cal D}_{n}},}&{n=k};\\\
{0,}&{else}.\end{array}},\right.$ (19)

where $\delta_{n,k}$ indicates Kronecker delta function,
$\left\langle{g\left(x\right),h\left(x\right)}\right\rangle$ denotes an inner
product defined on $2$-norm Lebesgue space $L^{2}(\mathbb{R},\mathcal{F},\mu)$
444Herein, $(\mathbb{R},\mathcal{F},\mu)$ is a measure space, where
$\mathcal{F}$ is $\delta$-algebra over $\mathbb{R}$. with the measure
$d\mu(x)$, that is,

$\left\langle{g\left(x\right),h\left(x\right)}\right\rangle=\int_{-\infty}^{+\infty}{g\left(x\right)h(x)d\mu(x)},$
(20)

and ${\cal D}_{n}$ is a non-zero parameter which will be specified later. With
the monic orthogonal polynomials $\boldsymbol{\mathcal{P}}(x)$,
$\hat{\psi}_{N}(x)$ can be written as

$\hat{\psi}_{N}(x){\rm{=}}\sum\limits_{i=0}^{N}{{\eta_{i}}{\mathcal{P}_{i}}\left(x\right)}={\boldsymbol{\eta}}^{\rm
T}\boldsymbol{\mathcal{P}}(x)$ (21)

where the column vector
${\boldsymbol{\eta}}=[\eta_{0},\eta_{1},\cdots,\eta_{N}]^{\rm{T}}$ can be
regarded as the corresponding coordinate vector of $\hat{\psi}_{N}(x)$ in the
space $\mathbb{P}_{N}$. Substituting (21) into (14), the PDF
${f_{I_{K}^{IR}}}(x)$ can be approximated as

${f_{I_{K}^{IR}}}(x){\approx}f_{{{\boldsymbol{\gamma}_{1:K}}},N}^{IR}(x)=\varphi(x)\hat{\psi}_{N}\left(x\right)=\varphi(x){\boldsymbol{\eta}}^{\rm
T}\boldsymbol{\mathcal{P}}(x).$ (22)

To guarantee the approximation accuracy, we need to find the optimal
polynomial correction term $\hat{\psi}_{N}(x)$ (i.e., the optimal
${\boldsymbol{\eta}}$ and $\boldsymbol{\mathcal{P}}(x)$) which can minimise
the MSE ${\rm{E}}\\{{e\left(x\right)^{2}}\\}$. It is noteworthy that the
resulting approximated PDF $f_{{{\boldsymbol{\gamma}_{1:K}}},N}^{IR}(x)$
should be normalized, such that

$\int_{0}^{\infty}f_{{{\boldsymbol{\gamma}_{1:K}}},N}^{IR}(x)dx=1.$ (23)

Considering this constraint and the fact that the monic orthogonal polynomials
$\boldsymbol{\mathcal{P}}(x)$ are uniquely determined by the measure
$d\mu(x)$, given the measure $d\mu(x)$ and $N$, the approximation problem can
be formulated as

$\displaystyle\underset{{\boldsymbol{\eta}}}{\text{min}}$ $\displaystyle{{\cal
S}_{mse}}({\boldsymbol{\eta}}|d\mu(x),N)$ (24)
$\displaystyle=\int\nolimits_{0}^{\infty}{\varphi\left(x\right){{\left({\sum\limits_{i=0}^{N}{{\eta_{i}}{{\cal
P}_{i}}\left(x\right)}-\psi\left(x\right)}\right)}^{2}}dx}$ s.t.
$\displaystyle\int_{0}^{\infty}f_{{{\boldsymbol{\gamma}_{1:K}}},N}^{IR}(x)dx=1$

This minimization problem can be solved by adopting the method of Lagrange
multiplier. The Lagrangian function corresponding to the minimization problem
can be written in matrix form as

$\displaystyle\Lambda({\boldsymbol{\eta}},\varsigma|d\mu(x),N)$
$\displaystyle=\int_{0}^{\infty}{\varphi\left(x\right){{\left({\sum\limits_{i=0}^{N}{{\eta_{i}}{{\cal
P}_{i}}\left(x\right)}-\psi\left(x\right)}\right)}^{2}}dx}$
$\displaystyle+\varsigma\left({\int_{0}^{\infty}{f_{{{\boldsymbol{\gamma}_{1:K}}},N}^{IR}(x)dx}-1}\right)$
$\displaystyle={{\boldsymbol{\eta}}^{\rm{T}}}{\bf{A\boldsymbol{\eta}}}-2{{\boldsymbol{\eta}}^{\rm{T}}}{\bf{b}}+\varsigma{{\boldsymbol{\eta}}^{\rm{T}}}{\bf{d}}+c-\varsigma$
(25)

where $\varsigma$ represents the Lagrange multiplier, $\bf A$, $\bf b$ and
$\bf d$ are given by (30)-(34) as shown on the top of next page, respectively,

$\displaystyle{{\bf{A}}=\left[{\begin{array}[]{*{20}{c}}{\int_{-\infty}^{\infty}{\varphi\left(x\right){{\mathcal{P}}_{0}}^{2}\left(x\right)dx}}&{\int_{0}^{\infty}{\varphi\left(x\right){{\mathcal{P}}_{0}}\left(x\right){{\mathcal{P}}_{1}}\left(x\right)dx}}&\cdots&{\int_{0}^{\infty}{\varphi\left(x\right){{\mathcal{P}}_{0}}\left(x\right){{\mathcal{P}}_{N}}\left(x\right)dx}}\\\
{\int_{0}^{\infty}{\varphi\left(x\right){{\mathcal{P}}_{1}}\left(x\right){{\mathcal{P}}_{0}}\left(x\right)dx}}&{\int_{0}^{\infty}{\varphi\left(x\right){{\mathcal{P}}_{1}}^{2}\left(x\right)dx}}&\cdots&{\int_{0}^{\infty}{\varphi\left(x\right){{\mathcal{P}}_{1}}\left(x\right){{\mathcal{P}}_{N}}\left(x\right)dx}}\\\
\vdots&\vdots&\ddots&\vdots\\\
{\int_{0}^{\infty}{\varphi\left(x\right){{\mathcal{P}}_{N}}\left(x\right){{\mathcal{P}}_{0}}\left(x\right)dx}}&{\int_{0}^{\infty}{\varphi\left(x\right){{\mathcal{P}}_{N}}\left(x\right){{\mathcal{P}}_{1}}\left(x\right)dx}}&\cdots&{\int_{0}^{\infty}{\varphi\left(x\right){{\mathcal{P}}_{N}}^{2}\left(x\right)dx}}\end{array}}\right]},$
(30)

$\displaystyle{\bf{b}}={\left[{\begin{array}[]{*{20}{c}}{\int_{0}^{\infty}{{f_{I_{K}^{IR}}}(x){{\cal
P}_{0}}\left(x\right)dx}}&{\int_{0}^{\infty}{{f_{I_{K}^{IR}}}(x){{\cal
P}_{1}}\left(x\right)dx}}&\cdots&{\int_{0}^{\infty}{{f_{I_{K}^{IR}}}(x){{\cal
P}_{N}}\left(x\right)dx}}\end{array}}\right]^{\rm{T}}},$ (32)

$\displaystyle{\bf{d}}={\left[{\begin{array}[]{*{20}{c}}{\int_{0}^{\infty}{\varphi(x){\mathcal{P}_{0}}\left(x\right)dx}}&{\int_{0}^{\infty}{\varphi(x){\mathcal{P}_{1}}\left(x\right)dx}}&\cdots&{\int_{0}^{\infty}{\varphi(x){\mathcal{P}_{N}}\left(x\right)dx}}\end{array}}\right]^{\rm{T}}},$
(34)

and $c$ is

$c=\int_{0}^{\infty}{\varphi\left(x\right){\psi^{2}}\left(x\right)dx}.$ (35)

According to the Karush-Kuhn-Tucker (KKT) conditions, the optimal solutions
${\boldsymbol{\eta}}$ and $\varsigma$ should satisfy the following conditions

$\left\\{{\begin{array}[]{*{20}{l}}{\frac{{\partial\Lambda\left({{\boldsymbol{\eta}},\varsigma|d\mu(x),N}\right)}}{{\partial{{\boldsymbol{\eta}}}}}{\rm{=2}}{\bf{A{\boldsymbol{\eta}}}}-2{\bf{b}}+\varsigma{\bf{d}}=0}\\\
{\frac{{\partial\Lambda\left({{\boldsymbol{\eta}},\varsigma|d\mu(x),N}\right)}}{{\partial{\varsigma}}}{={\boldsymbol{\eta}}^{\rm{T}}}{{\bf{d}}-1}=0}\end{array}}\right..$
(36)

As proved in Appendix B, the matrix $\bf A$ is invertible. Then the solution
to (36) is unique and follows as

$\left\\{{\begin{array}[]{*{20}{l}}{{\boldsymbol{\eta}}={{\bf{A}}^{-1}}\left({{\bf{b}}-\frac{\varsigma}{2}{\bf{d}}}\right)}\\\
{\varsigma=\frac{{2\left({{{\bf{b}}^{\rm{T}}}{{\bf{A}}^{-1}}{\bf{d}}-1}\right)}}{{{{\bf{d}}^{\rm{T}}}{{\bf{A}}^{-1}}{\bf{d}}}}}\end{array}}\right.$
(37)

Clearly from (30), (32), (34) and (37), the monic orthogonal polynomials
${{\mathcal{P}}_{n}}\left(x\right)$ need to be determined before the
calculation of the coefficient ${\boldsymbol{\eta}}$. In other words, the
measure $d\mu(x)$ should be determined first. In what follows, the selection
of the measure and the orthogonal polynomials and the calculation of the
coefficient ${\boldsymbol{\eta}}$ are then discussed in detail.

####  III-B1 Selection of $d\mu(x)$ and $\boldsymbol{\mathcal{P}}(x)$

After analyzing the MSE of the fitting error
${\rm{E}}\\{{e\left(x\right)^{2}}\\}$, we have the following property.

######  Property 1.

For any measure $d\mu(x)$, the same minimal MSE can be attained. Moreover, the
optimal polynomial $\hat{\psi}_{N}(x)$ is unique irrespective of the choice of
$d\mu(x)$.

###### Proof:

Please see Appendix C. ∎

Although the choice of the measure $d\mu(x)$ would not affect the solution of
the optimal polynomial $\hat{\psi}_{N}(x)$ as shown in Property 1, it does
affect the computational complexity in deriving the optimal polynomial. More
specifically, the major computation comes from the inverse operation of the
matrix $\bf A$ as seen from (37). It is clear from the definition of $\bf A$
(30) that the complexity of the inverse operation depends on the choice of the
orthogonal polynomials ${{\mathcal{P}}_{n}}\left(x\right)$ and thus depends on
the measure $d\mu(x)$. It is generally known that the complexity in computing
the inverse of a $N\times N$ matrix could be significantly reduced from
$\mathcal{O}(N^{3})$ to $\mathcal{O}(N)$ when the matrix is diagonal.
Therefore, to reduce the computational complexity, the measure $d\mu(x)$ is
suggested to be chosen such that the matrix $\bf A$ is diagonal. Clearly from
the structure of the matrix $\bf A$ (30), when the measure is chosen as
$d\mu(x)={\varphi\left(x\right)}dx$, with the orthogonality of the polynomials
in (19), the matrix $\bf A$ will reduce to a diagonal matrix as

${\bf{A}}=\rm{diag}\left({{\mathcal{D}_{0}},{\mathcal{D}_{1}},\cdots,{{\mathcal{D}}_{N}}}\right).$
(38)

Now with the measure $d\mu(x)={\varphi\left(x\right)}dx$, the monic orthogonal
polynomials
${{\mathcal{P}}_{n}}\left(x\right)=\sum\nolimits_{k=0}^{n}{{C_{n,k}}{x^{k}}}$
with $C_{n,n}=1$ can be uniquely determined by the method introduced in [34].
Specifically, the polynomial coefficients $C_{n,k}$ can be determined as
follows. Following the three-term recurrence relation (TTRR) in [34, Theorem
1.27], the monic orthogonal polynomials should satisfy

$\displaystyle{{{\mathcal{P}}_{n+1}}\left(x\right)=\left({x-{\alpha_{n}}}\right){{\mathcal{P}}_{n}}\left(x\right)-{\beta_{n}}{{\mathcal{P}}_{n-1}}\left(x\right),\,n=0,1,2,\cdots,}$
$\displaystyle{{{\mathcal{P}}_{-1}}\left(x\right)=0,\;\;\;{\kern
1.0pt}{{\mathcal{P}}_{0}}\left(x\right)=1,}$ (39)

where

$\displaystyle{\alpha_{n}}$
$\displaystyle=\frac{\left<x{\mathcal{P}}_{n}(x),{\mathcal{P}}_{n}(x)\right>}{\left<{\mathcal{P}}_{n}(x),{\mathcal{P}}_{n}(x)\right>}$
$\displaystyle=\frac{{\sum\limits_{i=0}^{n}{\sum\limits_{j=0}^{n}{{C_{n,i}}{C_{n,j}}{\nu_{i+j+1}}}}}}{{\sum\limits_{i=0}^{n}{\sum\limits_{j=0}^{n}{{C_{n,i}}{C_{n,j}}{\nu_{i+j}}}}}},\,n=0,1,2,\cdots,$
(40) $\displaystyle{\beta_{n}}$
$\displaystyle=\frac{\left<{\mathcal{P}}_{n}(x),{\mathcal{P}}_{n}(x)\right>}{\left<{\mathcal{P}}_{n-1}(x),{\mathcal{P}}_{n-1}(x)\right>}$
$\displaystyle=\frac{{\sum\limits_{i=0}^{n}{\sum\limits_{j=0}^{n}{{C_{n,i}}{C_{n,j}}{\nu_{i+j}}}}}}{{\sum\limits_{i=0}^{n-1}{\sum\limits_{j=0}^{n-1}{{C_{n-1,i}}{C_{n-1,j}}{\nu_{i+j}}}}}},\,n=1,2,\cdots.$
(41)

and $\nu_{n}$ denotes the $n$th moment with respect to the cumulative
distribution function (CDF) $\mu(x)$, that is,

$\displaystyle{\nu_{n}}$
$\displaystyle=\int_{0}^{\infty}{{x^{n}}d\mu\left(x\right)}=\int_{0}^{\infty}{{x^{n}}\varphi\left(x\right)dx}$
$\displaystyle=\int_{0}^{\infty}{\frac{{{x^{n+\zeta-1}}{e^{-\frac{x}{\theta}}}}}{{{\theta^{\zeta}}\Gamma\left(\zeta\right)}}dx}=\frac{{{\theta^{n}}\Gamma\left({n+\zeta}\right)}}{{\Gamma\left(\zeta\right)}}$
(42)

With
${{\mathcal{P}}_{n}}\left(x\right)=\sum\nolimits_{k=0}^{n}{{C_{n,k}}{x^{k}}}$
and the relation in (III-B1), the polynomial coefficients $C_{n,k}$ can be
obtained recursively as

${C_{n+1,k}}=\left\\{{\begin{array}[]{*{20}{c}}{{C_{n,k-1}}-{\alpha_{n}}{C_{n,k}}-{\beta_{n}}{C_{n-1,k}},}&{0\leq
k\leq n+1;}\\\ {0,}&{else.}\end{array}}\right.$ (43)

with ${C_{0,0}}=1$. Then the parameter ${\mathcal{D}}_{n}$ in (19) can be
determined as

${{\mathcal{D}}_{n}}=\sum\limits_{i=0}^{n}{\sum\limits_{j=0}^{n}{{C_{n,i}}{C_{n,j}}{\nu_{i+j}}}}.$
(44)

Specifically, ${{\mathcal{D}}_{0}}=1$.

####  III-B2 Calculation of $\boldsymbol{\eta}$

To compute $\boldsymbol{\eta}$ in (37), the vectors $\bf b$, $\bf d$ and
$\varsigma$ should be determined first. By the definition of $\mathcal{M}(i)$
and integrating out $x$ for (32), the vector $\bf{b}$ can be consequently
written as

${\bf{b}}=\left[{\begin{array}[]{*{20}{c}}{\sum\limits_{k=0}^{0}{{C_{0,k}}}{\cal
M}\left(k\right)}&{\sum\limits_{k=0}^{1}{{C_{0,k}}}{\cal
M}\left(k\right)}\end{array}}\right.\\\
{\left.{\begin{array}[]{*{20}{c}}\cdots&{\sum\limits_{k=0}^{N}{{C_{N,k}}}{\cal
M}\left(k\right)}\end{array}}\right]^{\rm{T}}}.$ (45)

Considering $\mathcal{P}_{0}(x)=1$ and the orthogonality among the polynomials
$\mathcal{P}_{n}(x)$, the vector d given in (34) is directly reduced as

${\bf{d}}={\left[{1,\quad\overbrace{0,\quad\cdots\quad,0}^{N-terms}}\right]^{\rm
T}}.$ (46)

On the other hand, by putting (38), (45) and (46) into (37), the Lagrange
multiplier $\varsigma$ is calculated as

$\varsigma=\frac{{{C_{0,0}}{\mathcal{M}}\left(0\right){D_{0}}-1}}{{{D_{0}}}}=0$
(47)

by recalling that $C_{0,0}=1$, ${\mathcal{M}}(0)=1$ and $\mathcal{D}_{0}=1$.
Accordingly, the coefficients $\boldsymbol{\eta}$ can be finally computed as

${\boldsymbol{\eta}}={{\bf{A}}^{-1}}{\bf{b}}\\\
={\left[{\begin{array}[]{*{20}{c}}{\frac{{\sum\limits_{k=0}^{0}{{C_{0,k}}}{\mathcal{M}}\left(k\right)}}{{{{\mathcal{D}}_{0}}}}}&{\frac{{\sum\limits_{k=0}^{1}{{C_{1,k}}}{\mathcal{M}}\left(k\right)}}{{{{\mathcal{D}}_{1}}}}}&\cdots&{\frac{{\sum\limits_{k=0}^{N}{{C_{N,k}}}{\mathcal{M}}\left(k\right)}}{{{{\mathcal{D}}_{N}}}}}\end{array}}\right]^{\rm
T}}.$ (48)

###  III-C Discussions

With the approximated PDF of the accumulated mutual information in (22), the
outage probability which is equivalent to the CDF of the accumulated mutual
information can be easily obtained. From the expression of the outage
probability, some interesting insights could be found in the following.
Moreover, as shown in (24), the MSE of the fitting error also depends on the
degree of the polynomials $N$. The selection of the degree will also be
briefly discussed here.

####  III-C1 Insights of Outage Probability

After determining $\boldsymbol{\eta}$, the approximated PDF
$f_{{{\boldsymbol{\gamma}_{1:K}}},N}^{IR}(x)$ is expressed with the expansion
of ${{{\mathcal{P}}_{n}}\left(x\right)}$ as

$f_{{{\boldsymbol{\gamma}_{1:K}}},N}^{IR}(x)=\varphi(x)\sum\limits_{n=0}^{N}{{\eta_{n}}\sum\limits_{i=0}^{n}{{C_{n,i}}{x^{i}}}}.$
(49)

Henceforth, the approximated CDF for $I_{K}^{IR}$ can be obtained as

$\displaystyle{F_{I_{K}^{IR}}}(x)$ $\displaystyle\approx
F_{{{\boldsymbol{\gamma}_{1:K}}},N}^{IR}(x)=\int_{0}^{x}{f_{{{\boldsymbol{\gamma}_{1:K}}},N}^{IR}(t)dt}$
$\displaystyle=\sum\limits_{n=0}^{N}{{\eta_{n}}\sum\limits_{i=0}^{n}{{C_{n,i}}\int_{0}^{x}{{t^{i}}\varphi(t)dt}}}.$
(50)

To facilitate the analysis, a family of functions $W_{i}(x)$ is defined as

$\displaystyle{W_{i}}\left(x\right)$
$\displaystyle=\frac{1}{{{\nu_{i}}}}\int_{0}^{x}{{t^{i}}\varphi\left(t\right)dt}=\frac{1}{{{\theta^{i+\zeta}}\Gamma\left({i+\zeta}\right)}}\int_{0}^{x}{{t^{i+\zeta-1}}{e^{-\frac{t}{\theta}}}dt}$
$\displaystyle=\frac{{\gamma\left({i+\zeta,\frac{x}{\theta}}\right)}}{{\Gamma\left({i+\zeta}\right)}}$
(51)

where $\gamma(\cdot)$ represents the lower incomplete Gamma function. It is
clear that $W_{i}(x)$ is the CDF of a Gamma RV with parameters
$(i+\zeta,\theta)$. Consequently, by exchanging the order of summations, the
CDF in (III-C1) can be rewritten as

$\displaystyle F_{{{\boldsymbol{\gamma}_{1:K}}},N}^{IR}(x)$
$\displaystyle=\sum\limits_{n=0}^{N}{{\eta_{n}}\sum\limits_{i=0}^{n}{{C_{n,i}}{\nu_{i}}{W_{i}}\left(x\right)}}$
$\displaystyle=\sum\limits_{i=0}^{N}{{W_{i}}\left(x\right){\nu_{i}}\sum\limits_{n=i}^{N}{{\eta_{n}}}{C_{n,i}}}$
$\displaystyle=\sum\limits_{i=0}^{N}{{\kappa_{i}}}{W_{i}}\left(x\right)$ (52)

where

${\kappa_{i}}={\nu_{i}}\sum\limits_{n=i}^{N}{{\eta_{n}}}{C_{n,i}}.$ (53)

Meanwhile, the parameters $\\{\kappa_{i}\\}$ satisfy

$\mathop{\lim}\limits_{x\to\infty}F_{{{\boldsymbol{\gamma}_{1:K}}},N}^{IR}(x)=\mathop{\lim}\limits_{x\to\infty}\sum\limits_{i=0}^{N}{{\kappa_{i}}}{W_{i}}\left(x\right)=1\\\
\Rightarrow\sum\limits_{i=0}^{N}{{\kappa_{i}}}=1.$ (54)

Therefore, under time-correlated Rayleigh fading channels, the CDF of the
accumulated mutual information $F_{I_{K}^{IR}}(x)$ can be written as the
weighted sum of the CDFs of Gamma RVs, such that

${F_{I_{K}^{IR}}}(x)\approx\sum\limits_{i=0}^{N}{{\kappa_{i}}}{W_{i}}\left(x\right).$
(55)

It means that the outage probability
$P_{out}^{IR}\left(K\right)=F_{I_{K}^{IR}}({\cal R})$ can be written as the
weighted sum of a number of outage probabilities, each associated with one
Gamma RV555The source code of our approximation is available at
https://sourceforge.net/projects/matlabgammaapproximation/files/Gamma%20Approximation/..
Specifically, if $N=0$, the approximation reduces to the ordinary Gamma
approximation which is valid for the case with independent fading channels
[21]. The approximation in (55) can ease the analysis of the system behaviors
with respect to various system parameters and facilitate system design to
achieve various objectives, e.g., the optimal rate design to maximize the long
term average throughput. This will be further illustrated in Section V.

####  III-C2 Choice of $N$

As shown in Property 1, the same minimal MSE ${{\cal
S}_{mse}}({\boldsymbol{\eta}}|d\mu(x),N)$ can be attained whatever the measure
$d\mu(x)$ is. Hereby, we define
${{\mathcal{S}}_{min\\_mse}}({\boldsymbol{\eta}}|N)$ as the minimal MSE given
$N$. By substituting (37) into the objective function of MSE in (24), it
yields the minimal MSE as

$\displaystyle{{\mathcal{S}}_{min\\_mse}}({\boldsymbol{\eta}}|N)$
$\displaystyle={{\boldsymbol{\eta}}^{\rm{T}}}\left({{\bf{A}\boldsymbol{\eta}}-2{\bf{b}}}\right)+c$
$\displaystyle=c-{{\boldsymbol{\eta}}^{\rm{T}}}\left({{\bf{b}}+\frac{\varsigma}{2}{\bf{d}}}\right)$
$\displaystyle=c-\left({{{\bf{b}}^{\rm
T}}-\frac{\varsigma}{2}{{\bf{d}}^{\rm{T}}}}\right){{\bf{A}}^{-1}}\left({{\bf{b}}+\frac{\varsigma}{2}{\bf{d}}}\right).$
(56)

Putting (38), (45), (46) and (47) into (III-C2), the minimal MSE is rewritten
as

${{\mathcal{S}}_{min\\_mse}}({\boldsymbol{\eta}}|N)=\int_{0}^{\infty}{\varphi\left(x\right){\psi^{2}}\left(x\right)dx}\\\
-\sum\limits_{n=0}^{N}{{{\mathcal{D}}_{n}}^{-1}{{\left({\sum\limits_{k=0}^{n}{{C_{n,k}}}{\mathcal{M}}\left(k\right)}\right)}^{2}}}\geq
0.$ (57)

Clearly, the minimum MSE decreases as the degree $N$ increases, i.e.
${{\mathcal{S}}_{min\\_mse}}({\boldsymbol{\eta}}|N)\geq{{\mathcal{S}}_{min\\_mse}}({\boldsymbol{\eta}}|N+1)$,
since
${\mathcal{D}_{n}}=\left\langle{{{\mathcal{P}}_{n}}\left(x\right),{{\mathcal{P}}_{n}}\left(x\right)}\right\rangle>0$.
It implies that the accuracy of the PDF approximation is limited by the degree
$N$ and would be improved as the degree $N$ increases. This result can be
further demonstrated by Fig. 2 where the approximated CDFs of $I_{K}^{IR}$
using different degrees are compared with the true CDFs obtained from Monte-
Carlo simulations, by taking a system with the following setting as an
example: the channel correlation coefficient $\rho_{k,l}=0.5$ and the mean of
the SNR ${\rm{E}}(\gamma_{l})=2{\sigma_{l}^{\prime}}^{2}=5$, where $1\leq
k\neq l\leq K$.

Figure 2: Comparison between the approximated CDFs ${F_{I_{K}^{IR}}}(x)$ with
different $N$ and the true CDFs obtained from Monte-Carlo simulations.

However, as shown in Fig. 2, the improvement of the approximation accuracy
becomes minor when the degree becomes relatively large. Moreover, the increase
of the degree $N$ would also cause the increase of the computational
complexity. It is thus necessary to properly choose the degree to balance the
approximation accuracy and computational complexity. To quantify the
contribution of increasing the degree in terms of the approximation accuracy,
a metric of MSE reduction is defined as

$\displaystyle{\Delta_{N}}$ $\displaystyle\triangleq{{\cal
S}_{min\\_mse}}({\boldsymbol{\eta}}|0)-{{\cal
S}_{min\\_mse}}({\boldsymbol{\eta}}|N)$
$\displaystyle=\sum\limits_{n=1}^{N}{{{\mathcal{D}}_{n}}^{-1}{{\left({\sum\limits_{k=0}^{n}{{C_{n,k}}}{\mathcal{M}}\left(k\right)}\right)}^{2}}}.$
(58)

For a good balance between approximation accuracy and computational
complexity, the degree is suggested to be chosen as

${N}=\min\left\\{{\min\left\\{{n\left|{{r_{n}}\leq\epsilon}\right.}\right\\},\hat{N}}\right\\}$
(59)

where $\hat{N}$ is a pre-determined upper bound of the degree to limit the
computational complexity, and $\epsilon$ denotes the tolerance for normalized
MSE reduction defined as
${r_{n}}=\frac{{{\Delta_{n}-{\Delta_{n-1}}}}}{{{\Delta_{n}}}}$. Clearly,
${{r_{N}}\leq\epsilon}$ holds. It roughly indicates that no significant
improvement on MSE can be expected if the degree is larger than $N$ in (59).
Therefore the degree of $N$ in (59) is sufficient for a good approximation.
Notice that $\Delta_{1}={\Delta_{2}}=0$ by using the orthogonality of the
polynomials. Thus the upper bound $\hat{N}$ and the degree $N$ should be set
greater than $2$ without doubt.

##  IV Diversity Order

Basically, HARQ-IR schemes exploit not only coding gain but also time
diversity to improve the transmission reliability. To better understand the
behavior of HARQ-IR schemes, diversity order which is another important
performance metric is also analyzed in this paper. To facilitate the analysis,
the transmit SNR in each HARQ round is set equal, i.e.
$P_{1}/\mathfrak{N}_{1}=P_{2}/\mathfrak{N}_{2}=\cdots=P_{M}/\mathfrak{N}_{M}=\gamma_{T}$.
The diversity order $d$ for the HARQ scheme is defined as [21, 35]

$d=-\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\log\left({P_{e}}\right)}}{{\log\left({{\gamma_{T}}}\right)}},$
(60)

where $P_{e}$ denotes the error probability. As shown in [21, 35], the error
probability can be well approximated as the outage probability
${P_{out}^{IR}\left(M\right)}$ when a capacity achieving code is applied.
Therefore the diversity order $d$ can be well approximated as

$d=-\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\log\left({P_{out}^{IR}\left(M\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}.$
(61)

Namely, diversity order quantifies the slope of
$-\log\left({P_{out}^{IR}\left(M\right)}\right)$ with respect to
$\log\left({{\gamma_{T}}}\right)$ when ${{\gamma_{T}}\to\infty}$.

Owing to the following inequalities

${\log_{2}}\left({1+\sum\limits_{l=1}^{M}{\gamma_{l}}}\right)\leq{\log_{2}}\left({\prod\limits_{l=1}^{M}{\left({1+\gamma_{l}}\right)}}\right)\\\
\leq M{\log_{2}}\left({1+{M^{-1}}\sum\limits_{l=1}^{M}{\gamma_{l}}}\right),$
(62)

the outage probability can be bounded as

$\Pr\left({M{{\log}_{2}}\left({1+{M^{-1}}\sum\limits_{l=1}^{M}{\gamma_{l}}}\right)<{\mathcal{R}}}\right)\leq
P_{out}^{IR}\left(M\right)\\\
\leq\Pr\left({{{\log}_{2}}\left({1+\sum\limits_{l=1}^{M}{\gamma_{l}}}\right)<{\mathcal{R}}}\right),$
(63)

where the right inequality in (62) holds by using the Jensen’s inequality.
Defining $Y=\sum\nolimits_{l=1}^{M}{\gamma_{l}}$, the bounds in (63) can be
rewritten as

$\Pr\left({Y<M({{2^{{M^{-1}}{\cal R}}}-1})}\right)\leq
P_{out}^{IR}\left(M\right)\leq\Pr\left({Y<{2^{\cal R}}-1}\right).$ (64)

Clearly, $Y$ represents a sum of correlated exponential RVs
$\\{\gamma_{l}\\}_{l=1}^{M}$. The CDF of a sum of correlated exponential RVs
has been derived in [36] and the result is summarized as the following
theorem.

######  Theorem 1.

[36] Given the joint PDF regarding to exponential RVs
$\\{\gamma_{l}\\}_{l=1}^{M}$ in (III-B), the CDF of
$Y=\sum\nolimits_{l=1}^{M}{\gamma_{l}}$ can be obtained as

${F_{Y}}\left(y\right)=\Pr\left({Y<y}\right)=\frac{{{y^{M}}}}{{\det\left({\bf{B}}\right)\Gamma\left({M+1}\right)}}\times\\\
\Phi_{2}^{\left(M\right)}\left({1,\cdots,1;M+1;-{\delta_{1}}^{-1}y,\cdots,-{\delta_{M}}^{-1}y}\right)$
(65)

where the notation $\rm det(\cdot)$ represents the determinant operation,
$\Phi_{2}^{\left(M\right)}\left(\cdot\right)$ denotes the confluent Lauricella
function [31, Def. A.19], $\\{\delta_{k}\\}_{k=1}^{M}$ are defined as the
eigenvalues of the matrix $\bf B=FE$, where $\bf F$ is an $M\times M$ diagonal
matrix with diagonal entries as $\\{2{\sigma_{k}^{\prime}}^{2}\\}_{k=1}^{M}$,
and $\bf E$ is an $M\times M$ positive definite matrix given by

${\bf{E}}=\left[{\begin{array}[]{*{20}{c}}1&{\sqrt{{\rho_{1,2}}}}&\cdots&{\sqrt{{\rho_{1,M}}}}\\\
{\sqrt{{\rho_{2,1}}}}&1&\cdots&{\sqrt{{\rho_{2,M}}}}\\\
\vdots&\vdots&\ddots&\vdots\\\
{\sqrt{{\rho_{M,1}}}}&{\sqrt{{\rho_{M,2}}}}&\cdots&1\end{array}}\right],\,0\leq\rho_{k,l}<1.$
(66)

Substituting (64) into (61), the bounds of the diversity order can be found as

$-\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\log\left({{F_{Y}}\left({{2^{\mathcal{R}}}-1}\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}\leq
d\\\
\leq-\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\log\left({{F_{Y}}\left({M\left({{2^{{M^{-1}}{\cal
R}}}-1}\right)}\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}.$ (67)

By using (65), the left inequality in (67) can be rewritten as

$\displaystyle
d\geq\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\left[{\frac{{\log\left({\det\left({\bf{B}}\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}}\right.-$
$\displaystyle{\left.{\frac{{\log\left({\Phi_{2}^{\left(M\right)}\left({1,\cdots,1;M+1;-{\delta_{1}}^{-1}y,\cdots,-{\delta_{M}}^{-1}y}\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}}\right]_{y={2^{R}}-1}}.$
(68)

The first term on the right hand side of (IV) can be simplified by using the
property of determinants as

$\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\log\left({\det\left({\bf{B}}\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}=\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\log\left({\det\left({\bf{F}}\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}+\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\log\left({\det\left({\bf{E}}\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}.$
(69)

Recalling that
${\bf{F}}={\gamma_{T}}{\rm{diag}}\left({2{\sigma_{1}}^{2},\cdots,2{\sigma_{M}}^{2}}\right)$
and ${\bf{E}}$ is irrelevant with ${\gamma_{T}}$, (69) can be reduced to

$\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\log\left({\det\left({\bf{B}}\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}=M.$
(70)

To derive the second term on the right hand side of the inequality (IV), we
define ${\bf
G}={\rm{diag}}\left({2{\sigma_{1}}^{2},\cdots,2{\sigma_{M}}^{2}}\right)\bf{E}$
for notational simplicity, such that ${\bf{B}}=\gamma_{T}\bf{G}$. By using the
series representation of the confluent Lauricella function [37], the confluent
Lauricella function in the second term on the right hand side of the
inequality (IV) can be rewritten as

$\Phi_{2}^{\left(M\right)}\left({1,\cdots,1;M+1;-{\delta_{1}}^{-1}y,\cdots,-{\delta_{M}}^{-1}y}\right)\\\
=1+\sum_{{m_{1}+\cdots+m_{M}}>0}{{\frac{{\prod\limits_{k=1}^{M}{{{\left({-{\delta_{k}}^{-1}y}\right)}^{{m_{k}}}}}}}{{{{\left(M+1\right)}_{{m_{1}}+\cdots+{m_{M}}}}}}}},$
(71)

where $(\cdot)_{n}$ denotes Pochhammer symbol. Since

$\displaystyle\left|\sum_{{m_{1}+\cdots+m_{M}}>0}{{\frac{{\prod\limits_{k=1}^{M}{{{\left({-{\delta_{k}}^{-1}y}\right)}^{{m_{k}}}}}}}{{{{\left(M+1\right)}_{{m_{1}}+\cdots+{m_{M}}}}}}}}\right|$
$\displaystyle\leq\sum_{{m_{1}+\cdots+m_{M}}>0}{{\frac{{{m_{1}}!\cdots{m_{M}}!}}{{{{\left({M+1}\right)}_{{m_{1}}+\cdots+{m_{M}}}}}}\frac{{\prod\limits_{k=1}^{M}{{{\left|{-{\delta_{k}}^{-1}y}\right|}^{{m_{k}}}}}}}{{{m_{1}}!\cdots{m_{M}}!}}}}$
$\displaystyle\underset{(a)}{\leq}\sum_{{m_{1}+\cdots+m_{M}}>0}{\frac{{\prod\limits_{k=1}^{M}{{{\left|{-{\delta_{k}}^{-1}y}\right|}^{{m_{k}}}}}}}{{{m_{1}}!\cdots{m_{M}}!}}}$
$\displaystyle=\sum\limits_{L=1}^{\infty}{\frac{1}{{L!}}\sum_{\sum\limits_{k=1}^{M}{{m_{k}}=L}}{\frac{{L!}}{{{m_{1}}!\cdots{m_{M}}!}}\prod\limits_{k=1}^{M}{{{\left|{-{\delta_{k}}^{-1}y}\right|}^{{m_{k}}}}}}}$
$\displaystyle=\sum\limits_{L=1}^{\infty}{\frac{1}{{L!}}{{\left({y\sum\limits_{k=1}^{M}{\left|{\delta_{k}}\right|^{-1}}}\right)}^{L}}}$
$\displaystyle={e^{y\sum\limits_{k=1}^{M}{\left|{\delta_{k}}\right|^{-1}}}}-1\underset{(b)}{=}{e^{y\sum\limits_{k=1}^{M}{\left|{\gamma_{T}\beta_{k}}\right|^{-1}}}}-1$
$\displaystyle={e^{y\gamma_{T}^{-1}\sum\limits_{k=1}^{M}{\left|{\beta_{k}}\right|^{-1}}}}-1$
(72)

where $(a)$ follows from
${m_{1}}!\cdots{m_{M}}!\leq{\left({M+1}\right)_{{m_{1}}+\cdots+{m_{M}}}}$,
$\beta_{k}$ denotes the eigenvalues of $\bf G$ and $(b)$ comes from
${\bf{B}}=\gamma_{T}\bf{G}$, together with
$\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}{e^{y\gamma_{T}^{-1}\sum\nolimits_{k=1}^{M}{\left|{\beta_{k}}\right|^{-1}}}}=1$,
we have

$\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\sum_{{m_{1}+m_{2}+\cdots+m_{M}}>0}{{\frac{{\prod\limits_{k=1}^{M}{{{\left({-{\delta_{k}}^{-1}y}\right)}^{{m_{k}}}}}}}{{{{\left(M+1\right)}_{{m_{1}}+\cdots+{m_{M}}}}}}}}=0.$
(73)

It follows the limit of the confluent Lauricella function as

$\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\Phi_{2}^{\left(M\right)}\left({1,\cdots,1;M+1;-{\delta_{1}}^{-1}y,\cdots,-{\delta_{M}}^{-1}y}\right)=1.$
(74)

and then the second term on the right hand side of the inequality (IV) reduces
as

$\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\log\left({\Phi_{2}^{\left(M\right)}\left({1,\cdots,1;M+1;-{\delta_{1}}^{-1}y,\cdots,-{\delta_{M}}^{-1}y}\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}\\\
=0.$ (75)

Substituting (70) and (75) into (IV) leads to $d\geq M$. Following the same
approach, the second inequality of (67) can also be derived as $d\leq M$. As a
result, under time correlated fading channels with $0\leq\rho_{k,l}<1$, the
diversity order $d$ is equal to the number of transmissions $M$, i.e., $d=M$.
Equivalently,

${P_{out}^{IR}\left(M\right)}\propto\frac{1}{{\gamma_{T}}^{M}}$ (76)

for large SNR. More precisely, the outage probability
${P_{out}^{IR}\left(M\right)}$ can be expressed as [1, 3.158]

$P_{out}^{IR}\left(M\right)=c\left({{\gamma_{T}},\sigma_{k},{\rho_{k,l}}},\mathcal{R},M\right){{{\gamma_{T}}^{-M}}},$
(77)

where the coefficient $c\left({{\gamma_{T}},{\sigma_{k}},{\rho_{k,l}},{\cal
R},M}\right)$ satisfies

$\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\log\left({c\left({{\gamma_{T}},\sigma_{k},{\rho_{k,l}},{\cal
R},M}\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}=0.$ (78)

From (77) and (78), we can find that when $\gamma_{T}$ is large, the slope of
$-\log\left({P_{out}^{IR}\left(M\right)}\right)$ with respect to
$\log\left({{\gamma_{T}}}\right)$ tends to be constant as $M$. In other words,
when $\rho_{k,l}<1$, the time correlation ${\rho_{k,l}}$ would not affect the
diversity order and the diversity order is constant as $M$. However the time
correlation would influence the coefficient
$c\left({{\gamma_{T}},{\sigma_{k}},{\rho_{k,l}},{\cal R},M}\right)$ and thus
affect the outage performance. This effect will be further investigated in
Section V.

The above analysis indicates that full diversity can be achieved by HARQ-IR
schemes even under time-correlated fading channels with $0\leq\rho_{k,l}<1$,
which further justifies the benefit of HARQ-IR.

Remark 1: The result of the diversity order is not applicable to the case with
fully correlated fading channels. Under fully correlated Rayleigh fading
channels, $|h_{1}|=|h_{2}|=\cdots=|h_{M}|\sim Rayleigh(2{\sigma}^{2})$ and
$\rho_{k,l}=1$. The outage probability $P_{out}^{IR}\left(M\right)$ can be
easily derived as

$P_{out}^{IR}\left(M\right)=1-\exp\left({-\frac{{{2^{\mathcal{R}/M}}-1}}{{{2\sigma^{2}\gamma_{T}}}}}\right).$
(79)

Putting (79) into (61), and by applying L’Hôpital’s rule and the method of
replacement with equivalent infinitesimal, the diversity order follows as

$\displaystyle d$
$\displaystyle=-\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\log\left({1-\exp\left({-\frac{{{2^{{\cal
R}/M}}-1}}{{{2\sigma^{2}\gamma_{T}}}}}\right)}\right)}}{{\log\left({{\gamma_{T}}}\right)}}$
$\displaystyle=\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{\left({{2^{{\cal
R}/M}}-1}\right)\exp\left({-\frac{{{2^{{\cal
R}/M}}-1}}{{{2\sigma^{2}\gamma_{T}}}}}\right)}}{{{2\sigma^{2}\gamma_{T}}\left({1-\exp\left({-\frac{{{2^{{\cal
R}/M}}-1}}{{{2\sigma^{2}\gamma_{T}}}}}\right)}\right)}}$
$\displaystyle=\mathop{\lim}\limits_{{\gamma_{T}}\to\infty}\frac{{{2^{{\cal
R}/M}}-1}}{{{2^{{\cal R}/M}}-1}}=1.$ (80)

It means that under fully correlated fading channels, the diversity order is
reduced to one.

##  V Numerical Results and Discussions

With the above analytical results, the performance of HARQ-IR over time-
correlated fading channels can be evaluated and optimal design of transmission
scheme is also enabled. In the following, we take systems with
${\rm{E}}(|h_{l}|^{2})=2\sigma_{l}^{2}=1$ and $\rho_{k,l}=\rho$ for $1\leq
k\neq l\leq M$ as examples for performance evaluation and optimal design.

###  V-A Verification of Analytical Results

To verify our analytical expressions for the performance metrics of HARQ-IR
over time-correlated Rayleigh fading channels, Monte-Carlo simulations are
conducted for comparison. For illustration, we set $\mathcal{R}=2\rm{bps/Hz}$
and $\rho=0.5$, and the outage probability versus transmit SNR $\gamma_{T}$ is
shown in Fig. 3. Apparently, there is a perfect match between the analytical
results and simulation results, which demonstrates the correctness of our
analytical results. Moreover, as expected, the outage probability
$P_{out}^{IR}(M)$ decreases as $M$ increases. For example, the outage
probability for $M=1$ is about $2.6*10^{-1}$ for a transmit SNR of
$10~{}\rm{dB}$. When $M$ is increased to $4$, the outage probability
significantly drops to $3*10^{-4}$. It demonstrates a notable performance gain
of HARQ-IR schemes.

Figure 3: Outage Probability $P_{out}^{IR}(M)$ versus transmit SNR
$\gamma_{T}$.

###  V-B Impact of Time Correlation

The impact of time correlation on the performance of HARQ-IR including the
outage probability, the average number of transmissions and the LTAT is now
investigated, and the results are shown in Fig. 4-6, respectively. Notice that
for fully correlated fading channels with $\rho=1$, the outage probability is
obtained as (79). Correspondingly, the average number of transmissions
$\bar{\mathcal{N}}$ and the LTAT $\bar{\mathcal{R}}$ can be obtained by
putting (79) into (9) and (10), respectively.

The outage probability versus the time correlation coefficient under various
$\gamma_{T}$ and $M$ is shown in Fig. 4. It is readily seen that the outage
probability increases with time correlation coefficient $\rho$. Specifically,
for the case with $M=4$ and $\gamma_{T}=7~{}\rm{dB}$, $P_{out}^{IR}(M)$
increases from $3*10^{-4}$ to $8*10^{-2}$ when $\rho$ increases from $0$ to
$1$. It indicates that time correlation does cause negative impact on the
outage performance. Additionally, it is observed that the gap between the
outage probabilities for two different $M$ becomes narrower when $\rho$ gets
higher. In other words, under highly correlated channels, further increase of
the number of transmissions will only lead to slight improvement on the outage
probability.

Figure 4: Outage Probability $P_{out}^{IR}\left(M\right)$ versus correlation
coefficient $\rho$.

In Fig. 5, the average number of transmissions $\bar{\mathcal{N}}$ versus
correlation coefficient $\rho$ is plotted. It can be seen that
$\bar{\mathcal{N}}$ increases with $\rho$ for $M>2$. This is also because of
the negative impact of the correlation in the channels. When the time
correlation increases, outage probability will increase, thus more HARQ rounds
are required to successfully deliver a message. Notice that when $M=1,2$, the
average number of transmissions $\bar{\mathcal{N}}$ is irrelevant to the time
correlation, which directly follows from (9).

Figure 5: Average number of transmissions $\bar{\mathcal{N}}$ versus
correlation coefficient $\rho$. ($\gamma_{T}=7$dB)

With respect to the effect of time correlation on the LTAT, similarly we can
find that the LTAT $\bar{\mathcal{R}}$ decreases with the increase of $\rho$
as shown in Fig. 6. For instance, the LTAT $\bar{\mathcal{R}}$ for the case
with $M=4$ decreases from $1.30~{}\rm{bps/Hz}$ to $1.05~{}\rm{bps/Hz}$ as
$\rho$ increases from $0$ to $1$. Interestingly, it can also be easily
observed that the LTAT shows opposite trends with the increase of the number
of transmissions under different correlation regions. More specifically, under
low-to-median correlation, the LTAT is improved when the number of
transmissions increases. However, when the time correlation is high,
additional transmission causes degradation of the LTAT when $M\geq 2$. This is
due to the twofold impact of increasing the number of transmissions. On one
hand, the increase of the number of transmissions would decrease the outage
probability. On the other hand, It would also increase the average number of
transmissions. The first impact dominates under low-to-median correlation,
while the second one dominates under high correlation. Therefore, from the
LTAT’s point of view, more transmissions may not be better when time
correlation is high.

Figure 6: LTAT $\bar{\mathcal{R}}$ versus correlation coefficient.
($\gamma_{T}=7$dB)

###  V-C Optimum Rate Design

As defined in (10), the LTAT is a complicated function of the transmission
rate $\mathcal{R}$ due to the implicit involvement of $\mathcal{R}$ in outage
probability and the average number of transmissions. To maximize the LTAT, the
rate should be properly designed. Mathematically, the problem of optimal rate
design can be formulated as

$\displaystyle\underset{\mathcal{R}}{\text{min}}$
$\displaystyle\bar{\mathcal{R}}=\frac{{\mathcal{R}\left({1-P_{out}^{IR}\left(M\right)}\right)}}{{\bar{\mathcal{N}}}}$
(81) s.t. $\displaystyle P_{out}^{IR}\left(M\right)\leq\varepsilon.$

With the approximation in (55) where $\kappa_{i}$, ${{\zeta}}$, and
${{\theta}}$ are irrelevant to the rate $\mathcal{R}$, the optimal rate can be
easily solved by using optimization tools.

Fig. 7 shows the optimal rate given various outage constraints $\varepsilon$
under different correlation scenarios. It is clear that the optimal
transmission rate increases with the transmit SNR and the allowable outage
probability $\varepsilon$. For examples, the optimum rate $\mathcal{R}$
increases by $3.2~{}\rm{bps/Hz}$ when the transmit SNR $\gamma_{T}$ is
increased from $0~{}\rm{dB}$ to $10~{}\rm{dB}$ for the case with
$\varepsilon=10^{-2}$ and $\rho=0.5$. It also increases from
$3.87~{}\rm{bps/Hz}$ to $6.57~{}\rm{bps/Hz}$ when the allowable outage
probability $\varepsilon$ increases from $0.01$ to $0.1$ for the case with
$\rho=0.5$ and $\gamma_{T}=10~{}\rm{dB}$. On the other hand, time correlation
of the channels has negative effect on the optimal rate.

Figure 7: Optimum transmission rate $R$ versus transmit SNR $\gamma_{T}$ for
different target outage probability $\varepsilon$ for $M=4$.

To further investigate the improvement of LTAT through optimal rate design,
the LTAT $\bar{\mathcal{R}}$ versus transmit SNR $\gamma_{T}$ for the schemes
with optimal rate design and a constant rate (which is set as the optimal rate
corresponding to the case with $\gamma_{T}=0~{}\rm{dB}$), is depicted in Fig.
8 by setting $\varepsilon=0.01$ and $M=4$. Apparently, notable improvement of
LTAT can be observed through optimal rate design and the contribution of the
optimal rate design becomes more significant when the transmit SNR gets
higher.

Figure 8: LTAT $\bar{\mathcal{R}}$ for the schemes with optimum rate design
and a constant rate.

##  VI Conclusions

Performance of HARQ-IR scheme operating over time-correlated Rayleigh fading
channels has been analyzed in this paper. By using polynomial fitting
technique, the PDF of the accumulated mutual information has been derived,
which enables the derivation of outage probability, average number of
transmission and LTAT in closed-forms. It has been found that the outage
probability can be written as a weighted sum of outage probabilities
corresponding to a number of Gamma RVs. Moveover, diversity order has been
analyzed and it has been revealed that full diversity can be achieved even
under time-correlated fading channels. The impact of time correlation in the
channels has also been investigated. It has been demonstrated that time
correlation has negative effect on the performance. Under highly correlated
channels, more transmissions would not necessarily lead to a higher LTAT and a
few transmissions may be sufficient. Finally, the analytical results have
enabled the optimal design of HARQ-IR scheme and optimal rate design has been
particularly discussed to demonstrate the significance of our analytical
results on HARQ-IR.

##  Appendix A Derivation of moments ${\mathcal{M}}(i)$

By definition, the $i$th moment of ${I_{K}^{IR}}$ is expressed as

${\cal
M}\left(i\right)=\int_{{x_{1}}=0}^{\infty}\cdots\int_{{x_{K}}=0}^{\infty}{{\left({\sum\limits_{k=1}^{K}{{{\log}_{2}}\left({1+{x_{k}}}\right)}}\right)}^{i}}\\\
\times{f_{{{\bf{\gamma}}_{1:K}}}}\left({{x_{1}},\cdots,{x_{K}}}\right)d{x_{1}}\cdots
d{x_{K}}.$ (82)

By substituting (III-B) into (82), it follows that

${\cal
M}\left(i\right)=\prod\limits_{k=1}^{K}{\frac{1}{{2{\sigma_{k}^{\prime}}{{}^{2}}\left({1-{\lambda_{k}}^{2}}\right)}}}\int\limits_{t=0}^{\infty}{{\rm{e}}^{-\left({1+\sum\limits_{k=1}^{K}{\frac{{{\lambda_{k}}^{2}}}{{1-{\lambda_{k}}^{2}}}}}\right)t}}\\\
\times{\int\limits_{0}^{\infty}{\cdots\int\limits_{0}^{\infty}{{{\left({\sum\limits_{k=1}^{K}{{{\log}_{2}}\left({1+{x_{k}}}\right)}}\right)}^{i}}\prod\limits_{k=1}^{K}{e^{-\frac{{{x_{k}}}}{{2{\sigma_{k}^{\prime}}{{}^{2}}\left({1-{\lambda_{k}}^{2}}\right)}}}}}}}\\\
\times{}_{0}{F_{1}}\left({;1;\frac{{{x_{k}}{\lambda_{k}}^{2}t}}{{2{\sigma_{k}^{\prime}}{{}^{2}}\left({1-{\lambda_{k}}^{2}}\right)^{2}}}}\right)d{x_{1}}\cdots
d{x_{k}}dt.$ (83)

Using binomial expansion, (83) can be rewritten as

${{\cal
M}\left(i\right)=\prod\limits_{k=1}^{K}{\frac{1}{{2{\sigma_{k}^{\prime}}{{}^{2}}\left({1-{\lambda_{k}}^{2}}\right)}}}\sum\limits_{{i_{1}}=0}^{i}{\sum\limits_{{i_{2}}=0}^{i-{i_{1}}}{\cdots\sum\limits_{{i_{K-1}}=0}^{i-\sum\limits_{l=1}^{K-2}{{i_{l}}}}}C_{i}^{{i_{1}}}}}\\\
\times C_{i-{i_{1}}}^{{i_{2}}}\cdots
C_{i-\sum\limits_{l=1}^{K-1}{{i_{l}}}}^{{i_{K}}}\int\limits_{t=0}^{\infty}{{\rm{e}}^{-\left({1+\sum\limits_{k=1}^{K}{\frac{{{\lambda_{k}}^{2}}}{{1-{\lambda_{k}}^{2}}}}}\right)t}}\\\
\times\prod\limits_{k=1}^{K}\int\limits_{0}^{\infty}\left({{\log}_{2}}\left({1+{x_{k}}}\right)\right)^{{i_{k}}}{e^{{-\frac{{{x_{k}}}}{{2{\sigma_{k}^{\prime}}{{}^{2}}{{\left({1-{\lambda_{k}}^{2}}\right)}}}}}}}\\\
\times{{{}_{0}{F_{1}}\left({;1;\frac{{{x_{k}}{\lambda_{k}}^{2}t}}{{2{\sigma_{k}^{\prime}}{{}^{2}}\left({1-{\lambda_{k}}^{2}}\right)^{2}}}}\right)}d{x_{k}}}dt.$
(84)

Then making change of variables yields

${\cal
M}\left(i\right)=\frac{{i!}}{{1+\sum\limits_{k=1}^{K}{\frac{{{\lambda_{k}}^{2}}}{{1-{\lambda_{k}}^{2}}}}}}\sum\limits_{\sum\limits_{l=1}^{K}{i_{l}}=i,{i_{l}}\geq
0}{\frac{1}{{{i_{1}}!{i_{2}}!\cdots{i_{K}}!}}}\\\
\times\int\limits_{t=0}^{\infty}{{{\rm{e}}^{-t}}\prod\limits_{k=1}^{K}{\int\limits_{0}^{\infty}{{e^{-y}}{{\log}_{2}}^{{i_{k}}}\left({1+{w_{k}}y}\right){}_{0}{F_{1}}\left({;1;{\varpi_{k}}yt}\right)}dy}}dt$
(85)

where
${\varpi_{k}}=\frac{{\frac{{{\lambda_{k}}^{2}}}{{1-{\lambda_{k}}^{2}}}}}{{1+\sum\limits_{l=1}^{K}{\frac{{{\lambda_{l}}^{2}}}{{1-{\lambda_{l}}^{2}}}}}}$,
${w_{k}}=2{\sigma_{k}^{\prime}}^{2}\left({1-{\lambda_{k}}^{2}}\right)$.

By applying Gaussian quadrature into (84) [33], it produces

${\mathcal{M}}\left(i\right)\approx\frac{{i!}}{{1+\sum\limits_{k=1}^{K}{\frac{{{\lambda_{k}}^{2}}}{{1-{\lambda_{k}}^{2}}}}}}\sum\limits_{\sum\limits_{l=1}^{K}{i_{l}}=i,{i_{l}}\geq
0}{\frac{1}{{{i_{1}}!{i_{2}}!\cdots{i_{K}}!}}}\\\
\times\int\limits_{t=0}^{\infty}{{{\mathop{\rm
e}\nolimits}^{-t}}\prod\limits_{k=1}^{K}{\sum\limits_{{q_{k}}=1}^{{N_{Q}}}{{\varrho_{{q_{k}}}}}{{\log}_{2}}^{{i_{k}}}\left({1+{w_{k}}{\xi_{{q_{k}}}}}\right){}_{0}{F_{1}}\left({;1;{\varpi_{k}}{\xi_{{q_{k}}}}t}\right)}}dt\\\
=\frac{{i!}}{{1+\sum\limits_{k=1}^{K}{\frac{{{\lambda_{k}}^{2}}}{{1-{\lambda_{k}}^{2}}}}}}\sum\limits_{\sum\limits_{l=1}^{K}{i_{l}}=i,{i_{l}}\geq
0}{\frac{1}{{{i_{1}}!{i_{2}}!\cdots{i_{K}}!}}\sum\limits_{{q_{k}}\in\left[{1,{N_{Q}}}\right],k\in\left[{1,K}\right]}}\\\
{\prod\limits_{k=1}^{K}{{\varrho_{{q_{k}}}}{{\log}_{2}}^{{i_{k}}}\left({1+{w_{k}}{\xi_{{q_{k}}}}}\right)}\int\limits_{t=0}^{\infty}{{{\mathop{\rm
e}\nolimits}^{-t}}\prod\limits_{k=1}^{K}{{}_{0}{F_{1}}\left({;1;{\varpi_{k}}{\xi_{{q_{k}}}}t}\right)}}dt}$
(86)

where $N_{Q}$ is the quadrature order, and the weights $\varrho_{q_{k}}$ and
abscissas ${{\xi_{q_{k}}}}$ are tabulated in [30, Table 25.9]. The residue
error becomes negligible if $N_{Q}$ is sufficiently large. By using the
following formula

$\displaystyle{\int\limits_{z=0}^{\infty}{{e^{-z}}\prod\limits_{l=1}^{K}{{}_{0}{F_{1}}\left({;1;{a_{l}}z}\right)}}dz}$
$\displaystyle=\int\limits_{z=0}^{\infty}{{e^{-z}}\prod\limits_{l=1}^{K}{\frac{1}{{2\pi{\rm
j}}}\int\limits_{{{\cal
C}_{l}}}{\frac{{\Gamma\left(s_{l}\right)}}{{\Gamma\left({1-s_{l}}\right)}}{{\left({-{a_{l}}z}\right)}^{-{s_{l}}}}d{s_{l}}}}}dz$
$\displaystyle={{\left({\frac{1}{{2\pi{\rm
j}}}}\right)}^{K}}\int\limits_{{{\cal C}_{1}}}{{\cdots\int\limits_{{{\cal
C}_{K}}}{\frac{{\Gamma\left({1-\sum\limits_{l=1}^{K}{{s_{l}}}}\right)\prod\limits_{l=1}^{K}{\Gamma\left({{s_{l}}}\right)}}}{{\prod\limits_{l=1}^{K}{\Gamma\left({1-{s_{l}}}\right)}}}}}}$
$\displaystyle\times{{\left({-{a_{1}}}\right)}^{-{s_{1}}}}\cdots{{\left({-{a_{K}}}\right)}^{-{s_{K}}}}d{s_{1}}\cdots
d{s_{K}}$
$\displaystyle=\Psi_{2}^{\left(K\right)}\left({1;\underbrace{1,1,\cdots,1}_{K-terms};{a_{1}},{a_{2}},\cdots,{a_{K}}}\right)$
(87)

where ${\rm j}=\sqrt{-1}$, and $\Psi_{2}^{\left(K\right)}(;;)$ denotes
confluent form of Lauricella hypergeometric function [31, Definition A.20],
[32], the final expression for (86) then follows as (III-B).

##  Appendix B Proof of Invertibility of Matrix $\bf{A}$

To complete the proof, it suffices to show that $\bf{A}$ is a positive
definite matrix. For arbitrary $(N+1)\times 1$ real vector $\bf u$, a function
$\bf{u^{\rm T}}{\bf A}{{\bf{u}}}$ can be derived as

${\bf{u^{\rm T}}}{\bf
A}{{\bf{u}}}=\int_{0}^{\infty}{\varphi\left(x\right)\left({{\bf{u}}^{\rm{T}}\boldsymbol{\mathcal{P}}(x)}\right)^{2}dx}.$
(88)

It is straightforward that
$\int_{0}^{\infty}{\varphi\left(x\right)\left({{\bf{u}}^{\rm{T}}\boldsymbol{\mathcal{P}}(x)}\right)^{2}dx}\geq
0$ and the equality holds if and only if $\bf{u}=\bf{0}$. It follows that
$\bf{A}\succ 0$, thus $\bf{A}$ is non-singular, that is, $\bf{A}$ is
invertible.

##  Appendix C Proof of Property 1

The proof is completed by contradiction. First, denote two measures as
$d\mu_{1}(x)$ and $d\mu_{2}(x)$, and their corresponding monic orthogonal
polynomials as $\boldsymbol{\mathcal{P}}_{1}(x)$ and
$\boldsymbol{\mathcal{P}}_{2}(x)$, respectively. As aforementioned, given the
measure $d\mu(x)$, the optimal solution to the problem of (24) is unique. We
then represent the unique optimal solutions corresponding to the two measures
$d\mu_{1}(x)$ and $d\mu_{2}(x)$ as $\boldsymbol{\eta}_{1}$ and
$\boldsymbol{\eta}_{2}$, respectively. Meanwhile, the corresponding optimal
polynomials are denoted as
$\hat{\psi}_{N}^{(1)}(x)={\boldsymbol{\eta}_{1}}^{\rm
T}\boldsymbol{\mathcal{P}}_{1}(x)$ and
$\hat{\psi}_{N}^{(2)}(x)={\boldsymbol{\eta}_{2}}^{\rm
T}\boldsymbol{\mathcal{P}}_{2}(x)$, while the corresponding minimal MSEs are
${{\mathcal{S}}_{mse}}({\boldsymbol{\eta}}_{1}|d\mu_{1}(x),N)$ and
${{\mathcal{S}}_{mse}}({\boldsymbol{\eta}}_{2}|d\mu_{2}(x),N)$.

Assume that
${{\mathcal{S}}_{mse}}({\boldsymbol{\eta}}_{1}|d\mu_{1}(x),N)\neq{{\mathcal{S}}_{mse}}({\boldsymbol{\eta}}_{2}|d\mu_{2}(x),N)$,
which implies that $\hat{\psi}_{N}^{(1)}(x)\neq\hat{\psi}_{N}^{(2)}(x)$.
Without loss of generality, we consider the case with
${{\mathcal{S}}_{mse}}({\boldsymbol{\eta}}_{1}|d\mu_{1}(x),N)>{{\mathcal{S}}_{mse}}({\boldsymbol{\eta}}_{2}|d\mu_{2}(x),N)$.
It is known that any polynomial $\hat{\psi}_{N}(x)\in\mathbb{P}_{N}$ has a
unique coordinate $\boldsymbol{\eta}$ regarding to a given basis of orthogonal
polynomials. Thus by taking $\boldsymbol{\mathcal{P}}_{1}(x)$ as a basis, we
have $\hat{\psi}_{N}^{(2)}(x)={\boldsymbol{\eta}_{3}}^{\rm
T}\boldsymbol{\mathcal{P}}_{1}(x)$, where ${\boldsymbol{\eta}}_{3}$ is the
unique coordinate associated with $\hat{\psi}_{N}^{(2)}(x)$ on the basis of
$\boldsymbol{\mathcal{P}}_{1}(x)$. Considering the unique optimality of
${\boldsymbol{\eta}}_{1}$ given the measure $d\mu_{1}(x)$, we have
${{\mathcal{S}}_{mse}}({\boldsymbol{\eta}}_{1}|d\mu_{1}(x),N)<{{\mathcal{S}}_{mse}}({\boldsymbol{\eta}}_{3}|d\mu_{1}(x),N)={{\mathcal{S}}_{mse}}({\boldsymbol{\eta}}_{2}|d\mu_{2}(x),N)$.
This is contradictory with our assumption. Therefore,
${{\mathcal{S}}_{mse}}({\boldsymbol{\eta}}_{1}|d\mu_{1}(x),N)={{\mathcal{S}}_{mse}}({\boldsymbol{\eta}}_{2}|d\mu_{2}(x),N)$
and $\hat{\psi}_{N}^{(1)}(x)=\hat{\psi}_{N}^{(2)}(x)$. Then the proof follows.

## References

  * [1] D. Tse and P. Viswanath, _Fundamentals of wireless communication_. Cambridge university press, 2005.
  * [2] E. Dahlman, S. Parkvall, J. Sköld, and P. Beming, _3G evolution: HSPA and LTE for mobile broadband_. Academic press, 2010.
  * [3] P. Wu and N. Jindal, “Performance of hybrid-ARQ in block-fading channels: a fixed outage probability analysis,” _IEEE Trans. Commun._ , vol. 58, no. 4, pp. 1129–1141, Apr. 2010.
  * [4] J. Wang, S. Park, D. J. Love, and M. D. Zoltowski, “Throughput delay tradeoff for wireless multicast using hybrid-ARQ protocols,” _IEEE Trans. Commun._ , vol. 58, no. 9, pp. 2741–2751, Sept. 2010.
  * [5] L. Szczecinski, S. R. Khosravirad, P. Duhamel, and M. Rahman, “Rate allocation and adaptation for incremental redundancy truncated HARQ,” _IEEE Trans. Commun._ , vol. 61, no. 6, pp. 2580–2590, Jun. 2013.
  * [6] B. Makki, A. Graell i Amat, and T. Eriksson, “Green communication via power-optimized HARQ protocols,” _IEEE Trans. Veh. Technol._ , vol. 63, no. 1, pp. 161–177, Jan. 2014.
  * [7] B. Makki and T. Eriksson, “On the average rate of HARQ-based quasi-static spectrum sharing networks,” _IEEE Trans. Wireless Commun._ , vol. 11, no. 1, pp. 65–77, Jan. 2012.
  * [8] P. Larsson, B. Smida, T. Koike-Akino, and V. Tarokh, “Analysis of network coded HARQ for multiple unicast flows,” in _Proc. IEEE Int. Conf. Commun. (ICC’10)_ , May 2010, pp. 1–6.
  * [9] S. Khosravirad, L. Szczecinski, and F. Labeau, “Rate adaptation for cooperative HARQ,” _IEEE Trans. Commun._ , vol. 62, no. 5, pp. 1469–1479, May 2014.
  * [10] J. Choi, D. To, Y. Wu, and S. Xu, “Energy-delay tradeoff for wireless relay systems using HARQ with incremental redundancy,” _IEEE Trans. Wireless Commun._ , vol. 12, no. 2, pp. 561–573, Feb. 2013.
  * [11] A. Chelli and M.-S. Alouini, “On the performance of hybrid-ARQ with incremental redundancy and with code combining over relay channels,” _IEEE Trans. Wireless Commun._ , vol. 12, no. 8, pp. 3860–3871, Aug. 2013\.
  * [12] S. M. Kim, W. Choi, T. W. Ban, and D. K. Sung, “Optimal rate adaptation for hybrid ARQ in time-correlated Rayleigh fading channels,” _IEEE Trans. Wireless Commun._ , vol. 10, no. 3, pp. 968–979, Mar. 2011.
  * [13] H. Jin, C. Cho, N.-O. Song, and D. K. Sung, “Optimal rate selection for persistent scheduling with HARQ in time-correlated Nakagami-m fading channels,” _IEEE Trans. Wireless Commun._ , vol. 10, no. 2, pp. 637–647, Feb. 2011.
  * [14] W. C. Jakes and D. C. Cox, _Microwave mobile communications_. Wiley-IEEE Press, 1994.
  * [15] _LTE; Evolved Universal Terrestrial Radio Access (E-UTRA); Medium Access Control (MAC) protocol specification (3GPP TS 36.321 version 8.12.0 Release 8)_ , European Telecommunications Standards Institute Std., Nov. 2011. [Online]. Available: http://www.etsi.org/deliver/etsi_ts/136300_136399/136321/09.04.00_60/ts_136321v090400p.pdf
  * [16] S. Sesia, I. Toufik, and M. Baker, _LTE: the UMTS long term evolution_ , 2nd ed. Wiley Online Library, 2011.
  * [17] J. Choi, W. Xing, D. To, Y. Wu, and S. Xu, “On the energy efficiency of a relaying protocol with HARQ-IR and distributed cooperative beamforming,” _IEEE Trans. Wireless Commun._ , vol. 12, no. 2, pp. 769–781, Feb. 2013\.
  * [18] N. C. Beaulieu and K. T. Hemachandra, “Novel simple representations for Gaussian class multivariate distributions with generalized correlation,” _IEEE Trans. Inf. Theory_ , vol. 57, no. 12, pp. 8072–8083, Dec. 2011.
  * [19] Y. Hu and A. Ribeiro, “Optimal wireless communications with imperfect channel state information,” _IEEE Trans. Signal Process._ , vol. 61, no. 11, pp. 2751–2766, Jun. 2013.
  * [20] A. Chelli and M. Pätzold, “On the performance of hybrid-ARQ with code combining over double Rayleigh fading channels,” in _Proc. IEEE International Symposium on Personal Indoor and Mobile Radio Communications (PIMRC’11)_ , Sept. 2011, pp. 2014–2019.
  * [21] A. Chelli, E. Zedini, M.-S. Alouini, J. Barry, and M. Patzold, “Performance and delay analysis of hybrid ARQ with incremental redundancy over double Rayleigh fading channels,” _IEEE Trans. Wireless Commun._ , vol. 13, no. 11, pp. 6245–6258, Nov. 2014.
  * [22] B. Makki and T. Eriksson, “On the performance of MIMO-ARQ systems with channel state information at the receiver,” _IEEE Trans. Commun._ , vol. 62, no. 5, pp. 1588–1603, May 2014.
  * [23] M. D. Springer and W. E. Thompson, “The distribution of products of independent random variables,” _SIAM Journal on Applied Mathematics_ , vol. 14, no. 3, pp. 511–526, May 1966.
  * [24] Z. A. Lomnicki, “On the distribution of products of random variables,” _Journal of the Royal Statistical Society. Series B (Methodological)_ , vol. 29, no. 3, pp. 513–524, 1967.
  * [25] M. D. Springer and W. E. Thompson, “The distribution of products of Beta, Gamma and Gaussian random variables,” _SIAM Journal on Applied Mathematics_ , vol. 18, no. 4, pp. 721–737, Jun. 1970.
  * [26] J. Salo, H. El-Sallabi, and P. Vainikainen, “The distribution of the product of independent Rayleigh random variables,” _IEEE Trans. Antennas Propag._ , vol. 54, no. 2, pp. 639–643, Feb. 2006.
  * [27] F. Yilmaz and M.-S. Alouini, “Product of shifted exponential variates and outage capacity of multicarrier systems,” in _Proc. European Wireless Conference (EW’09)_ , May 2009, pp. 282–286.
  * [28] ——, “Outage capacity of multicarrier systems,” in _Proc. IEEE International Conference on Telecommunications (ICT’10)_ , Apr. 2010, pp. 260–265.
  * [29] S. Primak, V. Kontorovitch, and V. Lyandres, _Stochastic methods and their applications to communications: stochastic differential equations approach_. John Wiley & Sons, 2005.
  * [30] M. Abramowitz and I. A. Stegun, _Handbook of mathematical functions: with formulas, graphs, and mathematical tables_. Dover Publications, 1965.
  * [31] A. Mathai, R. K. Saxena, and H. J. Haubold, _The H-function_. Springer, 2009.
  * [32] M. Saigo and V. K. Tuan, “Some integral representations of multivariable hypergeometric functions,” _Rendiconti del Circolo Matematico di Palermo_ , vol. 41, no. 1, pp. 69–80, Jan. 1992.
  * [33] G. Dahlquist and Å. Björck, _Numerical Methods in Scientific Computing, Volume I_. Society for Industrial and Applied Mathematics, 2008.
  * [34] W. Gautschi, _Orthogonal Polynomials: Computation and Approximation_ , ser. Numerical mathematics and scientific computation. Oxford University Press, 2004.
  * [35] L. Zheng and D. N. C. Tse, “Diversity and multiplexing: a fundamental tradeoff in multiple-antenna channels,” _IEEE Trans. Inf. Theory_ , vol. 49, no. 5, pp. 1073–1096, May 2003.
  * [36] S. Kalyani and R. Karthik, “The asymptotic distribution of maxima of independent and identically distributed sums of correlated or non-identical Gamma random variables and its applications,” _IEEE Trans. Commun._ , vol. 60, no. 9, pp. 2747–2758, Spet. 2012.
  * [37] V. Aalo, T. Piboongungon, and G. Efthymoglou, “Another look at the performance of MRC schemes in Nakagami-m fading channels with arbitrary parameters,” _IEEE Trans. Commun._ , vol. 53, no. 12, pp. 2002–2005, Dec. 2005.

|  Zheng Shi received the B.S. degree in communication engineering from Anhui
Normal University, China, in 2010 and the M.S. degree in communication and
information system from Nanjing University of Posts and Telecommunications
(NUPT), China, in 2013. Since Sep. 2013, he has been a Ph.D. student in
Department of Electrical and Computer Engineering, University of Macau, Macao.
His research interests include hybrid automatic repeat request (HARQ)
protocols, cooperative communications, full-duplex communications, massive
MIMO, and heterogeneous wireless networks.
---|---
|  Haichuan Ding is currently a Ph.D. student at the University of Florida. He
received the B.Eng. and M.S. degrees in electrical engineering from Beijing
Institute of Technology (BIT), Beijing, China, in 2011 and 2014, respectively.
From 2012 to 2014, he was with the Department of Electrical and Computer
Engineering, University of Macau, as a visiting student. During his M.S.
studies, he mainly worked on the analysis of HARQ techniques using the tools
of stochastic geometry. His current research is focused on cognitive radio
networks and security and privacy in distributed systems.
---|---
|  Shaodan Ma received her double Bachelor degrees in Science and Economics,
and her Master degree of Engineering, from Nankai University, Tianjin, China.
She obtained her Ph. D. degree in electrical and electronic engineering from
the University of Hong Kong, Hong Kong, in 2006. After graduation, she joined
the University of Hong Kong as a Postdoctoral Fellow. Since August 2011, she
has been with the University of Macau and is now an Associate Professor there.
She was a visiting scholar in Princeton University in 2010 and is currently an
Honorary Assistant Professor in the University of Hong Kong. Her research
interests are in the general areas of signal processing and communications,
particularly, transceiver design, resource allocation and performance
analysis.
---|---
|  Kam-Weng Tam (S’91-M’01-SM’05) received the B.Sc. and joint Ph.D. degrees
in electrical and electronics engineering from the University of Macau, Taipa,
Macao, China, and the University of Macau and Instituto Superior Técnico
(IST), Technical University of Lisbon, Lisbon, Portugal, in 1993 and 2000,
respectively. From 1993 to 1996, he was with the Instituto de Engenharia de
Sistemas e Computadores (INESC), Lisbon, Portugal, where he participated in
research and development on a broad range of applied microwave technologies
for satellite communication systems. From July 2000 to December 2001, he was
the Director of the Instituto de Engenharia de Sistemas e Computadores
(INESC)-Macau. In 2001, he cofounded the microelectronic design house Chipidea
Microelectrónica, Macau, China, where until 2003 he was the General Manager.
Since 1996, he has been with the University of Macau, where he is currently a
Professor and the Associate Dean (Research and Graduate Studies) with the
Faculty of Science and Technology. He has authored or coauthored over 100
journal and conference papers. His research interests have concerned
multifunctional microwave circuits, RFID, UWB for material analysis and
terahertz technology. Dr. Tam was interim secretary for the establishment of
the Macau Section in 2003. He supervised two IEEE Microwave Theory and
Techniques Society (IEEEMTT-S) Undergraduate Scholarship recipients in 2002
and 2003. He was founder of the IEEE Macau AP/MTT Joint Chapter in 2010 and
was chair in 2011–2012. He was a member of the organizing committees of 21
international and local conferences including co-chair of APMC2008, co-chair
of the Technical Program, IEEE MTT-S International Microwave Workshop Series
on Art of Miniaturizing RF and Microwave Passive Components (2008), and co-
chair of ISAP2010.
---|---