Data-Aided SNR Estimation for Bandlimited Optical Intensity Channels

Not only for radio frequency but also for optical communication systems, knowledge of the signal-to-noise ratio (SNR) is essential, e.g., for an adaptive network, where modulation schemes and/or error correction methods should be selected according to the varying channel states. In the current paper, this topic is discussed for a bandlimited optical intensity link under the assumption that the data symbols are known to the receiver unit in form of pilot sequences. This requires a unipolar signal design regarding the symbol constellation, but also a non-negative pulse shape satisfying the Nyquist criterion is necessary. Focusing on this kind of scenario, the modified Cramer–Rao lower bound is derived, representing the theoretical limit of the error performance of the data-aided SNR estimator developed in this context. Furthermore, we derive and analyze a maximum likelihood algorithm for SNR estimation, which turns out to be particularly simple for specific values of the excess bandwidth, among them the most attractive case of minimum bandwidth occupation. Numerical results confirming the analytical work conclude the paper.


Introduction
In a series of papers recently published by the author [1][2][3], parameter estimation and synchronization for a bandlimited optical intensity link have been discussed. In this context, a unipolar waveform design is indispensable with respect to pulse shaping and symbol constellation [4,5]. Furthermore, it is most helpful that pulse shapes satisfy the Nyquist criterion, which allows for a simple detection process in the receiver unit [6].
However, not only for radio frequency (RF) but also for optical wireless communication (OWC) solutions [7][8][9][10], the relevant transmission parameters have to be recovered reliably, because otherwise subsequent receiver stages, such as error correction algorithms, cannot be operated in an efficient way [11,12]. In particular, recovery of the symbol timing is of paramount importance in this respect, since this is a prerequisite for many other estimation and detection procedures. In [1][2][3], it has been shown how this might be achieved for a bandlimited optical intensity link under different conditions, e.g., whether data are known to the receiver unit or not in the form of pilot sequences, or if the estimator or synchronizer module is to be implemented in a feedforward or feedback manner.
Usually, the estimation of the signal-to-noise ratio (SNR) requires that the symbol timing has been established before by a properly selected algorithm. It is to be recalled that knowledge of the SNR is normally needed for adaptive communication systems to select modulation and coding schemes according to the given channel conditions [13], but also powerful error correction methods-such as turbo or LDPC algorithms-need this sort of information [14]. Scanning the open literature, numerous papers are available about SNR estimation in RF channels, e.g., the frequently cited overview by Pauluzzi and Beaulieu [15], but little or no information is published for OWC systems. This has been the main motivation of the current contribution addressing data-aided SNR estimation for a bandlimited optical intensity channel. Finally, it is to be noticed that the article was prepared for a Special Issue of Feature Papers 2022 in the Communications Section of MDPI Sensors.
The rest of the paper is organized as follows: the signal and channel model for analytical and simulation work is introduced in Section 2, whereas Section 3 focuses on the derivation of the Cramer-Rao lower bound (CRLB) as the theoretical limit of the jitter performance for any algorithm discussed in the context of SNR estimation. In Section 4, we derive a maximum likelihood (ML) algorithm and analyze it in terms of mean value and variance. Numerical results are shown in Section 5, and Section 6 concludes the paper.

Signal and Channel Model
As already mentioned in the introductory section, properly selected pulse shapes and modulation schemes are necessary to satisfy the non-negativity as well as the Nyquist constraints required for a bandlimited optical intensity channel. In this respect, it is assumed that the real-valued data symbols a k , k ∈ Z, are independent and identically distributed (i.i.d.) elements of an M-ary PAM alphabet A. It makes sense to organize the alphabet such that the symbols are normalized to unit energy, i.e., E[a 2 k ] = 1, where E[·] denotes the expectation operator. Then, with η M = 1 This means that the average value is given by On the other hand, the signal at the output of the opto-electrical receiver module is obtained as where A > 0 is the channel gain, h(t) describes the pulse shape, T and τ symbolize the symbol period and the propagation delay between receiver and transmitter station, respectively. We assume that A is a constant regarding the observation interval used for estimation purposes, because variations of the channel state are normally slow enough so that fading effects need not be taken into account. As already required previously, h(t) must satisfy the non-negativity as well as the Nyquist criterion, e.g., achieved by a squared raised cosine function [1,6]. In line with the investigations carried out in [1][2][3][4][5][6], the receiver signal in (2) is also assumed to be distorted by additive white Gaussian noise (AWGN), in the following expressed by w(t), with zero mean and variance σ 2 w . In addition, we introduce the average optical power as P 0 = µ a h, where so that the average electrical SNR at the receiver can be defined as However, before being treated in further stages of operation, the signal in (2) has to pass the receiver filter q(t), whose output is given by z(t) = q(t) ⊗ r(t), where ⊗ denotes the convolutional operator. For convenient reasons, this is summarized in Figure 1. Since there exists no simple solution for a matched filter structure in the context of a bandlimited optical intensity link [6], it is suggested that q(t) exhibits a flat behavior over the spectrum occupied by the user component in (2). This straightforward approach guar- Since there exists no simple solution for a matched filter structure in the context of a bandlimited optical intensity link [6], it is suggested that q(t) exhibits a flat behavior over the spectrum occupied by the user component in (2). This straightforward approach guarantees that the waveform will not be distorted, but the price to be paid is an increased amount of noise the subsequent receiver stages have to cope with. In particular, this means that the transfer function of the filter performs a rectangular shape in the frequency domain, i.e., Q( f ) = F [q(t)] = √ T for | f | ≤ (1 + α)/T and Q( f ) = 0 elsewhere, with α as the roll-off factor (excess bandwidth) of the selected pulse shape; recall that α = 0 indicates the minimum bandwidth scenario. The related impulse response is then furnished by application of the inverse Fourier transform [16], i.e., with sinc(x) = sin(πx)/(πx). Of course, the noise signal at this flat filter output develops as n(t) = w(t) ⊗ q(t) representing a zero-mean non-white Gaussian process. Assuming in the next step that the symbol timing has been reliably recovered and corrected, e.g., by the algorithm proposed in [1], the T-spaced samples at the output of the receiver filter are obtained as

Derivation of the Log-Likelihood Function
The Cramer-Rao lower bound (CRLB) is a major figure of merit when it comes to the estimation of a parameter [17]. The reason behind this is the fact that the bound represents the theoretical limit of the jitter (error) variance of any estimator developed in this context.
According to the signal model specified previously, we have to consider the parameter vector u = (u 1 , u 2 ) = (A, σ w ). The CRLB for u i is determined by where [·] i indicates the i-th diagonal entry of the inverted Fisher information matrix (FIM) expressed by J(u). In the case that no nuisance parameter needs to be taken into account, the FIM entries are computed as with z as the given vector of observables, Λ(z;u) denotes the log-likelihood function (LLF) characterizing the communication link, and E[·] symbolizes expectation with respect to the noise model. By inspection of (8), it is clear that the computation of the CRLB requires the knowledge of the LLF describing the subject of investigation. To this end, we assume that a sequence of L receiver samples (6) forms the vector z expressed by The vector a of known data symbols specifies the pilot sequence, which is to be used for estimation purposes in the sequel, and n denotes the noise vector with covariance matrix where the entries for line i and column k of Ω are given by ω ik = sinc[2(1 + α)(i − k)] = ω ki forming this way a symmetric Toeplitz matrix [18]. As a result, the likelihood function for our estimation problem can be written as [19,20], However, instead of using u = (A, σ w ), it is easier to concentrate in the following on the average electrical SNR normalized by P 2 0 , i.e., ρ s = γ s /P 2 0 = A 2 /σ 2 w . Then, by introduction of P n = σ 2 w , we have that u = (ρ s , P n ) and the related LLF is furnished by which has been achieved by Ψ = Ω −1 as well as omitting all immaterial constants and factors not depending on u.

Modified Cramer-Rao Lower Bound
In the next step, the FIM entries are obtained by computing the second-order derivatives according to (8), the results of which have then to be averaged with respect to the noise vector n. However, this approach means that the CRLB will be a function of the selected pilot sequence a. Therefore, it is suggested to extend the averaging procedure to a as well, which creates the so-called modified Cramer-Rao lower bound (MCRLB) [21][22][23]. Doing so, we get after some algebra: Evaluating (7) for u 1 = ρ s , the corresponding MCRLB is given by Substituting (13)-(15) into (16) and scaling the result with respect to ρ 2 s , we obtain the normalized MCRLB expressed as Introducing the auxiliary terms where ψ ik is the entry of Ψ indicating line i and column k, the expected operations in (17) can be written as and Finally, by plugging (19) and (20) into (17), we have that Nevertheless, the relationship might be simplified for α ∈ 0, 1 2 , 1 , because in this case ω ik = sin c[2(1 + α)(i − k)] = 1 for i = k and zero elsewhere. This means that Ω = Ω −1 = Ψ = I L , with I L as the L-dimensional identity matrix, which means also that Ψ 0 = Ψ 2 = 1 and Ψ 1 = 0. Hence, the normalized bound boils down to

Derivation of the Estimator Algorithm
By means of the LLF in (12), we are basically in the position to derive a maximum likelihood (ML) algorithm for SNR estimation. However, the SNR parameter is composed of two ingredients-channel gain A and the noise power P n -the estimates of which are needed to compute the SNR estimate. This is simply achieved by substituting ρ s = A 2 /P n into (12), deriving the resulting LLF with respect to A and P n , equating both relationships to zero and solving them analytically. Doing this, we obtain for u = (A, P n ) and Then, by introduction of M aa = a T Ψ a, M az = z T Ψ a, and M zz = z T Ψ z, we find the estimates for channel gain and noise power in closed form: By inspection of (25) and (26), it is clear that a viable solution is only achievable for M aa > 0, i.e., a pilot sequence consisting of zero symbols only would not work. Finally, according to the invariance principle for ML estimates [24], the SNR solution is given bŷ

Probability Analysis
In this subsection, we want to derive the probability density function (PDF) of the SNR estimate in (27) and analyze it in terms of bias and variance. It turns out that this is possible in closed form for Ψ = I L , i.e., α ∈ 0, 1 2 , 1 , whereas for other values of α, it is verified in Section 5 that the analytical results achieved with Ψ = I L are very close to the true ones obtained by numerical means.
By plugging Ψ = I L into (25), the estimate for the channel gain develops aŝ which means thatÂ is a zero-mean Gaussian variate. Computing the variance of the latter, we have to consider that the noise samples n k are zero-mean and i.i.d. in case that Ψ = I L . Therefore, where σ 2 n = E[n 2 k ] = 2(1 + α) σ 2 w . The related PDF is then straightforwardly given by On the other hand, Y =Â 2 corresponds to a non-central Gamma variate [19] characterized by the distribution If we consider in the next step the estimate of the noise power for Ψ = I L , we havê With z = A · a + n andÂ determined by (28), the numerator in (32) simplifies to which represents a central Gamma variate with variance σ 2 n and L-1 degrees of freedom [20,25]. Hence, by introduction of X =P n , the PDF of (32) can be written as where σ 2 x = σ 2 n 2(1+α)L = σ 2 w L . Employing in the following the PDFs in (31) and (34), the distribution of the SNR estimate is determined by (A2) derived in the Appendix A, i.e., Regarding (A3) and (A4), the parameters β, λ, and K 0 are functions of the true SNR value denoted by ρ s , the roll-off factor α, the observation length L, as well as the selected pilot sequence a.
By means of the relationships (A9) and (A10) detailed in the Appendix A, we can specify the first-and second-order moments of (35) as follows: Therefore, bias and variance ofρ s , normalized by ρ s and ρ 2 s , respectively, are given by and Via M aa , it is obvious that (38) and (39) depend on the selected pilot sequence. In order to avoid this drawback, we could average the relationships with respect to a. The problem in this context is that there exists no closed form solution. A way out of this dilemma is Jensen's inequality [26] (Appendix 1B), which provides us with where κ a denotes the symbol kurtosis of the PAM alphabet, which is given by By taking into account the auxiliary results in (40) and (41), we finally obtain and as lower bounds for the relationships in (38) and (39), respectively.

Numerical Results
The analytical results achieved for SNR estimation in Sections 3 and 4 will be verified by Monte Carlo (MC) simulations. In the following, the former are indicated by lines, whereas the latter are shown by markers. Each point in the diagrams below has been obtained by averaging a number of 10 5 estimates, which turned out to be large enough to verify the analytical results with sufficient accuracy.
Assuming a 4-PAM constellation operated with ρ s = 0 dB and a roll-off factor α ∈ {0.0, 1.0}, Figure 2 illustrates the evolution of the normalized bias as a function of the observation length L. It is to be noticed that the lines in different colors represent the lower bound given by (43); verified by simulation results, we observe that the lower limit is very tight over the full range of L. We observe that the bias decreases rapidly with increasing values of L, which is also confirmed by (43). In addition, the diagram depicts the results obtained for 16-PAM, ρ s = 10 dB, and α ∈ {0.2, 0.8}. In the strict sense, the relationship in (43) applies only to values of α ∈ {0.0, 0.5, 1.0}, but the 16-PAM scenario in Figure 2 demonstrates that it represents also a very good approximation for other values of the excess bandwidth. and α α ρ ρ κ ρ as lower bounds for the relationships in (38) and (39), respectively.

Numerical Results
The analytical results achieved for SNR estimation in Sections 3 and 4 will be verified by Monte Carlo (MC) simulations. In the following, the former are indicated by lines, whereas the latter are shown by markers. Each point in the diagrams below has been obtained by averaging a number of 10 5 estimates, which turned out to be large enough to verify the analytical results with sufficient accuracy.
Assuming a 4-PAM constellation operated with ρs = 0 dB and a roll-off factor α ∈ {0.0, 1.0}, Figure 2 illustrates the evolution of the normalized bias as a function of the observation length L. It is to be noticed that the lines in different colors represent the lower bound given by (43); verified by simulation results, we observe that the lower limit is very tight over the full range of L. We observe that the bias decreases rapidly with increasing values of L, which is also confirmed by (43). In addition, the diagram depicts the results obtained   The evolution of the normalized bias has been simulated and verified for modulation schemes other than 4-PAM and 16-PAM, e.g., 2-PAM and 8-PAM, as well as for roll-off factors different to those exemplified in Figure 2. It turned out that the bias of the estimator algorithm is reflected accurately enough by the formula in (43), disappearing for very large values of L irrespective of the selected values of M or α.
Using a 4-PAM scheme with L = 10 and the same roll-off factors as before, Figure 3 illustrates the evolution of the normalized variance as a function of the true SNR value in dB. For comparison purposes, the normalized MCRLB expressed by (22) is shown in dashed style. We observe that the latter is fairly loose for such small observation windows, whereas the lower bound of the variance in (44) appears to be very tight as confirmed by simulation results, in particular at larger SNR values. However, the diagram illustrates also that the MCRLB is more and more approximated by the jitter variance of the related estimator algorithm, when we increase the observation length in Figure 3, verified for 16-PAM, L = 100, and α ∈ {0.2, 0.8}; it is to be recalled that for α / ∈ {0.0, 0.5, 1.0}, the NMCRLB is furnished by (21). Finally, we see that the MC output is very close to (44) over the full SNR range, although the relationship is, in a strict sense, only applicable to roll-off factors α ∈ {0.0, 0.5, 1.0}. These observations also hold true for modulation schemes and roll-off factors other than those used in Figure 3; especially, one can see that for L 1 and ρ s 1, the normalized variance is simply given by 2/L. is furnished by (21). Finally, we see that the MC output is very close to (44) over the full SNR range, although the relationship is, in a strict sense, only applicable to roll-off factors   {0.0, 0.5, 1.0}. These observations also hold true for modulation schemes and roll-off factors other than those used in Figure 3; especially, one can see that for L  1 and s  1, the normalized variance is simply given by 2/L.

Concluding Remarks
Assuming a data-aided situation, i.e., data symbols are known to the receiver in the form of a pilot sequence, SNR estimation for a bandlimited optical intensity link has been investigated in the current paper. This requires a signal design achieved by an M-ary PAM scheme and a non-negative pulse shape also satisfying the Nyquist criterion. By means of a flat receiver filter, it is avoided that the waveforms of the user signal are distorted, but the price to be paid is an additional amount of noise which the subsequent receiver stages are suffering from.
Conditioned on reliable recovery and correction of the symbol timing, the modified CRLB could be derived as the theoretical limit of the jitter variance produced by the SNR estimator developed in the context of this paper. With respect to the latter, an ML solution has been obtained in closed form, which turned out to be particularly simple from a computational point of view for specific values of excess bandwidth, among them being the minimum bandwidth scenario. For these values, the analytical relationships for bias and jitter variance have been obtained in closed form as well.

Concluding Remarks
Assuming a data-aided situation, i.e., data symbols are known to the receiver in the form of a pilot sequence, SNR estimation for a bandlimited optical intensity link has been investigated in the current paper. This requires a signal design achieved by an M-ary PAM scheme and a non-negative pulse shape also satisfying the Nyquist criterion. By means of a flat receiver filter, it is avoided that the waveforms of the user signal are distorted, but the price to be paid is an additional amount of noise which the subsequent receiver stages are suffering from.
Conditioned on reliable recovery and correction of the symbol timing, the modified CRLB could be derived as the theoretical limit of the jitter variance produced by the SNR estimator developed in the context of this paper. With respect to the latter, an ML solution has been obtained in closed form, which turned out to be particularly simple from a computational point of view for specific values of excess bandwidth, among them being the minimum bandwidth scenario. For these values, the analytical relationships for bias and jitter variance have been obtained in closed form as well.
Verified by simulation results, it could be shown that-irrespective of the chosen PAM constellation and the value of the excess bandwidth-the bias effect vanishes more and more with increasing values of the true SNR value and the observation length L the link is operated with. This is also confirmed in view of jitter performance insofar as the CRLB is successively approached by increasing values of L.
Substituting in the sequel the PDFs given by (31) and (34), we obtain by means of [28] (3.462/1) and (9.240) after some lengthy but straightforward manipulations, where 1 F 1 (·) denotes the confluent hypergeometric (Kummer) function [29] and the parameters β, λ, and K 0 are specified as follows: Deriving the m-th order moment of (A2), we first express the confluent hypergeometric function by its Meijer G-equivalent [30]