Adaptive Coded Modulation for IM/DD Free-Space Optical Backhauling: A Probabilistic Shaping Approach

In this paper, we propose a practical adaptive coding modulation scheme to approach the capacity of free-space optical (FSO) channels with intensity modulation/direct detection based on probabilistic shaping. The encoder efficiently adapts the transmission rate to the signal-to-noise ratio, accounting for the fading induced by the atmospheric turbulence. The transponder can support an arbitrarily large number of transmission modes using a low complexity channel encoder with a small set of supported rates. Hence, it can provide a solution for FSO backhauling in terrestrial and satellite communication systems to achieve higher spectral efficiency. We propose two algorithms to determine the capacity and capacity-achieving distribution of the scheme with unipolar M-ary pulse amplitude modulation (M-PAM) signaling. Then, the signal constellation is probabilistically shaped according to the optimal distribution, and the shaped signal is channel encoded by an efficient binary forward error correction scheme. Extensive numerical results and simulations are provided to evaluate the performance. The proposed scheme yields a rate close to the tightest lower bound on the capacity of FSO channels. For instance, the coded modulator operates within 0.2 dB from the M-PAM capacity, and it outperforms uniform signaling with more than 1.7 dB, at a transmission rate of 3 bits per channel use.

In fact, the optimal capacity-achieving distribution should maximize mutual information, rather than the source entropy. For the implementation of the scheme in [21], multilevel coding (MLC) with multistage decoding (MSD) is considered. In this scheme, the encoder necessitates MLC with multiple FEC encoders to generate the probabilistically shaped symbols with the desired distribution. At the receiver, MSD is required, leading to error propagation and long latency due to the successive decoding of bit levels [22]. 1 In [23], a sub-optimal coded modulation scheme for IM/DD channels is proposed with unipolar M -PAM. In this approach, only the even-indexed symbols are freely probabilistically shaped, while the probability of the odd-indexed symbols is forced to equal the probability of the preceding even symbol. However, it is a sub-optimal approach, as the input distribution can not be fully optimized to match the capacity-achieving distribution of the channel, leading to an increased shaping gap and a rate loss. To the best of the authors' knowledge, efficient and practical adaptive coded modulation schemes with fine granularity and their capacity-achieving input distributions have not been well investigated for IM/DD FSO channels.
In fiber-optical communications, PS has recently gained increased interest, followed by the introduction of the probabilistic amplitude shaping (PAS) scheme to approach the capacity of fiber-optical channels in [17]. For the PAS architecture, the capacity-achieving distribution should be symmetric around zero. In this case, the uniformly distributed parity bits from the FEC encoder can modulate the sign of the symbols. Therefore, it is only suitable for bipolar input signals [17]- [19], [24]- [28]. Unfortunately, the PAS can not be directly extended to IM/DD in FSO channels, as the constellation symbols are constrained to be non-negative, i.e., unipolar signaling.

B. Contributions
In this work, we propose an adaptive coded modulation scheme with fine granularity to approach the capacity of IM/DD FSO channels through M -PAM signaling with non-negativity and optical power constraints. Realistic models for the atmospheric turbulence are considered, i.e., Gamma-Gamma and Lognormal distributions [29], [30]. The proposed encoder considers probabilistic shaping of a unipolar M -PAM constellation with a constant composition distribution matcher (CCDM), followed by an efficient FEC encoder. Therefore, the information symbols can be probabilistically shaped, while the parity check bits, generated by the FEC encoder, are uniformly distributed. In the decoder, the FEC decoding is performed before the distribution dematching (i.e., reverse concatenation architecture). We compute the capacity of the proposed scheme (i.e., the maximum achievable rate) with both optimal symbol-metric decoders (SMDs) and practical low complexity bit-metric decoders (BMDs), for various SNRs. The ergodic capacity of the coding scheme is calculated for Gamma-Gamma and Lognormal induced-turbulence fading. The encoder can operate when the channel state information (CSI) is known at the encoder and decoder or only at the decoder. In this regard, we derive the outage probability due to the turbulence-induced fading, when the CSI is not available at the transmitter. The contributions of the manuscript can be summarized as follows.
• We propose a coded modulation scheme to achieve the capacity of FSO channels with probabilistically shaped unipolar M -PAM signaling.
• An algorithm is provided to compute the input distribution that achieves the capacity of the coded modulation scheme, utilizing SMDs and FECs with optimized code rate.
• For practical channel encoders with a finite set of coding rates and BMDs, an algorithm is proposed to estimate the distribution which approaches the capacity of the proposed scheme.
• The ergodic capacity of the proposed coding scheme is analyzed for different SNRs and turbulence conditions.
• The outage probability is derived for the blind encoder, i.e., the CSI is not known at the transmitter.
• An approach to design the blind encoder with an arbitrary outage probability is proposed.
The proposed transponder has several distinguishing features. First, the IM/DD with M -PAM signaling has lower computational complexity and it can be efficiently implemented compared to coherent modulation [2]. Moreover, the architecture considers a CCDM, which is asymptotically optimal in the frame length, and a single binary FEC encoder and decoder at the transmitter and receiver, respectively. Also, the rate can be adapted to the SNR for various turbulence conditions with fine granularity. Finally, the encoder can operate within a pre-designed maximum outage probability in the absence of the CSI at the transmitter. These features for the proposed encoder allow it to become a potential candidate for FSO long-distance high-speed links in both terrestrial and satellite systems [31], [32].

C. Paper Organization and Notations
The work is organized as follows. In Section II, the signal model is described, and the capacity of M -PAM signals is calculated. In Section III, the proposed scheme is introduced, the capacity of the sparse-dense signaling is obtained for SMD and BMD. Also, the optimal operating rate using the proposed encoder is analyzed. In Section IV, an algorithm to estimate the capacity of the coding scheme is provided, assuming optimal SMD with a fully-controllable channel coding rate. For channel encoders with a finite set of possible rates, the capacity and the computational complexity of the proposed scheme are analyzed in SectionV. The outage probability is derived in Section VI, while numerical results are provided in Section VII to attest the scheme performance.
Finally, we conclude the work in Section VIII.
Throughout this paper, we denote random variables (r.v.s) with capital letters and their realizations with small letters. We use P{X = x} to denote the probability that a discrete r.v. X equals x. The expression E{X} denotes the expected value of the r.v. X. For vectors, we use bold letters, e.g., x = [x 1 , x 2 , · · · , x n ]. All logarithmic functions used throughout the paper are of base 2.

II. SIGNAL MODEL
Let us consider n transmissions (channel uses) over a discrete-time FSO channel, also known as optical direct-detection channel with Gaussian post-detection noise [4], [21], [33]. In this channel, the input signal modulates the light intensity, while a photo-detector at the receiver produces a noisy signal that is proportional to the intensity. The dominant noise sources are thermal noise, intensity fluctuation noise by the laser source, and shot noise induced by ambient light. The contributions from all the noise sources can be modeled as AWGN [4]. Hence, the received signal at time instant i can be written as where X i is the channel input, W i is a Gaussian noise with zero mean and variance σ 2 , and G is a r.v. representing the fading due to the atmospheric turbulence. In FSO systems, the channel changes slowly with respect to the bit rate. Hence, G is considered as a block fading process, and it is assumed to be fixed over the entire frame of n symbols. For moderate and strong turbulence, the irradiance fluctuations can be modeled as a Gamma-Gamma distribution [29].
The probability density function (PDF) of the Gamma-Gamma r.v. is where K a (·) is the modified Bessel function of the second kind of order a, and the parameters α and β are the effective number of small scale and large scale cells of the scattering environment, which can be expressed as a function of the Rytov variance, σ 2 R [29]. For example, considering plane waves from [29, (14)- (19)], we have A common format for the channel input X i that can achieve a near capacity performance for IM/DD systems is unipolar M -PAM [3], [21]. Considering that the r.v.s {X i } n i=1 are independent, identically distributed (i.i.d.), the probability mass function (PMF) of X ∈ {a 0 , a 1 , · · · , a M −1 } can be written as p [p 0 , p 1 , · · · , p M −1 ], where a j is the jth element of a [0, ∆, · · · , (M − 1)∆], ∆ > 0 is the spacing between the symbols, and p j is the probability assigned to the constellation symbol a j . Let us define the set that includes all possible symbol distributions as The signal is also subject to an average optical power constraint, i.e., where P is the average optical power limit. In fact, the optical power constraint is more relevant for FSO signals than the electrical power limit (i.e., E {X 2 } ≤ P ), widely adopted in RF systems [4], [33]. Also, the instantaneous optical SNR, defined as gP/σ, is usually adopted in FSO communication, rather than the electrical SNR [3]- [5], [33].
In order to guarantee reliable communication with an arbitrarily low probability of error, the transmission rate should be less than the achievable rate of the coded modulation scheme. 2 The achievable rate depends on the distribution of the signal at the channel input. Hence, the input distribution should be optimized to maximize the achievable rate. In this regard, the capacity of unipolar M -PAM, C(g), can be found for a given g as the optimal value of the following optimization problem maximize ∆>0, p∈S is the mutual information between X and Y given are differential conditional entropy functions, and is the distribution of the received signal conditioned on the channel gain. 3 For notation simplicity, let I ∆ (p|g) I(X; Y |G = g) to emphasize its dependence on the parameters of the input distribution, i.e., p and ∆. For a fixed ∆, the problem (6) is a convex optimization problem in p [34, Theorem 2.7.4]. In fact, I ∆ (p|g) is a concave function in p from the concavity of the conditional entropy function in (7) and the composition with an affine mapping property [35, Sec. 3.2.2]. Hence, it can be efficiently solved using any suitable convex optimization algorithm, e.g., interior-point method, and the optimal p at this ∆, p ∆ , can be computed. Now, it is required to find the optimal value of the constellation spacing that maximizes the mutual information.
Since I ∆ (p ∆ |g) is unimodal in ∆, any efficient optimization algorithm in one dimension can be employed to compute the optimal constellation spacing, ∆ * , e.g., golden section search and successive parabolic interpolation [17].
In order to achieve the maximum capacity, the channel input X should be probabilistically shaped to have a distribution with the optimal parameters obtained from (6). Unfortunately, there  is no known computationally efficient and practical scheme in the literature which achieves this rate for unipolar M -PAM signals. For example, the efficient PAS encoder, usually adopted for coherent fiber optical communications, can not be used due to the asymmetric signaling around zero in unipolar M -PAM [17].

III. PROPOSED CODED MODULATION SCHEME WITH M -PAM FOR FSO
We propose two practical adaptive coded modulation schemes to increase the bandwidth efficiency of FSO communications by probabilistically shaping the input distribution. The first method considers that the CSI is available at both the encoder and decoder, as shown in Fig. 1.
For the second scheme, the CSI is assumed to be known only at the receiver (blind encoder). In this section, we start by describing the encoder, and the maximum mutual information between the unipolar M -PAM signal at the channel input and the received noisy signal. Then, the decoder, its achievable rate, and the optimal feasible operating rate are discussed.

A. Sparse-Dense Encoder
The target of the encoder is to convert the input binary string of uniformly distributed bits into probabilistically shaped unipolar M -PAM symbols and to perform channel coding. Hence, the decoder can reliably recover the original data from noisy measurements at the receiver. The proposed scheme is described in detail below.
1) Capacity-achieving Distribution: First, we compute the capacity-achieving distribution of the proposed scheme, parameterized with p * and ∆ * , which permits the highest reliable communication rate for a given SNR. Section V proposes an algorithm to obtain the optimal distribution when CSI is known at the encoder, or when only channel statistics are available, as in Section VI.
2) Distribution Matching: The distribution matching transforms the uniformly distributed input bit string, u ∈ {0, 1} kp , into unipolar M -PAM symbols, x p ∈ {0, ∆ * , · · · , (M − 1)∆ * } np with the target distribution, p * . Several distribution matching (DM) techniques have been proposed in the literature with various computational complexity, rate loss, and parallelization ability [36]- [40]. The CCDM is an invertible mapping from a fixed-length vector of uniformly distributed bits to a fixed-length sequence of shaped symbols (i.e., amplitudes) [36]. The empirical distributions of all possible output sequences are identical, i.e., they have a constant composition.
Therefore, every output sequence follows, to some extent, the target distribution. The target distribution should be quantized such that the probability of each symbol can be represented as a rational number, where the denominator is the frame length length, n p . In other words, the PMF of X p , p * , is approximated by what is called n p -type distribution in the form of where z j is an integer representing the number of times at which the jth symbol appears and The discrepancy between the target and n p -type distributions decreases with the output sequence length. Hence, for asymptotically large n p , the quantization error for p * is negligible, i.e., lim np→∞p = p * .
In order to quantify the CCDM rate, R DM k p /n p , the number of input bits should be computed. The number of bits, k p , that are required to be transformed to n p shaped symbols depends on the number of possible configurations (i.e., permutations) of the output symbols that have empirical distributionp. More precisely, we have where z j ! is the factorial of z j and · ·, ·, · is the multinomial coefficient that determines the number of permutations [41, 24.1.2 ]. The rate of the CCDM, i.e., the number of bits per output symbol, converges to the entropy of the source, for asymptotically large number of output symbols [36], i.e., where H(X p ) − M −1 j=0 p j log(p j ) is the entropy of the discrete r.v.. On the other hand, for finite block lengths, the CCDM exhibits a rate lower than the source entropy. The rate loss can be upper bounded from [36] as where H(Xp) is a entropy of the n p -type distributed r.v. Xp. For example, R loss < 7.5 × 10 −4 bits/symbol for 8-PAM with frame length n p = 64800, adopted in the second generation digital video broadcasting over satellite (DVB-S2) [42].
For the invertible mapping between the bits and the symbols, large lookup table (LUT) can be used, in principle. However, for long block lengths, the size of the table is too large to be useful, i.e., 2 R DM np . Alternatively, the mapping can be achieved in an algorithmic manner, e.g., arithmetic coding is considered for the CCDM implementation in [36].
3) Channel Coding: In order to achieve reliable communication with high spectral efficiency close to the channel capacity, an FEC scheme should be employed. Since one of the main targets is to keep the computational complexity low, we opt for binary FEC encoders, as they have low complexity compared to non-binary methods. In this regard, the probabilistically shaped M -PAM signal is first mapped into a binary string using a mapper B, where each element of x p is labeled where b i, is the th bit level of the ith symbol. The vector b b 1, , b 2, , · · · , b np, contains all the bits of level ∈ {1, 2, · · · , m}. A single binary string, z p ∈ {0, 1} m np , is formed from the concatenation of the mapped bits for all the n p symbols, where The proper choice of the binary mapper B improves the performance of the scheme, e.g., the reflected binary mapping (Gray code) yields good performance [17].
In this scheme, any systematic binary FEC encoder with rate c, dimensionk m n p , and block lengthñ k /c can be employed. The redundant information in terms of the parity bits can be generated by where the multiplication is in the Galois field of two elements, and P ∈ {0, 1}k ×(1−c)ñ can be found by putting the code generator matrix in the standard form, Ik|P , with Ik denoting thẽ k ×k identity matrix.
Although the vector z p at the input of the FEC is probabilistically shaped, the parity bits, z u , tend to be uniformly distributed. This is attributed to the fact that each redundancy bit results from a modulo-2 sum of a large number of bits [17], [43]. These parity bits are mapped to the corresponding unipolar M -PAM symbols by binary demapper B −1 , i.e., which are also uniformly distributed.

B. Capacity of the Sparse-Dense Transmission
The capacity of sparse-dense signaling can be considered as an upper bound on the achievable rate of our system, for a given FEC rate. Hence, it is beneficial to compute the maximum mutual information (i.e., capacity) of the SDT, regardless if such a rate is achievable or not by the proposed SpaDCoM scheme. First, let us define X p as the r.v. representing the probabilistically , and X u as the r.v. representing the uniformly distributed symbols, {X i } n i=n c+1 . The capacity of the SDT, C SDT (g), can be found from (1) and (6) as the optimal value for the following optimization problem maximize ∆>0,p∈S where p j P{X p = a j } and u j P{X u = a j } = 1/M , for j ∈ {0, 1, · · · , M − 1}. Similar to (6), the optimization problem is convex in the probability vector p, for a fixed ∆, where (1 − c)I ∆ (u|g) and 0.5 ∆(1 − c)(M − 1) are constants that do not depend on p [34]. Hence, the interior-point algorithm can be used to find the optimal probabilities. Regarding ∆, the golden section method can be adopted to obtain the value of ∆ that maximizes the achievable rate within its feasibility range. The range of ∆ that satisfies the power constraint can be found from where the second inequality holds with equality if p = [1, 0, · · · , 0]. The SDT capacity can now be expressed, for a given g, as C SDT (g) = R SDT (g, ∆ * , p * ), where ∆ * and p * are the optimal symbol spacing and probabilities obtained as the solution of optimization problem (17).

C. Bit Metric Decoder
In the decoder, the FEC decoding is performed before the distribution dematching (i.e., reverse concatenation architecture). The reverse concatenation method prevents the common problem of the burst of errors after the distribution dematching, due to the receipt of erroneous symbols from the channel [22]. For the FEC decoding, the capacity of the SDT, computed in the Section III-B, could be achieved using an optimal symbol-metric decoder (SMD). However, the computational complexity of SMD is high. On the contrary, a bit-metric decoder (BMD) can yield a rate close to C SDT (g), while having lower complexity. For BMD with soft decisions, one real number is computed for each bit level of each received symbol, representing the likelihood of this bit. The log-likelihood ratio (LLR) of the th bit level can be written given the received symbol y i as where g is the value of the channel gain which is correctly estimated at the receiver, B i, is a r.v. representing the th level bit of the ith symbol, and the sets X 0 and X 1 include all the values of x such that the th level of their binary mapping equals 0 and 1, respectively. Since the distribution of the signal X i depends on the time instant i for the SDT scheme, the LLR can be written as , Since the LLRs are the sufficient statistics to recover the transmitted bits; they are provided to the FEC decoder for soft-decision decoding of the binary bits. Finally, the estimated bits are mapped into the corresponding symbols. The distribution dematching maps the first n p estimated symbols into their k p associated bits.
The achievable rate of the BMD for PAS has been investigated in [17], [44]. For SDT, let Then, an achievable rate with BMD for the SDT can be written as where [x] + max(x, 0) gives the maximum between x and zero. Equation (21) is due to the one-to-one mapping between X and its binary vector representation B(X). The maximum achievable rate using BMD can be found by maximizing R BMD (g, ∆, p) subject to an average power constraint. Nevertheless, the problem is not convex [45]; hence, we propose to obtain an achievable rate (not the maximum) by considering the distribution obtained by solving (17).
The achievable rate of the SDT with BMD can now be written as C BMD (g) R BMD (g, ∆ * , p * ). 4 Although the achievable rate of BMD is less than that provided by SMD, the rate loss is small, as illustrated in the numerical results.
Another metric to quantify an achievable information rate is the generalized mutual information (GMI), defined in [22,Eq. (12)]. The GMI quantifies the number of transmitted bits per symbol in a way similar to what mutual information does [46]- [48]. However, the GMI considers a mismatched decoding metric, e.g., BMD, in contract to the implied optimal decoder to achieve the rate indicated by the mutual information. For PAS, it has been shown that the GMI equals the achievable rate with BMD [22]. 5 For the proposed PS scheme with SDT, the GMI has the same expression as the achievable rate of the proposed scheme under BMD, R BMD (g, ∆, p). This can be proved starting from [22,Eq. (13)] and by noting that the transmitted bits are probabilistically shaped only for 100 c % of the channel uses.
In order to achieve reliable communication, the transmission rate should be less than the maximum achievable rate for the proposed SpaDCoM with BMD. Hence, it is essential to determine the transmission rate of the SpaDCoM. Since the information bits need to be transmitted are k p bits in n channel uses, the overall transmission rate can be written from (11) when n → ∞ as It is clear from (22) that the transmission rate depends on the input distribution and the FEC rate. Therefore, it is essential in the system design to quantify the maximum FEC rate such that the transmission rate is achievable, i.e., R(p) ≤ R BMD (g, ∆, p). In the literature, the achievable binary code rate and normalized generalized mutual information (NGMI) metrics are usually adopted to quantify the number of information bits per transmitted bits [22], [46]- [48]. In our scheme, the achievable binary code rate is represented from (21) and (22) for a given p as f The NGMI has the same expression as the achievable binary code rate for PAS, under some conditions on the decoding metric, as illustrated in [49].

D. Optimal Operating Point
For a fixed FEC, the optimal distribution p * that maximizes the achievable rate in (17) can lead to an unachievable transmission rate, i.e., R(p * ) > R BMD (g, ∆ * , p * ). In this regard, we investigate the optimal operating SNR that permits the maximum transmission rate such that reliable communication is still feasible. First, we describe the case when the CSI is known at both the transmitter and receiver. For a given instantaneous optical SNR, gP/σ, the optimal distribution that maximizes (17) is calculated. Accordingly, the capacity C SDT (g), the achievable rate with BMD, and the transmission rate are computed from (17), (21), and (22), respectively. Then, the optimal operating point is the intersection between the achievable and transmission rates, i.e., R(p * ) = R BMD (g, ∆ * , p * ). In order to achieve reliable communication at the optimal operating point, an asymptotically long FEC code should be used. For practical finite-rate codes, lower rates should be considered, e.g., R = R BMD − R backoff , where R backoff > 0 is a back-off rate to account for the non-optimal FEC codes.
In Fig. 2, the capacity, achievable rate, and transmission rate of BMD, SDT, and uniform signaling are depicted for various SNRs. Several upper and lower bounds on the capacity of intensity channels, i.e., [3, (3) and (11)], [4, (26) and (28) We can see that for SNRs lower than the optimal point, the transmission rate is unachievable, while at higher SNRs, the rate gap increases.
For Fig. 2b with c = 0.8, the optimal operating point is at SNR of 2.8 dB and transmission rate of 1.115 bpcu, which are lower than their corresponding values for c = 0.9. Also, we can notice that the difference between the achievable rate of M -PAM, i.e., all the symbols can be probabilistically shaped, and SDT is 0.125 bpcu for c = 0.8, which is larger than the corresponding value for c = 0.9, i.e., 0.075 bpcu. The reason is that the number of uniformly distributed symbols (not probabilistically shaped) is inversely proportional to the coding rate.
Regarding the blind encoder, i.e., the CSI is available only at the decoder, we can not guarantee that the transmission rate is always less than the achievable rate. This is attributed to the fact that the achievable rate is a r.v., as it is a function of the channel irradiance fluctuation. Hence, there is a non-zero probability that the transmission rate is not achievable, leading to an outage probability. In this case, the encoder can be designed with the worst-case channel condition to guarantee that the outage probability is less than a predefined threshold. This can be achieved by considering a fixed channel gain g =ḡ such that the outage probability is upper bounded by the target level, as described in Section VI.

IV. THE CAPACITY OF SPADCOM WITH SDT AND OPTIMAL FEC RATE
In the previous section, we obtain a single optimal point for each FEC rate, i.e., the rate at which the SDT capacity equals the transmission rate. In this section, we provide a technique to obtain the maximum achievable transmission rate of the proposed scheme and the corresponding input distribution for any given SNR. This can be done by also optimizing the rate of the channel encoder. Additionally, in order to guarantee that the transmission rate is achievable for the considered SNR, an additional constraint is added such that the input distribution yields an achievable transmission rate. More precisely, the capacity of the proposed scheme can be formulated as maximize In order to simplify (24), the rate constraint can be eliminated by choosing the channel coding rate such that the achievable rate equals the transmission rate. In particular, the FEC rate is written for a fixed p from (17) and (22) as Since the transmission rate equals the achievable rate, it is equivalent to maximize R(p) or R SDT (g, ∆, p). Then, by substituting the channel coding rate into the (24), the optimization problem becomes maximize ∆>0, p∈S where β ∆ (0.5∆ (M − 1) − P )/I ∆ (u|g) and the equality constraint is eliminated, as it is always active because of the proper formulation of the coding rate in (25). Let us define the optimal value of the maximization problem by R SDOR (g) R SDT (g, ∆, p * ), where the PMF p * maximizes (26) for a given ∆. Since R SDOR (g) is the maximum possible rate after the PMF optimization, we have with equality if p = p * . By multiplying both sides by H(X p ) − I ∆ (p|g) + I ∆ (u|g) ≥ 0, equation (27) can be written as with equality if p = p * . Hence, we need to minimize (28) to obtain near optimal p, i.e., Generally, the capacity R SDOR (g) is not known a-priori. Hence, we can substitute it with a preliminary estimate, r R SDT (g, ∆,ṗ), around an initial pointṗ that is iteratively updated, as suggested in [50]. After some manipulations and by dividing the objective function by the constant I ∆ (u|g), the problem becomes minimize p∈S f 0 (p,r) The optimization problem (30) is not convex, but it can be reformulated as a difference of convex (DC) problem. In fact, both the objective function and power constraint are represented as DC functions, for r > I ∆ (u|g) and β ∆ ≥ 0 (i.e., ∆ ≥ 2 P/[M −1]). This is because −I ∆ (p|g) and −H(X p ) are convex functions, while a T p is an affine.
A local minimum for the DC problem can be obtained through many iterative algorithms, e.g., convex-concave procedure [51]. In this method, the subtracted convex function is approximated by a Taylor expansion of the first-order around a feasible initial point,ṗ, which is successively updated till convergence. More precisely, ψ 0 (p, r) is replaced with Similarly, the constraint f 1 (p) − ψ 1 (p) ≤ 0 can be convexified by substituting ψ 1 (p) with Algorithm 1 SpaDCoM with Optimal FEC Rate 1: Input ∆,ṗ, δ 1 , δ 2 %δ 1 and δ 1 are the stopping criteria tolerances 2: repeat 3: r := R SDT (g, ∆,ṗ) 4: repeat 5: Solve the following convex problem updateṗ :=p 8: until |δ f | < δ 1 9: until |R SDT (g, ∆,ṗ) − r| < δ 2 10: p ∆ :=ṗ 11: Output p ∆ , R SDT (g, ∆, p ∆ ), R BMD (g, ∆, p ∆ ) The convex-concave procedure is described in Algorithm 1 for a fixed ∆. Then, the golden section method can be used to search for the optimal constellation spacing, ∆ * , and the associated probabilities p * p ∆ * that lead to the maximum achievable rate. Additionally, it can be proved that the solution of the convexified problem is still subject to the average power constraint in (30b). From the convexity of ψ 1 (p), we have ψ 1 (p) ≥ψ 1 (p,ṗ), leading to Although the solution is feasible, some of the feasibility range is lost because of the power constraint relaxation. Hence, we can obtain a feasible sub-optimal input distribution that yields an achievable rate of the scheme with SMD metric.
The capacity of the proposed scheme can be obtained through the above-indicated procedure; however, optimal FEC encoders are required. More precisely, the optimal channel coding rate, c(p * ), can take any real value between zero and one; however, most of the encoders allow only several operational modes with predefined FEC rates. In this case, the system can sub-optimally operate on the maximum allowable FEC rate, which is slightly less than c(p * ). This leads to a transmission rate that is lower than the achievable rate. Nevertheless, there is no guarantee that the average power constraint is not being violated when operating on lower FEC rate. In the following section, we propose an approach to calculate the achievable rate of SpaDCoM with practical BMD with a finite set of coding rates.

ENCODERS
In this section, we describe how the proposed adaptive scheme can be implemented using practical off-the-shelf channel encoders with predefined FEC rates and BMDs. For example, the DVB-S2 considers LDPC channel encoders with coding rate while the newer standard DVB-S2X permits more rates [42], [52]. Then, the computational complexity of the proposed scheme is analyzed.

A. Capacity-Achieving Distribution
We develop a method to obtain the input distribution that maximizes the spectral efficiency of the proposed SpaDCoM with BMD. For a given FEC rate, we would like to compute the maximum possible transmission rate that is achievable by the SpaDCoM scheme. The transmission rate is maximized, while satisfying both the power and rate constraints. The maximum achievable transmission rate can be found by solving the following optimization problem maximize ∆>0, p∈S where R backoff ≤ (1 − c)I ∆ (u|g) is a back-off rate to account for the reduced achievable rate with BMD, which can be iteratively set as indicated in Algorithm 2. 6 The three main differences between the optimization problem in (34) when compared to (24) is that the transmission rate is maximized rather than the mutual information, the rate is forced to be less than the mutual information with a margin allowing the use of a low complexity BMD, and the FEC rate is fixed.
These modifications allow obtaining the maximum achievable rate of the SpaDCoM scheme with practical FEC encoders and BMD. Unfortunately, the optimization problem (34) is not convex, because of the second constraint. Nevertheless, it can be convexified by noting that the constraint can be reformulated as a DC problem. Then, the convex-concave procedure is applied as before.
More precisely, we iteratively solve the following convex problem whereṗ is an initial feasible point, as shown in Algorithm 2. The optimal constellation spacing is obtained also using golden section search method. The procedure is repeated for each supported FEC rate, and the optimal channel coding rate c * ∈ R c is the one that yields the maximum transmission rate for the proposed scheme, denoted by R SpaDCoM (g).
The same procedure can be repeated for all the modulation orders M supported by the encoder.
Then, the M -ary modulation that yields the largest transmission rate, for the considered SNR, is selected. Alternatively, in order to reduce the modulation complexity, we can opt for the minimum M -ary modulation that yields the target rate for the considered SNR.
The proposed scheme adapts the rate according to the channel condition. Hence, the ergodic transmission rate for the SpaDCoM scheme can be written as where R SpaDCoM (g) is the maximum transmission rate for a given instantaneous SNR, i.e., gP/σ, and f G (g) is the PDF of the irradiance in (2).

B. Computational Complexity of the Proposed Scheme
In the following, the computational complexity of the proposed scheme is analyzed. The complexity is mainly due to the numerical optimization, required for computing the input distribution, and distribution matching/dematching through the CCDM. 7 For computing the capacity-achieving distribution of the proposed scheme, Algorithm 2 should be run for each value of ∆. Let us define z δ as the number of iterations till rate convergence in Algorithm 2,  Solve the convex problem (35) with interior-point method 6: Assign the optimal PMF top 7: Updateṗ :=p 9: until |δ R | < δ 10: p * ∆ :=ṗ 11: δ BMDRp := R BMD (g, ∆, p * ∆ ) − R(p * ∆ ) % The transmission rate back-off with respect to the BMD 12: if δ BMDRp < 0 then 13: which depends on the stopping criteria δ, z GS as the number of iterations till convergence for the golden section algorithm to find the optimal ∆, and z c as the number of code rates supported by the FEC coder. Hence, in order to find the optimal distribution, we need to solve z δ z GS z c convex optimization problems. The computation complexity for solving a single convex optimization problem using interior-point method is on the order of M 3 , i.e., O(M 3 ) [53], [54].
The optimization can either be performed online or offline, depending on the available computational capability at the transmitter. For instance, considering offline optimization, the optimal modulation order, probabilities of symbols, constellation spacing, and FEC rate can be obtained for a predefined set of SNRs with fine arbitrary quantization and stored in the memory, reducing the computational complexity at the expense of some rate loss due to the quantization.
For the DM, the CCDM, implemented using arithmetic coding [36], requires k p iterations for the matching, while it needs n p for dematching. Each iteration involves M additions, multiplications, and comparisons [55]. Unfortunately, the algorithms for arithmetic coding are sequential in nature; hence, it is a challenging task to parallelize the implementation [39].
In the literature, several methods have been proposed to reduce the computational complexity and increase the parallelization capability of DM methods, including CCDM [37]- [40]. In particular, efficient implementation of the CCDM using finite-precision algorithms has been proposed in [38]. Also, a parallel architecture for a CCDM has been proposed with a subset ranking algorithm rather than arithmetic codes [40]. Alternatively, the Multiset-Partition Distribution Matching (MPDM) is a non-constant composition DM that can achieve lower rate loss compared to CCDM for a given n p [39]. However, the output sequences have different empirical distributions that match the target distribution only on average, i.e., the ensemble average of the output sequences imitates the desired distribution. Hence, the empirical distribution of a specific sequence may deviate from the desired distribution.
Besides DM algorithms, there are other methods for PS with indirect signal shaping algorithms (e.g., sphere shaping and shell mapping) that opt for the most energy-efficient codewords [56], [57]. On the contrary to DM methods that target a specific pre-designed distribution in the output, the indirect methods are designed to achieve a target rate, i.e., a constant k p for a given n p . For example, in shell mapping, 2 kp codewords are chosen from all possible sequences that fulfill the energy budget, and the others are neglected. The potential codewords are on the surface or inside an n p -dimensional sphere [10]. The computational complexity of enumerative sphere shaping [56] and shell mapping [57,Algorithm 2] is O(n p ) and O(n p 3 ), respectively [55]. The rate loss for indirect signal shaping can be smaller than CCDM; however, the rate loss as a function of the sequence length is not the only parameter to judge the performance of the PS algorithms. In fact, the CCDM can afford a longer frame length with lower computational complexity compared to shell mapping [22]. Therefore, it could be preferable to use CCDM over sphere mapping for longer frame lengths, even if the CCDM rate loss is higher for a fixed frame length. Another issue for sphere shaping and shell mapping is that they are not optimized for IM/DD optical communications. This is attributed to the fact that shell and sphere mapping are asymptotically optimal for AWGN channels with electrical power constraint (restricted signal variance) [10]. On the other hand, the input distribution for IM/DD is subject to non-negativity and optical power constraints (restricted signal mean); hence, a sphere is not necessarily the optimal shape for the shell in this case.
It may be of interest also to investigate the power consumption of the proposed scheme.
In fact, the power dissipation can be attributed to two primary sources: the communication and computation power consumption. The proposed scheme reduces the communication power consumption, as it can achieve the same rate compared to uniform signaling with less transmitted power (up to 2 dB). The reason is that the CCDM is asymptotically optimal from a power efficiency perspective [22]. Nevertheless, the power dissipated in computations increases due to the additional operations required for the CCDM. The savings in the communication power increases significantly with the distance, while the computational power dissipation is distanceindependent. Hence, the proposed scheme is energy efficient for backhauling applications. However, a careful analysis is required to precisely judge the total power consumption, accounting for the considered hardware components, e.g., the microprocessor and power amplifier, their impact on the actual dissipated energy for each arithmetic operation, and the savings in the transmitted power [58], [59].

VI. OUTAGE PROBABILITY DUE TO THE IRRADIANCE FLUCTUATIONS
The outage probability of the proposed SpaDCoM is the probability that a given transmission rate is not achievable, because of the irradiance fluctuations. The outage probability can be written as where R is the transmission rate. In the proposed SpaDCoM when the CSI is available at the transmitter, the encoder adapts the transmission rate according to the channel condition g. From the rate constraint (34), the rate is always achievable with SMD and optimal FEC, even with zero back-off rate. For BMD and practical FEC, an appropriate R backoff in (34) should be considered to achieve an acceptable error performance.
Regarding the proposed blind SpaDCoM when the CSI is known only at the receiver, there is a non-zero outage probability. In order to calculate it we start by defining the threshold for the irradiance where R + is the set of all positive real numbers. The achievable rate R BMD (g o , ∆, p) is nondecreasing in g. Hence, from (37), (38), and [60,Eq. 11], the outage probability can be rewritten in closed-form as where G 2,1 1,3 (·|·) is the Meijer G-function [61].  [5] and [3].
For the blind SpaDCoM, we propose to design the encoder such that the outage probability is less than a predefined thresholdγ. In this regard, the transmission rate is calculated for a fixed irradianceḡ(γ) defined asḡ which can be found from the inverse of the turbulence fading cumulative distribution function.
Finally, the transmission rate is obtained for a given P/σ by substituting g withḡ in the procedures indicated in section V. Consequently, the outage probability is bounded belowγ as required, i.e., P outage (R SpaDCoM (ḡ)) ≤γ.

VII. NUMERICAL RESULTS
In this section, Monte Carlo simulations and numerical results are depicted to evaluate the performance of the proposed SpaDCoM schemes. The uniform signaling, M -PAM capacity, capacity bounds in [5, (1)] and [3, (3)], and coded modulation scheme in [23] are used as benchmarks for the performance. In all numerical results, the optimal distribution for the proposed scheme is obtained according to the procedures in Section V. The back-off rate, required input parameter for Algorithm 2, is set as R backoff = 0.05. For the channel coding, we consider the LDPC DVB-S2 with c ∈ R c , defined in Section V.
In Fig. 3 In Fig. 4, the optimal PMF obtained following the procedures in Section V is shown for various SNRs and modulation orders M ∈ {4, 8}. It can be noticed that the symbols with low amplitudes tend to have higher probabilities than the symbols with larger amplitudes at low SNRs. This permits large constellation spacing ∆ * without violating the average optical power constraint. For instance, the maximum value of ∆ in OOK increases with the probability of the zero symbol, i.e., ∆ ≤ P/(1 − p 0 ). For high instantaneous SNRs and a fixed M , the distribution is almost uniform, and the distance between the symbols is small. This is attributed to the fact that for asymptotically high SNR the mutual information approaches the source entropy that is maximized using equiprobable symbols. 8 In Fig. 5, the performance of the proposed scheme is compared with both uniform signaling and pairwise coded modulation [23], in terms of the frame error rate (FER) using Monte Carlo simulation. The transmission rates of the schemes are kept constant at 1.5 bpcu. The minimum SNR that permits transmission at this rate can be found from Fig. 2a for the capacity lower  SpaDCoM, the minimum SNR that is required for an achievable rate of 1.5 bpcu is 5 dB, as shown in Fig. 3. The optimal symbol probabilities, constellation spacing, and FEC rate are p * = [0.53, 0.25, 0.14, 0.08], ∆ * = 1.18, and c = 9/10, respectively. These parameters are  shown. The FEC rate and the pairwise distribution, i.e., two consecutive symbols have the same probability, are adjusted to yield the target rate and power. In particular, we set the pairwise PMF as [0.405, 0.405, 0.095, 0.095] with 1.14 constellation spacing, while c = 0.9. The LDPC code adopted for DVB-S2 is considered for all the schemes with word length 64800. For the decoding of the proposed scheme, the LLR is computed from (20) assuming perfect estimation of the channel gain g. We can see that the FER exhibits a phase transition phenomenon around the SNR which is associated with the transmission rate. Also, the proposed scheme achieves about 1.3 dB and 0.25 dB reduction in transmitted power compared to uniform and pairwise signaling, respectively, for FER= 10 −3 and R = 1.5.
The ergodic rates of the SpaDCoM CSI aware scheme, (36), and the uniform approach are As a benchmark, we calculate ergodic upper and lower bounds on the capacity from the bounds in [4]. For comparison, we also consider the lognormal model for the atmospheric turbulence induced fading in [62,Eq. (34)]. The proposed scheme achieves better ergodic performance compared to the uniform signaling based method, e.g., 2.5 dB at R = 0.56. It can be seen that as the turbulence increases (i.e, the Rytov variance σ 2 R ), the ergodic rate decreases. Also, for small turbulence with σ R = 0.1, the ergodic rate approaches the transmission rate for AWGN channels, while the gap between the rate and the ergodic upper bound of the capacity in [4] is about 2 dB at R = 1.2 bpcu.
Finally, the performance of the blind SpaDCoM scheme when CSI is known only at the decoder is illustrated in Fig. 7. The transmission rate of the proposed scheme versus P/σ is compared with the uniform signaling for different turbulence conditions with Gamma-Gamma distributed turbulence. The design criteria is according to Section VI with upper bounded outage probability such that P outage ≤γ = 10 −4 . For the proposed scheme, we set the transmission rate as c = 0.9.
The FEC rate for uniform signaling is chosen such that the transmission rate, c log(M ), is less than the achievable rate with uniform signaling, R BMD (ḡ(γ), ∆ u , u). The SpaDCoM encoder is superior to the uniform method with about 2 dB for R = 0.5 bpcu, and with 1 dB for R = 1.5 and σ R = 0.5.

VIII. CONCLUSION
In this paper, we propose a coded modulation scheme for free-space optical based backhaul applications. The encoder is adaptive to the atmospheric turbulence-induced fading with arbitrary fine granularity. In particular, the signal constellation is probabilistically shaped by a low complexity fixed-to-fixed length distribution matcher to approach the capacity of FSO channels with IM/DD. The proposed method can employ any efficient off-the-shelf FEC encoder, and it can also operate when the CSI is known only at the receiver. The proposed scheme approaches the capacity of unipolar M -PAM. Moreover, it outperforms the uniform signaling based encoders.
For instance, the probabilistic based scheme can achieve a reduction in the transmitted power up to 2 dB compared to the uniform signaling, at an ergodic rate of 0.5 bpcu.