Parameter Estimation of the Lognormal-Rician Channel Model Using Saddlepoint Approximation

In this article, the challenges of effective channel estimation for the lognormal-Rician turbulence model are addressed. We present a novel maximum likelihood estimation algorithm involving a saddlepoint approximation (SAP) method to estimate the shaping parameters of the lognormal-Rician distribution. An additional parameter <inline-formula> <tex-math notation="LaTeX">$k$ </tex-math></inline-formula> needs to be estimated in addition to <inline-formula> <tex-math notation="LaTeX">$r$ </tex-math></inline-formula> and <inline-formula> <tex-math notation="LaTeX">$\sigma _{z}^{2}$ </tex-math></inline-formula> under the SAP representation. The accuracy of the proposed estimator is investigated by using the mean square error and normalized mean square error. The simulated results show that the proposed estimator exhibits satisfactory performance over a wide range of turbulence conditions, and it can be easily applied to both noiseless and noisy situations. The effect of the estimated shaping parameters errors on the bit error rate (BER) for the on-off key (OOK) modulation is also investigated; it is shown that the BER performance derived with the SAP estimator becomes closer to the system performance with perfect shaping parameters as <inline-formula> <tex-math notation="LaTeX">$r$ </tex-math></inline-formula> increases. In particular, we present a qualitative comparison between the proposed SAP estimator and other algorithms available in the literature.


I. INTRODUCTION
Free space optical (FSO) communication, also known as outdoor optical wireless communication, has recently experienced extraordinary advances. High-speed data transmission can be provided by FSO communication in terms of its large available bandwidth. In addition, there are additional advantages of FSO communication over radio frequency (RF) communication including easy deployment, greater security, and unregulated spectrum resources [1]. Moreover, FSO communication has been intensively considered as a complementary technology to the RF communication in many application scenarios, which is known as hybrid (RF-FSO) communication [2]- [5]. Nevertheless, the performance of FSO links can be severely degraded by environment effects such as aerosol scattering, building-sway, and turbulence-induced scintillation. Hence, in order to improve the reliable operation of FSO systems under a wide range of weather conditions, various types of processing techniques have been proposed to mitigate the adverse effects of the atmosphere (e.g., [6]- [9]).
Many statistical models have been proposed to characterize FSO fading channels, and they are used with a specific range The associate editor coordinating the review of this manuscript and approving it for publication was Zilong Liu . of turbulence conditions. For weak turbulence, the lognormal distribution is generally accepted. It is derived from the method of first Rytov approximation [10]. The negative exponential distribution is used to characterize the limiting case of saturated scintillation [11]. The K-distribution provides a good fit with experimental results in the very strong refractive turbulence region [12]. Furthermore, there also exist several universal models for matching a broad range of turbulence conditions, including the I-K, Gamma-Gamma, and lognormal-Rician distributions. Among these models, the lognormal-Rician distribution provides an excellent agreement with experimental data [13] and is more accurate than the Gamma-Gamma distribution, especially under a spherical wave assumption [14]. The I-K distribution is a generalized form of the K distribution to cover weak fluctuation regimes, in which the K distribution is not theoretically applicable [15], [16]. However, both the K and the I-K models have limitations as probability density function (PDF) models for the irradiance in extended turbulence [17]. In this correspondence, it is crucial to study some properties about the lognormal-Rician distribution and the problem of estimating the shaping parameters of the fading channel is considered. Furthermore, we are interested in obtaining the bit error rate (BER) performance of a wireless communication system under a specified set of channel conditions. The estimated shaping parameters can be used to demodulate signals when the noncoherent detection is employed at the receiver side. Interested readers can refer to [18]- [20] for more details.
Given the complicated integral form of the lognormal-Rician distribution, one can expect that it is infeasible to estimate the shaping parameters by the method of traditional maximum likelihood estimation (MLE). To address this concern, alternative approaches have been proposed. An estimation of the shaping parameters for this distribution has first been introduced by Churnside and Clifford [13]. They obtain parameter values based on a physical model of the turbulenceinduced scattering. However, the accuracy of this approach depends heavily on the scattering physical model, which may not be readily available. In [21], the authors apply the Hansen two-step generalized method of moments (GMM) method to estimate the shaping parameters. However, this approach requires 10 6 data samples to have satisfactory performance; this indicates that the system would exhibit latency on the order 10 6 milliseconds = 1000 seconds since quasi-static turbulence is typical millisecond timescale, which is unacceptable for practical FSO communication systems. More recently, the authors use the expectation-maximization (EM) algorithm to estimate the parameters of the lognormal-Rician distribution [22]. This estimation approach is highly accurate and it operates with relatively low quantities of data samples. However, its success requires the computation of complicated integrals, and this impedes its hardware implementation in FSO communication. In this correspondence, one should consider a more efficient estimator with a reasonable data size for practical FSO communication.
In this article, we propose the use of the MLE with a simple Newton-Raphson iterative method to estimate the parameters of the lognormal-Rician distribution. It adopts an alternative likelihood function derived from a saddlepoint approximation (SAP). The performance of such an estimator is investigated in terms of the NMSE (normalized mean squared error) and MSE (mean squared error) by Monte Carlo simulation. The simulation results show that the proposed estimation approach provides satisfactory performance over a wide range of turbulence conditions. A qualitative analysis between the proposed SAP and other methods is also shown. Also, the MSE and NMSE performance of the SAP estimator in a noisy situation is studied.

II. PRELIMINARIES
Throughout this article, bold uppercase letters (e.g, X) denote a random length-N sequence, and each element of this sequence is noted as X [·]. The symbol · denotes an estimate of a parameter.

A. LOGNORMAL-RICIAN PROBABILITY DENSITY FUNCTION
For the resulting lognormal-Rician fading channel model, the optical irradiance I can be obtained by [13] where U C is a real deterministic quantity, U G is a circular complex Gaussian random variable (RV), and χ is a real Gaussian RV. Then, we can deduce that |U C +U G | is a Rician RV, and |U C + U G | 2 follows a noncentral chi-square RV with a degree of freedom of two and exp(2χ ) is a lognormal RV. Equation (1) can be partitioned into factors that behave like a modulation process where the Rician and lognormal aspects arise respectively from small-scale and large-scale atmospheric effects. The PDF of the lognormal-Rician distribution is given by [13] f (I ; r, where z is a lognormal RV that represents exp(2χ), r denotes the coherence parameter defined as r is the variance of the logarithm of the irradiance modulation factor z, and I 0 is the zero-order modified Bessel function of the first kind. The moments of the lognormal-Rician distribution are given by The lognormal-Rician distribution is a universal model for scintillation, and it comprises other well-known distributions. For example, when the coherence parameter r → ∞, the lognormal-Rician distribution becomes the lognormal distribution, which is given by In [22], the authors relate these two empirical parameters to atmospheric conditions in terms of the Rytov variance σ 2 R indirectly.

B. SADDLEPOINT APPROXIMATION REPRESENTATION
The SAP method, initially proposed by Daniels (1954), is a valuable tool for approximating a density or mass function from its associated moment generating function or cumulant generating function [23]. Due to the highly accurate approximation it achieves, the SAP method has been widely used with great success by many authors. In [24], the authors employ the SAP to evaluate the outage probability of generalized wireless channels provided their moment generating function (MGF) M X (t) = E(e tX ) exists. Also, using the SAP technique, [25] consider the outage probability performance of maximum ratio combining (MRC) receiver in a variety of fading scenarios. Another application of SAP is shown in [26], where the bounds of the error probability in channel coding are derived. Here, we apply the SAP method to approximate complicated integrals. VOLUME 8, 2020 According to the SAP theory described in [27], we can approximate the positive function h(x) as is a cumulant generating function (CGF), and t(x) is the saddle point. For each x,t(x) satisfies ∂κ(x,t) ∂t = 0 and ∂ 2 κ(x,t) ∂t 2 < 0. As can be seen, such an approximation is exact if κ(x, t) is a quadratic function of t for each x, which indicates that a more Gaussian-like distribution results in a highly accurate approximation.
Using the representation in (5), the lognormal-Rician distribution (2) can be rewritten in the following way where k denotes a normalized factor and must also be estimated. κ(z, r, σ 2 z , I ) is as follows κ z, r, σ 2 z , I = −2 ln z + ln I 0 2 Given the CGF κ z, r, σ 2 z , I in (7), we then find the saddle point z 0 which solves the following saddle point equation using the Newton-Raphson method: By combining (6), (7), and (8), a complicated integral is transformed into an algebraic computation by saddlepoint approximation. Note that the distribution of the saddle point z 0 derived from (8) does not follow the lognormal distribution unless r → ∞; this is shown in the appendix. Furthermore, for the Gaussian noisy channel, the observed data Y is then expressed as when the all-one signals are transmitted, where I , N g represent the lognormal-Rician sample and additive zero-mean Gaussian noise sample with variance σ 2 g . For convenience, the signal-to-noise ratio (SNR) is defined as E[I ] 2 /E[N 2 ] = 1/σ 2 g . According to (9), the PDF of Y is derived in (10), as shown at the bottom of the next page.
Correspondingly, we have with κ (z, r, σ z , I , y) defined as

C. EFFECT OF THE ESTIMATION ERRORS ON THE BER
Note that such all-one signals in (9) can be considered as training symbols to estimate the shaping parameters of a fading channel. For a noncoherent detection with perfect shaping parameters, the maximal likelihood detector with onoff keying (OOK) modulation scheme can be formulated aŝ Accordingly, the maximal likelihood detector for the imperfect estimates of shaping parameters can be defined as follows where r m , σ 2 z m represent the estimated shaping parameters of r and σ 2 z for the m-th trial, respectively. In such cases, the corresponding BER averaged the number of trials M is In the next section, we give a detailed description of parameters estimation using SAP under noiseless conditions.

III. SAP FOR THE LOGNORMAL-RICIAN PARAMETER ESTIMATION
Based on the above information, three types of parameters need to be estimated under SAP representation: the normalized factor k, coherence parameter r, and variance of the lognormal distribution σ . (16) The right-hand side of (16) can be solved by the multivariate Newton-Raphson algorithm or BFGS Quasi-Newton algorithm [28]. However, the iteration of parameters can be carried out due to the independence of r and σ 2 z . Furthermore, the first (second) derivative of L approx (I; θ ) with respect to r or σ 2 z can be individually approximated by the finite difference.
The initial value of the coherence parameter r can be obtained by solving the polynomial equation [13, (9)] .
The initial value of σ 2 z can be given by [13, (10)] Recalling that f (I ) is actually a density function, we can renormalize the approximation by calculating the constant of the right side of (6) so that f (I )dI = 1. By doing this, the normalized factor k is theoretically carried out by computing in each iteration, where the Interpolation function is adopted for any two adjacent discrete values of C, and (19) can be easily solved by Mathematica software. In Fig. 1, we present the simulated normalized factor k, the absolute maximal difference between empirical cumulative distribution function (AbsMaxDiffCDF) of the random variable I , F I (·), the approximate distribution f appro (I ), F appro (·)) of r, and three kinds of σ 2 z . As can be seen in the figure, the parameter normalized factor k is around the constant 1 over a wide range channel conditions. It converges to 1 as r goes to infinity by the following theorem Theorem 1: The normalized factor k converges to 1 as r goes to infinity regardless of σ 2 z . Proof: A proof is given in the appendix. Thus, it is reasonable to set k 0 = 1 as the initial value of k. Also, we employ the Kolmogorov-Smirnov (KS) goodnessof-fit statistical tests to investigate the validity of the proposed approximation. According to [29, (8.320)], the KS test statistic is defined as  For a given significance level α = 5% and sample size N = 10 4 , the critical value T max = 0.0136 is derived [29, (8.322)] and is also included (dash-dotted line) in Fig. 1 as a benchmark. It can be observed in the figure that the approximate distribution f appro (I ) is a highly efficient approximation of the distribution of f (I ) over a wide range of turbulence conditions. The same conclusion can also be drawn based on Fig. 2, Fig. 3 and Fig. 4, which shows the approximation accuracy in noisy situations that the SNR ranges from low to high. However, by comparing the AbsMaxDiffCDF curves in these figures, we find a higher approximation accuracy can be provided in the situation of low SNR because a more Gaussian-like shaped PDF is achieved in this case. Furthermore, we show a qualitative comparison of the three types of  estimation algorithms for the noiseless situation in Table 1. Note that evaluating moments also involves the computation of the complicated integrals that are time-consuming to accomplish for hardware implementation [30], [31].

IV. NUMERICAL SIMULATION AND ANALYSIS
We investigate the estimator performance of using the saddlepoint approximation by using the MSE [32] are defined for the componentwise operation. The size of the data samples for estimating the lognormal-Rician parameters in both noiseless and noisy situations are N = 10 4 and N = 10 3 , respectively. M = 100 trials are employed to calculate the MSE and NMSE of the estimator.
In Fig. 5, we present the simulated MSE and NMSE performance of r and σ 2 z when σ 2 z = 0.25 and r ranges from 2 to 9. It can be seen that the MSE of r increases with the value of r while the MSE of σ 2 z is not sensitive. In other words, changes to the value of r have minimal effects on the estimation performance of σ 2 z . These results are consistent with the observation shown in [21] and [22]. The same conclusion can be drawn for the NMSE of r and σ 2 z as they both stay relatively flat as r increases; this conclusion is also drawn in [21]. In addition, we find that the MSE and NMSE performance of the normalized factor k both decrease with increasing values of r.  As for the MSE and NMSE performance of the normalized factor k, they both decrease gradually as σ 2 z approaches 0.3, and then tend to increase gradually. Fig. 7 presents the simulated MSE and NMSE performance for σ 2 z = 0.25 and SNR = 10 dB. It can be observed that as r increases, the performance behavior of the parameters r, σ 2 z , and k are consistent with the results shown for noiseless situations. Table. 2 compares the BER performance with OOK modulation under the perfect and imperfect shaping parameters. We demonstrate that the relative error decreases slightly with increasing r; this indicates that the BER performance derived with the SAP estimator becomes closer to the system performance with perfect shaping parameters as r increases.
In conclusion, using the Monte Carlo simulation, the derived results for the MSE and NMSE in Fig. 5 and Fig. 6 are indicative of accurate estimations of r and σ 2 z . To see the benefits of our proposed method in noiseless situation, the performance of the MLE with the EM algorithm is included as a benchmark in both figures. As can be seen, the performance of the estimation method proposed here is comparable with the EM algorithm. Fig. 7 confirms that it is feasible to use the SAP estimator for lognormal-Rician r and σ 2 z parameters in a noisy environment.

V. CONCLUSION
In this article, we present a novel method for the estimation of the parameters of the lognormal-Rician distribution. With the help of the SAP method, the proposed technique used the MLE to estimate the shaping parameters of the lognormal-Rician distribution in both noiseless and noisy conditions. The performance has been simulated in terms of the MSE and NMSE. The numerical results show that the SAP estimator can achieve satisfactory estimates over a wide range of turbulence conditions even though an additional parameter k needs to be estimated. A qualitative comparison between the SAP estimator and other algorithms is also presented which shows that the proposed estimator outperforms the others in terms of the hardware implementation.