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

Authors: Maoke Miao; Weiwei Cai; Xiaofeng Li
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

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 valuewz, 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.


---

## Body

## 1. 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 weather-induced 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 andK 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.


## 2. 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].(1)y=hx+n,where y is the output of the channel, x is the transmitted symbol and is real valued, i.e., x∈0,1, and n is the additive white Gaussian noise (AWGN) with zero mean and variance σn2=N0/2, i.e., n∼N0,N0/2. The real-valued channel gain (irradiance) h is considered to be a product of three factors: h=hlhahp, where hl denotes the deterministic path loss that is described by the exponential Beers-Lambert law [21], ha is the fading induced by the atmospheric turbulence, and hp is the fading due to geometric spread and pointing errors. Without loss of generality, hl is set to 1. Specifically, the irradiance ha is normalized in the sense that Eha=1. The average electrical SNR in the absence of pointing errors and path loss is defined by γ¯=Eha2/N0, or simply γ¯=1/N0.To cover a wide range of turbulence conditions (weak-to-strong), the Gamma-Gamma distribution model is adopted here. The PDF ofha is given by [5](2)fhaha=2αβα+β/2ΓαΓβhaα+β/2−1Kα−β2αβha,ha>0,where α and β are the shaping parameters related to the atmospheric conditions, Γ· is the gamma function, and Kv· 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 fhaha and experimental data. Assuming a spherical wave propagation, the expressions corresponding to the Gamma-Gamma parameters are given by [22–24](3)α=exp0.49σR21+0.18d2+0.56σR12/57/6−1−1,(4)β=exp0.51σR21+0.69σR12/5−5/61+0.9d2+0.62d2σR12/55/6−1−1,where d=kD2/4L and σR2=0.5Cn2k7/6L11/6. Here, k=2π/λ is the optical wave number, D=2Ra is the diameter of the receiver aperture, L denotes the link distance, and Cn2 is the refractive index structure parameter, which plays a vital role in determining the turbulence strength. Also, it is acknowledged that Cn2 typically varies from 10−17m−2/3 for weak turbulence to 10−13m−2/3 for strong turbulence. Specifically, when D≥ρsp, where ρsp=0.55Cn2k2L−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.Figure 1
The ratio of the Gamma-Gamma shaping parametersβ/α as function of receiver aperture D and Cn2.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](5)hpr,L≈A0exp−2r2wzeq2,where r denotes the radical distance, v=πRa/2wz, A0=erfv is the fraction of the collected power at r=0, Ra is the radius of a circular detection aperture, wzeq=wz2πerfv/2vexp−v2 is the equivalent beam width, and wz 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 wz/Ra>6. Furthermore, two useful parameters are defined: φx=wzeq/2σx and φy=wzeq/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=x2+y2, 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 x0 and y0 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](6)frr=r2πσxσx∫02πexp−rcosθ−μx22σx2−rsinθ−μy22σy2dθ.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(7)fhh=∫0A0∫02π2αβα+β/2ΓαΓβ1hphhpα+β/2−1Kα−β2hαβhpwzeq24hp12πσxσyexp−−1/2wzeq2lnhp/A0cosθ−μx22σx2−−1/2wzeq2lnhp/A0sinθ−μy22σy2dhpdθ.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.


## 3. 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, A0, and wzeq, are assumed to be known perfectly in the receiver. In what follows, y=y0,y1,⋯,yK denotes the independent and identically distributed K observations of the output signals, and C=C0,C1,⋯,CK is a discrete sequence after sorting y in ascending order. The symbol ∙^ denotes the estimation of the shaping parameter in each iteration while ∙~ stands for the saddle point.


### 3.1. 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 gx that is written in p-dimension integral form can be approximated by(8)gx=∫D⊆Rpcxmx,t1,t2,⋯,tpdt1dt2⋯dtp=∫D⊆Rpcxexpκx,t1,t2,⋯,tpdt1dt2⋯dtp≈g^x=cx2πp/2κ″x,t∧1,t∧21/2expκx,t~1,t~2,⋯,t~p,where D denotes the domain of integration, R represents the field of real numbers, cx is some function that only depends on x, mx,t1,t2,⋯,tp is a positive function, κx,t1,t2,⋯,tp=lnmx,t1,t2,⋯,tp, 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,t1,t2,⋯,tp for each x. Note that g^x is called the saddlepoint density when function gx 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 of t~1,t~2,⋯,t~p for each x.Note that an all-one training sequence of lengthK 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(9)fyx=1=fhy∗fNy=∫0∞∫0A0∫02π2αβα+β/2ΓαΓβ1hphhpα+β/2−1Kα−β2hαβhpwzeq24hp12πσxσyexp−−1/2wzeq2lnhp/A0cosθ−μx22σx2−−1/2wzeq2lnhp/A0sinθ−μy22σy2γ¯πexp−γ¯y−h2dhdhpdθ.Also, whenγ¯ approaches ∞, or equivalently when σN2=N0/2 approaches zero, the Gaussian distribution approaches a Dirac delta function δ·. Hence, the PDF of y in a noiseless case is given by(10)fyx=1=limγ¯⟶∞fhy∗fNy=fhy.Therefore, using the SAP representation in Eq. (8), the PDF of fyx=1 in the noisy and noiseless cases can be approximated by fappro1yx=1 and fappro2yx=1, respectively. They are written in the following way(11)fyx=1≈fappro1yx=1=k1c1y2π3/2κ1″y,h~,θ~,h~p1/2expκ1y,h~,θ~,h~p,fyx=1≈fappro2yx=1=k2c2y2πκ2″y,θ~,h~p1/2expκ2y,θ~,h~p,where their corresponding κ· and c· are derived as(12)κ1y,h,θ,hp=α+β2−1lnh−α+β2+1lnhp+lnKα−β2hαβhp−−1/2wzeq2lnhp/A0cosθ−μx22σx2−−1/2wzeq2lnhp/A0sinθ−μy22σy2−γ¯y−h2,c1y=αβα+β/2ΓαΓβwzeq24πσxσyγ¯πy0,(13)κ2y,θ,hp=−α+β2+1lnhp+lnKα−β2yαβhp−−1/2wzeq2lnhp/A0cosθ−μx22σx2−−1/2wzeq2lnhp/A0sinθ−μy22σy2,c2y=αβα+β/2ΓαΓβwzeq24πσxσyyα+β/2−1.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 k1 and k2 are introduced in Eq. (11) since fyx=1 is actually a density function, i.e., ∫fyx=1=1.The details of the derivation of these parameters are provided in Section3.2. Furthermore, according to Eqs. (12) and (13), the obtained saddle point θ~ that maximizes a function κ1y,h,θ,hp or κ2y,θ,hp satisfies the following relation(14)−12wzeq2lnhpA0cosθ~=μx,−12wzeq2lnhpA0sinθ~=μy,and this leads to(15)µxsinθ~=µycosθ~.Equation (15) indicates that the saddle point θ~ only depends on µx,µy.


### 3.2. KS Goodness-of-Fit Tests

In this subsection, we investigate the validity of the proposed approximation in Eq. (12) by employing the KS goodness-of-fit statistical tests. It is shown in [32] that the KS goodness-of-fit tests measure the maximum value of the absolute difference between the empirical cumulative density function (CDF) of random variable (RV) y, Fyx=1 and the CDF Fapproy∣x=1 of the approximate PDF fapproyx=1. Note that we here omit the subscript (1 and 2) for notational simplicity. Hence, the KS test statistic is defined as(16)T≜maxFyx=1−Fapproyx=1.The approximate CDFFapproyx=1 derived in the previous subsection is not expressed in closed-form. Note that fapproyx=1 is a discrete sequence with K elements, and we can use interpolation functions to achieve continuity.Intuitively, the larger the number of the observationsK, 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, Fapproyx=1 is derived as(17)Fapproyx=1≈k∫C0yInterpolationfapprocx=1dc.The normalized factork is then carried out by(18)k≈1∫C0CKInterpolationfapprocx=1dc.In Figure2, we present the KS test statistic results for different combinations of the parameters Cn2,wz,γ¯, 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=104 data samples of y. The critical value Tmax=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.Figure 2
KS goodness-of-fit tests for the distributionFapproxy to approximate the distribution Fy with a 5% significance level. In (a), (c), and (e), the effects of the different jitter variances and nonzero boresight displacements in the noiseless case are considered. In (b), (d), and (f), the effects of the different γ¯ values in the noisy case are considered, where μx=μy=1 m and σx=σy=0.55 m.
(a)
μx=μy=0.3 m, σx=σy=0.55 m(b)
Low SNRγ¯=−10 dB(c)
μx=μy=1 m, σx=σy=0.55 m(d)
Medium SNRγ¯=35 dB(e)
μx=μy=1 m, σx=σy=1 m(f)
High SNRγ¯=65 dBTable 1
Normalized factork2 for the different parameters of pointing errors at the receiver when wz/Ra=4000 and Cn2=1×10−14 m−2/3.
μx (m)μy (m)σx (m)σy (m)Results of Eq. (22)Simulated results0.230.30.40.551.047731.047770.650.80.450.551.019351.019360.90.91.81.21.083331.083290.20.31.80.81.301711.30172110.550.5510.9999481.81.82.12.110.999509Figure 3
KS goodness-of-fit tests for the distributionFapproxy to approximate the distribution Fy with a 5% significance level, where μx=μy=1 m and σx=σy=0.55 m.
(a)
D=0.3 m(b)
D=0.3m, γ¯=65 dBIn Figures2(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 wz, smaller jitter variance, and nonzero boresight displacement lead to a better approximation of fapproyx=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 wz is lower than the threshold and is below the critical when wz is larger than such a threshold.In Figures2(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 wz≤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 wz with increasing γ¯.The same conclusion can be also seen in Figure3, 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.


### 3.3. SAP Estimator

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


#### 3.3.1. Estimator in a Noisy Situation

For theK observations y, the MLE of the unknown parameters θ=k1,α,β in a noisy situation is obtained by maximizing the approximate log-likelihood function as(19)Ly;α,β≈Lappro1y;θ≈Klnk1+1K∑k=0K−1c1yk+32ln2π+1K∑k=0K−1κ1yk,h~,θ~,h~p−12K∑k=0K−1lnκ1″yk,h~,θ~,h~p,where h~,θ~,h~p denotes the saddle point with respect to yk and the current estimated parameters α^ and β^.


#### 3.3.2. 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-likelihood function(20)Ly;α,β≈Lappro2y;θ≈Klnk2+1K∑k=0K−1C2yk+ln2π+1K∑k=0K−1K2yk,h~,θ~,h~p−12K∑k=0K−1lnk2″yk,h~,θ~,h~p,where θ~,h~p represents the saddle point given yk;α^ and β^ are the current parameter estimates.The normalized factorsk1 and k2 are derived according to Eq. (18) given the estimated parameters α^ and β^. It must be emphasized that the θ~ 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 approximation 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 Lappro1θ or Lappro2θ 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.Proposition 1.
AsA0⟶1 and wz⟶∞, h~p closes to 1.Proof.
A proof is presented in the Appendix.☐Corollary 2.
WhenA0⟶1 and wz is large enough, we have(21)h~p=1−ε=1−2μx2+μy2wzeq2.Proof.
A proof is given in the Appendix.☐Theorem 3.
WhenA0⟶1 and wz is large enough, the normalized factor k2 only depends on the parameters μx, μy, σx, and σy. This is carried out by(22)k2=μy2σy2+μx2σx2μx2σy2+μy2σx2μx2+μy2σxσy.
A proof is given in the Appendix.Table1 presents the results for the normalized factor k2 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×104. 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).


#### 3.3.3. 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(23)Ehn=μn≈φmod2Amodnφmod2+n1αβΓα+nΓβ+nΓαΓβn,where Ehn denotes the expectation of hn. The notations φmod2 and Amod are, respectively, defined as(24)φmod2=wzeq22σmod,Amod=A0exp1φmod2−12φx2−12φy2−μx22σx2φx2−μy22σy2φy2.Hence, the initial values for theα and β in a noiseless situation are computed by ([19], Eq. (9))(25)α∧0≈cd−12c2d2−4d,β∧0≈cd+12c2d2−4d,where c and d are given by(26)c=n2Ey2φmod2+2/φmod2Amod2−Eyn+1/Eynφmod2+n+1/φmod2+nAmod−n2−1n2−n,d=nEy2φmod2+2/φmod2Amod2−Eyn+1/Eynφmod2+n+1/φmod2+nAmod−n−1n−n2.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(27)c=4Ey2−σn2φmod2+2/φmod2Amod2−Ey3−3σn2φmod2Amod/φmod2+1φmod2+3/Amod/Ey2−σn2φmod2+2−32d=2Ey2−σn2φmod2+2/φmod2Amod2−Ey3−3σn2φmod2Amod/φmod2+1φmod2+3/Amod/Ey2−σn2φmod2+2−1−2Substituting 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(28)α∧0=β∧0=1Ey2−σ2φmod2+2/φmod2Amod2−1,which is reasonable for some specified circumstances, such as FSO systems with small received apertures that operate in weak turbulence conditions according to Figure1.The estimation process of the SAP estimator is summarized in Figure4. 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.Figure 4
Flow chart of the SAP estimator for the Gamma-Gamma shaping parametersα and β.


## 3.1. 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 gx that is written in p-dimension integral form can be approximated by(8)gx=∫D⊆Rpcxmx,t1,t2,⋯,tpdt1dt2⋯dtp=∫D⊆Rpcxexpκx,t1,t2,⋯,tpdt1dt2⋯dtp≈g^x=cx2πp/2κ″x,t∧1,t∧21/2expκx,t~1,t~2,⋯,t~p,where D denotes the domain of integration, R represents the field of real numbers, cx is some function that only depends on x, mx,t1,t2,⋯,tp is a positive function, κx,t1,t2,⋯,tp=lnmx,t1,t2,⋯,tp, 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,t1,t2,⋯,tp for each x. Note that g^x is called the saddlepoint density when function gx 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 of t~1,t~2,⋯,t~p for each x.Note that an all-one training sequence of lengthK 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(9)fyx=1=fhy∗fNy=∫0∞∫0A0∫02π2αβα+β/2ΓαΓβ1hphhpα+β/2−1Kα−β2hαβhpwzeq24hp12πσxσyexp−−1/2wzeq2lnhp/A0cosθ−μx22σx2−−1/2wzeq2lnhp/A0sinθ−μy22σy2γ¯πexp−γ¯y−h2dhdhpdθ.Also, whenγ¯ approaches ∞, or equivalently when σN2=N0/2 approaches zero, the Gaussian distribution approaches a Dirac delta function δ·. Hence, the PDF of y in a noiseless case is given by(10)fyx=1=limγ¯⟶∞fhy∗fNy=fhy.Therefore, using the SAP representation in Eq. (8), the PDF of fyx=1 in the noisy and noiseless cases can be approximated by fappro1yx=1 and fappro2yx=1, respectively. They are written in the following way(11)fyx=1≈fappro1yx=1=k1c1y2π3/2κ1″y,h~,θ~,h~p1/2expκ1y,h~,θ~,h~p,fyx=1≈fappro2yx=1=k2c2y2πκ2″y,θ~,h~p1/2expκ2y,θ~,h~p,where their corresponding κ· and c· are derived as(12)κ1y,h,θ,hp=α+β2−1lnh−α+β2+1lnhp+lnKα−β2hαβhp−−1/2wzeq2lnhp/A0cosθ−μx22σx2−−1/2wzeq2lnhp/A0sinθ−μy22σy2−γ¯y−h2,c1y=αβα+β/2ΓαΓβwzeq24πσxσyγ¯πy0,(13)κ2y,θ,hp=−α+β2+1lnhp+lnKα−β2yαβhp−−1/2wzeq2lnhp/A0cosθ−μx22σx2−−1/2wzeq2lnhp/A0sinθ−μy22σy2,c2y=αβα+β/2ΓαΓβwzeq24πσxσyyα+β/2−1.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 k1 and k2 are introduced in Eq. (11) since fyx=1 is actually a density function, i.e., ∫fyx=1=1.The details of the derivation of these parameters are provided in Section3.2. Furthermore, according to Eqs. (12) and (13), the obtained saddle point θ~ that maximizes a function κ1y,h,θ,hp or κ2y,θ,hp satisfies the following relation(14)−12wzeq2lnhpA0cosθ~=μx,−12wzeq2lnhpA0sinθ~=μy,and this leads to(15)µxsinθ~=µycosθ~.Equation (15) indicates that the saddle point θ~ only depends on µx,µy.


## 3.2. KS Goodness-of-Fit Tests

In this subsection, we investigate the validity of the proposed approximation in Eq. (12) by employing the KS goodness-of-fit statistical tests. It is shown in [32] that the KS goodness-of-fit tests measure the maximum value of the absolute difference between the empirical cumulative density function (CDF) of random variable (RV) y, Fyx=1 and the CDF Fapproy∣x=1 of the approximate PDF fapproyx=1. Note that we here omit the subscript (1 and 2) for notational simplicity. Hence, the KS test statistic is defined as(16)T≜maxFyx=1−Fapproyx=1.The approximate CDFFapproyx=1 derived in the previous subsection is not expressed in closed-form. Note that fapproyx=1 is a discrete sequence with K elements, and we can use interpolation functions to achieve continuity.Intuitively, the larger the number of the observationsK, 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, Fapproyx=1 is derived as(17)Fapproyx=1≈k∫C0yInterpolationfapprocx=1dc.The normalized factork is then carried out by(18)k≈1∫C0CKInterpolationfapprocx=1dc.In Figure2, we present the KS test statistic results for different combinations of the parameters Cn2,wz,γ¯, 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=104 data samples of y. The critical value Tmax=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.Figure 2
KS goodness-of-fit tests for the distributionFapproxy to approximate the distribution Fy with a 5% significance level. In (a), (c), and (e), the effects of the different jitter variances and nonzero boresight displacements in the noiseless case are considered. In (b), (d), and (f), the effects of the different γ¯ values in the noisy case are considered, where μx=μy=1 m and σx=σy=0.55 m.
(a)
μx=μy=0.3 m, σx=σy=0.55 m(b)
Low SNRγ¯=−10 dB(c)
μx=μy=1 m, σx=σy=0.55 m(d)
Medium SNRγ¯=35 dB(e)
μx=μy=1 m, σx=σy=1 m(f)
High SNRγ¯=65 dBTable 1
Normalized factork2 for the different parameters of pointing errors at the receiver when wz/Ra=4000 and Cn2=1×10−14 m−2/3.
μx (m)μy (m)σx (m)σy (m)Results of Eq. (22)Simulated results0.230.30.40.551.047731.047770.650.80.450.551.019351.019360.90.91.81.21.083331.083290.20.31.80.81.301711.30172110.550.5510.9999481.81.82.12.110.999509Figure 3
KS goodness-of-fit tests for the distributionFapproxy to approximate the distribution Fy with a 5% significance level, where μx=μy=1 m and σx=σy=0.55 m.
(a)
D=0.3 m(b)
D=0.3m, γ¯=65 dBIn Figures2(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 wz, smaller jitter variance, and nonzero boresight displacement lead to a better approximation of fapproyx=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 wz is lower than the threshold and is below the critical when wz is larger than such a threshold.In Figures2(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 wz≤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 wz with increasing γ¯.The same conclusion can be also seen in Figure3, 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.


## 3.3. SAP Estimator

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


### 3.3.1. Estimator in a Noisy Situation

For theK observations y, the MLE of the unknown parameters θ=k1,α,β in a noisy situation is obtained by maximizing the approximate log-likelihood function as(19)Ly;α,β≈Lappro1y;θ≈Klnk1+1K∑k=0K−1c1yk+32ln2π+1K∑k=0K−1κ1yk,h~,θ~,h~p−12K∑k=0K−1lnκ1″yk,h~,θ~,h~p,where h~,θ~,h~p denotes the saddle point with respect to yk and the current estimated parameters α^ and β^.


### 3.3.2. 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-likelihood function(20)Ly;α,β≈Lappro2y;θ≈Klnk2+1K∑k=0K−1C2yk+ln2π+1K∑k=0K−1K2yk,h~,θ~,h~p−12K∑k=0K−1lnk2″yk,h~,θ~,h~p,where θ~,h~p represents the saddle point given yk;α^ and β^ are the current parameter estimates.The normalized factorsk1 and k2 are derived according to Eq. (18) given the estimated parameters α^ and β^. It must be emphasized that the θ~ 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 approximation 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 Lappro1θ or Lappro2θ 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.Proposition 1.
AsA0⟶1 and wz⟶∞, h~p closes to 1.Proof.
A proof is presented in the Appendix.☐Corollary 2.
WhenA0⟶1 and wz is large enough, we have(21)h~p=1−ε=1−2μx2+μy2wzeq2.Proof.
A proof is given in the Appendix.☐Theorem 3.
WhenA0⟶1 and wz is large enough, the normalized factor k2 only depends on the parameters μx, μy, σx, and σy. This is carried out by(22)k2=μy2σy2+μx2σx2μx2σy2+μy2σx2μx2+μy2σxσy.
A proof is given in the Appendix.Table1 presents the results for the normalized factor k2 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×104. 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).


### 3.3.3. 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(23)Ehn=μn≈φmod2Amodnφmod2+n1αβΓα+nΓβ+nΓαΓβn,where Ehn denotes the expectation of hn. The notations φmod2 and Amod are, respectively, defined as(24)φmod2=wzeq22σmod,Amod=A0exp1φmod2−12φx2−12φy2−μx22σx2φx2−μy22σy2φy2.Hence, the initial values for theα and β in a noiseless situation are computed by ([19], Eq. (9))(25)α∧0≈cd−12c2d2−4d,β∧0≈cd+12c2d2−4d,where c and d are given by(26)c=n2Ey2φmod2+2/φmod2Amod2−Eyn+1/Eynφmod2+n+1/φmod2+nAmod−n2−1n2−n,d=nEy2φmod2+2/φmod2Amod2−Eyn+1/Eynφmod2+n+1/φmod2+nAmod−n−1n−n2.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(27)c=4Ey2−σn2φmod2+2/φmod2Amod2−Ey3−3σn2φmod2Amod/φmod2+1φmod2+3/Amod/Ey2−σn2φmod2+2−32d=2Ey2−σn2φmod2+2/φmod2Amod2−Ey3−3σn2φmod2Amod/φmod2+1φmod2+3/Amod/Ey2−σn2φmod2+2−1−2Substituting 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(28)α∧0=β∧0=1Ey2−σ2φmod2+2/φmod2Amod2−1,which is reasonable for some specified circumstances, such as FSO systems with small received apertures that operate in weak turbulence conditions according to Figure1.The estimation process of the SAP estimator is summarized in Figure4. 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.Figure 4
Flow chart of the SAP estimator for the Gamma-Gamma shaping parametersα and β.


## 3.3.1. Estimator in a Noisy Situation

For theK observations y, the MLE of the unknown parameters θ=k1,α,β in a noisy situation is obtained by maximizing the approximate log-likelihood function as(19)Ly;α,β≈Lappro1y;θ≈Klnk1+1K∑k=0K−1c1yk+32ln2π+1K∑k=0K−1κ1yk,h~,θ~,h~p−12K∑k=0K−1lnκ1″yk,h~,θ~,h~p,where h~,θ~,h~p denotes the saddle point with respect to yk and the current estimated parameters α^ and β^.


## 3.3.2. 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-likelihood function(20)Ly;α,β≈Lappro2y;θ≈Klnk2+1K∑k=0K−1C2yk+ln2π+1K∑k=0K−1K2yk,h~,θ~,h~p−12K∑k=0K−1lnk2″yk,h~,θ~,h~p,where θ~,h~p represents the saddle point given yk;α^ and β^ are the current parameter estimates.The normalized factorsk1 and k2 are derived according to Eq. (18) given the estimated parameters α^ and β^. It must be emphasized that the θ~ 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 approximation 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 Lappro1θ or Lappro2θ 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.Proposition 1.
AsA0⟶1 and wz⟶∞, h~p closes to 1.Proof.
A proof is presented in the Appendix.☐Corollary 2.
WhenA0⟶1 and wz is large enough, we have(21)h~p=1−ε=1−2μx2+μy2wzeq2.Proof.
A proof is given in the Appendix.☐Theorem 3.
WhenA0⟶1 and wz is large enough, the normalized factor k2 only depends on the parameters μx, μy, σx, and σy. This is carried out by(22)k2=μy2σy2+μx2σx2μx2σy2+μy2σx2μx2+μy2σxσy.
A proof is given in the Appendix.Table1 presents the results for the normalized factor k2 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×104. 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).


## 3.3.3. 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(23)Ehn=μn≈φmod2Amodnφmod2+n1αβΓα+nΓβ+nΓαΓβn,where Ehn denotes the expectation of hn. The notations φmod2 and Amod are, respectively, defined as(24)φmod2=wzeq22σmod,Amod=A0exp1φmod2−12φx2−12φy2−μx22σx2φx2−μy22σy2φy2.Hence, the initial values for theα and β in a noiseless situation are computed by ([19], Eq. (9))(25)α∧0≈cd−12c2d2−4d,β∧0≈cd+12c2d2−4d,where c and d are given by(26)c=n2Ey2φmod2+2/φmod2Amod2−Eyn+1/Eynφmod2+n+1/φmod2+nAmod−n2−1n2−n,d=nEy2φmod2+2/φmod2Amod2−Eyn+1/Eynφmod2+n+1/φmod2+nAmod−n−1n−n2.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(27)c=4Ey2−σn2φmod2+2/φmod2Amod2−Ey3−3σn2φmod2Amod/φmod2+1φmod2+3/Amod/Ey2−σn2φmod2+2−32d=2Ey2−σn2φmod2+2/φmod2Amod2−Ey3−3σn2φmod2Amod/φmod2+1φmod2+3/Amod/Ey2−σn2φmod2+2−1−2Substituting 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(28)α∧0=β∧0=1Ey2−σ2φmod2+2/φmod2Amod2−1,which is reasonable for some specified circumstances, such as FSO systems with small received apertures that operate in weak turbulence conditions according to Figure1.The estimation process of the SAP estimator is summarized in Figure4. 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.Figure 4
Flow chart of the SAP estimator for the Gamma-Gamma shaping parametersα and β.


## 4. Numerical Results

To investigate the SAP estimator performance, the MSE ofθ^ is studied using a Monte Carlo simulation, where MSE θ^=Varθ^+Eθ∧−θ2 are defined for the component-wise operation. Varθ^ denotes the variance, i.e., Varθ^=1/K−1∑k=0K−1θ∧k−θ∧¯2,θ^¯ is the sample mean of the estimator. The NMSE is an important performance metric to measure the quality of an estimator and is defined as NMSE=MSEθ^/θ2 [38]. In the simulation, the size of data samples K for D=0.1 m is 2×103 while K is set to 7×103 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 Figure5, we present the simulated MSE and NMSE performance of α and β for different Cn2 values. It can be observed that increasing the strength of atmospheric turbulence Cn2 decreases the MSE performance of α while it does not substantially change the NMSE of α. The MSE and NMSE performances of β and k2 are not sensitive to the value of Cn2, and they remain relatively flat as Cn2 increases. In order to investigate the impacts of wz on the estimator performance, we also show the MSE performance for different wz, i.e., wz=5 m and wz=10 m in Figure 5. By comparing these curves, we demonstrate there is no significant difference in their MSE performance for α, β, and k2. This may be due to their nearly equal maximum CDF errors under the SAP representation according to Figure 2(c).Figure 5
MSE and NMSE performance of the SAP estimator in noiseless channel conditions whenμx=μy=0.3 m, σx=σy=0.55 m, and D=0.1 m.In Figure6, 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.Figure 6
MSE and NMSE performance of the SAP estimator in noisy channel conditions whenμx=μy=1 m, σx=σy=0.55 m, wz=10 m, D=0.1 m, and γ¯=95 dB, where dashed lines represent the performance of SAP estimator for noiseless cases.Also, the estimation performance results for the noiseless case are included in Figure6 as a benchmark. We find that the MSE performances for α and β are close when Cn2≥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 k2 can be achieved for jitters with a larger mean at the receiver, and this is.Figure7 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 Cn2≤1×10−14 m−2/3. These results suggest that the performance of the SAP estimators is sensitive to the value of Cn2 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 k2.Figure 7
MSE and NMSE performance of the SAP estimator in noiseless channel conditions whenμx=μy=1 m, σx=σy=0.55 m, wz=10 m, and D=0.2 m.In conclusion, using the Monte Carlo simulation, the derived results for the MSE and NMSE performance in Figures5–7 are indicative of producing accurate estimations for α and β.


## 5. 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 conditions. 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.



---

*Source: 1301878-2021-07-17.xml*