Parameter Estimation of Gamma-Gamma Fading with Generalized Pointing Errors in FSO Systems

In this paper, the challenges of parameter estimation for the Gamma-Gamma turbulence channels with generalized pointing errors are addressed. The Kolmogorov-Smirnov goodness-of-fit statistical test results indicate that the approximate probability density function obtained by the saddlepoint approximation (SAP) method provides a better approximation for a larger value wz , and this means that the proposed method is more efficient for the FSO links over long distances when the transmit divergence angle at the transmitter side is fixed. Also, an additional parameter k needs to be estimated in addition to the shaping parameters α and β under the SAP method. An estimation scheme for the shaping parameters is proposed based on the SAP method. The performance of the proposed estimation is investigated by using the mean square error (MSE) and normalized mean square error (NMSE). The results indicate the proposed estimator exhibits satisfactory performance in both noiseless and noisy environments. The effects of the receiver aperture on the estimation performance are also discussed.


Introduction
Free-space optical (FSO) communication systems offer strong alternatives to radio-frequency (RF) communication systems because of their merits including an unregulated spectrum, ease of deployment, enormous available bandwidth, and high level of security. These kinds of systems are well suited for a wide range of applications including wireless access points, wireless local area networks, and metropolitan area networks [1]. However, FSO system reliability can be easily deteriorated in the presence of various weatherinduced disruptions. For example, the authors investigated the quality of FSO systems in Malaysia under various weather conditions and demonstrated that heavy rains cause very high BER and thus limit maximum link length [2]. The same conclusion can be also seen in [3], where the FSO link availability under the heavy-rain condition in Changsha, China, was studied.
The inhomogeneities in the temperature and pressure of the atmosphere result in variations of the refractive index; these variations cause atmospheric turbulence [4]. As shown in [5], the major limiting performance for FSO links is the atmospheric turbulence-induced fading (also called the scintillation). Numerous statistical models have been proposed to describe the scintillation in the literature; these models are effective with specific ranges of turbulence conditions. The lognormal distribution is generally accepted for weak turbulence and is suitable for characterizing FSO communication in clear sky links over several hundred meters [6]. The K distribution [7] and I-K distribution [8] are proposed as the probability density function (PDF) of irradiance fluctuations in strong and weak-to-strong refractive turbulence regions. The lognormal-Rician distribution [9] is shown to have excellent fit with both simulation and experimental data over a wide range of weak-to-strong turbulence conditions. However, the applications of lognormal-Rician fading distribution in the literature are currently limited since it does not have a tractable closed-form PDF. The adoption of the Gamma-Gamma distribution has been popular because it fits the measured data with high accuracy, and it has a simple mathematical form from a performance analysis point of view [10]. Apart from the scintillation effects, the performance of FSO links is also vulnerable to pointing errors caused by building sway. With respect to the models for pointing errors, numerous proposals have been put forward, including Rayleigh, Rician, Hoyt, and Beckmann distributions. Note that the Beckmann distribution is a versatile model because it encompasses the three other models as special cases [11,12].
To apply turbulence models to the analyses of practical FSO systems, an estimator is needed to effectively estimate the corresponding unknown parameters. The parameter estimation methods for the lognormal distribution and K distribution have been discussed in [13,14]. The maximum likelihood estimation (MLE) with expectation-maximization (EM) or the saddlepoint approximation algorithm is applied to characterize the lognormal-Rician turbulence model parameters [15,16]. The iterative EM algorithm based on the generalized Newton method using a nonquadratic approximation for the MLE of Gamma-Gamma parameters is investigated in [17]. A computationally efficient estimator using the method of moments is proposed in [18,19]. However, all of the above results are derived in the absence of pointing errors, which cannot be ignored in practical FSO systems. To the best of the authors' knowledge, an estimator for the parameters of the Gamma-Gamma distribution with a generalized model addressing pointing errors (Beckmann distribution) has not been reported in literature at this time.
In this paper, we study the problem of parameter estimation over Gamma-Gamma turbulence channels with generalized pointing errors. To avoid calculating the complicated integral, a saddlepoint approximation (SAP) method is introduced. With this approach, a highly accurate approximation can be achieved over a wide range of channel conditions according to the results of Kolmogorov-Smirnov (KS) goodness-of-fit statistical tests. Based on this method, an approximate log-likelihood function is developed. The performance of the proposed SAP estimator is investigated in terms of mean square error (MSE) and normalized mean square error (NMSE) by Monte Carlo simulation. The simulation results demonstrate that the proposed estimator provides satisfactory performance in both noiseless and noisy scenarios.

System and Channel Model
For a single-input single-output FSO system using on-off keying (OOK) modulation with intensity modulation/direct detection (IMDD), the received electrical signal is given by [20].
where y is the output of the channel, x is the transmitted symbol and is real valued, i.e., x ∈ f0, 1g, and n is the additive white Gaussian noise (AWGN) with zero mean and variance σ 2 n = N 0 /2, i.e., n ∼ N ð0, N 0 /2Þ . The real-valued channel gain (irradiance) h is considered to be a product of three factors: h = h l h a h p , where h l denotes the deterministic path loss that is described by the exponential Beers-Lambert law [21], h a is the fading induced by the atmospheric turbulence, and h p is the fading due to geometric spread and pointing errors. Without loss of generality, h l is set to 1. Specifically, the irradiance h a is normalized in the sense that E½h a = 1. The average electrical SNR in the absence of pointing errors and path loss is defined by γ = ðE½h a Þ2/N0, or simply γ = 1/N 0 .
To cover a wide range of turbulence conditions (weak-tostrong), the Gamma-Gamma distribution model is adopted here. The PDF of h a is given by [5] where α and β are the shaping parameters related to the atmospheric conditions, Γð·Þ is the gamma function, and K v ð·Þ is the modified Bessel function of the second order v. Note that the parameters α and β are flexible; they are chosen to achieve an excellent fit between f h a ðh a Þ and experimental data. Assuming a spherical wave propagation, the expressions corresponding to the Gamma-Gamma parameters are given by [22][23][24] α = exp 0:49σ 2 R 1 + 0:18d 2 + 0:56σ 12/5 β = exp 0:51σ 2 R 1 + 0:69σ 12/5 R À Á −5/6 1 + 0:9d 2 + 0:62d 2 σ 12/5 where d = ffiffiffiffiffiffiffiffiffiffiffiffiffiffi kD 2 /4L p and σ 2 R = 0:5C 2 n k 7/6 L 11/6 . Here, k = 2π/λ is the optical wave number, D = 2R a is the diameter of the receiver aperture, L denotes the link distance, and C 2 n is the refractive index structure parameter, which plays a vital role in determining the turbulence strength. Also, it is acknowledged that C 2 n typically varies from 10 −17 m −2/3 for weak turbulence to 10 −13 m −2/3 for strong turbulence. Specifically, when D ≥ ρ sp , where ρ sp = ð0:55C 2 n k 2 LÞ −3/5 is the spherical wave coherence radius [25], we need to consider the effect of aperture averaging on the FSO system performance. Furthermore, the relationship β > α consistently holds according to Figure 1. Due to the lack of perfect alignment of the transmitter and receiver, the effect of pointing errors should be considered. The attenuation due to geometric spread and pointing errors can be approximated by [26] where r denotes the radical distance, v = Wireless Communications and Mobile Computing ½erf ðvÞ is the fraction of the collected power at r = 0, R a is the radius of a circular detection aperture, w zeq = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi is the equivalent beam width, and w z represents the received beam width. Also, the results presented in [26] indicate that Eq. (5) provides a good agreement with the exact value when the normalized beam width w z /R a > 6. Furthermore, two useful parameters are defined: φ x = w zeq /2σ x and φ y = w zeq /2σ y . These parameters are the ratios between the equivalent beam radius at the receiver and the corresponding pointing error displacement standard deviation at the receiver. Note that the radical displacement r at the receiver plane is expressed as r = ffiffiffiffiffiffiffiffiffiffiffiffiffi , where x and y represent the horizontal and elevation displacements, respectively. Both x and y are modeled as Gaussian random variables, i.e., x ∼ N ðµ x , σ x Þ and y ∼ N ð µ y , σ y Þ, and they are assumed to be independent of each other. It must be emphasized that the correlated variables x and y can be transformed into equivalent x 0 and y 0 that are independent when employing the method proposed in [24]. Then, the radial displacement r is distributed according to the well-known Beckmann distribution, which is given by [11] Although the effects of turbulence and pointing errors are not strictly independent, they can be approximated as independent for small jitter values [27]. Then, Eqs. (2), (5), and (6) are combined and reorganized using algebra; the resulting PDF of irradiance h takes the integral form as It is rather challenging to directly deal with the Eq. (7) due to the presence of complicated integral [28]. Thus, the estimation algorithms, such as MLE and EM with MLE, are time-consuming to accomplish for engineering applications in practical. However, it is demonstrated that the complicated integral in Eq. (7) can be efficiently computed by the SAP method.

SAP for the Shaping Parameter Estimation
In this section, we first present a fundamental description of the SAP method. Based on this technique, the SAP estimator for estimating parameters α and β is established in both noiseless and noisy cases. For convenience, the parameters included in the pointing error model, such as μ x , μ y , σ x , σ y , A 0 , and w zeq , are assumed to be known perfectly in the receiver. In what follows, y = ½y½0, y½1, ⋯, y½K denotes the independent and identically distributed K observations of the output signals, and C = ½C½0, C½1, ⋯, C½K is a discrete sequence after sorting y in ascending order. The symbol b • denotes the estimation of the shaping parameter in each iteration while e • stands for the saddle point.

Saddlepoint Approximation Representation.
The SAP method initially proposed by Daniels [29] is a valuable tool for approximating a density for the sum of random variables. More specifically, it can be used to obtain the asymptotic approximations to the computation of integrals. A basic description of the SAP method is illustrated in [29,30], and it is presented below with minor modifications for convenience. A function gðxÞ that is written in p-dimension integral form can be approximated by where D denotes the domain of integration, R represents the

Wireless Communications and Mobile Computing
field of real numbers, cðxÞ is some function that only depends on x, mðx, t 1 , t 2 , ⋯, t p Þ is a positive function, κðx, t 1 , t 2 , ⋯, t p Þ = ln ðmðx, t 1 , t 2 , ⋯, t p ÞÞ, and κ 00 is the p × p Hessian matrix of second derivatives for κ.t 1 ,t 2 , ⋯,t p are the p -dimension saddle points that maximize κðx, t 1 , t 2 , ⋯, t p Þ for each x. Note thatĝðxÞ is called the saddlepoint density when function gðxÞ is a density function. It can be easily deduced that such an approximation (8) provides high accuracy if the degree of κð·Þ is a "quadratic-looking" in a neighborhood oft 1 ,t 2 , ⋯,t p for each x.
Note that an all-one training sequence of length K is transmitted when the parameter estimation is performed at the receiver. Hence, according to Eqs. (1) and (7), the PDF y in a noisy case is obtained as Also, when γ approaches ∞, or equivalently when σ 2 N = N 0 /2 approaches zero, the Gaussian distribution approaches a Dirac delta function δð·Þ. Hence, the PDF of y in a noiseless case is given by Therefore, using the SAP representation in Eq. (8), the PDF of f ðyjx = 1Þ in the noisy and noiseless cases can be approximated by f appro1 ðyjx = 1Þ and f appro2 ðyjx = 1Þ, respectively. They are written in the following way where their corresponding κð·Þ and cð·Þ are derived as Since traditional gradient-based methods always converge to a local maximum for nonlinear and nondifferentiable functions [31], a more robust method is needed. Here, the saddle point in κð·Þ is found using the differential evolution (DE) method. Note the normalized factors k 1 and k 2 are introduced in Eq. (11) since f ðyjx = 1Þ is actually a density function, i.e., Ð f ðyjx = 1Þ = 1. The details of the derivation of these parameters are provided in Section 3.2. Furthermore, according to Eqs. (12) and (13), the obtained saddle point e θ that maximizes a function κ 1 ðy, h, θ, h p Þ or κ 2 ðy, θ, h p Þ satisfies the following relation and this leads to Equation (15) indicates that the saddle point e θ only depends on µ x , µ y .    Wireless Communications and Mobile Computing

KS
The approximate CDF F appro ðyjx = 1Þ derived in the previous subsection is not expressed in closed-form. Note that f appro ðyjx = 1Þ is a discrete sequence with K elements, and we can use interpolation functions to achieve continuity.
Intuitively, the larger the number of the observations K, the smaller the residual error caused by interpolation. Note that the interpolation functions can be readily evaluated and efficiently programmed in most standard software packages (e.g., MATLAB and MATHEMATICA). In all of the Monte Carlo simulations, we consider the first-order polynomial interpolation, namely, a linear interpolation for computational simplicity. Therefore, F appro ðyjx = 1Þ is derived as The normalized factor k is then carried out by In Figure 2, we present the KS test statistic results for different combinations of the parameters C 2 n , w z , γ, jitter variance, and boresight displacements in each axis of receiver plane. Note that D = 0:1 m for all of the figures in Figure 2 and Table 1. Other system configuration parameters, such as L and λ, have been used in most practical terrestrial FSO systems [25,33]. These are shown in Section 4. The present results are obtained by averaging the results of 100 simulation runs, and each run uses at least K = 10 4 data samples of y. The critical value T max = 0:0136 corresponds to a significance level of α = 5% ( [32], Eq. (9.73)). It is also included in Figures 2 and 3 as a benchmark.
In Figures 2(a) and 2(c) and Figure 2(e), the effects of the jitter variance and the nonzero boresight error on the approximation accuracy are investigated. These figures demonstrate that a larger value w z , smaller jitter variance, and nonzero boresight displacement lead to a better approximation of f appro ðyjx = 1Þ. In other words, the proposed method is more efficient for the FSO links over long distances when the transmit divergence angle at the transmitter side is fixed. Also, we find an interesting phenomenon that is also shown in [34]: the KS statistic significantly exceeds the critical value when w z is lower than the threshold and is below the critical when w z is larger than such a threshold.
In Figures 2(b), 2(d), and 2(f), we present the influence of the average SNR γ on the approximation accuracy when μ x = μ y = 1 m and σ z = σ y = 0:55 m. In comparison to the results presented in Figure 2(c), a higher approximation accuracy can be provided in the low SNR regime even if w z ≤ 4 m. This is because a more Gaussian-like shaped PDF is achieved at low SNRs. In addition, when comparing the curves in the three figures, we demonstrate that the KS test statistics gradually lose accuracy for small w z with increasing γ.
The same conclusion can be also seen in Figure 3, where the effects of the diameter of the receiver aperture D on the approximation accuracy are considered. Also, the results shown in Figure 3 indicate that an efficient approximation can be achieved for large aperture receivers in both noiseless and noisy scenarios.
Furthermore, we note that the KS goodness-of-fit test results play an important role in the estimation performance. More specifically, a better estimate can be achieved for the smaller KS goodness-of-fit test results.

SAP Estimator.
To reveal the importance of the proposed approximation, an estimator based on the SAP is established. As shown in Section 3.2, the proposed SAP estimator needs to estimate three types of parameters: the normalized factor k in addition to the shaping parameters α and β.

Estimator in a Noisy
Situation. For the K observations y , the MLE of the unknown parameters θ = fk 1 , α, βg in a noisy situation is obtained by maximizing the approximate log-likelihood function as where fh, e θ,h p g denotes the saddle point with respect to y½k and the current estimated parameters b α and b β.

Estimator in a Noiseless
Situation. The SAP estimator in a noiseless situation can be developed similarly, which is given by maximizing the following approximate log- where f e θ,h p g represents the saddle point given y½k ; b α and b β are the current parameter estimates.
The normalized factors k 1 and k 2 are derived according to Eq. (18) given the estimated parameters b α and b β. It must be emphasized that the e θ can be obtained according to Eq. (14). Note that the right-hand side of Eqs. (19) and (20) can be considered to have only one global maximum providing that the number of observations is large enough and the approx-imation is highly accurate. Thus, they can be solved efficiently by using a conventional iterative method such as the Newton-Raphson algorithm or the BFGS Quasi-Newton algorithm [35]. The first (second) derivative of L appro1 ðθÞ or L appro2 ðθÞ with respect to α and β can be approximated by the finite difference individually in each iteration.
Furthermore, we find that the obtained saddle pointh p in a noiseless situation satisfies the following proposition.
Proof. A proof is presented in the Appendix. ☐ Corollary 2. When A 0 ⟶ 1 and w z is large enough, we havẽ Proof. A proof is given in the Appendix. ☐

Wireless Communications and Mobile Computing
Theorem 3. When A 0 ⟶ 1 and w z is large enough, the normalized factor k 2 only depends on the parameters μ x , μ y , σ x , and σ y . This is carried out by A proof is given in the Appendix. Table 1 presents the results for the normalized factor k 2 for different combinations of the parameters for pointing errors. The simulated results presented here are the average values derived from 100 repeated simulations and K = 1 × 10 4 . From the table, we observe the theoretical results are almost the same as the simulated results over the range of the parameters examined (their absolute differences are not greater than 5 × 10 −4 ).

Initial Estimation for
Parameters α and β. The initial estimation for the parameters α and β can be derived by using a moment-based estimation. This has been used in [19] for Gamma-Gamma turbulence channels in the absence of pointing errors. It must be noted that obtaining the closed-form expression for the n-th moment of Eq. (7) is remarkably difficult. However, using the method proposed in [12], the n-th moment of Eq. (7) can be efficiently approximated by where E½h n denotes the expectation of h n . The notations φ 2 mod and A mod are, respectively, defined as Hence, the initial values for the α and β in a noiseless situation are computed by ( [19], Eq. (9)) where c and d are given by

Wireless Communications and Mobile Computing
Here, we use the fractional moments, i.e., n = 1/2, instead of positive moments to achieve better performance. Note that the successful application of fractional moments in the study of atmospheric scintillation has been discussed in [36].
However, fractional moments are not suitable for noisy situations. The relations defined in ( [37], Eqs. (21), (22), and (23)) hold in medium and high SNR regimes. Taking the first two equations, namely, ( [37], Eqs. (21) and (22)), we obtain the following by using algebraic manipulations Substituting Eq. (27) into Eq. (25) yields the initial estimation for parameters α and β in the noisy case. When the SNR is low, parameters α and β can be initialized as which is reasonable for some specified circumstances, such as FSO systems with small received apertures that operate in weak turbulence conditions according to Figure 1.
The estimation process of the SAP estimator is summarized in Figure 4. It should be noted that the spacing h should not be too small since it may cause the derived approximate log-likelihood function to vibrate in a small range and thus fail to reach the maximum value. Furthermore, the initial estimation corresponding to the different cases contributes to speed-up algorithm convergence.

Numerical Results
To investigate the SAP estimator performance, the MSE of b θ is studied using a Monte Carlo simulation, where MSE ½ b θ = Var½ b θ + ðE½θ∧ − θÞ 2 are defined for the component-wise    tor. The NMSE is an important performance metric to measure the quality of an estimator and is defined as [38]. In the simulation, the size of data samples K for D = 0:1 m is 2 × 10 3 while K is set to 7 × 10 3 for D = 0:2 m. M = 50 trials are employed to calculate the MSE and NMSE of the proposed SAP estimator. Note that a wavelength λ = 1550 nm, a link distance L = 5 km, and the spacing h = 10 −4 are chosen for all of the figures.
In Figure 5, we present the simulated MSE and NMSE performance of α and β for different C 2 n values. It can be observed that increasing the strength of atmospheric turbulence C 2 n decreases the MSE performance of α while it does not substantially change the NMSE of α. The MSE and NMSE performances of β and k 2 are not sensitive to the value of C 2 n , and they remain relatively flat as C 2 n increases. In order to investigate the impacts of w z on the estimator performance, we also show the MSE performance for different w z , i.e., w z = 5 m and w z = 10 m in Figure 5. By comparing these curves, we demonstrate there is no significant difference in their MSE performance for α, β, and k 2 . This may be due to their nearly equal maximum CDF errors under the SAP representation according to Figure 2(c).
In Figure 6, the simulated MSE and NMSE performance results in a noisy case are depicted. Performances similar to Figure 5 can be observed. Thus, we conclude that the proposed SAP estimator is feasible at the target SNR level.
Also, the estimation performance results for the noiseless case are included in Figure 6 as a benchmark. We find that the MSE performances for α and β are close when C 2 n ≥ 2 × 10 −14 m −2/3 . Furthermore, when comparing the curves in Figures 5 and 6, we draw the conclusion that a slight performance improvement for the estimation of the parameters α and k 2 can be achieved for jitters with a larger mean at the receiver, and this is. Figure 7 depicts the MSE and NMSE performance for larger receiver apertures, i.e., D = 0:2 m, in a noiseless case. Similar observations can be made 330 in this figure. However, we note that the MSE and NMSE of α improve when C 2 n ≤ 1 × 10 −14 m −2/3 . These results suggest that the performance of the SAP estimators is sensitive to the value of C 2 n in this case. Moreover, when comparing the curves in Figures 6 and 7, we find that a large aperture at the receiver results in better estimations for the normalized factor k 2 .
In conclusion, using the Monte Carlo simulation, the derived results for the MSE and NMSE performance in Figures 5-7 are indicative of producing accurate estimations for α and β.

Conclusion
In this work, we have proposed to use the SAP technique for the estimation of shaping parameters over Gamma-Gamma turbulence channel with generalized pointing errors. With the help of the SAP method, the approximate log-likelihood functions are developed in both noiseless and noisy condi-tions. The performance of proposed estimation scheme is investigated in terms of MSE and NMSE. The simulated results show that the SAP estimator provides a satisfactory performance over a wide range of channel conditions. Also, we find that larger receiver aperture and mean of jitters result to a better estimation for the normalized factor in a noiseless situation. x + μ 2  ðA:5Þ It is easy to see that the first line of the right-hand side of Eq. (A.6) is the Gamma-Gamma distribution. Hence, the normalized factor k 2 is derived as

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.