Design and Performance Analysis of Probabilistically Shaped QAM Signals for Coherent FSO Systems with Gamma-Gamma Turbulence Channels

: Probabilistic shaping (PS) is a promising technique to approach the Shannon limit. In this paper, we design a practical coded modulation scheme based on PS to improve the capacity of coherent free-space optical (FSO) links with quadrature amplitude modulation (QAM), where the fading channel follows the Gamma-Gamma distribution. The aim of this paper is to optimize the probability mass function (PMF) of the QAM signal points to achieve the maximum channel capacity. Due to the complexity of the objective function, the heuristic algorithm was employed to solve the optimization problem. To the best of the authors’ knowledge, the closed-form pairwise error probability (PEP) is ﬁrst derived with the non-uniform signals under the turbulence channel. In addition, we measure the average symbol error rate (SER) and post-FEC bit error rate (BER) by the Monte Carlo simulation method. The numerical simulation results of both capacity and BER show that the proposed PS scheme is better than the uniform distribution. The post-FEC BER results show that the proposed PS scheme provides signiﬁcant gains compared with the uniform scheme.


Introduction
Free-space optical (FSO) communication has attracted more and more attention due to its large capacity, high security, high data rates and license-free operation [1][2][3]. Nevertheless, the demand for data transmission has shown an explosive growth in the past few years, as a result of the emergence of big data-related technologies. In order to confront this trend, the FSO communication system needs to continuously increase its capacity to provide higher performance. Compared with the intensity modulation and direct detection (IM/DD) systems, the coherent free space optical communication provides impressive performance enhancements, for instance, the quadrature amplitude modulation (QAM), even though their implementation is more complex.
There are several technologies that can significantly improve the channel capacity in an FSO system, for instance, high-level modulation formats [4], wavelength division multiplexing [5], space division multiplexing [6], and angular momentum multiplexing [7]. However, the gap between the capacity of an FSO system and Shannon limit still exists. Recently, probabilistic shaping (PS) has been widely investigated as a means to approach the Shannon limit. The emergence of PS provides enhanced flexibility for FSO systems, without increasing their complexity. To be specific, PS changes the probabilities of the points based on off-the-shelf constellations, hence, incurring no additional complexity in system design and implementation [8][9][10].
With the deepening of research on PS, researchers mainly focus on four aspects. Firstly, under different channels, the input distributions are designed to be suitable for various conditions [11,12]. Yao et al. [11] optimized the input distributions by the modified alternating optimized method under the Rayleigh fading channel. In [12], Wang et al. proposed two succinct PS schemes over 16QAM FSO systems, which are based on mapping rules. Secondly, various distribution matching schemes applied to PS have been presented [13][14][15][16]. For instance, constant composition distribution matching (CCDM) [13], hierarchical distribution matching (HiDM) [14], multiset-partition distribution matching (MPDM) [15] and multi-level distribution matching (MLDM) [16]. Thirdly, the design of forward error correction (FEC) codes suitable for PS schemes were proposed [17,18]. Meanwhile, Sheikh et al. [17] proposed a coded modulation scheme based on binary staircase codes as the FEC code, and achieved 0.56-1.55 dB gains of generalized mutual information (GMI) using only a single staircase code. In [18], Git et al. designed protograph-based low-density parity-check (LDPC) codes for PS techniques.The tailored protograph-based LDPC codes have a significant improvement over the standard schemes. Lastly, the structure of the PS implementation was discussed [19]. In [19], Böcherer et al. proposed the probabilistic amplitude shaping (PAS) scheme, which concatenates the distribution matching and FEC code, without distorting the shaped distribution and has no additional error bursts. In this paper, we mainly focus on designing the input distributions suitable for FSO system to improve their capacity.
To the best of the authors' knowledge, there are two main kinds of methods to find the optimal distribution that maximizes capacity. First, the researchers choose the objective function (the mutual information or the achievable information rate) and optimize the input distributions by convex optimization [11,20]. The mutual information is discussed as the objective function in [11] and the alternating optimization is developed. In order to achieve the capacity and capacity-achieving distribution by the interior-point algorithm, the achievable information rate under the FSO channel is utilized in [20]. Secondly, the optimization algorithms based on machine learning are proposed [21,22]. Fallahpour et al. adopted a known sequence, solving the optimization problem based on error counting at the receiver [21]. While Stark et al. introduced a machine learning method to solve the optimization problem of PS in [22]. The Gumbel-Softmax trick is used to train a sampling mechanism for symbols drawn from the finite set.
Motivated by the above methods, we propose a coded modulation scheme based on the probabilistic shaped QAM signal to achieve the capacity of the FSO channel. We choose the heuristic algorithm to optimize the probability mass function (PMF) of the discrete constellations for maximizing the achievable information rate (AIR) of bit-metric decoding (BMD) in a coherent FSO system, which is modeled as a Gamma-Gamma channel for the first time. Considering the reliable communication, the transmission rate should be less than the AIR of the coded modulation scheme and the AIR depends on the input distribution of the QAM signal, hence the input distribution should be optimized to maximize the AIR. Furthermore, to measure the performance of the PS scheme under a Gamma-Gamma channel, we also simulate the average symbol error rate (SER) and the bit error rate after the FEC decoder (post-FEC BER). Different from [11], we consider the PS scheme under the Gamma-Gamma channel instead of the Rayleigh fading channel. Moreover, in [20], Elzanaty et al. also discussed the capacity of the PS scheme for an FSO channel, which is modeled as a Gamma-Gamma channel. They adopted the pulse amplitude modulation (PAM) signals under the assumption of IM/DD. Distinguishing this with [20], we consider the optimal input distribution of the capacity under the coherent FSO system.
The main contributions of this manuscript are twofold. First, under the FSO system, the heuristic algorithm is proposed to optimize the PMF of the constellation signal points and maximize the AIR. Secondly, we theoretically analyze the error performance of the PS scheme for the turbulence fading channel. We derive a closed-form pairwise error probability (PEP) expression for the non-uniform signal sets under the Gamma-Gamma channel for the first time.
The structure of this paper is as follows. In Section 2, the system model is described. We present the achievable information rate of bit-metric decoding and optimize the input distributions of the turbulence fading channel in Section 3. The error performance analysis is provided in Section 4. The simulation results are demonstrated in Section 5, and Section 6 gives the conclusion.
Notations: Random variables are denoted by capital italic letters and their realizations with small italic letters. Scalars and vectors are distinguished using normal and bold fonts. E[·] denotes expectation. The entropy of random variables is denoted by H[·].

System Structure
The FSO communication system based on a rectangular QAM with probabilistic shaping is adopted in this paper, as shown in Figure 1. In the transmitter section, the uniform bitstream is reshaped into the desired distribution by the PS encoder, which will be illustrated in detail in Section 5. After that, the shaped bits encoded by the LDPC codes are mapped to QAM symbols. A Mach-Zehnder modulator (MZM) is used to modulate the laser. The in-phase and quadrature branches of the MZM are driven by the QAM electric signals. The modulated optical signal is amplified by an erbium-doped fiber amplifier (EDFA) and transmitted in the turbulence channel. To simulate the turbulence channel, the Gamma-Gamma turbulence model is adopted in our work for emulating the turbulence conditions from weak to strong. Its probability distribution function (PDF) is given by [23]: where Γ(•) denotes the Gamma function and K n (•) is the modified Bessel function of the second kind of order n. The parameters α and β represent the effective numbers of large scale and small scale turbulence eddies, respectively, which can be expressed as a function of Rytov variance σ 2 R . Considering plane waves, we have [23]: In the coherent receiver, the received optical signals are mixed with a local oscillator. The output signals are detected by two balanced photodetectors (BPD). The corresponding electrical in-phase and quadrature signals are further processed by high-speed digital signal processing (DSP). Eventually, the bit streams are obtained by an LPDC decoder and a PS decoder. The coherent detection currents can be expressed as [2]: where η represents the responsivity of the photodiode and P L is the optical power of the local oscillator beams. P S stands for the received optical power in the present turbulence. θ sig denotes the modulation phase of the received signal. θ IF is the phase difference between the local oscillator and transmitter laser. θ h represents the phase fluctuation caused by atmospheric turbulence. n I and n Q are the compound additive white Gaussian noise. Perfect channel state information (CSI) is assumed at the receiver and the system is not affected by intersymbol interference. Considering the equivalent baseband model, we define the input X to be complex QAM symbols after a bit mapper block. The value of the non-equiprobable QAM symbols are taken from where P X (x i ) stands for the probability when the signal x i is transmitted. The corresponding electrical signal Y before demapper block is given by: where N is the complex zero mean Gaussian noise with variance σ 2 and G is intensity attenuation related to atmospheric turbulence.

Problem Formulation and Solution
In this section, we aim to search for the optimal PMF of the discrete constellations under a Gamma-Gamma channel. We compute the AIR of the proposed scheme with bit-metric decoding, for various symbol to noise ratios (SNRs).
In the following, we review the basic principle of probabilistic shaping. We mainly focus on the AIR of the BMD, rather than symbol-metric decoding (SMD) due to the higher computational complexity of SMD. In this scheme, the input signal X is composed of m bit levels B = B 1 , . . . , B m . The BMD rate is defined as follows [24]: where [•] + means the largest number of • and 0, namely max(•, 0). H(B i | Y, G) calculates the uncertainty between B i at i-th position and the channel output Y under the turbulence condition. These log-likelihood ratios (LLRs) are computed with the auxiliary channel as: where X 0 i and X 1 i denote the sets of PS-QAM symbols whose i-th bit is 0 and 1. In an auxiliary channel, the condition probability q Y|X (Y | X, g) in Equation (7) can be given by: The AIR depends on the distribution of the input signal X, for various SNRs. Hence, the objective function at the specified SNR can be written as: where • represents the summation formula. The value of PMF is a non-negative number less than one with a sum of one. As mentioned above, the objective function R BMD is a lack of the closed-form expression. Furthermore, we have to calculate the BMD rate by Monte Carlo simulations. Thus, finding an efficient algorithm that finds the optimal PMF is difficult. Under these circumstances, we adopt the heuristic algorithm to maximize the BMD rate to derive the optimal PMF of PS-QAM signals.
In the following, we will introduce the heuristic algorithm. The first stage of the algorithm is initializing the population. We produce individuals with the size of P pop . In this algorithm, each individual contains two vectors, namely position and velocity. Considering the symmetric property of the 2-dim 16-ary QAM constellation, we can initialize the population with 2-dim searching space instead of 16-dim. That is, each position represents the 2dim probability p 2 dim , which can correspond to the PMF of 16QAM. Generalized p 2 dim = [p 1 p 2 ], then P X can be expressed as P X = p 1 p 2 p 1 p 2 × p 1 p 2 p 1 p 2 , where × denotes the matrix product. We take p 2 dim = [0.25, 0.25] as an example, then P X = [0.625, . . . , 0.625]. The velocity is also a 2-dim vector formed by the change of each position p 2 dim in the previous and subsequent iterations. Besides, we can initialize the position of the population individuals near the Gaussian distribution in order to achieve faster convergence. We set f lag = 0 as the identification of the iteration. Conversely, the iteration terminates when f lag = 1.
For each SNR, there is an optimal PMF to maximize the BMD rate R BMD . In this situation, the position needs to be measured with the fitness function in each iteration at the specified SNR. This means that the fitness function can be used to evaluate the performance of each individual. Obviously, individuals who can make the BMD rate higher should have a larger fitness value and greater possibility to be retained in the subsequent iteration process. Considering the Monte Carlo simulation of N samples, the fitness function n of the n-th iteration can be calculated by 1 The individual optimum fitness and population optimum fitness is determined in each individual fitness value. Specifically, if the fitness value of the new p 2 dim is better than the historical individual optimal fitness, the p 2 dim is updated to the latest position; if there is a individual in the entire population whose position is better than the global optimal position, the vector is updated to the global optimal position.
Since the iterative process leads to a stochastic manipulation of velocities and positions according to the best experiences of the swarm to search for the global optimum in solution space, it is of vital importance to update the position and velocity, which will directly affect the convergence speed of this heuristic algorithm. Hence, the stochastic inertia weight is employed in this algorithm, which can avoid the lack of local search ability in the early stage of iteration. It can also avoid the lack of global search ability in the late iteration, and its convergence speed and global convergence are significantly improved compared with the normal inertia weight. We adjust the weight through normal random variables to achieve the best balance between the global and local searches.
We introduce the mutation δ in the population so that it may escape from a local optimum. In this algorithm, we choose the number of iteration G as a sign of whether to terminate the algorithm. The detailed algorithm is shown in Table 1.  Initialize: 2 Produce initial population with P pop individuals. 3 Initialize the position and velocity of the individuals. 4 Set the iteration identification f lag = 0. 5 Set mutation probability δ = 0.85. 6 Set the genetic iteration number G = 0. 7 while ( f lag == 0). 8 Reproduction:. 9 Update the velocity and position. 10 Calculate individual fitness function n . 11 Determine the individual optimum fitness and the population optimum fitness. 12 Mutation: 13 if rand > δ 14 Update the velocity and position. 15 if G = 100 update f lag = 1. 16 Output: return the optimum value p 2 dim .

Error Performance Analysis
After depicting the optimization in Section 3, the average SER under the desired PMF, which is optimized by the heuristic algorithm, needs to be discussed. Therefore, in this section, we will evaluate the error performance of the PS scheme under the Gamma-Gamma fading channel. For a BICM system, PEP is applied to derive the union bound in SER. The upper bound of the average SER can be expressed as [25]: where PEP P o ij represents the probability of choosing x j as a transmitted signal, when in fact x i was transmitted. Under the assumption of the perfect CSI available at the receiver, the PEP concerning turbulence attenuation, is given by: where • denotes the norm operation of a vector.
3 , we obtained: For the convenience of representation, we define three variables, namely ζ = 1 2 ln . The PDF of the intensity attenuation H is given in Equation (13) by means of rewriting the K α−β (2 αβg).

Simulation
In this section, we design a practical LDPC coded modulation system based on the PS scheme. The system block diagram is shown in Figure 2. At the transmitter, in order to generate the 16-ary PS-QAM symbols, the pseudo-random binary sequence (PRBS) is generated. Therefore, we can see that the uniform distribution binary sequence is divided into two parts, one of which is fed to the CCDM. The CCDM converts the uniformly distributed bits into the desired distributed 16-ary QAM symbols. After binary converting, they are encoded together with the remaining binary bits. We utilize the DVB-S2 LDPC code, whose code rate is 2 3 . At the receiver, the LLRs are calculated for bit-metric decoding. We then evaluate the achievable information rate of the signal. We also measure the average SER and the post-FEC BER.

Simulation Results
In the following, we will discuss the simulation results of the proposed system. It has been shown that for each SNR, there is an optimum PMF, which maximizes the R BMD over the Gamma-Gamma channel. Hence, the distribution of constellation points with different SNRs and the uniform distributions are illustrated in Figure 3. The specific probabilities are listed in Table 2. When the probability distribution deviates from the uniform distribution, the entropy decreases. Moreover, it can be seen that in the case of lower SNR, the probabilities of the constellation points near the origin are higher than that far away from the origin. Figure 3 illustrates that the PMF begins as a Gaussian-like shape; with an increase of the SNR, PMF eventually tends to a uniform distribution. This is because when the SNR is lower, more probabilistic shaping is required. Specifically, there is a great difference between the constellation points with high probabilities and those with low probabilities.  Theoretically, each SNR has a corresponding optimal PMF. The effect of optimizing the PMF at a corresponding SNR is shown in Figure 4. The gap to the Shannon limit is closed by the optimized PS. We compare the performances of the optimized solutions and uniform distributions under different turbulence conditions. Under week turbulence, the corresponding optimized PS scheme is closer to the Shannon limit. With the increase of the SNR, the shaping gain between the optimized PS and the uniform distribution decreases gradually. Eventually, the PMF approaches the uniform distribution. These phenomena are consistent with the tendency of the constellation with the input distribution, shown in Figure 3. The PMF tends to a uniform distribution, as the effect of the noise is decreasing. Hence, the shaping gain obtained by the optimized PS becomes smaller.  The result of applying PS with a fixed entropy is shown in Figure 5. It can be seen that the PS with H = 3.7964 (resp. PS with H = 3.9451) outperforms the uniform distribution by 0.4 dB (resp. 0.3 dB) at the achievable information rate of about 1.8 bit/symbol. Furthermore, the PS with H = 3.9451 is superior to others at the medium SNR. When the E s/ N 0 is about 8 dB, the performance of the PS with H = 3.9451 is better than the PS with H = 3.7964. Besides, after 16 dB, the uniform distribution has more outstanding performance than others. In addition, we compare the performance of the Maxwell-Boltzmann (MB) distribution with fixed entropy. As we can see, when the E s/ N 0 is less than 6 dB , there is little difference between the PS scheme and MB distribution, however when E s/ N 0 is greater than 6 dB , it can be seen that the PS scheme performance is gradually better than MB distribution. As we can see, the PS with H = 3.7964 is worst at a higher SNR. This is because its PMF has a marked difference with the uniform distribution. To make a tradeoff between the complexity and performance, we can use only three PMFs under the unique SNR with no need to adjust the input distribution under the conditions that a small reduced shaping gain can be accepted.  As discussed in Equation (14), the input distribution has a significant impact on the error performance. The average SER upper bound from Equation (10) and from the Monte Carlo simulation are shown in Figure 6. As the entropies of the probability distributions increase, the performance of the SER gradually deteriorates. This is because increasing entropies increase the probabilities of the outermost constellation points, which are dominant sources of errors. Moreover, the gain between the two PS schemes and the uniform distribution increases with the increase of the SER. As we can see, the PS scheme with H = 3.7964 achieves about 1.5 dB gains over the uniform distribution at an SER of 4 × 10 −1 . It is shown that the theoretical curves and simulation curves are basically consistent, because PEP is calculated to approximate the average SER. This result shows that the PS schemes outperform the uniform distribution with the performance of the average SER. 6  To further investigate the performance of PS schemes, we also evaluate the post-FEC BER under the turbulence channel in Figure 7, which is consistent with the trend of the average SER. Compared with the average SER, post-FEC BER falls sharply when its value is about 4 × 10 −2 . The post-FEC BER would converge asymptotically as the BERs approach zero. As we can see, the PS with H = 3.7964 firstly approaches the falling region. It can also be concluded that the smaller the entropy of the PS, the earlier the curves converge. At the post-FEC BER of 1 × 10 −2 , the proposed PS scheme with H = 3.7964 (resp. H = 3.9451) obtains nearly a 1.3 dB (resp. 0.5 dB) gain over the uniform QAM. The reason is that the PS with H = 3.7964 increases the probability of the constellation points with better BER performance to a greater extent, and reduces that with poorer BER performance.

Conclusions
In this paper, the PS technique is investigated for a coherent FSO communication system. We aim to optimize the PMF of the input QAM signal over the Gamma-Gamma channel. To begin with, we first present the BMD rate R BMD . After that, the optimal PMF was obtained by the heuristic algorithm, in order to maximize the R BMD . Correspondingly, the achievable information rate is simulated by a Monte Carlo simulation. Compared with the uniform distribution, PS has a larger gain in the low SNR than that in the high SNR. The results show that in the low SNR regime, the PS with H = 3.7964 could offer the best performance, which has an approximate 0.4 dB gain, while the optimal PMF tends to be a uniform distribution in the high SNR regime. In addition, it can be seen from the results that in each case of the SNR, the desired optimal PMFs are different. In addition, we derived a closed-form expression of the PEP for a non-uniform signal under the condition of the Gamma-Gamma channel. Finally, we evaluate the performance of the average SER and post-FEC BER. An approximately 1.3 dB shaping gain was achieved by PS with H = 3.7964, which proves the effectiveness of the proposed scheme in a coherent FSO communication system. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The study did not report any data.