Non-Data-Aided SNR Estimation for Bandlimited Optical Intensity Channels

Powerful and reliable estimation of transmission parameters is an indispensable task in each receiver unit—not only for radio frequency, but also for optical wireless communication systems. In this context, the signal-to-noise ratio (SNR) plays an eminent role, especially for adaptive scenarios. Assuming a bandlimited optical intensity channel, which requires a unipolar waveform design, an algorithm for SNR estimation is developed in this paper, which requires no knowledge of the transmitted data. This non-data-aided approach benefits to a great extent from the fact that very long observation windows of payload symbols might be used for the estimation process to increase the accuracy of the result; this is in striking contrast to a data-aided approach based on pilot symbols reducing the spectral efficiency of a communication link. Since maximum likelihood, moment-based or decision-directed algorithms are not considered for complexity and performance reasons, an expectation-maximization solution is introduced whose error performance is close to the Cramer-Rao lower bound as the theoretical limit, which has been derived as well.


Introduction
When comparing radio frequency (RF) to optical wireless communication (OWC) techniques, the advantages of the latter are well known: there are no regulatory and license issues, they are rather inexpensive and easy to deploy, have extremely high throughput and cause no problems with data security, just to mention the most significant aspects in this context [1][2][3][4].
However, not only for RF but also for OWC solutions, the relevant transmission parameters have to be recovered by powerful algorithms, because otherwise subsequent receiver stages, e.g., detectors or error correction algorithms cannot be reliably operated [5,6]. Of course, in case that bandlimited optical intensity solutions are envisaged, a unipolar waveform design is indispensable with respect to pulse shaping and symbol constellation. Investigated in [7,8] for a PAM scheme and root-raised cosines as typical pulse shapes used for RF solutions, this is simply achieved by a suitably selected bias or offset signal. Unfortunately, such concepts are not very efficient in terms of power and energy if no harvesting is implemented. Hence, squared raised cosine and double jump functions have been suggested in [9] as viable alternatives.
Nevertheless, focusing on the pulse shapes proposed in [9] also means that recovery methods developed in the RF context are not automatically applicable. Of course, synchronization of carrier frequency and phase need not be considered in case of intensity modulation, whereas the recovery of the symbol timing is still of paramount importance since this is a pre-requisite for many other estimation and detection procedures. This problem has been tackled in a couple of papers recently published by the author in [10][11][12].
Apart from symbol timing and clock recovery, some knowledge about the signal-tonoise ratio (SNR) is equally important for the reliable transmission of data, e.g., for adaptive communication systems to select modulation and coding schemes according to the given

Signal and Channel Model
Of course, for the NDA scenario to be investigated in this contribution, we are working with the same signal and channel model used in the companion paper addressing a DA situation [17]. For clarity and readability reasons, this model is briefly recapitulated in the sequel.
Due to the unipolar waveform design mentioned previously, it is assumed that the data symbols a k , k ∈ Z, are independent and identically distributed (i.i.d.) elements of an M-ary PAM alphabet A. In this context, it makes sense to normalize them to unit energy, i.e., E[a 2 k ] = 1, where E[·] denotes the expectation operator. Therefore, by definition of η M = 1 6 (M − 1)(2M − 1), we obtain a k ∈ A = 1 √ η M {0, 1, . . . , M − 1}. As a consequence, the average value is given by If h(t) describes the pulse shape satisfying the non-negativity as well as the Nyquist constraint, the signal at the output of the opto-electrical receiver unit can be expressed by where T is the symbol period and τ indicates the propagation delay between transmitter and receiver. The channel gain A > 0 is assumed to be a constant, which is justified by the fact that the coherence time of fading events is usually much larger than the observation window needed for estimation purposes. In line with the previous publication about this topic in [17], the noise component in (2) is assumed to be a zero-mean white Gaussian process with variance σ 2 w . In addition, if we introduce the average optical power as P 0 = µ a h, where Nevertheless, before r(t) is processed in further receiver stages, e.g., in the SNR estimator to be developed in the sequel, it must be filtered appropriately. Assuming an impulse response q(t), the filter output is determined by z(t) = q(t) ⊗ r(t), where ⊗ denotes the convolutional operator. For convenient reasons, the signal model used for SNR estimation is illustrated in Figure 1.
the average electrical SNR at the receiver can be defined as Nevertheless, before r(t) is processed in further receiver stages, e.g., in the SNR estimator to be developed in the sequel, it must be filtered appropriately. Assuming an impulse response q(t), the filter output is determined by = ⊗ ( ) ( ) ( ) z t q t r t , where ⊗ denotes the convolutional operator. For convenient reasons, the signal model used for SNR estimation is illustrated in Figure 1. It has been proved in [9] that there exists no simple solution for a matched receiver filter in case of a bandlimited optical intensity link. Hence, it is suggested to focus on a solution performing a rectangular shape over the complete spectrum occupied by the user component in (2). By application of the Fourier-transform [18], we have that ( elsewhere, with α as the roll-off factor (excess bandwidth) of the selected pulse shape. In this case, the signal parts of r(t) and z(t) are the same, whereas the noise component is determined by n(t) = w(t) ⊗ q(t) representing a zero-mean non-white Gaussian process. Under the assumption that the symbol timing has been reliably recovered and corrected, e.g., by one of the algorithms proposed in [10][11][12], the T-spaced samples at the output of the receiver filter are obtained as i k w n n i k , and sinc(x) = π π sin( ) ( ) x x .

Log-Likelihood Function and Fisher Information Matrix
For parameter estimation, in general, the Cramer-Rao lower bound (CRLB) is a major figure of merit [19]. It turns out to be most helpful for comparison purposes, since the bound represents the theoretical limit of the error performance of any algorithm developed in this context. By detailed inspection of (4), it is obvious that the average electrical SNR is a function of the channel gain and the standard deviation of the noise component in (2), A and σw, respectively, whereas P0 might be considered as a constant factor depending on the PAM scheme and the selected pulse shape. Therefore, it makes sense to focus on the SNR . On top of that, it is assumed that L observables given by (5) are available for the estimation procedure, which might be elegantly expressed in vector form: It has been proved in [9] that there exists no simple solution for a matched receiver filter in case of a bandlimited optical intensity link. Hence, it is suggested to focus on a solution performing a rectangular shape over the complete spectrum occupied by the user component in (2). By application of the Fourier-transform [18], we have that and Q( f ) = 0 elsewhere, with α as the roll-off factor (excess bandwidth) of the selected pulse shape. In this case, the signal parts of r(t) and z(t) are the same, whereas the noise component is determined by n(t) = w(t) ⊗ q(t) representing a zero-mean non-white Gaussian process. Under the assumption that the symbol timing has been reliably recovered and corrected, e.g., by one of the algorithms proposed in [10][11][12], the T-spaced samples at the output of the receiver filter are obtained as , and sinc(x) = sin(πx)/(πx).

Log-Likelihood Function and Fisher Information Matrix
For parameter estimation, in general, the Cramer-Rao lower bound (CRLB) is a major figure of merit [19]. It turns out to be most helpful for comparison purposes, since the bound represents the theoretical limit of the error performance of any algorithm developed in this context.
By detailed inspection of (4), it is obvious that the average electrical SNR is a function of the channel gain and the standard deviation of the noise component in (2), A and σ w , respectively, whereas P 0 might be considered as a constant factor depending on the PAM scheme and the selected pulse shape. Therefore, it makes sense to focus on the SNR normalized by P 2 0 , henceforth denoted by ρ s = γ s /P 2 0 , and to employ ρ s and P n = σ 2 w , instead of A and σ w , as elements constituting the parameter vector u to be estimated in the sequel, i.e., u = (u 1 , u 2 ) = (ρ s , P n ). On top of that, it is assumed that L observables given by (5) are available for the estimation procedure, which might be elegantly expressed in vector form: It is to be recalled that the noise samples in n are not independent. The related covariance matrix is given by E[n · n T ] = 2(1 + α) σ 2 w Ω, where Ω describes a symmetric Toeplitz matrix [20] with entries ω ik = sinc[2(1 + α)(i − k)] for line i and column k. Conditioned on the knowledge of the data sequence a and the unknown but deterministic parameter vector u, the log-likelihood function (LLF) characterizing the estimation problem has been derived in [17] as where Ψ = Ω −1 . However, the computation of the CRLB for NDA estimation requires that the LLF does not depend on a, which is achieved by averaging the related likelihood function, i.e., Pr(z|a; u) = e Λ(z|a;u) , with respect to a ∈ A L , where A L denotes the L-dimensional symbol space spanned by L i.i.d. elements of A. Therefore, we have that As a consequence, the entries of the Fisher information matrix (FIM) are obtained as where E z [·] denotes expectation with respect to z, i.e., an averaging procedure with respect to a and n, which is only possible by numerical means. According to theory, the CRLB for parameter u i is formally given by the i-th diagonal entry of the inverted FIM. Since we are only interested in the estimation of ρ s , representing the first element of u in our definition introduced before, the corresponding CRLB develops as (10)

Low-Complexity Solution
From the complexity point of view, it is clear that the LLF in (8) is the most demanding ingredient for the computation of the CRLB in (10). This is mainly due to the averaging procedure, which is in the order of M L operations. Even smaller values of M and L are challenging in this respect, but for values of L between 100 and 1000, which are typical for an NDA scenario, the computational load would be intractable. Luckily, for some values of the excess bandwidth, in particular for α ∈ 0, 1 2 , 1 , it turns out that Ω boils down to an L-dimensional identity matrix I L , which means that Ψ = Ω −1 = I L . In consequence, the likelihood function in (8) can be re-organized as Because of Ψ = I L the elements of z might be considered as statistically independent entries. By taking into account that the symbols a i ∈ A are i.i.d., we just obtain where Finally, the averaging of Pr(z|a; u) with respect to a ∈ A L is simply achieved, if we exchange in (12) the order of sum and product, i.e., resulting in a computational complexity in the order of M × L, which is much less compared to M L . In the next step, with Λ(z; u) = log Pr(z; u), the first-order derivatives in (9) are expressed by and where and For the simplified scenario, the computation of FIM entries and CRLB does not differ from the general case such that the relationships in (9) and (10) might be used in the same way.

Asymptotic Scenario
For increasing values of M, the density of the PAM symbols a i will increase accordingly, when we assume that their average energy is normalized to unity, i.e., E[a 2 k ] = 1. Hence, in case that M → ∞, the symbol alphabet A turns out to be equally distributed between 0 and √ 3. Applying the framework developed previously to compute the CRLB for such a scenario, it is clear that the sums over a i ∈ A in (15) and (16) have to be replaced by the related integrals, i.e., where m ∈ {0, 1, 2}. By taking into account the solutions for (19) computed in Appendix A, the first-order derivatives for the FIM entries in (9) are after some lengthy but straightforward manipulations given by It is to be recalled that the simplified computation of the CRLB applies in the strict sense only to values of α ∈ 0, 1 2 , 1 , where Ψ = Ω −1 = I L . However, since the entries of Ω are given by ω ik = sinc[2(1 + α)(i − k)], it is clear that the off-diagonal elements of the matrix are rapidly decaying for α / ∈ 0, 1 2 , 1 , which means that Ω and Ψ might be approximated by an identity matrix of the same dimension such that the simplification would be applicable. This results in a tight approximation of the true bound confirmed by numerical results in Section 5.

Expectation-Maximization Estimator
Since the normalized SNR value is given by ρ s = A 2 /P n , it is clear that any estimator algorithm must provide the estimates of channel gain as well as noise power, which means that we are focusing in the sequel on a parameter vector u = (A, P n ). Formally, a maximum likelihood solution for u is easily obtained by using the LLF in (8), i.e.,û = argmaxũΛ(z;ũ). However, this problem cannot be solved in closed form so that we must resort to numerical methods, e.g., the iterative Newton-Raphson procedure [21]. Apart from troubles in terms of initialization, convergence and stability, it is the computational complexity which prohibits this approach, even if the simplified variant with Ψ = I L would be envisaged.
Alternatively, a rather simple solution based on first-and second-order moments given by E[z k ] and E[z 2 k ], respectively, delivered reliable estimates only in the very low SNR range. On the other hand, for an algorithm based on symbol decisions [15], useful results were solely achievable for very large SNRs. As a consequence, an expectation-maximization (EM) estimator [22][23][24] is proposed, whose error variance turned out to be close to the CRLB over a wide SNR range as will be demonstrated in Section 5.
The EM algorithm is an iterative procedure using in step l the parameter estimates calculated in step l -1. The conditional LLF for this approach is for Ψ = I L expressed by This is very similar to (7), but u is replaced by the trial valueũ used for optimization purposes; a as well as a T a are substituted by the corresponding soft decisions in step l -1, henceforth denoted by whereas .. a (l−1) , representing a scalar, is a finite double sum given by ..
By inspection of (23) and (24), we observe that the soft decisions are characterized by an averaging of the symbols a i ∈ A via the a posteriori probabilities [25]: ] develops as Since the denominator in (25) does not depend on a i , it might be replaced by a constant including also the a priori probability Pr[a i ∈ A|û Finally, computing the first-order derivatives of (22) with respect to A as well as P n and equating them to zero for A =Â (l) and P n =P (l) n , the parameter estimates for step l are achieved in closed form byÂ Nevertheless, the algorithm has to be initialized by appropriately selected values for the probabilities in (25). Since the symbols a i are assumed to be i.i.d., it makes sense to start the EM algorithm with P It is not difficult to see that the iterative step of the EM algorithm has a computational complexity in the order of M × L additive and multiplicative operations, only relationship (26) requires the evaluation of square root and exponential functions, which might be elegantly handled via look-up tables. For the numerical results in Section 5, the iterative procedure is organized such that it stops as soon as the relative error between two successive estimates achieves a predefined value of 10 −3 or when a maximum number of 10 3 iterations is achieved; by extensive tests it turned out that this would be a good compromise between complexity and accuracy. Furthermore, it is to be remembered that the algorithm applies in the strict sense only to Ψ = I L , i.e., α ∈ 0, 1 2 , 1 . However, since good results are obtained for other values of α as well, it makes sense to employ the EM algorithm in these cases as well.

Numerical Results
In the following, the EM algorithm developed previously will be verified by numerical means in terms of jitter (error) performance and bias. For comparison purposes, the CRLB is included to the related diagrams, which show also the modified Cramer-Rao lower bound (MCRLB). Derived in closed form in [17], the MCRLB is much simpler than the CRLB, but it is in general less tight [26][27][28][29], i.e., MCRLB(ρ s ) ≤ CRLB(ρ s ), in particular at lower SNRs as demonstrated subsequently.
For 2-PAM and 4-PAM signals operated with α = 0 and two different observation lengths, L = 100 and 1000, Figure 3 illustrates the evolution of the error performance as a function of the SNR value; for convenient reasons, error performance and theoretical limits are normalized by ρ 2 s . By detailed inspection, we observe that • For larger SNRs, the CRLB (dashed line) approaches the corresponding MCRLB (solid line) irrespective of the selected modulation scheme and the value of L. • For very low SNRs, the ratio between CRLB and MCRLB seems to approach a small but non-negligible constant, which decreases somewhat by increasing values of M.

•
In the medium SNR range, we see a significant difference between MCRLB and CRLB whose maximum grows with increasing values of M and which moves to larger SNR values.

•
For medium-to-low SNRs and L = 100, the error performance of the EM estimator, indicated by markers in different style, is characterized by a considerable difference to the CRLB, which shrinks more and more with increasing values of the SNR. This degradation is basically explained by the fact that the algorithm performs a bias effect evolving in the same way, which is depicted in Figure 2 (in this case, the dashed lines do not correspond to an analytical relationship; they are due to an interpolation procedure in order to achieve a better readability of these numerical results). This drawback might be circumvented with larger observation windows, in   Figure 4 includes also the evolution of the CRLB for M → ∞ as it has been derived in Section 4. One can see that the theoretical limit for M → ∞ is close to that computed for M = 16 as long as the latter does not start to approach the MCRLB, whereas for M → ∞ it continues to increase with increasing SNRs. This property suggests that NDA estimation of the SNR becomes more and more problematic in case we increase the order of the selected PAM scheme. degradation is basically explained by the fact that the algorithm performs a bias effect evolving in the same way, which is depicted in Figure 3 (in this case, the dashed lines do not correspond to an analytical relationship; they are due to an interpolation procedure in order to achieve a better readability of these numerical results). This drawback might be circumvented with larger observation windows, in Figures 2 and 3 exemplified by L = 1000.   Figure 4 includes also the evolution of the CRLB for M → ∞ as it has been derived in Section 4. One can see that the theoretical limit for M → ∞ is close to that computed for M = 16 as long as the latter does not start to approach the MCRLB, whereas for M → ∞ it continues to increase with increasing SNRs. This property suggests that NDA estimation of the SNR becomes more and more problematic in case we increase the order of the selected PAM scheme.   Figure 4 includes also the evolution of the CRLB for M → ∞ as it has been derived in Section 4. One can see that the theoretical limit for M → ∞ is close to that computed for M = 16 as long as the latter does not start to approach the MCRLB, whereas for M → ∞ it continues to increase with increasing SNRs. This property suggests that NDA estimation of the SNR becomes more and more problematic in case we increase the order of the selected PAM scheme.   The diagrams above visualize the error and bias performance for different PAM constellations operated with α = 0.0, which is perhaps most interesting in practice, since it represents the scenario with minimum bandwidth. Nevertheless, in order to complete the portrait, Figure 6 shows the normalized error performance for a 4-PAM signal operated with L = 1000 and α ∈ {0.0, 0.3, 1.0}. According to the exact relationship derived in [17], the MCRLB is proportional to α + 2 (1 ) L for very low SNRs, whereas for ρs → ∞ it approaches 2 L regardless of the selected α, or in other words: with respect to α = 0, the MCRLB for α > 0 appears as shifted to the right depending on the chosen value.
Exactly the same behavior is reflected by the CRLB, although the results are in the very low SNR range somewhat higher than those achieved with the MCRLB, but for rather large values, say ρs > 25 dB, the CRLB approaches more and more the horizontal floor characterizing the MCRLB performance. Of course, in the medium SNR range, the CRLB deviates significantly from the MCRLB, which applies also to the simplified computation for α = 0.3.  Figure 6 includes also the normalized jitter variance of the EM algorithm developed in the previous section. By detailed inspection, we observe that the performance differs The diagrams above visualize the error and bias performance for different PAM constellations operated with α = 0.0, which is perhaps most interesting in practice, since it represents the scenario with minimum bandwidth. Nevertheless, in order to complete the portrait, Figure 6 shows the normalized error performance for a 4-PAM signal operated with L = 1000 and α ∈ {0.0, 0.3, 1.0}. According to the exact relationship derived in [17], the MCRLB is proportional to 2(1 + α)/L for very low SNRs, whereas for ρ s → ∞ it approaches 2/L regardless of the selected α, or in other words: with respect to α = 0, the MCRLB for α > 0 appears as shifted to the right depending on the chosen value.
Exactly the same behavior is reflected by the CRLB, although the results are in the very low SNR range somewhat higher than those achieved with the MCRLB, but for rather large values, say ρ s > 25 dB, the CRLB approaches more and more the horizontal floor characterizing the MCRLB performance. Of course, in the medium SNR range, the CRLB deviates significantly from the MCRLB, which applies also to the simplified computation for α = 0.3. The diagrams above visualize the error and bias performance for different PAM constellations operated with α = 0.0, which is perhaps most interesting in practice, since it represents the scenario with minimum bandwidth. Nevertheless, in order to complete the portrait, Figure 6 shows the normalized error performance for a 4-PAM signal operated with L = 1000 and α ∈ {0.0, 0.3, 1.0}. According to the exact relationship derived in [17], the MCRLB is proportional to α + 2 (1 ) L for very low SNRs, whereas for ρs → ∞ it approaches 2 L regardless of the selected α, or in other words: with respect to α = 0, the MCRLB for α > 0 appears as shifted to the right depending on the chosen value.
Exactly the same behavior is reflected by the CRLB, although the results are in the very low SNR range somewhat higher than those achieved with the MCRLB, but for rather large values, say ρs > 25 dB, the CRLB approaches more and more the horizontal floor characterizing the MCRLB performance. Of course, in the medium SNR range, the CRLB deviates significantly from the MCRLB, which applies also to the simplified computation for α = 0.3.  Figure 6 includes also the normalized jitter variance of the EM algorithm developed in the previous section. By detailed inspection, we observe that the performance differs for α = 0 in the lower SNR domain somewhat from the CRLB, which is mainly due to a  Figure 6 includes also the normalized jitter variance of the EM algorithm developed in the previous section. By detailed inspection, we observe that the performance differs for α = 0 in the lower SNR domain somewhat from the CRLB, which is mainly due to a residual bias effect, whereas for medium-to-high SNR values the jitter variances are very close to the corresponding CRLBs regardless of the selected excess bandwidth.

Concluding Remarks
The availability of reliable SNR estimates is most helpful in many communication systems, particularly when adaptive solutions have to be considered in terms of modulation and coding schemes. This is not only true for radio frequency, but also for optical wireless links. Assuming a bandlimited optical intensity channel, an algorithm for SNR estimation has been developed, which does not require any knowledge about the transmitted data symbols. Such an NDA approach is very appreciated, because the larger observation lengths do not adversely affect the spectral efficiency as it would happen with a DA solution. Maximum likelihood, moment-based and decision-directed methods were out of scope because of complexity and/or performance reasons, but it turned out that the developed expectation-maximization algorithm exhibits an error performance close to the CRLB as the theoretical limit, which is mainly true for longer observation windows where bias effects are negligible.