**2. Atmospheric turbulence**

2 Will-be-set-by-IN-TECH

for mitigation of the degrading effects of atmospheric turbulence are error control coding in conjunction with interleaving [70], [83] and maximum likelihood sequence detection (MLSD) [82]. However, for most scenarios the first one requires large-size interleavers to achieve the promised coding gains. On the other hand, MLSD requires complicated multidimensional integrations and suffers from excessive computational complexity. Some sub-optimal temporal-domain fading mitigation techniques are further explored in [83] and

Another promising solution is the use of diversity techniques and the most popular scheme is the spatial diversity, i.e., the employment of multiple transmit/receive apertures, a well known diversity technique in Radio-Frequency (RF) systems [46, 48, 49, 61, 69, 72, 78, 80]. By using multiple apertures at the transmitter and/or the receiver, the inherent redundancy of spatial diversity has the potential to significantly enhance the performance. Moreover, the possibility for temporal blockage of the laser beams by obstructions is further reduced and longer distances can be covered through heavier weather conditions [46]. Concerning the performance analysis of FSO systems employing spatial diversity, the technical literature is rather rich. Representative past examples can be found in [14, 19, 21, 26, 27, 29, 46, 54–56, 62,

On the other hand, various statistical models, e.g. the log normal, the gamma gamma (G − G), the I-K, the K, the negative exponential, and the Rician log normal distribution, have been used in order to describe the optical channel characteristics with respect to the atmospheric turbulence strength [4, 9, 10, 12, 23, 30, 32, 35, 38, 44, 47, 51–53, 59, 65, 71]. Recently, Al-Habash et al. proposed the G-G distribution [4] as a tractable mathematical model for atmospheric turbulence. This model is a two parameter distribution which is based on a doubly stochastic theory of scintillation and assumes that small-scale irradiance fluctuations are modulated by large-scale irradiance fluctuations of the propagating wave, both governed by independent gamma distributions. This distribution has become the dominant fading channel model for FSO links due to its excellent agreement with measurement data for a wide range of

For many practical FSO applications, however, irradiance is temporally correlated. Thus the derivation of a correlated G − G model is of significant theoretical and practical interest. It is noted that multivariate distributions have recently attracted the interest within the research community due to their importance in studying the performance of diversity systems operating over a multipath fading channels, see, e.g., [1, 2, 5, 15, 33, 45, 49, 58, 64, 68] and references therein. These distributions are particularly useful in the performance analysis of practical systems configurations where antenna branches are closely spaced and the correlation between diversity signals cannot be ignored [68, Chapt. 9]. In the past, several spatial correlation models have been proposed [5, 68] and used for the performance evaluation of wireless communication systems over correlated fading channels. Among them, the exponential correlation model has gained particular interest [2, 33, 45, 64]. This model corresponds to the scenario of multichannel reception by equispaced diversity antennas, in which the correlation between pairs of combined signals decays as the spacing between

In this chapter, we investigate the performance of multiple-input-multiple output (MIMO) FSO links over both independent and identically distributed (i.i.d.) G − G turbulence channels as well as for exponentially correlated G − G turbulent channels. A closed-form expression

[84].

63, 69, 73–75].

turbulence conditions [4].

antenna branches increases [68].

The refractive index of the atmosphere in the area between the transmitter and the receiver of a wireless optical link fluctuates randomly due to the atmospheric turbulence [12]. These fluctuations are induced mainly due to temperature oscillations among the atmosphere, the ground and the oceans [8, 12, 13, 24, 39]. More specifically atmospheric turbulence is a phenomenon belonging to different spatial and temporal scales. In the Planetary Boundary Layer (PBL), where human activities take place, it is generated by the wind's interaction with the earth's surface, which is said to be in a state of turbulent motion [18, 20, 37, 76, 77].

Turbulence is responsible for the transfer of heat, matter, and momentum within the PBL. However, it is random in nature and remains a complicated phenomenon with many unsolved aspects. The scientific community relies on the combination of experiments, theory, and computer models to understand it. Turbulence is created by thermal convection, wind shear and by the wind flowing over ground obstacles. Within the PBL, turbulence has a diurnal variation reaching a maximum about midday when the solar radiation is at a maximum. Solar radiation heats the surface which, in turn radiates heat to the air above it that becomes warmer and more buoyant and rises while cooler, denser air descends to displace it. The resulting vertical movement of air, together with flow disturbances around surface obstacles, makes low-level winds extremely irregular thus turbulent. The turbulence intensity depends primarily on the temperature lapse rate, i.e. *λ* = *dT*/*dH*, which is essentially the rate of temperature increase or decrease with increasing height. Under unstable conditions the temperature decreases with height and a hypothetical parcel of air which is warmer than its surrounding air would tend to rise [18, 20, 37, 76, 77].

Turbulence in the atmosphere is carried by rotational-like motions called eddies, which exist in different length scales characterized by different velocity and time scales. The larger eddies are unstable and break up into smaller ones with the subsequent transfer of the kinetic energy of the initial large eddy. The smaller eddies in turn break up into even smaller eddies and the energy is passed on to the new smaller ones. The energy is passed down from the larger scales to the smaller ones until the viscosity of the fluid (in this case air) can dissipate the

#### 4 Will-be-set-by-IN-TECH 418 Optical Communication

kinetic energy into internal energy (or thermal energy). This is called the energy cascade and it is one of the main characteristics of turbulent motion. Other important features of turbulence are irregularity (or chaotic), diffusivity (mixing), turbulent diffusion (molecular diffusivity), rotationality (always three dimensional, the mechanism that aids the energy cascade), dissipation (the transfer of energy from larger to smaller eddies), length scales of turbulent eddies. The size of eddies spans from the order of a few millimetres to meters, namely the inner and outer scales, respectively [18, 20, 37, 76, 77]. This inner and outer scales is the main phenomenon for the choice of the appropriate distribution for turbulence's mathematical representation.

It is commonly known that in turbulent flow the actual flow velocity is broken into the mean velocity *U* plus the fluctuating turbulent velocity component *u*, in the three directions, respectively. Thus, in the *x*-axis the instantaneous flow velocity is *Ux* + *ux*. The values of the turbulent components namely *ux*, *uy* and *uz* may not be expressed as functions of time, but a statistical description is only possible. Turbulence within the PBL can degrade the performance of free-space optical links, particularly over ranges of the order of 1 km or longer. This phenomenon is not caused by the turbulent eddies which possess different velocities but only by the parcels of air with different temperatures (and thus different densities) rising or descending that cross the path of the FSO links. They cause inhomogeneities in the temperature and pressure profiles of the atmosphere and lead to variations of the refractive index along the transmission path. Based on the fluctuations of the air density, the scientific community has developed the refractive index structure parameter *C*<sup>2</sup> *<sup>n</sup>* (related to the temperature structure one, i.e. *C*<sup>2</sup> *<sup>T</sup>*) which takes different values depending on the strength of turbulence. The parameter *C*<sup>2</sup> *<sup>n</sup>* may be measured experimentally of computed theoretically if one knows the outer scale of turbulence (i.e., the large eddy size scale) and the potential refractive index [18, 20, 37, 76, 77]. Additionally, as we will demonstrate below, there are many mathematical models for the estimation of the *C*<sup>2</sup> *<sup>n</sup>* value [6, 7, 44, 52]

Thus, these variation of the refraction index in the free space area that the beam of the optical link propagates and causes deflections of the light beam into and mostly out of the transmit path [17, 52]. This random radiation of the laser beam results in fluctuations of the optical signal's irradiance at the receiver's side. This phenomenon is the so-called scintillation [4, 12, 17, 31, 42–44, 51]. The influence of scintillation in the performance of the wireless optical communication systems is very strong for the terrestrial links because induces fading of the signal arriving at the receiver in a random way. Thus, in order to estimate the optical signal arriving at the receiver it is necessary to study the appropriate statistical distribution which describes the fading statistics of each area.

Many statistical models have been proposed for the simulation of these fading statistics caused by the atmospheric turbulence effect. Some of them have been arising from experimental results, while, others, from theoretical studies. It is obvious that each location has irregularities, depending on the ground's morphology, the weather conditions, the time of the day [7, 24] and the turbulence strength. The proposed statistical models concern weak, moderate, strong or very strong turbulence conditions and the turbulence strength can be estimated through the turbulence parameter *C*<sup>2</sup> *<sup>n</sup>*, which depends on many parameters of the weather conditions.

One of the statistical parameters that we are using for the practical estimation of scintillation's influence at the wireless optical links' performance,is the scintillation index which is given by

the following mathematical expression [44]:

4 Will-be-set-by-IN-TECH

kinetic energy into internal energy (or thermal energy). This is called the energy cascade and it is one of the main characteristics of turbulent motion. Other important features of turbulence are irregularity (or chaotic), diffusivity (mixing), turbulent diffusion (molecular diffusivity), rotationality (always three dimensional, the mechanism that aids the energy cascade), dissipation (the transfer of energy from larger to smaller eddies), length scales of turbulent eddies. The size of eddies spans from the order of a few millimetres to meters, namely the inner and outer scales, respectively [18, 20, 37, 76, 77]. This inner and outer scales is the main phenomenon for the choice of the appropriate distribution for turbulence's

It is commonly known that in turbulent flow the actual flow velocity is broken into the mean velocity *U* plus the fluctuating turbulent velocity component *u*, in the three directions, respectively. Thus, in the *x*-axis the instantaneous flow velocity is *Ux* + *ux*. The values of the turbulent components namely *ux*, *uy* and *uz* may not be expressed as functions of time, but a statistical description is only possible. Turbulence within the PBL can degrade the performance of free-space optical links, particularly over ranges of the order of 1 km or longer. This phenomenon is not caused by the turbulent eddies which possess different velocities but only by the parcels of air with different temperatures (and thus different densities) rising or descending that cross the path of the FSO links. They cause inhomogeneities in the temperature and pressure profiles of the atmosphere and lead to variations of the refractive index along the transmission path. Based on the fluctuations of the air density, the

scientific community has developed the refractive index structure parameter *C*<sup>2</sup>

if one knows the outer scale of turbulence (i.e., the large eddy size scale) and the potential refractive index [18, 20, 37, 76, 77]. Additionally, as we will demonstrate below, there are

Thus, these variation of the refraction index in the free space area that the beam of the optical link propagates and causes deflections of the light beam into and mostly out of the transmit path [17, 52]. This random radiation of the laser beam results in fluctuations of the optical signal's irradiance at the receiver's side. This phenomenon is the so-called scintillation [4, 12, 17, 31, 42–44, 51]. The influence of scintillation in the performance of the wireless optical communication systems is very strong for the terrestrial links because induces fading of the signal arriving at the receiver in a random way. Thus, in order to estimate the optical signal arriving at the receiver it is necessary to study the appropriate statistical distribution which

Many statistical models have been proposed for the simulation of these fading statistics caused by the atmospheric turbulence effect. Some of them have been arising from experimental results, while, others, from theoretical studies. It is obvious that each location has irregularities, depending on the ground's morphology, the weather conditions, the time of the day [7, 24] and the turbulence strength. The proposed statistical models concern weak, moderate, strong or very strong turbulence conditions and the turbulence strength can be

One of the statistical parameters that we are using for the practical estimation of scintillation's influence at the wireless optical links' performance,is the scintillation index which is given by

*<sup>n</sup>* (related to

*<sup>T</sup>*) which takes different values depending on the strength

*<sup>n</sup>* value [6, 7, 44, 52]

*<sup>n</sup>*, which depends on many parameters of the

*<sup>n</sup>* may be measured experimentally of computed theoretically

mathematical representation.

the temperature structure one, i.e. *C*<sup>2</sup>

describes the fading statistics of each area.

estimated through the turbulence parameter *C*<sup>2</sup>

weather conditions.

many mathematical models for the estimation of the *C*<sup>2</sup>

of turbulence. The parameter *C*<sup>2</sup>

$$
\sigma\_I^2 = \frac{ - *^2}{*^2} \tag{1}**$$

with *I* being the optical signal's irradiance at the receiver and <> represents the ensemble average value. In the weak scintillation theory, under the assumption of plane wave propagation, the scintillation index is proportional to the Rytov variance and is given as [44]:

$$
\sigma\_{I,R}^2 = 1.23 \mathcal{C}\_n^2 k^{7/6} L^{11/6} \tag{2}
$$

where *C*<sup>2</sup> *<sup>n</sup>* is the parameter of turbulence, *k* = 2*π*/*λ* is the optical wavenumber, while *L* is the link's length [44, 49, 52].

If we assume spherical wave propagation, *σI*,*<sup>R</sup>* can be expressed as [44]:

$$
\sigma\_{I,R}^2 = 0.5 \mathbf{C}\_n^2 \mathbf{k}^{7/6} \mathbf{L}^{11/6} \tag{3}
$$

In order to estimate the influence of the atmospheric turbulence conditions in the wireless optical communication system's performance, we assume a horizontal propagation path of up to a few kilometers (usually the FSO link is smaller than 4-5 km), where the turbulence strength value is considered as constant [50, 52]. Many models has been proposed for the estimation of the turbulence parameter *C*<sup>2</sup> *<sup>n</sup>* [6, 12, 34, 44]. The most widely employed in research literature, are the so-called Hufnagel-Valley model and the Hufnagel and Stanley model, the SLC Day model and the SLC Night model [6, 7, 12, 34, 44, 49, 52].

The estimation of the parameter *C*<sup>2</sup> *<sup>n</sup>* through the Hufnagel-Valley model (or HV5/7 model) [7, 44], depends on the wind speed, *v*, and the altitude, *h*, where the wireless optical link operates. The values 5 and 7 of the abbreviation of the Hufnagel-Valley model (i.e. HV5/7) correspond to the atmospheric coherence length (*r*0) in cm and the isoplantic angle (*θ*0) in *μrad* respectively, for *λ* = 0.55*μm* while for *λ* = 1.315*μm* the above values are 14cm and 20 *μrad* respectively [7, 44, 52]. The mathematical expression of the HV5/7 is the following [7, 44, 51],

$$\mathcal{C}\_{n}^{2}(h) = 0.00594(\mu/27)^{2}(10^{-5}h)^{10}e^{h/1000} + 2.7 \times 10^{-6}e^{-h/1500} + \mathcal{C}\_{n}^{2}(0)e^{-h/1000} \tag{4}$$

where *C*<sup>2</sup> *<sup>n</sup>*(0) is the value of *C*<sup>2</sup> *<sup>n</sup>* at the sea level in *m*−2/3. In general, the turbulence parameter *C*2 *<sup>n</sup>* varies between the values 10−<sup>17</sup> *m*−2/3 and 10−<sup>13</sup> *m*−2/3 for weak up to very strong turbulence cases, respectively, [44, 51].

Another significant model for the estimation of the turbulence parameter is proposed by Hufnagel and Stanley [12, 34, 52]. Its mathematical expression is the following,

$$\mathbf{C}\_{\rm n}^{2}(h) = \mathbf{K}\_{0}h^{-1/3}e^{-h/h\_{0}} \tag{5}$$

with *K*<sup>0</sup> being the turbulence strength parameter and depends on the characteristics of the specific location where the wireless optical link has been installed and *h* is the altitude.

Another model for the estimation of the atmospheric turbulence parameter, *C*<sup>2</sup> *<sup>n</sup>*, is the SLC Day model which depends only on the height where the FSO link operates [7, 52]. The mathematical expression for this model, which is accurate enough for daytime hours, is given as:

$$C\_n^2(h) = \begin{cases} 1.7 \times 10^{-14} \text{ if } 0 < h < 18.5 \text{m;}\\ \frac{3.13 \times 10^{-13}}{h^{1.05}} & \text{if } 18.5 < h < 240 \text{m;}\\ 1.3 \times 10^{-15} \text{ if } 240 < h < 880 \text{m;}\\ \frac{8.87 \times 10^{-7}}{h^3} & \text{if } 880 < h < 7200 \text{m;}\\ \frac{2.0 \times 10^{-16}}{h^{0.5}} & \text{if } 7200 < h < 20000 \text{m;} \end{cases} \tag{6}$$

while for accurate results for night hours are given by the SLC Night model [7, 52], depends only on the height where the optical link operates like the previous model are given by the following expression:

$$C\_n^2(h) = \begin{cases} 8.4 \times 10^{-15} \text{ if } 0 < h < 18.5 \text{m};\\ \frac{2.87 \times 10^{-12}}{h^2} & \text{if } 18.5 < h < 110 \text{m};\\ 2.5 \times 10^{-16} \text{ if } 110 < h < 1500 \text{m};\\ \frac{8.87 \times 10^{-7}}{h^3} & \text{if } 1500 < h < 7200 \text{m};\\ \frac{2.0 \times 10^{-16}}{h^{0.5}} & \text{if } 7200 < h < 20000 \text{m}; \end{cases} \tag{7}$$

The above mentioned models for the estimation of the turbulence parameter *C*<sup>2</sup> *<sup>n</sup>*, (i.e HV5/7, Hufnagel and Stanley, SLC Day and SLC Night), are valid mainly for wireless optical links which have been installed over terrestrial area. Thus, the estimation of the turbulence parameter for paths over maritime area can be done through the following accurate approximation which is valid for low altitude and for specific constant parameters *c*1, *c*2, *c*3, *c*<sup>4</sup> and *c*5. The values of the constant parameters are given in [41, 44] and the corresponding mathematical form is the following:

$$\mathbf{C}\_{n}^{2}(h) = c\_{1} + c\_{2}e^{-h/c\_{3}} + c\_{4}e^{-h/c\_{5}} \tag{8}$$

Obviously, as we mentioned above, for the estimation of the turbulence strength through the above Equations, we assumed that the value of *C*<sup>2</sup> *<sup>n</sup>* remains constant for relatively long time interval and for the whole, horizontal, propagation path. This assumption is not always very accurate and thus, in some cases, it is necessary to handle this parameter as a random variable (RV), which follows a specific distribution [35, 52].

#### **3. The Gamma-Gamma wireless optical channel model revisited**

In the G − G channel model for atmospheric turbulence channels the PDF of the irradiance *I* can be derived from the product of two independent Gamma-distributed RVs *x* and *y* with suitably defined parameters [4]. The PDF of the G − G distribution is given by [44]

$$f\_I(I) = \frac{2(km)^{(k+m)/2}}{\Gamma(k)\Gamma(m)\overline{I}} \left(\frac{I}{\overline{I}}\right)^{\frac{k+m}{2}-1} K\_{k-m} \left[2\sqrt{km\left(\frac{I}{\overline{I}}\right)}\right] \tag{9}$$

where *Kα*(·) is the modified Bessel function of the second kind and order *α*, Γ(·) is the gamma function and **E**�*I*� with **E**�·�. denoting expectation. The parameters *k* > 0 and *m* > 0 can be properly adjusted to provide good agreement between theoretical and experimental data. Assuming spherical wave propagation, *k*, *m* can be directly related to atmospheric conditions through the following expressions [14]

$$k = \left[ \exp\left(\frac{0.49\sigma\_{I,R}^2}{\left(1 + 0.18d^2 + 0.56\sigma\_{I,R}^{12/5}\right)^{7/6}}\right) - 1\right]^{-1} \tag{10}$$

$$m = \left[ \exp\left(\frac{0.51\sigma\_{I,R}^2 \left(1 + 0.69\sigma\_{I,R}^{12/5}\right)^{-5/6}}{1 + 0.90d^2 + 0.62d^2\sigma\_{I,R}^{12/5}}\right) - 1\right]^{-1} \tag{11}$$

where *<sup>σ</sup>I*,*<sup>R</sup>* is the Rytov variance defined above and *<sup>d</sup>* <sup>=</sup> <sup>√</sup>*kD*2/4*L*, with *<sup>D</sup>* being the aperture diameter of the receiver.

The corresponding CDF can be expressed as [21]

6 Will-be-set-by-IN-TECH

mathematical expression for this model, which is accurate enough for daytime hours, is given

while for accurate results for night hours are given by the SLC Night model [7, 52], depends only on the height where the optical link operates like the previous model are given by the

8.4 <sup>×</sup> <sup>10</sup>−<sup>15</sup> if 0 <sup>&</sup>lt; *<sup>h</sup>* <sup>&</sup>lt; 18.5m;

*<sup>h</sup>*<sup>2</sup> if 18.5 < *h* < 110m; 2.5 <sup>×</sup> <sup>10</sup>−<sup>16</sup> if 110 <sup>&</sup>lt; *<sup>h</sup>* <sup>&</sup>lt; 1500m;

*<sup>h</sup>*<sup>3</sup> if 1500 < *h* < 7200m;

<sup>−</sup>*h*/*c*<sup>3</sup> + *c*4*e*

*<sup>h</sup>*0.5 if 7200 < *h* < 20000m; .

1.7 <sup>×</sup> <sup>10</sup>−<sup>14</sup> if 0 <sup>&</sup>lt; *<sup>h</sup>* <sup>&</sup>lt; 18.5m;

*<sup>h</sup>*1.05 if 18.5 < *h* < 240m; 1.3 <sup>×</sup> <sup>10</sup>−<sup>15</sup> if 240 <sup>&</sup>lt; *<sup>h</sup>* <sup>&</sup>lt; 880m;

*<sup>h</sup>*<sup>3</sup> if 880 < *h* < 7200m;

*<sup>h</sup>*0.5 if 7200 < *h* < 20000m; .

(6)

(7)

*<sup>n</sup>*, (i.e

(9)

<sup>−</sup>*h*/*c*<sup>5</sup> (8)

*<sup>n</sup>* remains constant for relatively long time

��

as:

following expression:

*C*2 *<sup>n</sup>*(*h*) =

*C*2 *<sup>n</sup>*(*h*) =

mathematical form is the following:

⎧ ⎪⎪⎪⎪⎪⎨

3.13×10−<sup>13</sup>

8.87×10−<sup>7</sup>

2.0×10−<sup>16</sup>

2.87×10−<sup>12</sup>

8.87×10−<sup>7</sup>

2.0×10−<sup>16</sup>

*<sup>n</sup>*(*h*) = *c*<sup>1</sup> + *c*2*e*

**3. The Gamma-Gamma wireless optical channel model revisited**

suitably defined parameters [4]. The PDF of the G − G distribution is given by [44]

� *I I*

The above mentioned models for the estimation of the turbulence parameter *C*<sup>2</sup>

HV5/7, Hufnagel and Stanley, SLC Day and SLC Night), are valid mainly for wireless optical links which have been installed over terrestrial area. Thus, the estimation of the turbulence parameter for paths over maritime area can be done through the following accurate approximation which is valid for low altitude and for specific constant parameters *c*1, *c*2, *c*3, *c*<sup>4</sup> and *c*5. The values of the constant parameters are given in [41, 44] and the corresponding

Obviously, as we mentioned above, for the estimation of the turbulence strength through the

interval and for the whole, horizontal, propagation path. This assumption is not always very accurate and thus, in some cases, it is necessary to handle this parameter as a random variable

In the G − G channel model for atmospheric turbulence channels the PDF of the irradiance *I* can be derived from the product of two independent Gamma-distributed RVs *x* and *y* with

> � *<sup>k</sup>*+*<sup>m</sup>* <sup>2</sup> −1

where *Kα*(·) is the modified Bessel function of the second kind and order *α*, Γ(·) is the gamma function and **E**�*I*� with **E**�·�. denoting expectation. The parameters *k* > 0 and *m* > 0 can be properly adjusted to provide good agreement between theoretical and experimental data.

*Kk*−*<sup>m</sup>*

� 2 � *km* � *I I*

⎪⎪⎪⎪⎪⎩

⎧ ⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎩

*C*2

*fI*(*I*) = <sup>2</sup>(*km*)(*k*+*m*)/2

Γ(*k*)Γ(*m*)*I*

above Equations, we assumed that the value of *C*<sup>2</sup>

(RV), which follows a specific distribution [35, 52].

$$F\_I(I) = \frac{1}{\Gamma(k)\,\Gamma(m)} \, G\_{1,3}^{2,1} \left[ \frac{km}{\overline{I}} \, I \Big|\_{k,m,0}^{-1} \right] \tag{12}$$

with *<sup>G</sup> <sup>m</sup>*,*<sup>n</sup> <sup>p</sup>*,*<sup>q</sup>* [·] being the Meijer-G function [60, Eq. (8.2.1.1)] 1.

The *<sup>ν</sup>*−th moment of *<sup>I</sup>* defined as **<sup>E</sup>**{*Xν*} <sup>=</sup> � <sup>∞</sup> <sup>0</sup> *fI*(*I*)d*I* can be obtained using [28, Eq. (6.561/16)] as

$$\mathbb{E}\{I^{\nu}\} = \left(\frac{k\,m}{\overline{I}}\right)^{-\nu} \frac{\Gamma(k+\nu)\Gamma(m+\nu)}{\Gamma(k)\Gamma(m)}\tag{13}$$

For the correlated test case, to obtain an analytical expression for the multivariate G-G distribution, we first define **N** independent Gamma-distributed RVs *Wn*, *n* = 1, 2, . . . , **N**, with marginal PDFs given by [70, Eq. (2)]

$$f\_{\mathcal{W}\_{\mathbb{R}}}(w\_{\mathbb{N}}) = \frac{k^k w\_{\mathbb{N}}^{k-1}}{\Gamma(k)} \exp\left(-kw\_{\mathbb{N}}\right) H\left(w\_{\mathbb{N}}\right) \tag{14}$$

where *k* ≥ 1/2 is the distribution's shaping parameter and *H* (·) the unit step function [28, p. xliv]. Also, let *Yn*'s be correlated Gamma-distributed RVs with correlation matrix given by **Σ***i*,*<sup>j</sup>* = *ρ*|*i*−*j*<sup>|</sup> , where 0 ≤ *ρ* < 1 is the power correlation coefficient [68, Eq. (9.195)]. Performing **N** RVs transformations in [33, Eq. (3)] and using [28, Eq. (8.445)], the joint PDF of **Y** = [*Y*<sup>1</sup> *Y*<sup>2</sup> ··· *Y***N**] is obtained as

$$f\_{\mathbf{Y}}(\mathbf{y}) = \frac{(1-\rho^2)^m}{\Gamma(m)} \exp\left[ -\frac{m}{\overline{\pi}(1-\rho^2)} \left( y\_1 + y\_\mathbf{N} \right) - \frac{m}{\overline{\pi}} \frac{(1+\rho^2)}{(1-\rho^2)} \sum\_{j=2}^{\mathbf{N}-1} y\_j \right] \tag{15}$$

$$\times \sum\_{i\_1, i\_2, \dots, i\_{\mathbf{N}-1} = 0}^{\infty} \left[ \frac{m}{\overline{\pi} \left( 1-\rho^2 \right)} \right]^{\mathbf{N} \cdot m + 2\sum\_{j=1}^{\mathbf{N}-1} i\_j} y\_1^{q\_1 - 1} y\_{\mathbf{N}}^{q\_2 - 1} \prod\_{j=2}^{\mathbf{N}-1} y\_j^{q\_j - 1} \prod\_{n=1}^{\mathbf{N}-1} \left[ \frac{1}{\Gamma \left( m + i\_n \right) i\_n!} \right] \tag{15}$$

<sup>1</sup> The Meijer-G function is a standard built-in function available in the most popular mathematical software packages, such as Maple or Mathematica

where **y** = [*y*<sup>1</sup> *y*<sup>2</sup> ··· *y***N**], with *yn* > 0, ∀ *n*, *m* ≥ 1/2 is the distribution's shaping parameter, and

$$q\_j = \begin{cases} m + i\_1 & \text{if } j = 1; \\ m + i\_{\mathbf{N} - 1} & \text{if } j = \mathbf{N}; \\ m + i\_{j - 1} + i\_j & \text{if } j = 2, 3, \dots, \mathbf{N} - 1. \end{cases} \tag{16}$$

Then, following the approach presented in [70] and [40], the multivariate G-G distribution with exponential correlation can be derived as [70, Eq. (5)]

$$f\_{\mathbf{I}}(\mathbf{I}) = \int\_0^\infty \int\_0^\infty \cdots \int\_0^\infty \prod\_{n=1}^\mathbf{N} \left[ y\_n^{-1} \, f \mathbf{w}\_n \left( \frac{I\_n}{y\_n} \right) \right] \, f \mathbf{y} \left( \mathbf{y} \right) \, d\mathbf{y} \tag{17}$$

where **I** = [*I*<sup>1</sup> *I*<sup>2</sup> ··· *I***N**], with *In* = *Wn Yn*, ∀ *n*, Substituting (14) and (15) to (17) and using [28, Eq. (3.471/9)], the joint PDF of **I** is obtained as

$$\begin{split} f\_{\mathbf{I}}(\mathbf{I}) &= \frac{2^{\mathbf{N}} \left(1 - \rho^{2}\right)^{\mathbf{m}}}{\left[\Gamma(k)\right]^{\mathbf{N}} \Gamma(m)} \sum\_{i\_{1}, i\_{2}, \dots, i\_{\mathbf{N}-1} = 0}^{\infty} \Xi^{\frac{\mathbf{N}(k+m)}{2} + \sum\_{j=1}^{\mathbf{N}-1} i\_{j}} \rho^{2} \Sigma\_{\boldsymbol{\gamma} - 1}^{\mathbf{N}-1} i\_{j} \\ & \quad \times I\_{1}^{\omega\_{1}} \, \mathbf{K}\_{q\_{1} - k} \left[2\sqrt{\boldsymbol{\varpi} \, \boldsymbol{I}\_{1}}\right] u\_{\mathbf{N}}^{\omega\_{\mathbf{N}}} \, \mathbf{K}\_{q\_{\mathbf{N}} - k} \left[2\sqrt{\boldsymbol{\varpi} \, \boldsymbol{I}\_{\mathbf{N}}}\right] \\ & \quad \times \prod\_{j=2}^{\mathbf{N}-1} \left\{ \frac{I\_{j}^{\omega\_{j}}}{(1 + \rho^{2})^{\frac{q\_{j} - k}{2}}} \, \mathbf{K}\_{q\_{j} - k} \left[2\sqrt{\boldsymbol{\varpi} \, (1 + \rho^{2}) \, \boldsymbol{I}\_{j}}\right] \right\} \prod\_{n=1}^{\mathbf{N}-1} \left[\frac{1}{\Gamma\left(m + i\_{n}\right) \, i\_{n}!} \right] \end{split} \tag{18}$$

with Ξ = (*k m*)/[*I* � <sup>1</sup> <sup>−</sup> *<sup>ρ</sup>*2� ] , *ω<sup>n</sup>* = (*k* + *qn*) /2 − 1, ∀ *n* = 1, 2, . . . , **N**. Note that for *ρ* = 0, all infinite series terms that appear in (18) with *ij* �= 0, ∀ *j*, vanish and (18) simplifies to the product of **N** independent G-G distributions [16, Eq. (2)].

Next, important statistical properties of the exponentially correlated multivariate G-G distribution, namely the joint G-G CDF and MGF will be presented.

#### **3.1. Joint CDF**

The joint CDF of **I** is given by

$$F\_{\mathbf{I}}\left(\mathbf{I}\right) = \int\_{0}^{I\_{1}} \int\_{0}^{I\_{2}} \cdots \int\_{0}^{I\_{N}} f\_{\mathbf{I}}\left(\mathbf{t}\right) \,\mathrm{d}\mathbf{t} \tag{19}$$

with **t** = [*t*<sup>1</sup> *t*<sup>2</sup> ··· *t***N**]. In order to obtain an analytical expression for (19), the Bessel functions that appear in (18) are expressed in terms of Meijer-G functions [60, Eq. (8.4.3.1)]. Then, using [60, Eq. (8.4.23.1)], [60, Eq. (2.24.1.1)], and [60, Eq. (8.2.2.15)], an analytical expression for the joint G-G CDF with exponential correlation is derived in infinite series form as

$$\begin{split} F\_{\mathbf{I}}(\mathbf{I}) &= \frac{\left(1-\rho^{2}\right)^{m}}{\left[\Gamma(k)\right]^{\mathbf{N}}\Gamma(m)} \sum\_{i\_{1},i\_{2},\ldots,i\_{\mathbf{N}-1}=0}^{\infty} \prod\_{n=1}^{\mathbf{N}-1} \left[\frac{\rho^{2i\_{n}}}{\Gamma\left(m+i\_{n}\right)} \, {}\_{i\_{1}}!\right] \mathcal{G}\left(q\_{1},k,\boldsymbol{\varpi}\,I\_{1}\right) \\ & \times \mathcal{G}\left(q\_{\mathbf{N}},k,\boldsymbol{\varpi}\,I\_{\mathbf{N}}\right) \prod\_{j=2}^{\mathbf{N}-1} \left\{\mathcal{G}\left[q\_{j},k,\boldsymbol{\varpi}\left(1+\rho^{2}\right)\, {}\_{j}\right] \, {}\_{i\_{j}}\right] \left(1+\rho^{2}\right)^{-q\_{j}} \end{split} \tag{20}$$

where

$$\mathcal{G}(A,B,\mathbf{x}) = \, \, \, \mathcal{G}^{2,1}\_{1,3} \left[ \mathbf{x} \Big| \, \Big| \, \begin{matrix} 1 \\ A,B,0 \end{matrix} \right]. \tag{21}$$

#### **3.2. Joint MGF**

(16)

(18)

(20)

*f***<sup>Y</sup>** (**y**) d**y** (17)

8 Will-be-set-by-IN-TECH

where **y** = [*y*<sup>1</sup> *y*<sup>2</sup> ··· *y***N**], with *yn* > 0, ∀ *n*, *m* ≥ 1/2 is the distribution's shaping parameter,

Then, following the approach presented in [70] and [40], the multivariate G-G distribution

where **I** = [*I*<sup>1</sup> *I*<sup>2</sup> ··· *I***N**], with *In* = *Wn Yn*, ∀ *n*, Substituting (14) and (15) to (17) and using [28,

all infinite series terms that appear in (18) with *ij* �= 0, ∀ *j*, vanish and (18) simplifies to the

Next, important statistical properties of the exponentially correlated multivariate G-G

� 2 �<sup>Ξ</sup> *<sup>I</sup>***<sup>N</sup>** �

Ξ (1 + *ρ*2) *Ij*

� ⎫ ⎬ ⎭

] , *ω<sup>n</sup>* = (*k* + *qn*) /2 − 1, ∀ *n* = 1, 2, . . . , **N**. Note that for *ρ* = 0,

**N**−1 ∏*n*=1

� 1

*f***<sup>I</sup>** (**t**) d**t** (19)

G (*q*1, *k*, Ξ *I*1)

�−*qj* �

. (21)

Γ (*m* + *in*) *in*!

�

Ξ **N**(*k*+*m*) <sup>2</sup> <sup>+</sup>∑**N**−<sup>1</sup> *<sup>j</sup>*=<sup>1</sup> *ij ρ* 2 ∑**N**−<sup>1</sup> *<sup>j</sup>*=<sup>1</sup> *ij*

*Kqj*−*<sup>k</sup>* � 2 �

**N** ∏*n*=1 � *y*−<sup>1</sup> *<sup>n</sup> fWn*

*<sup>m</sup>* + *ij*−<sup>1</sup> + *ij* if *<sup>j</sup>* = 2, 3, . . . , **<sup>N</sup>** − 1.

� *In yn* ��

*m* + *i*<sup>1</sup> if *j* = 1; *<sup>m</sup>* + *<sup>i</sup>***N**−<sup>1</sup> if *<sup>j</sup>* = **<sup>N</sup>**;

*qj* =

with exponential correlation can be derived as [70, Eq. (5)]

� ∞ 0

� ∞ 0 ··· � ∞ 0

∞ ∑ *i*1,*i*2,...,*i***N**−<sup>1</sup>=0

> *qj*−*k* 2

distribution, namely the joint G-G CDF and MGF will be presented.

*F***<sup>I</sup>** (**I**) =

� *I*<sup>1</sup> 0

joint G-G CDF with exponential correlation is derived in infinite series form as

∞ ∑ *i*1,*i*2,...,*i***N**−<sup>1</sup>=0

> **N**−1 ∏ *j*=2

<sup>G</sup>(*A*, *<sup>B</sup>*, *<sup>x</sup>*) = *<sup>G</sup>* 2,1

� G � *qj*, *k*, Ξ � 1 + *ρ*<sup>2</sup> � *Ij* � � 1 + *ρ*<sup>2</sup>

� *I*<sup>2</sup> 0 ··· � *I***<sup>N</sup>** 0

with **t** = [*t*<sup>1</sup> *t*<sup>2</sup> ··· *t***N**]. In order to obtain an analytical expression for (19), the Bessel functions that appear in (18) are expressed in terms of Meijer-G functions [60, Eq. (8.4.3.1)]. Then, using [60, Eq. (8.4.23.1)], [60, Eq. (2.24.1.1)], and [60, Eq. (8.2.2.15)], an analytical expression for the

> **N**−1 ∏*n*=1

> > 1,3 � *x* � � � 1 *A*, *B*, 0 �

�

*ρ*2*in* Γ (*m* + *in*) *in*! �

*f***<sup>I</sup>** (**I**) =

Eq. (3.471/9)], the joint PDF of **I** is obtained as

<sup>1</sup> <sup>−</sup> *<sup>ρ</sup>*2�*<sup>m</sup>*

**<sup>N</sup>** Γ(*m*)

⎧ ⎨ ⎩ � 2 �<sup>Ξ</sup> *<sup>I</sup>*<sup>1</sup> � *uω***<sup>N</sup> <sup>N</sup>** *Kq***N**−*<sup>k</sup>*

*I ωj j* (1 + *ρ*2)

product of **N** independent G-G distributions [16, Eq. (2)].

*<sup>f</sup>***<sup>I</sup>** (**I**) <sup>=</sup>2**<sup>N</sup>** �

[Γ(*k*)]

× *I ω*<sup>1</sup> <sup>1</sup> *Kq*1−*<sup>k</sup>*

× **N**−1 ∏ *j*=2

> � <sup>1</sup> <sup>−</sup> *<sup>ρ</sup>*2�

with Ξ = (*k m*)/[*I*

**3.1. Joint CDF**

where

The joint CDF of **I** is given by

*F***<sup>I</sup>** (**I**) =

�

[Γ(*k*)]

<sup>1</sup> <sup>−</sup> *<sup>ρ</sup>*2�*<sup>m</sup>*

**<sup>N</sup>** Γ(*m*)

× G (*q***N**, *k*, Ξ *I***N**)

⎧ ⎨ ⎩

and

Using *f***<sup>I</sup>** (**I**), the joint MGF of **I** can be obtained as

$$\mathcal{M}\_{\mathbf{I}}\left(\mathbf{s}\right) = \int\_{0}^{\infty} \int\_{0}^{\infty} \cdots \int\_{0}^{\infty} \exp\left(-\sum\_{n=1}^{\mathbf{N}} s\_{n} I\_{\mathbf{I}}\right) f\_{\mathbf{I}}\left(\mathbf{I}\right) \,\mathrm{d}\mathbf{I} \tag{22}$$

with **s** = [*s*<sup>1</sup> *s*<sup>2</sup> ··· *s***N**]. Substituting (18) in (22) and using [60, Eq. (8.4.3.1)], [60, Eq. (8.4.23.1)], [60, Eq. (2.24.1.1)], and [60, Eq. (8.2.2.15)], an analytical expression for the joint G-G MGF is derived as

$$\begin{split} \mathcal{M}\_{\mathbf{I}}\left(\mathbf{s}\right) &= \frac{\left(1-\rho^{2}\right)^{m}}{\left[\Gamma\left(k\right)\right]^{\mathbf{N}}\Gamma\left(m\right)} \sum\_{\substack{i\_{1},i\_{2},\ldots,i\_{\mathbf{N}-1}=0 \ j=1}}^{\infty} \prod\_{j=1}^{\mathbf{N}-1} \left[\frac{\rho^{2i\_{j}}}{\Gamma\left(m+i\_{j}\right)} \right] \mathcal{H}\left(q\_{1},k\_{\star}\frac{\Xi}{s\_{1}}\right) \\ & \times \mathcal{H}\left(q\_{\mathbf{N}},k\_{\star}\frac{\Xi}{s\_{\mathbf{N}}}\right) \prod\_{j=2}^{\mathbf{N}-1} \left\{\mathcal{H}\left[q\_{j},k\_{\star}\frac{\Xi\left(1+\rho^{2}\right)}{s\_{j}}\right] \left(1+\rho^{2}\right)^{-q\_{j}}\right\} \end{split} \tag{23}$$

where

$$\mathcal{H}(A,B,\mathbf{x}) = \, \, G\_{1,2}^{2,1} \left[ \mathbf{x} \Big| \, \Big|\_{A,B}^{1} \right]. \tag{24}$$

#### **4. System model**

We consider a MIMO FSO system where the information signal is transmitted via *M* apertures and received by *N* apertures. The information bits are modulated using IM/DD with on-off keying (OOK) and transmitted through the *M* apertures using repetition coding [14]. A high-energy FSO system whose performance is limited by background radiation and thermal noise is assumed. Under this assumption, the use of the AWGN model as a good approximation of the Poisson photon counting detection model is applicable [46]. The received signal at the *n*-th receive aperture is expressed as

$$r\_{\text{fl}} = \eta s \sum\_{m=1}^{M} I\_{\text{min}} + v\_{\text{fl}} \cdot n = 1, \dots, N \tag{12}$$

where *η* is the optical-to-electrical conversion coefficient, *s* ∈ {0, 1} represents the information bits and *vn* is the AWGN with zero mean and variance *σ*<sup>2</sup> *<sup>v</sup>* = *N*0/2. In the following analysis, we assume that the *Imn*-s are either independent or exponentially correlated G-G random variables. It is noted that the assumption of independence can significantly simplify the underlying mathematical analysis and it is justified for link distances of the order of kilometers and for aperture separation distances of the order of centimeters [14]. Finally, we define the instantaneous electrical signal-to-noise ratio (SNR) as *μmn* (*ηImn*)2/*N*<sup>0</sup> and the corresponding average electrical SNR as *<sup>μ</sup>mn* (*η***E**�*Imn*�)2/*N*0.

#### **5. Performance analysis of SISO links**

The ABEP of the considered system in the presence of AWGN and under the assumption of perfect channel state information (CSI) at the receiver side is given by

$$P\_{\ell} = P(1)P(e|1) + P(0)P(e|0) \tag{13}$$

where *P*(1) and *P*(0) are the probabilities of sending 1 and 0 bits, respectively, and *P*(*e*|1), *P*(*e*|0) are the conditional bit-error probabilities bit 1 or 0 has been transmitted, respectively. Without loss of generality, we assume *P*(1) = *P*(0) = 0.5 and *P*(*e*|1) = *P*(*e*|0), a fact also justified by the symmetry of the problem. Using the analysis presented in [46] one obtains

$$P(e|I) = P(e|1, I) = P(e|0, I) = Q\left(\frac{\eta I}{\sqrt{2\text{N}\_0}}\right) \tag{14}$$

where *Q*(*x*) = √ 1 2*π* ∞ *<sup>x</sup>* exp − *t* 2 2 d*t* is the Gauss-Q function. The ABEP can be obtaining by averaging *P*(*e*|*I*) over the PDF of *I*, namely

$$\overline{P}\_{b\varepsilon} = \int\_0^\infty P(e|I) f\_I(I) dI \tag{15}$$

By expressing *P*(*e*|*I*) in terms of the electrical SNR *μ*, namely *Q ηI* <sup>√</sup>2*N*<sup>0</sup> = *Q <sup>μ</sup>* 2 = 0.5erfc <sup>√</sup>*<sup>μ</sup>* 2 , *Pbe* can be obtained by performing a change of variables and averaging over the PDF of the electrical SNR, *μ*, instead of the PDF of *I*. Assuming that *I* follows a G − G distribution with parameters *k* and *m* and by applying a simple power transformation of RVs, the PDF of *μ* can be expressed as

$$f\_{\mu}(\mu) = \frac{(km)^{(k+m)/2}}{\Gamma(k)\Gamma(m)\sqrt{\overline{\mu}\mu}} \left(\sqrt{\frac{\mu}{\overline{\mu}}}\right)^{\frac{k+m}{2}-1} K\_{k-m} \left[2\sqrt{km\sqrt{\frac{\mu}{\overline{\mu}}}}\right] \tag{16}$$

Finally *Pbe* can be expressed in closed form as

$$\overline{P}\_{b\varepsilon} = 0.5 \mathcal{F}\left(k, m, \overline{\mu}, \frac{1}{2}\right) \tag{17}$$

where <sup>F</sup>(*k*, *<sup>m</sup>*, *<sup>μ</sup>*,*s*) is given by [57, Eq. (27)] <sup>2</sup>

$$\mathcal{F}(k,m,\overline{\mu},s) = \frac{\Xi^{k+m}s^{-(k+m)/2}}{4\pi^{\frac{3}{2}}\Gamma(k)\Gamma(m)}G\_{2,5}^{4,2}\left[\frac{k^2m^2}{16\overline{\mu}s^2}\,\middle|\,^{a\_p}\_{b\_q}\right] \tag{16}$$

$$\text{with } a\_p = \{1 - \frac{k+m}{4}, \frac{1}{2} - \frac{k+m}{4}\}, b\_q = \{\frac{1}{2} + \frac{k-m}{4}, \frac{k-m}{4}, \frac{1}{2} + \frac{m-k}{4}, \frac{m-k}{4}, -\frac{k+m}{4}\} \text{ and } \Xi \stackrel{\Delta}{=} \sqrt{\frac{km}{\sqrt{\mu}}}.$$

Another important performance metric is the outage probability. The outage probability is defined as the probability that the SNR of the signal at the output of the receiver, *μ* , falls below a specified threshold, *μ*th. This metric is considered as an important parameter for FSO links to be operated as a part of a data network and is critical in the design of both transport and network layer [21]. With the help of (12) and assuming *I* = 1, the outage probability is obtained in closed form as

$$P\_{\rm out} = \Pr\{\mu < \mu\_{\rm th}\} = \Pr\left\{\frac{I^2 \eta^2}{N\_0} < \mu\_{\rm th}\right\} = \Pr\left\{I < \sqrt{\frac{\mu\_{\rm th}}{\overline{\mu}}}\right\} = F\_I\left(\sqrt{\frac{\mu\_{\rm th}}{\overline{\mu}}}\right) \tag{17}$$

<sup>2</sup> It is noted that Eq. (27) in [57] includes typos which we have corrected in (16)

**Figure 1.** Average Bit Error Probability of SISO OOK receivers employing intensity modulation and direct detection and operating over G − G channels, for different values of parameters *k* and *m*, as a function of the Average Electrical SNR, *μ*.

Using (16), in Figure 1 the ABEP of SISO OOK receivers employing intensity modulation and direct detection and operating over G − G channels is depicted, for different values of parameters *k* and *m*, as a function of the Average Electrical SNR, *μ*. Moreover, in Figure 2 the OP of the system under consideration is depicted for the same set of parameters as a function of the inverse normalized outage threshold, *<sup>μ</sup>*th *μ* .

As it can be observed, for a SISO FSO link, both the ABEP and OP performance is quite poor (i.e., higher than 10−<sup>3</sup> in the SNR range of 30-40 dB , especially over strong atmospheric turbulence conditions (that corresponds to small values of *k* or *m*) and therefore the use of diversity techniques is absolutely justified. The use of spatial diversity can be implemented either at the transmitter (MISO) or at the receiver (SIMO case) or at both of them (MIMO case). In both figures, our numerically evaluated results are accompanied with semi-analytical Monte Carlo simulations. In our simulations, more than 10<sup>6</sup> square G <sup>−</sup> G samples are generated to guarantee statistical convergence. The following listings, written in Matlab 2008a demonstrate the evaluation of the ABEP and OP for the SISO case using Monte Carlo simulations.

**function** y = BER\_OOK ( SNR\_dB, k , m, SIZE ) *%SNR\_dB : a v e r a g e e l e c t i c a l SNR %k , m: t h e d i s t r i b u t i o n p a r a m e t e r s* mu\_bar = 1 0^ ( 0. 1 \* SNR\_dB ) ; *% The f i r s t gamma p r o c e s s :* g = gamrnd ( k , **sqrt** (mu\_bar )/k , 1 , SIZE ) ; *% The s e c o n d gamma p r o c e s s :*

10 Will-be-set-by-IN-TECH

where *P*(1) and *P*(0) are the probabilities of sending 1 and 0 bits, respectively, and *P*(*e*|1), *P*(*e*|0) are the conditional bit-error probabilities bit 1 or 0 has been transmitted, respectively. Without loss of generality, we assume *P*(1) = *P*(0) = 0.5 and *P*(*e*|1) = *P*(*e*|0), a fact also justified by the symmetry of the problem. Using the analysis presented in [46] one obtains

> *ηI* <sup>√</sup>2*N*<sup>0</sup>

d*t* is the Gauss-Q function. The ABEP can be obtaining by

*P*(*e*|*I*)*fI*(*I*)d*I* (15)

 *ηI* <sup>√</sup>2*N*<sup>0</sup> = *Q* (14)

(16)

(17)

(16)

(17)

 √ *km μ* .

*<sup>μ</sup>* 2 =

*P*(*e*|*I*) = *P*(*e*|1, *I*) = *P*(*e*|0, *I*) = *Q*

 ∞ 0

, *Pbe* can be obtained by performing a change of variables and averaging over

 *<sup>k</sup>*+*<sup>m</sup>* <sup>2</sup> −1

*Kk*−*<sup>m</sup>*

1 2 

*G* 4,2 2,5

<sup>2</sup> <sup>+</sup> *<sup>m</sup>*−*<sup>k</sup>*

 *k*2*m*<sup>2</sup> 16*μs*<sup>2</sup>

<sup>4</sup> , *<sup>m</sup>*−*<sup>k</sup>*

*μ*th *μ* = *FI*

 *ap bq* 

<sup>4</sup> , <sup>−</sup> *<sup>k</sup>*+*<sup>m</sup>*

<sup>4</sup> } and Ξ

*μ*th *μ* 

 2 *km μ μ* 

the PDF of the electrical SNR, *μ*, instead of the PDF of *I*. Assuming that *I* follows a G − G distribution with parameters *k* and *m* and by applying a simple power transformation of RVs,

> *μ μ*

> > *k*, *m*, *μ*,

<sup>2</sup> Γ(*k*)Γ(*m*)

<sup>4</sup> , *<sup>k</sup>*−*<sup>m</sup>* <sup>4</sup> , <sup>1</sup>

Another important performance metric is the outage probability. The outage probability is defined as the probability that the SNR of the signal at the output of the receiver, *μ* , falls below a specified threshold, *μ*th. This metric is considered as an important parameter for FSO links to be operated as a part of a data network and is critical in the design of both transport and network layer [21]. With the help of (12) and assuming *I* = 1, the outage probability is

> < *μ*th = Pr *I* <

<sup>2</sup> <sup>+</sup> *<sup>k</sup>*−*<sup>m</sup>*

*Pbe* =

*μμ*

*Pbe* = 0.5F

<sup>F</sup>(*k*, *<sup>m</sup>*, *<sup>μ</sup>*,*s*) = <sup>Ξ</sup>*k*<sup>+</sup>*ms*−(*k*+*m*)/2 4*π* <sup>3</sup>

> *I*2*η*<sup>2</sup> *N*0

<sup>2</sup> It is noted that Eq. (27) in [57] includes typos which we have corrected in (16)

<sup>4</sup> }, *bq* <sup>=</sup> { <sup>1</sup>

By expressing *P*(*e*|*I*) in terms of the electrical SNR *μ*, namely *Q*

*<sup>f</sup>μ*(*μ*) = (*km*)(*k*+*m*)/2 Γ(*k*)Γ(*m*)

Finally *Pbe* can be expressed in closed form as

where <sup>F</sup>(*k*, *<sup>m</sup>*, *<sup>μ</sup>*,*s*) is given by [57, Eq. (27)] <sup>2</sup>

<sup>4</sup> , <sup>1</sup>

*P*out = Pr{*μ* < *μ*th} = Pr

<sup>2</sup> <sup>−</sup> *<sup>k</sup>*+*<sup>m</sup>*

where *Q*(*x*) = √

 √*<sup>μ</sup>* 2 

with *ap* <sup>=</sup> {<sup>1</sup> <sup>−</sup> *<sup>k</sup>*+*<sup>m</sup>*

obtained in closed form as

0.5erfc

1 2*π* ∞ *<sup>x</sup>* exp − *t* 2 2 

the PDF of *μ* can be expressed as

averaging *P*(*e*|*I*) over the PDF of *I*, namely

**Figure 2.** Outage Probability of SISO receivers employing intensity modulation and direct detection and operating over G − G channels, for different values of parameters *k* and *m*, as a function of the Inverse Normalized Outage Threshold, *μ*/*μth*.

```
da ta = gamrnd (m, g/m, 1 , SIZE ) ;
% The fi n al data
data = data .^2;
% ABEP o f OOK
y = sum( 0.5 * erfc (0.5 * sqrt ( data ) ) )/ SIZE ;
function y = outage ( k , m, inverse_ th_dB , SIZE )
th = 10^(0.1 * inverse_th_dB );
% f i r s t gamma p r o c e s s
f i r s t = gamrnd ( k , 1/k , 1 , SIZE ) ;
% final process
da ta = ( gamrnd (m, f i r s t /m, 1 , SIZE ) ) ;
y = length ( find ( data < 1/ sqrt ( th ) ) )/ SIZE ;
```
In the following analysis, analytical results for both SIMO and MIMO cases will be presented.

#### **5.1. SIMO case: Selection Diversity (SD)**

The SD combining scheme is the least complicated among the considered combining schemes because of the fact that it processes the aperture with the maximum received irradiance (or electrical SNR). Consequently, the selection is made according to

$$I\_{\rm SD} = \max\{I\_1, I\_2, \dots, I\_N\} \tag{18}$$

If *In* are i.i.d random variables, then using a similar to the previous section analysis, the OP of SD receivers is readily obtained as

$$P\_{\rm out} = \left[ F\_I \left( \sqrt{\frac{\mu\_{\rm th}}{\overline{\mu}}} \right) \right]^N \tag{19}$$

If exponentially correlated irradiance is considered, the CDF of *I*SD is readily obtained using (20) as

$$P\_{\rm out} \left( \gamma\_{\rm th} \right) = F\_{\rm I} \left( \sqrt{\frac{\mu\_{\rm th}}{\overline{\mu}}}, \dots, \sqrt{\frac{\mu\_{\rm th}}{\overline{\mu}}} \right) \,. \tag{20}$$

In Figure 3 the OP of three- and four-branch SIMO receivers employing intensity modulation and direct detection and operating over exponentially correlated G − G channels, for *ρ* = 0.25 and different values of parameters *k* and *m*, is depicted as a function of the Inverse Normalized Outage Threshold, *μ*/*μth*. As it can be observed, the employment of spatial diversity significantly enhances the outage performance of the considered system. To double-check the correctness of the proposed analysis, the analytical results are accompanied with numerical ones obtained using Monte-Carlo simulations. Because of the long computational time, inherent to Monte-Carlo methods, simulation results of up to 10−<sup>6</sup> are given. In order to generate exponentially correlated G − G samples, we make use of the fact that *In* can be written as *In* = *Wn Yn*, ∀ *n*. where *Wn* and *Yn* are samples of an uncorrelated and a correlated gamma process, respectively. Uncorrelated gamma samples with given parameters can be easily generated using the standard Matlab built-in function *gamrnd*(). The generation of correlated gamma random samples with given correlation matrix has been addressed in several past works. In the context of this work, we used the method presented in [67], which yields accurate results for the given system parameters, despite the fact that it imposes certain conditions or constraints on the PDF parameters. Finally, it is noted that although the analytical expression for the outage probability is given in terms of infinite series, it converges rapidly and steadily, requiring few terms to obtain sufficient numerical accuracy. As it was shown in [56], the number of terms depends on the values of the parameters *k* and *m*, the correlation coefficient as well as the SNR. In our terms, a truncation of the infinite series to 25 terms and 18 terms for the three- and four-branch case, respectively, resulted in an excellent match of the numerical results with the Monte-Carlo simulations.

#### **5.2. SIMO case: Optimal combining**

12 Will-be-set-by-IN-TECH

**Figure 2.** Outage Probability of SISO receivers employing intensity modulation and direct detection and operating over G − G channels, for different values of parameters *k* and *m*, as a function of the Inverse

In the following analysis, analytical results for both SIMO and MIMO cases will be presented.

The SD combining scheme is the least complicated among the considered combining schemes because of the fact that it processes the aperture with the maximum received irradiance (or

*I*SD = max{*I*1, *I*<sup>2</sup> ··· , *IN*} (18)

Normalized Outage Threshold, *μ*/*μth*.

*% The fi n al data* data = data .^2; *% ABEP o f OOK*

da ta = gamrnd (m, g/m, 1 , SIZE ) ;

th = 10^(0.1 \* inverse\_th\_dB );

f i r s t = gamrnd ( k , 1/k , 1 , SIZE ) ;

da ta = ( gamrnd (m, f i r s t /m, 1 , SIZE ) ) ; y = **length** ( **find** ( data < 1/ **sqrt** ( th ) ) )/ SIZE ;

**5.1. SIMO case: Selection Diversity (SD)**

*% f i r s t gamma p r o c e s s*

*% final process*

y = **sum**( 0.5 \* **erfc** (0.5 \* **sqrt** ( data ) ) )/ SIZE ;

**function** y = outage ( k , m, inverse\_ th\_dB , SIZE )

electrical SNR). Consequently, the selection is made according to

When receive diversity with optimal combining (OC) is applied, following a similar analysis as in [46], the ABEP is obtained as

$$\overline{P}\_{be} = \int\_0^\infty \int\_0^\infty \dots \int\_0^\infty f\_\mathbf{I}(\mathbf{I}) \, \mathbf{Q} \left( \frac{\eta}{\sqrt{2 \, N \, N\_0}} \sqrt{\sum\_{n=1}^N I\_n^2} \right) \, \mathbf{d} \mathbf{I} \tag{21}$$

where the average irradiance is considered to be normalized to one.

In the following analysis, we assume that *In* are exponentially correlated G-G random variables. In general, the *N*-fold integral in (21) is difficult, if not impossible to be obtained in closed form. To circumvent this problem, we utilize a simple and accurate exponential

**Figure 3.** Outage Probability of triple and quadruple branch SIMO receivers employing intensity modulation and direct detection with selection diversity and operating over exponentially correlated G − G channels, for *ρ* = 0.25 and different values of parameters *k* and *m*, as a function of the Inverse Normalized Outage Threshold, *μ*/*μth*.

approximation for the Gaussian-*Q* function using [22, eq. (2)] and [22, eq. (14)] namely

$$Q(\mathbf{x}) \cong \frac{1}{12} \exp\left(-\frac{\mathbf{x}^2}{2}\right) + \frac{1}{4} \exp\left(-\frac{2\mathbf{x}^2}{3}\right). \tag{22}$$

Substituting (22) in (21) and using [60, Eq. (8.4.3.1)], [60, Eq. (8.4.23.1)], and [60, Eq. (2.24.1.1)], yields

$$\overline{P}\_{b\epsilon} \cong \frac{1}{12} \mathcal{N}\left(k, m, \sqrt{\frac{\mu}{4\,\mathrm{N}}}\right) + \frac{1}{4} \mathcal{N}\left(k, m, \sqrt{\frac{\mu}{3\,\mathrm{N}}}\right) \tag{23}$$

where *μ* is the average electrical SNR and

$$\begin{split} \mathcal{N}\left(k,m,s\right) &= \frac{2^{N\left(m+k-2\right)}\left(1-\rho^{2}\right)^{m}}{\pi^{N}\left[\Gamma(k)\right]^{N}\Gamma(m)} \sum\_{i\_{1},i\_{2},\ldots,i\_{N-1}=0}^{\infty} \left(2\rho\right)^{2} \Sigma\_{i\_{1}-1}^{N-1} i\_{j} \prod\_{n=1}^{N} \left[\frac{1}{\Gamma\left(m+i\_{n}\right)i\_{n}!}\right] \\ & \quad \times \mathcal{Z}\left(q\_{1},k,\Sigma,s\right) \mathcal{Z}\left(q\_{N},k,\Sigma,s\right) \prod\_{j=2}^{N-1} \left\{\mathcal{Z}\left[q\_{j},k,\Sigma\left(1+\rho^{2}\right),s\right] \left(1+\rho^{2}\right)^{-q\_{j}}\right\} \end{split} \tag{24}$$

with

$$\mathcal{Z}(A,B,\mathbb{C},\mathbf{x}) = \mathbf{G}\_{4,1}^{1,4} \left[ \frac{16\,\mathrm{x}^2}{\mathbb{C}^2} \Big| \begin{array}{c} \frac{1-A}{2}, \frac{2-A}{2}, \frac{1-B}{2}, \frac{2-B}{2} \\ 0 \end{array} \right]. \tag{25}$$

Figure 4 illustrates the ABEP performance of OC FSO links, operating at *λ* = 1550nm with *N* = 3 and 4 receive apertures, versus *μ*. Similarly to [70], exponentially correlated atmospheric turbulence channels are considered with *ρ* = 0.175. It is assumed that *C*<sup>2</sup> *<sup>n</sup>* = 1.7 <sup>×</sup> <sup>10</sup>−14, which is a typical value of refractive index for FSO links near the ground during daytime. Furthermore, it is assumed that *D* � *L* leading to *d* = 0, and hence, no aperture averaging is possible. Three different link distances are considered: *L* = 3, 4, and 5 km. The resulting values for *k* and *m* are obtained via (10) and (11), respectively. As expected, the ABEP improves as *N* and/or *μ* increase and/or *L* decreases. Furthermore, as it is shown comparing numerically evaluated results for the ABEP with the equivalent exact ones obtained via Monte Carlo simulations, the proposed ABEP approximation given by (23) provides a tight upper bound for all test cases under consideration.

**Figure 4.** ABEP of triple and quadruple branch SIMO FSO links operating over exponentially correlated G − G atmospheric turbulent channels and employing optimal combining, versus average electrical SNR *μ* for *ρ* = 0.175 and different link distances.

#### **5.3. MIMO case: Equal gain combining**

14 Will-be-set-by-IN-TECH

**Figure 3.** Outage Probability of triple and quadruple branch SIMO receivers employing intensity modulation and direct detection with selection diversity and operating over exponentially correlated G − G channels, for *ρ* = 0.25 and different values of parameters *k* and *m*, as a function of the Inverse

approximation for the Gaussian-*Q* function using [22, eq. (2)] and [22, eq. (14)] namely

 *μ* 4 *N* + 1 4 N *k*, *m*,

∞ ∑ *i*1,*i*2,...,*iN*−<sup>1</sup>=0

> 4,1 16 *x*<sup>2</sup> *C*2 1−*A* <sup>2</sup> , <sup>2</sup>−*<sup>A</sup>* <sup>2</sup> , <sup>1</sup>−*<sup>B</sup>* <sup>2</sup> , <sup>2</sup>−*<sup>B</sup>* 2

*N*−1 ∏ *j*=2

Figure 4 illustrates the ABEP performance of OC FSO links, operating at *λ* = 1550nm with *N* = 3 and 4 receive apertures, versus *μ*. Similarly to [70], exponentially correlated

Substituting (22) in (21) and using [60, Eq. (8.4.3.1)], [60, Eq. (8.4.23.1)], and [60, Eq. (2.24.1.1)],

(2 *ρ*)

 I *qj*, *k*, Ξ

2 ∑*N*−<sup>1</sup> *<sup>j</sup>*=<sup>1</sup> *ij N* ∏*n*=1

> 1 + *ρ*<sup>2</sup> ,*s* 1 + *ρ*<sup>2</sup>

0

 −<sup>2</sup> *<sup>x</sup>*<sup>2</sup> 3 

> *μ* 3 *N*

> > 1

Γ (*m* + *in*)*in*!

. (22)

−*qj* 

. (25)

(23)

(24)

 − *x*2 2 + 1 <sup>4</sup> exp

*<sup>Q</sup>*(*x*) <sup>∼</sup><sup>=</sup> <sup>1</sup>

12 N

<sup>1</sup> <sup>−</sup> *<sup>ρ</sup>*2*<sup>m</sup>*

<sup>I</sup>(*A*, *<sup>B</sup>*, *<sup>C</sup>*, *<sup>x</sup>*) = *<sup>G</sup>* 1,4

*<sup>N</sup>* Γ(*m*)

× I (*q*1, *k*, Ξ,*s*) I (*qN*, *k*, Ξ,*s*)

*Pbe* <sup>∼</sup><sup>=</sup> <sup>1</sup>

where *μ* is the average electrical SNR and

*π<sup>N</sup>* [Γ(*k*)]

<sup>N</sup> (*k*, *<sup>m</sup>*,*s*) <sup>=</sup>2*<sup>N</sup>* (*m*+*k*−2)

<sup>12</sup> exp

 *k*, *m*,

Normalized Outage Threshold, *μ*/*μth*.

yields

with

Assuming perfect Channel State Information (CSI), the ABEP of the considered FSO system is given by [46]

$$\overline{P}\_{be} = \int\_0^\infty \int\_0^\infty \dots \int\_0^\infty f\_{\mathbf{I}}(\mathbf{I}) \, \mathbf{Q} \left( \frac{\eta}{N \, M \sqrt{2 \, N\_0}} \sum\_{n=1}^N I\_n \right) \, \mathbf{d} \mathbf{I} \tag{26}$$

which can be further expressed as [21]

$$\overline{P}\_{be} = \frac{1}{2} \int\_0^\infty f\_I(I) \operatorname{erfc} \left( \frac{\eta}{2 \, MN \sqrt{N\_0}} I \right) \, \mathrm{d}I \tag{27}$$

where *I* = ∑*<sup>N</sup> <sup>n</sup>*=<sup>1</sup> *In*. This integral is very difficult to be obtained in closed form since an analytical expression for the statistical distribution of *I* is required. In the following it will be shown that when *In* are i.i.d G-G random variables, the distribution of *I* can be accurately approximated with the so-called *α* − *μ* distribution [79].
