**3. Turbulent atmospheric channel model**

There is an extensive literature on the subject of the theory of line-of-sight propagation through the atmosphere (Andrews & Phillips, 1998; Andrews et al., 2000; Fante, 1975; Ishimaru, 1997; Strohbehn, 1978; Tatarskii, 1971). One of the most important works was developed by Tatarskii (Tatarskii, 1971). He supposed a plane wave that is incident upon the random medium (the atmosphere in this particular case). It is assumed an atmosphere having no free charges with a constant magnetic permeability. In addition, it is suppossed that the electromagnetic field has a sinusoidal time dependence (a monochromatic wave). Under these

Operating, assuming that |∇*ψ*1| � |∇*ψ*0|, due to *<sup>n</sup>*1(**r**) � 1, neglecting *<sup>n</sup>*<sup>2</sup>

<sup>∇</sup>2*ψ*<sup>0</sup> + (∇*ψ*0)

expressions are obtained:

exp (*jkx*), this solution can be written as:

order corrections (*χ*<sup>1</sup> and *S*1).

wave at distance *L* from the source is represented by:

*I* = |*U*0(*L*,**r**)|

the irradiance of the random field shown in Eq. (12) takes the form:

<sup>2</sup> exp (*ψ*<sup>1</sup> + *ψ*<sup>∗</sup>

where, from now onwards, we denote *χ*<sup>1</sup> as *χ* for simplicity in the notation. Hence,

that the fading does not attenuate or amplify the average power, i.e., *E*[*I*] = |*U*0|

*<sup>ψ</sup>*1(*L*,**r**) = *<sup>k</sup>*<sup>2</sup>

2*π*

Optical Communications Using an Alternative Matrix Decomposition

 *V n*1(**r** � ) exp *jk* |**r** − **r**�

to 2*n*1(**r**), and equating the terms with the same order of perturbation, then the following

<sup>433</sup> Upper Burst Error Bound for Atmospheric Correlated

The first one is the differential equation for ∇*ψ* in the absence of the fluctuation whereas turbulent atmosphere induced perturbation are found in the second expression. The resolution of Eq. (10b) is detailed in (Fante, 1975; Ishimaru, 1997). For the particular case of a monochromatic optical plane wave propagating along the positive x-axis, i.e., *U*0(*L*,**r**) =

where the position (*L*,**r**) denotes a position in the receiver plane (at *x* = *L*) whereas (*x*�

represents any position at an arbitrary plane along the propagation path. The mathematical development needed to solve Eq. (11) can be consulted in (Andrews & Phillips, 1998; Ishimaru, 1997). Furthermore, the statistical nature of *ψ*1(*L*,**r**) can be deduced in an easy way. Equation (11) has the physical interpretation that the first-order Rytov perturbation, *ψ*1(*L*,**r**) is a sum of spherical waves generated at various points **r**� throughout the scattering volume *V*, the strength of each sum wave being proportional to the product of the unperturbed field term *U*<sup>0</sup> and the refractive-index perturbation, *n*1, at the point **r**� (Andrews & Phillips, 1998). Thus it is possible to apply the central limit theorem. According to such a theorem, the distribution of a random variable which is a sum of *N* independent random variables approaches normal as *N* → ∞ regardless of the distribution of each random variable. Application of the central limit theorem to this integral equation leads to the prediction of a normal probability distribution for *ψ*. Since we can substitute Ψ = *χ* + *jS*, where *χ* and *S* are called the log-amplitude and phase, respectively, of the field, then application of the central limit theorem also leads to the prediction of a Gaussian (normal) probability distribution for both *χ* and *S*, at least up to first

Accordingly, under this first-order Rytov approximation, the field of a propagating optical

with *U*0(*L*,**r**) being the unperturbed portion of the field in the absence of turbulence. Hence,

In Eq. (13), operator ∗ denotes the complex conjugate, |*U*0| is the amplitude of the unperturbed field and *I*<sup>0</sup> is the level of irradiance fluctuation in the absence of air turbulence that ensures

thought of as a conservation of energy consideration and requires the choice of *<sup>E</sup>*[*χ*] = <sup>−</sup>*σ*<sup>2</sup>

as was explained in (Fried, 1967; Strohbehn, 1978), where *E*[*χ*] is the ensemble average of

<sup>2</sup> + *k*2*n*<sup>2</sup>

<sup>∇</sup>2*ψ*<sup>1</sup> <sup>+</sup> <sup>2</sup>∇*ψ*0∇*ψ*<sup>1</sup> <sup>+</sup> <sup>2</sup>*k*2*n*1(**r**) = 0. (10b)


*U* = exp (*ψ*) = *U*0(*L*,**r**) exp (*ψ*1), (12)

*I* = *I*<sup>0</sup> exp (2*χ*), [*w*/*m*2]. (14)

<sup>1</sup> ) = *<sup>I</sup>*<sup>0</sup> exp (2*χ*1), [*w*/*m*2] (13)

<sup>|</sup>**<sup>r</sup>** <sup>−</sup> **<sup>r</sup>**�| *<sup>d</sup>*3**r**�


, (11)

2. This may be

*χ*,

,**r**� )

<sup>0</sup>(**r**) = 0; (10a)

circumstances, the vector wave equation becomes

$$
\nabla^2 \mathbf{E} + k^2 n^2(\mathbf{r}) \mathbf{E} + 2\nabla \left( \mathbf{E} \cdot \nabla \log n(\mathbf{r}) \right) = 0,\tag{3}
$$

where **E** is the vector amplitude of the electric field, *k*=2*π*/*λ* is the wave number of the electromagnetic wave with *λ* being the optical wavelength; whereas *n* is the atmospheric refractive index whose time variations have been suppressed, and being a random function of position, **r** = (*x*, *y*, *z*). The ∇ operator is the well-known vector derivative (*∂*/*∂x*, *∂*/*∂y*, *∂*/*∂z*). Equation (3) can be simplified by imposing certain characteristics of the propagation wave. In particular, since the wavelength *λ* for optical radiation is much smaller than the smallest scale of turbulence, *<sup>l</sup>*0, (Strohbehn, 1968) the maximum scattering angle is roughly *<sup>λ</sup>*/*l*<sup>0</sup> <sup>≈</sup> <sup>10</sup>−<sup>4</sup> rad. As a consequence, the last term on the left-hand side of Eq. (3) is negligible. Such a term is related to the change in polarization of the wave as it propagates (Strohbehn, 1971; Strohbehn & Clifford, 1967). This conclusion permit us to drop the last term and Eq. (3) then reduces to

$$
\nabla^2 \mathbf{E} + k^2 n^2(\mathbf{r}) \mathbf{E} = 0. \tag{4}
$$

Because Eq. (4) is easily decomposed into three scalar equations, one for each component of the electric field, **E**, we may solve one scalar equation and ignore the vector character of the wave until the final solution. Therefore if we let *U*(**r**) denote one of the scalar components that is transverse to the direction of propagation along the positive x-axis (Andrews & Phillips, 1998), then Eq. (4) may be replaced by the scalar stochastic differential equation

$$
\nabla^2 \mathcal{U} + k^2 n^2(\mathbf{r}) \mathcal{U} = 0. \tag{5}
$$

The index of refraction, *n*(**r**)=*n*<sup>0</sup> + *n*1(**r**), fluctuates about the average value *n*<sup>0</sup> = *E*[*n*(**r**)] ∼= 1, whereas *n*1(**r**) � 1 is the fluctuation of the refractive index from its free space value. Thus

$$
\nabla^2 \mathcal{U} + k^2 (n\_0 + n\_1(\mathbf{r}))^2 \mathcal{U} = 0. \tag{6}
$$

For weak fluctuation, it is necessary to obtain an approximate solution of Eq. (6) for small *n*1. Most of the literature since 1960 has followed the approach of using the so-called Rytov method, which substitutes *U* in a series:

$$\mathcal{U} = \exp\left(\psi\_0 + \psi\_1 + \psi\_2 + \dots\right) = \exp\left(\psi\right). \tag{7}$$

In Eq. (7), *ψ*1, *ψ*<sup>2</sup> are the first and second order complex phase perturbations, respectively, whereas *ψ*<sup>0</sup> is the phase of the optical wave in free space. The Rytov solution is widely used in line-of-sight propagation problems because it simplifies the procedure of obtaining both amplitude and phase fluctuations. From the Rytov solution, the wave equation becomes:

$$(\nabla^2 \psi + (\nabla \psi)^2 + k^2 (n\_0 + n\_1(\mathbf{r}))^2 = 0. \tag{8}$$

This is a nonlinear first order differential equation for ∇*ψ* and is known as the Riccati equation. Consider now a first order perturbation, then

$$
\psi(L, \mathbf{r}) = \psi\_0(L, \mathbf{r}) + \psi\_1(L, \mathbf{r});\tag{9a}
$$

$$n(\mathbf{r}) = n\_0 + n\_1(\mathbf{r}); \qquad \qquad n\_0 \cong 1. \tag{9b}$$

4 Numerical Simulations / Book 1

where **E** is the vector amplitude of the electric field, *k*=2*π*/*λ* is the wave number of the electromagnetic wave with *λ* being the optical wavelength; whereas *n* is the atmospheric refractive index whose time variations have been suppressed, and being a random function of position, **r** = (*x*, *y*, *z*). The ∇ operator is the well-known vector derivative (*∂*/*∂x*, *∂*/*∂y*, *∂*/*∂z*). Equation (3) can be simplified by imposing certain characteristics of the propagation wave. In particular, since the wavelength *λ* for optical radiation is much smaller than the smallest scale of turbulence, *<sup>l</sup>*0, (Strohbehn, 1968) the maximum scattering angle is roughly *<sup>λ</sup>*/*l*<sup>0</sup> <sup>≈</sup> <sup>10</sup>−<sup>4</sup> rad. As a consequence, the last term on the left-hand side of Eq. (3) is negligible. Such a term is related to the change in polarization of the wave as it propagates (Strohbehn, 1971; Strohbehn & Clifford, 1967). This conclusion permit us to drop the last term and Eq. (3) then

Because Eq. (4) is easily decomposed into three scalar equations, one for each component of the electric field, **E**, we may solve one scalar equation and ignore the vector character of the wave until the final solution. Therefore if we let *U*(**r**) denote one of the scalar components that is transverse to the direction of propagation along the positive x-axis (Andrews & Phillips,

The index of refraction, *n*(**r**)=*n*<sup>0</sup> + *n*1(**r**), fluctuates about the average value *n*<sup>0</sup> = *E*[*n*(**r**)] ∼= 1, whereas *n*1(**r**) � 1 is the fluctuation of the refractive index from its free space value. Thus

For weak fluctuation, it is necessary to obtain an approximate solution of Eq. (6) for small *n*1. Most of the literature since 1960 has followed the approach of using the so-called Rytov

In Eq. (7), *ψ*1, *ψ*<sup>2</sup> are the first and second order complex phase perturbations, respectively, whereas *ψ*<sup>0</sup> is the phase of the optical wave in free space. The Rytov solution is widely used in line-of-sight propagation problems because it simplifies the procedure of obtaining both amplitude and phase fluctuations. From the Rytov solution, the wave equation becomes:

This is a nonlinear first order differential equation for ∇*ψ* and is known as the Riccati equation.

1998), then Eq. (4) may be replaced by the scalar stochastic differential equation

**E** · ∇ log *n*(**r**)

<sup>∇</sup>2**<sup>E</sup>** <sup>+</sup> *<sup>k</sup>*2*n*2(**r**)**<sup>E</sup>** <sup>=</sup> 0. (4)

<sup>∇</sup>2*<sup>U</sup>* <sup>+</sup> *<sup>k</sup>*2*n*2(**r**)*<sup>U</sup>* <sup>=</sup> 0. (5)

<sup>∇</sup>2*<sup>U</sup>* <sup>+</sup> *<sup>k</sup>*2(*n*<sup>0</sup> <sup>+</sup> *<sup>n</sup>*1(**r**))2*<sup>U</sup>* <sup>=</sup> 0. (6)

*U* = exp (*ψ*<sup>0</sup> + *ψ*<sup>1</sup> + *ψ*<sup>2</sup> + ...) = exp (*ψ*). (7)

<sup>∇</sup>2*<sup>ψ</sup>* + (∇*ψ*)<sup>2</sup> <sup>+</sup> *<sup>k</sup>*2(*n*<sup>0</sup> <sup>+</sup> *<sup>n</sup>*1(**r**))<sup>2</sup> <sup>=</sup> 0. (8)

*ψ*(*L*,**r**) = *ψ*0(*L*,**r**) + *ψ*1(*L*,**r**); (9a)

*n*(**r**) = *n*<sup>0</sup> + *n*1(**r**); *n*<sup>0</sup> ∼= 1. (9b)

= 0, (3)

<sup>∇</sup>2**<sup>E</sup>** <sup>+</sup> *<sup>k</sup>*2*n*2(**r**)**<sup>E</sup>** <sup>+</sup> <sup>2</sup><sup>∇</sup>

circumstances, the vector wave equation becomes

method, which substitutes *U* in a series:

Consider now a first order perturbation, then

reduces to

Operating, assuming that |∇*ψ*1| � |∇*ψ*0|, due to *<sup>n</sup>*1(**r**) � 1, neglecting *<sup>n</sup>*<sup>2</sup> <sup>1</sup>(**r**) in comparison to 2*n*1(**r**), and equating the terms with the same order of perturbation, then the following expressions are obtained:

$$
\nabla^2 \psi\_0 + (\nabla \psi\_0)^2 + k^2 n\_0^2(\mathbf{r}) = 0;\tag{10a}
$$

$$
\nabla^2 \psi\_1 + 2\nabla \psi\_0 \nabla \psi\_1 + 2k^2 n\_1(\mathbf{r}) = 0. \tag{10b}
$$

The first one is the differential equation for ∇*ψ* in the absence of the fluctuation whereas turbulent atmosphere induced perturbation are found in the second expression. The resolution of Eq. (10b) is detailed in (Fante, 1975; Ishimaru, 1997). For the particular case of a monochromatic optical plane wave propagating along the positive x-axis, i.e., *U*0(*L*,**r**) = exp (*jkx*), this solution can be written as:

$$\psi\_1(L, \mathbf{r}) = \frac{k^2}{2\pi} \iiint\_V n\_1(\mathbf{r}') \frac{\exp\left(jk\left[|\mathbf{r} - \mathbf{r}'| - |L - \mathbf{x}'|\right]\right)}{|\mathbf{r} - \mathbf{r}'|} d^3 \mathbf{r}' \tag{11}$$

where the position (*L*,**r**) denotes a position in the receiver plane (at *x* = *L*) whereas (*x*� ,**r**� ) represents any position at an arbitrary plane along the propagation path. The mathematical development needed to solve Eq. (11) can be consulted in (Andrews & Phillips, 1998; Ishimaru, 1997). Furthermore, the statistical nature of *ψ*1(*L*,**r**) can be deduced in an easy way. Equation (11) has the physical interpretation that the first-order Rytov perturbation, *ψ*1(*L*,**r**) is a sum of spherical waves generated at various points **r**� throughout the scattering volume *V*, the strength of each sum wave being proportional to the product of the unperturbed field term *U*<sup>0</sup> and the refractive-index perturbation, *n*1, at the point **r**� (Andrews & Phillips, 1998). Thus it is possible to apply the central limit theorem. According to such a theorem, the distribution of a random variable which is a sum of *N* independent random variables approaches normal as *N* → ∞ regardless of the distribution of each random variable. Application of the central limit theorem to this integral equation leads to the prediction of a normal probability distribution for *ψ*. Since we can substitute Ψ = *χ* + *jS*, where *χ* and *S* are called the log-amplitude and phase, respectively, of the field, then application of the central limit theorem also leads to the prediction of a Gaussian (normal) probability distribution for both *χ* and *S*, at least up to first order corrections (*χ*<sup>1</sup> and *S*1).

Accordingly, under this first-order Rytov approximation, the field of a propagating optical wave at distance *L* from the source is represented by:

$$\mathcal{U} = \exp\left(\psi\right) = \mathcal{U}\_0(L, \mathbf{r}) \exp\left(\psi\_1\right), \tag{12}$$

with *U*0(*L*,**r**) being the unperturbed portion of the field in the absence of turbulence. Hence, the irradiance of the random field shown in Eq. (12) takes the form:

$$I = |\mathcal{U}l\_0(L, \mathbf{r})|^2 \exp\left(\psi\_1 + \psi\_1^\*\right) = I\_0 \exp\left(2\chi\_1\right), \qquad \left[w/m^2\right] \tag{13}$$

where, from now onwards, we denote *χ*<sup>1</sup> as *χ* for simplicity in the notation. Hence,

$$I = I\_0 \exp\left(2\chi\right), \qquad \left[w/m^2\right]. \tag{14}$$

In Eq. (13), operator ∗ denotes the complex conjugate, |*U*0| is the amplitude of the unperturbed field and *I*<sup>0</sup> is the level of irradiance fluctuation in the absence of air turbulence that ensures that the fading does not attenuate or amplify the average power, i.e., *E*[*I*] = |*U*0| 2. This may be thought of as a conservation of energy consideration and requires the choice of *<sup>E</sup>*[*χ*] = <sup>−</sup>*σ*<sup>2</sup> *χ*, as was explained in (Fried, 1967; Strohbehn, 1978), where *E*[*χ*] is the ensemble average of

white Gaussian noise, *N*(*t*), is assumed to include any front-end receiver thermal noise as well as shot noise caused by ambient light much stronger than the desired signal. In the following section, we complete the scheme presented in (Jurado-Navas et al., 2007) to afford the inclusion of *m* correlated scintillation sequences approximating to the model introduced

<sup>435</sup> Upper Burst Error Bound for Atmospheric Correlated

As indicated in the last section, we can simplify the model proposed in (Jurado-Navas & Puerta-Notario, 2009) by an unrealistic but reasonably accurate space-time separable statistics model with a reduced computational load if it is compared with the complete model. Through this section, we develop the space-time separable statistics model.

When multiple receivers are considered, then, as shown in (Zhu & Kahn, 2002), the real

*χρ<sup>d</sup>*<sup>12</sup> ... *<sup>σ</sup>*<sup>2</sup>

... ... ... ...

where *dij* and *ρdij* are the distance and its normalized CC coefficient respectively between points *i* and *j* in the receiver plane. In (Jurado-Navas & Puerta-Notario, 2009; Zhu & Kahn, 2002), a Gaussian spatial covariance function for the log-amplitude fluctuations is employed

In order to satisfy the hypothesis of frozen turbulence, we can adopt an AR model as in (Baddour & Beaulieu, 2002; Jurado-Navas & Puerta-Notario, 2009) as a possible solution. By doing so, the multichannel Yule-Walker equations are expressed as

where **<sup>A</sup>***H*[*k*], *<sup>k</sup>* <sup>=</sup> 1, 2, ...*<sup>p</sup>* are *<sup>m</sup>* <sup>×</sup> *<sup>m</sup>* matrices containing the multichannel AR model coefficients; whereas **C**χ[*j*] is the covariance matrix evaluated in the '*j*'-time instant. Then, the system of equations in Eq. (22) can be solved efficiently via the Levinson-Wiggins-Robinson algorithm (Kay, 1988). Once the **A***H*[*k*] coefficient matrices have been determined, we can obtain the *m* × *m* covariance matrix of the driving noise vector process of the AR model from:

> *p* ∑ *k*=1

by a transpose operator due to all the samples are real; the driving noise process, ω[*n*],

**C**<sup>w</sup> = **C**χ[0] +

*<sup>χ</sup>* ... *<sup>σ</sup>*<sup>2</sup>

*χρdn*<sup>2</sup> ... *<sup>σ</sup>*<sup>2</sup>

*χρ<sup>d</sup>*1*<sup>m</sup>*

⎞

⎟⎟⎠

*m*×*m*

**<sup>A</sup>H**[1]*m*×*<sup>m</sup>* **<sup>A</sup>H**[2]*m*×*<sup>m</sup>* ... **<sup>A</sup>H**[*p*]*m*×*<sup>m</sup>*

⎞

⎟⎟⎠ <sup>=</sup> <sup>−</sup>

**<sup>C</sup>**χ[−*k*]**A***H*[*k*]. (23)

*<sup>T</sup>*}, where the hermitian operator has been substituted

⎛

**C**χ[1]*m*×*<sup>m</sup>* **C**χ[2]*m*×*<sup>m</sup>* ... **C**χ[*p*]*m*×*<sup>m</sup>*

⎞

⎟⎟⎠ .

(22)

⎜⎜⎝

*χρ<sup>d</sup>*2*<sup>m</sup>*

*χ*

⎞

⎛

⎜⎜⎝

⎟⎟⎠

*<sup>i</sup>*,*j*=1, at *m* receivers in

(21)

symmetric auto-covariance matrix of the log-amplitude, **<sup>C</sup>**χ={c*χ*(i,j)}*<sup>m</sup>*

⎛

*σ*2 *<sup>χ</sup> σ*<sup>2</sup>

*χρ<sup>d</sup>*<sup>21</sup> *<sup>σ</sup>*<sup>2</sup>

that approximates the theoretical covariance function resulting from Rytov theory.

*σ*2

*σ*2 *χρdm*<sup>1</sup> *<sup>σ</sup>*<sup>2</sup>

⎜⎜⎝

a plane transverse to the direction of propagation is given by:

Optical Communications Using an Alternative Matrix Decomposition

**C**<sup>χ</sup> =

**C**χ[0]*m*×*<sup>m</sup>* **C**χ[−1]*m*×*<sup>m</sup>* ... **C**χ[−*p* + 1]*m*×*<sup>m</sup>* **C**χ[1]*m*×*<sup>m</sup>* **C**χ[0]*m*×*<sup>m</sup>* ... **C**χ[−*p* + 2]*m*×*<sup>m</sup>* ... ... ... ... **C**χ[*p* − 1]*m*×*<sup>m</sup>* **C**χ[*p* − 2]*m*×*<sup>m</sup>* ... **C**χ[0]*m*×*<sup>m</sup>*

**4. Proposed approximation: space-time separable statistics channel model**

in (Jurado-Navas & Puerta-Notario, 2009).

**4.1 Spatial diversity reception**

(Jurado-Navas & Puerta-Notario, 2009):

Thus, after obtaining **C**<sup>w</sup> = *E*{ω[*n*]ω[*n*]

⎛

⎜⎜⎝

log-amplitude, whereas *σ*<sup>2</sup> *<sup>χ</sup>* is its variance depending on the structure parameter, *<sup>C</sup>*<sup>2</sup> *<sup>n</sup>*. With all of these expressions, we have modeled the irradiance of the random field, *I*, in the space at a single instant in time. Now, because the state of the atmospheric turbulence varies with time, the intensity fluctuations is also temporally correlated. Then, Eq. (14) can be expressed as:

$$I = \mathfrak{a}\_{\rm sc}(t) \cdot I\_{0\prime} \tag{15}$$

whereas *αsc*(*t*) = exp (2*χ*(*t*)) is the temporal behavior of the scintillation sequence and represents the effect of the intensity fluctuations on the transmitted signal. In Eq. (15), a space-to-time statistical conversion has been assumed by employing the well-known Taylor's hypothesis of frozen turbulence (Tatarskii, 1971; Taylor, 1938).

As analyzed before, and by the central limit theorem, the marginal distribution of the logamplitude, *χ*, is Gaussian. Thus,

$$f\_{\chi}(\chi) = \left(\frac{1}{2\pi\sigma\_{\chi}^2}\right)^{1/2} \exp\left[-\frac{(\chi - E[\chi])^2}{2\sigma\_{\chi}^2}\right].\tag{16}$$

Hence, from the Jacobian statistical transformation (Papoulis, 1991),

$$f\_I(I) = \frac{f\_\chi(\chi)}{|\frac{dI}{d\chi}|},\tag{17}$$

the probability density function of the intensity, *I*, can be identified to have a lognormal distribution typical of weak turbulence regime. Then:

$$f\_I(I) = \left(\frac{1}{2I}\right) \left(\frac{1}{2\pi\sigma\_\chi^2}\right)^{1/2} \exp\left[-\frac{(\ln I - \ln I\_0)^2}{8\sigma\_\chi^2}\right].\tag{18}$$

Theoretical and experimental studies of irradiance fluctuations generally center around the scintillation index. It was evaluated in (Mercier, 1962) and it is defined as the normalized variance of irradiance fluctuations:

$$
\sigma\_I^2 = \frac{E[I^2]}{\left(E[I]\right)^2} - 1.\tag{19}
$$

Hence it is possible to define the weak turbulence regimes as those regimes for which the scintillation index given in Eq. (19) is less than unity.

With all these considerations taken into account, an efficient channel model for FSO communications using intensity modulation and direct detection (IM/DD) was presented in (Jurado-Navas et al., 2007; 2011a) under the assumption of weak turbulence regime. For these systems, the received optical power, *Y*(*t*), can be written as

$$Y(t) = a\_{\rm sc}(t)X(t) + N(t),\tag{20}$$

being *X*(*t*) the received optical power without scintillation; whereas *αsc*(*t*) = exp 2*χ*(*t*) is the temporal behavior of the scintillation sequence and represent the effect of the intensity fluctuations on the transmitted signal. To generate *αsc*(*t*), a scheme based on a lowpass filtering of a random Gaussian signal, *z*(*t*), is implemented as in (Jurado-Navas et al., 2007; 2011a). *χ*(*t*) is, as was explained above, the log-amplitude of the optical wave governed by Gaussian statistics with ensemble average *E*[*χ*] and variance *σ*<sup>2</sup> *<sup>χ</sup>*. Finally, the additive white Gaussian noise, *N*(*t*), is assumed to include any front-end receiver thermal noise as well as shot noise caused by ambient light much stronger than the desired signal. In the following section, we complete the scheme presented in (Jurado-Navas et al., 2007) to afford the inclusion of *m* correlated scintillation sequences approximating to the model introduced in (Jurado-Navas & Puerta-Notario, 2009).

### **4. Proposed approximation: space-time separable statistics channel model**

As indicated in the last section, we can simplify the model proposed in (Jurado-Navas & Puerta-Notario, 2009) by an unrealistic but reasonably accurate space-time separable statistics model with a reduced computational load if it is compared with the complete model. Through this section, we develop the space-time separable statistics model.

#### **4.1 Spatial diversity reception**

6 Numerical Simulations / Book 1

of these expressions, we have modeled the irradiance of the random field, *I*, in the space at a single instant in time. Now, because the state of the atmospheric turbulence varies with time, the intensity fluctuations is also temporally correlated. Then, Eq. (14) can be expressed as:

whereas *αsc*(*t*) = exp (2*χ*(*t*)) is the temporal behavior of the scintillation sequence and represents the effect of the intensity fluctuations on the transmitted signal. In Eq. (15), a space-to-time statistical conversion has been assumed by employing the well-known Taylor's

As analyzed before, and by the central limit theorem, the marginal distribution of the log-

exp 

*fI*(*I*) = *<sup>f</sup>χ*(*χ*) | *d I dχ* |

the probability density function of the intensity, *I*, can be identified to have a lognormal

1/2

Theoretical and experimental studies of irradiance fluctuations generally center around the scintillation index. It was evaluated in (Mercier, 1962) and it is defined as the normalized

Hence it is possible to define the weak turbulence regimes as those regimes for which the

With all these considerations taken into account, an efficient channel model for FSO communications using intensity modulation and direct detection (IM/DD) was presented in (Jurado-Navas et al., 2007; 2011a) under the assumption of weak turbulence regime. For these

the temporal behavior of the scintillation sequence and represent the effect of the intensity fluctuations on the transmitted signal. To generate *αsc*(*t*), a scheme based on a lowpass filtering of a random Gaussian signal, *z*(*t*), is implemented as in (Jurado-Navas et al., 2007; 2011a). *χ*(*t*) is, as was explained above, the log-amplitude of the optical wave governed

being *X*(*t*) the received optical power without scintillation; whereas *αsc*(*t*) = exp

*<sup>I</sup>* <sup>=</sup> *<sup>E</sup>*[*I*2] *E*[*I*]

exp 

<sup>−</sup> (*<sup>χ</sup>* <sup>−</sup> *<sup>E</sup>*[*χ*])<sup>2</sup> 2*σ*<sup>2</sup> *χ*

> <sup>−</sup> (ln *<sup>I</sup>* <sup>−</sup> ln *<sup>I</sup>*0)<sup>2</sup> 8*σ*<sup>2</sup> *χ*

*Y*(*t*) = *αsc*(*t*)*X*(*t*) + *N*(*t*), (20)

1/2

hypothesis of frozen turbulence (Tatarskii, 1971; Taylor, 1938).

*fχ*(*χ*) =

distribution typical of weak turbulence regime. Then:

 1 2*I*  1 2*πσ*<sup>2</sup> *χ*

*σ*2

*fI*(*I*) =

scintillation index given in Eq. (19) is less than unity.

systems, the received optical power, *Y*(*t*), can be written as

by Gaussian statistics with ensemble average *E*[*χ*] and variance *σ*<sup>2</sup>

variance of irradiance fluctuations:

 1 2*πσ*<sup>2</sup> *χ*

Hence, from the Jacobian statistical transformation (Papoulis, 1991),

*<sup>χ</sup>* is its variance depending on the structure parameter, *<sup>C</sup>*<sup>2</sup>

*I* = *αsc*(*t*) · *I*0, (15)

. (16)

. (18)

2*χ*(*t*) is

*<sup>χ</sup>*. Finally, the additive

, (17)

<sup>2</sup> <sup>−</sup> 1. (19)

*<sup>n</sup>*. With all

log-amplitude, whereas *σ*<sup>2</sup>

amplitude, *χ*, is Gaussian. Thus,

When multiple receivers are considered, then, as shown in (Zhu & Kahn, 2002), the real symmetric auto-covariance matrix of the log-amplitude, **<sup>C</sup>**χ={c*χ*(i,j)}*<sup>m</sup> <sup>i</sup>*,*j*=1, at *m* receivers in a plane transverse to the direction of propagation is given by:

$$\mathbf{C}\_{\mathcal{X}} = \begin{pmatrix} \sigma\_{\mathcal{X}}^2 & \sigma\_{\mathcal{X}}^2 \rho\_{d\_{12}} & \dots \sigma\_{\mathcal{X}}^2 \rho\_{d\_{1m}} \\ \sigma\_{\mathcal{X}}^2 \rho\_{d\_{21}} & \sigma\_{\mathcal{X}}^2 & \dots \sigma\_{\mathcal{X}}^2 \rho\_{d\_{2m}} \\ \dots & \dots & \dots & \dots \\ \sigma\_{\mathcal{X}}^2 \rho\_{d\_{m1}} & \sigma\_{\mathcal{X}}^2 \rho\_{d\_{m2}} & \dots & \sigma\_{\mathcal{X}}^2 \end{pmatrix}\_{m \times m} \tag{21}$$

where *dij* and *ρdij* are the distance and its normalized CC coefficient respectively between points *i* and *j* in the receiver plane. In (Jurado-Navas & Puerta-Notario, 2009; Zhu & Kahn, 2002), a Gaussian spatial covariance function for the log-amplitude fluctuations is employed that approximates the theoretical covariance function resulting from Rytov theory.

In order to satisfy the hypothesis of frozen turbulence, we can adopt an AR model as in (Baddour & Beaulieu, 2002; Jurado-Navas & Puerta-Notario, 2009) as a possible solution. By doing so, the multichannel Yule-Walker equations are expressed as (Jurado-Navas & Puerta-Notario, 2009):

$$
\begin{pmatrix}
\mathbf{C\_{X}}[0]\_{m\times m} & \mathbf{C\_{X}}[-1]\_{m\times m} & \dots & \mathbf{C\_{X}}[-p+1]\_{m\times m} \\
\mathbf{C\_{X}}[1]\_{m\times m} & \mathbf{C\_{X}}[0]\_{m\times m} & \dots & \mathbf{C\_{X}}[-p+2]\_{m\times m} \\
\dots & \dots & \dots & \dots \\
\mathbf{C\_{X}}[p-1]\_{m\times m} & \mathbf{C\_{X}}[p-2]\_{m\times m} & \dots & \mathbf{C\_{X}}[0]\_{m\times m}
\end{pmatrix}
\begin{pmatrix}
\mathbf{A^{H}}[1]\_{m\times m} \\
\mathbf{A^{H}}[2]\_{m\times m} \\
\dots \\
\mathbf{A^{H}}[p]\_{m\times m}
\end{pmatrix} = -\begin{pmatrix}
\mathbf{C\_{X}}[1]\_{m\times m} \\
\mathbf{C\_{X}}[2]\_{m\times m} \\
\dots \\
\mathbf{C\_{X}}[p]\_{m\times m}
\end{pmatrix}.
\tag{22}
$$

where **<sup>A</sup>***H*[*k*], *<sup>k</sup>* <sup>=</sup> 1, 2, ...*<sup>p</sup>* are *<sup>m</sup>* <sup>×</sup> *<sup>m</sup>* matrices containing the multichannel AR model coefficients; whereas **C**χ[*j*] is the covariance matrix evaluated in the '*j*'-time instant. Then, the system of equations in Eq. (22) can be solved efficiently via the Levinson-Wiggins-Robinson algorithm (Kay, 1988). Once the **A***H*[*k*] coefficient matrices have been determined, we can obtain the *m* × *m* covariance matrix of the driving noise vector process of the AR model from:

$$\mathbf{C}\_{w} = \mathbf{C}\_{\mathbf{X}}[0] + \sum\_{k=1}^{p} \mathbf{C}\_{\mathbf{X}}[-k] \mathbf{A}^{H}[k]. \tag{23}$$

Thus, after obtaining **C**<sup>w</sup> = *E*{ω[*n*]ω[*n*] *<sup>T</sup>*}, where the hermitian operator has been substituted by a transpose operator due to all the samples are real; the driving noise process, ω[*n*],

(Holton, 2004). In (Abramovich et al., 2001), the particular problem of non-positive definite matrices that arise from characterizing a realistic system is investigated in detail, showing that in 988 outcomes of their 1000 Monte Carlo trials, gave rise to a non-positive definite covariance matrix, although the corresponding "contracted" matrix may well be positive definite. In addition, in most practical situations the covariance matrix must be estimated from scintillation measurements (Monserrat et al., 2007; Moore et al., 2005) and outliers and other special situations as urban canopies may induce that correlation matrices derived from measurements are not always guaranteed to be positive definite, although a theoretical

<sup>437</sup> Upper Burst Error Bound for Atmospheric Correlated

Optical Communications Using an Alternative Matrix Decomposition

In such cases, a real Schur decomposition is proposed (Golub & van Loan, 1996). Since the symmetric matrix **<sup>C</sup>**<sup>χ</sup> ∈ �*mxm*, then there exists an orthogonal matrix, **<sup>Q</sup>** ∈ �*mxm*, such that

as was shown in (Golub & van Loan, 1996), where superscript *T* denotes the transpose and *λi*, ∀*i* = 1...*m*, are the eigenvalues of **C**χ. In this case, we allow **C**<sup>χ</sup> to have zero-eigenvalues. But however, as we will explain later, all the eigenvalues must be either zeros or positives. Unluckily, there will exist particular cases where **C**<sup>χ</sup> might have any negative eigenvalue, as was indicated above. When this circumstance occurs, an adjustment of the Schur decomposition will be accomplished in order to obtain a positive semi-definite approximation, **<sup>C</sup>**<sup>χ</sup>={<sup>c</sup>*χ*(i,j)}, to the original auto-covariance matrix of the log-amplitude, **<sup>C</sup>**χ, by minimizing

*δ*(**C**χ) = min

semi-definite matrix if any eigenvalue of **C**<sup>χ</sup> is negative. In this case,

**Q** (**Λ** + **Υ**)

**C**<sup>χ</sup> = **KK***<sup>T</sup>* =

eigenvalue for 0. The resulting matrix may be written as

**<sup>C</sup>**<sup>χ</sup> <sup>=</sup> **KK***<sup>T</sup>* <sup>=</sup>

Let us see these two scenarios in detail:

**4.2.1 General scenario**

however, if all eigenvalues of **C**<sup>χ</sup> are positive, then

**C**<sup>χ</sup>=**C** *<sup>T</sup>* <sup>χ</sup>≥0

where �·� represents the norm of the matrix **C**<sup>χ</sup> − **C**<sup>χ</sup> . We distinguish two main situations (Jurado-Navas & Puerta-Notario, 2009) corresponding to minimize the Frobenius norm (Golub & van Loan, 1996) or the the *p*-norms (with *p* = 1, 2, ∞), respectively. In both cases, we can substitute the diagonal matrix **Λ** by a resulting matrix written as **R** = **Λ** + **Υ**, where **Υ** is a perturbation matrix obtained after minimizing Eq. (27) in any of the two situations mentioned above. After that, and for any of these two scenarios, we can form the coloring matrix, **K**, as **K** = **Q**(**R**)**1**/**<sup>2</sup>** in order to generate the correlated log-normal scintillation samples. Here it is shown why we need that **C**<sup>χ</sup> were approximated by a positive

1/2 **<sup>Q</sup>** (**<sup>Λ</sup>** <sup>+</sup> **<sup>Υ</sup>**)

**QΛ**1/2 **QΛ**1/2*<sup>T</sup>*

If any *λ<sup>i</sup>* in **Λ** is negative, then we can approximate the matrix **Λ** substituting every negative

1/2*<sup>T</sup>*

**R** = **Λ** + **Υ**, (30)

**Q***T***C**χ**Q** = **Λ** = diag(*λ*1, ...*λm*) (26)

�**C**<sup>χ</sup> − **C**<sup>χ</sup>�, (27)

= **QRQ***T*; (28)

= **QΛQ***T*. (29)

correlation matrix must always be.

the distance (Halmos, 1972):

can be accomplished as was indicated in (Baddour & Beaulieu, 2002). But first we need to compute the factorization **C**<sup>w</sup> = **LL***T*. If **C**<sup>w</sup> is positive definite, a Cholesky decomposition is performable using the interval Cholesky method proposed in (Alefeld & Mayer, 1993) in order to include the desired correlation among scintillation sequences, with **L** being a lower triangular matrix (Beaulieu & Merani, 2000; Ertel & Reed, 1998). So, the driving process is then generated by the product ω[*n*] = **Lz**[*n*], where **z**[*n*] is an *m* × 1 vector of independent zero mean Gaussian variates with unit variance and a autocorrelation function expressed as:

$$\mathbf{R}\_{ZZ} = E\{\mathbf{z}[n]\mathbf{z}^T[n]\} = \mathbf{I}\_{\mathbf{m}}.\tag{24}$$

where **Im** is the 'm'-element identity matrix. Finally,

$$\chi[n] = -\sum\_{k=1}^{p} \mathbf{A}[k] \chi[n-k] + \mathbf{w}[n]. \tag{25}$$

But if **C**<sup>w</sup> is not positive definite, then an alternative efficient decomposition algorithm that let the matrix be factorized was proposed in (Jurado-Navas & Puerta-Notario, 2009).

In addition, as said in the introduction, we can simplify the model proposed in (Jurado-Navas & Puerta-Notario, 2009) by an unrealistic but reasonably accurate space-time separable statistics model with a reduced computational load if it is compared with the complete model. Hence, we work directly with the covariance matrix, **C**χ, unlike the work proposed in (Jurado-Navas & Puerta-Notario, 2009), where the matrix to be decomposed is **C**w.

Thus, when **C**<sup>χ</sup> is positive definite, as in the Gaussian approximation employed in (Zhu & Kahn, 2002), a Cholesky decomposition is performable using the interval Cholesky method proposed in (Alefeld & Mayer, 1993) in order to include the desired correlation among scintillation sequences. Hence we find a lower triangular matrix **L** such that **C**<sup>χ</sup> = **LL***<sup>H</sup>* in a similar way as in (Beaulieu & Merani, 2000; Ertel & Reed, 1998).

#### **4.2 Constructing the coloring matrix**

As we have already said, the AR model presented in (Jurado-Navas & Puerta-Notario, 2009) can be simplified by directly coloring independent Gaussian sequences first between them and then in time (space-time separable statistical model), avoiding the realization of the AR model. Obviously, this procedure is unrealistic because it implies that Taylor's hypothesis is not considered. However, numerical results obtained (and included through this chapter) show that this simplification behaves as a good approximation in terms of burst error rate.

Thus, we can directly consider the auto-covariance matrix of the log-amplitude, **C**χ. As commented before, when **C**<sup>χ</sup> is positive definite, a Cholesky decomposition is performable by using the interval Cholesky method proposed (Alefeld & Mayer, 1993) in order to include the desired correlation among scintillation sequences. Hence we find a lower triangular matrix, **L**, as a coloring matrix such that **C**<sup>χ</sup> = **LL***<sup>H</sup>* in a similar way as in (Ertel & Reed, 1998).

Nevertheless, **C**<sup>χ</sup> may not be positive definite so the Cholesky decomposition is not always feasible. Although any theoretical covariance matrix must be positive or, at least, semipositive definite, when such matrices are processed on a computer, small deviations from theory are introduced by the limits of floating point computations, inducing numeric inconsistencies. In this respect, there are no guarantees that the finite precision representation of the matrix can ensure positive definiteness, even more when the Cholesky algorithm is unstable for positive definite matrices that have one or more eigenvalues close to 0, as it is explained in 8 Numerical Simulations / Book 1

can be accomplished as was indicated in (Baddour & Beaulieu, 2002). But first we need to compute the factorization **C**<sup>w</sup> = **LL***T*. If **C**<sup>w</sup> is positive definite, a Cholesky decomposition is performable using the interval Cholesky method proposed in (Alefeld & Mayer, 1993) in order to include the desired correlation among scintillation sequences, with **L** being a lower triangular matrix (Beaulieu & Merani, 2000; Ertel & Reed, 1998). So, the driving process is then generated by the product ω[*n*] = **Lz**[*n*], where **z**[*n*] is an *m* × 1 vector of independent zero mean Gaussian variates with unit variance and a autocorrelation function expressed as:

where **Im** is the 'm'-element identity matrix. Finally,

**C**w.

χ[*n*] = −

similar way as in (Beaulieu & Merani, 2000; Ertel & Reed, 1998).

**4.2 Constructing the coloring matrix**

*p* ∑ *k*=1

let the matrix be factorized was proposed in (Jurado-Navas & Puerta-Notario, 2009).

But if **C**<sup>w</sup> is not positive definite, then an alternative efficient decomposition algorithm that

In addition, as said in the introduction, we can simplify the model proposed in (Jurado-Navas & Puerta-Notario, 2009) by an unrealistic but reasonably accurate space-time separable statistics model with a reduced computational load if it is compared with the complete model. Hence, we work directly with the covariance matrix, **C**χ, unlike the work proposed in (Jurado-Navas & Puerta-Notario, 2009), where the matrix to be decomposed is

Thus, when **C**<sup>χ</sup> is positive definite, as in the Gaussian approximation employed in (Zhu & Kahn, 2002), a Cholesky decomposition is performable using the interval Cholesky method proposed in (Alefeld & Mayer, 1993) in order to include the desired correlation among scintillation sequences. Hence we find a lower triangular matrix **L** such that **C**<sup>χ</sup> = **LL***<sup>H</sup>* in a

As we have already said, the AR model presented in (Jurado-Navas & Puerta-Notario, 2009) can be simplified by directly coloring independent Gaussian sequences first between them and then in time (space-time separable statistical model), avoiding the realization of the AR model. Obviously, this procedure is unrealistic because it implies that Taylor's hypothesis is not considered. However, numerical results obtained (and included through this chapter) show that this simplification behaves as a good approximation in terms of burst error rate. Thus, we can directly consider the auto-covariance matrix of the log-amplitude, **C**χ. As commented before, when **C**<sup>χ</sup> is positive definite, a Cholesky decomposition is performable by using the interval Cholesky method proposed (Alefeld & Mayer, 1993) in order to include the desired correlation among scintillation sequences. Hence we find a lower triangular matrix, **L**, as a coloring matrix such that **C**<sup>χ</sup> = **LL***<sup>H</sup>* in a similar way as in (Ertel & Reed, 1998). Nevertheless, **C**<sup>χ</sup> may not be positive definite so the Cholesky decomposition is not always feasible. Although any theoretical covariance matrix must be positive or, at least, semipositive definite, when such matrices are processed on a computer, small deviations from theory are introduced by the limits of floating point computations, inducing numeric inconsistencies. In this respect, there are no guarantees that the finite precision representation of the matrix can ensure positive definiteness, even more when the Cholesky algorithm is unstable for positive definite matrices that have one or more eigenvalues close to 0, as it is explained in

**<sup>R</sup>**zz <sup>=</sup> *<sup>E</sup>*{**z**[*n*]**z***T*[*n*]} <sup>=</sup> **Im**, (24)

**A**[*k*]χ[*n* − *k*] + **w**[*n*]. (25)

(Holton, 2004). In (Abramovich et al., 2001), the particular problem of non-positive definite matrices that arise from characterizing a realistic system is investigated in detail, showing that in 988 outcomes of their 1000 Monte Carlo trials, gave rise to a non-positive definite covariance matrix, although the corresponding "contracted" matrix may well be positive definite. In addition, in most practical situations the covariance matrix must be estimated from scintillation measurements (Monserrat et al., 2007; Moore et al., 2005) and outliers and other special situations as urban canopies may induce that correlation matrices derived from measurements are not always guaranteed to be positive definite, although a theoretical correlation matrix must always be.

In such cases, a real Schur decomposition is proposed (Golub & van Loan, 1996). Since the symmetric matrix **<sup>C</sup>**<sup>χ</sup> ∈ �*mxm*, then there exists an orthogonal matrix, **<sup>Q</sup>** ∈ �*mxm*, such that

$$\mathbf{Q}^T \mathbf{C}\_\mathbf{X} \mathbf{Q} = \boldsymbol{\Lambda} = \text{diag}(\lambda\_1, \dots \lambda\_m) \tag{26}$$

as was shown in (Golub & van Loan, 1996), where superscript *T* denotes the transpose and *λi*, ∀*i* = 1...*m*, are the eigenvalues of **C**χ. In this case, we allow **C**<sup>χ</sup> to have zero-eigenvalues. But however, as we will explain later, all the eigenvalues must be either zeros or positives. Unluckily, there will exist particular cases where **C**<sup>χ</sup> might have any negative eigenvalue, as was indicated above. When this circumstance occurs, an adjustment of the Schur decomposition will be accomplished in order to obtain a positive semi-definite approximation, **<sup>C</sup>**<sup>χ</sup>={<sup>c</sup>*χ*(i,j)}, to the original auto-covariance matrix of the log-amplitude, **<sup>C</sup>**χ, by minimizing the distance (Halmos, 1972):

$$\delta(\mathbf{C}\_{\mathbf{X}}) = \min\_{\tilde{\mathbf{C}}\_{\mathbf{X}} = \tilde{\mathbf{C}}\_{\mathbf{X}}^{\top} \ge 0} \| \mathbf{C}\_{\mathbf{X}} - \tilde{\mathbf{C}}\_{\mathbf{X}} \|\_{\prime} \tag{27}$$

where �·� represents the norm of the matrix **C**<sup>χ</sup> − **C**<sup>χ</sup> . We distinguish two main situations (Jurado-Navas & Puerta-Notario, 2009) corresponding to minimize the Frobenius norm (Golub & van Loan, 1996) or the the *p*-norms (with *p* = 1, 2, ∞), respectively. In both cases, we can substitute the diagonal matrix **Λ** by a resulting matrix written as **R** = **Λ** + **Υ**, where **Υ** is a perturbation matrix obtained after minimizing Eq. (27) in any of the two situations mentioned above. After that, and for any of these two scenarios, we can form the coloring matrix, **K**, as **K** = **Q**(**R**)**1**/**<sup>2</sup>** in order to generate the correlated log-normal scintillation samples. Here it is shown why we need that **C**<sup>χ</sup> were approximated by a positive semi-definite matrix if any eigenvalue of **C**<sup>χ</sup> is negative. In this case,

$$\tilde{\mathbf{C}}\_{\mathbf{X}} = \mathbf{K} \mathbf{K}^T = \left( \mathbf{Q} \left( \mathbf{A} + \mathbf{Y} \right)^{1/2} \right) \left( \mathbf{Q} \left( \mathbf{A} + \mathbf{Y} \right)^{1/2} \right)^T = \mathbf{Q} \mathbf{R} \mathbf{Q}^T;\tag{28}$$

however, if all eigenvalues of **C**<sup>χ</sup> are positive, then

$$\mathbf{C}\_{\mathbf{X}} = \mathbf{K} \mathbf{K}^{T} = \left(\mathbf{Q} \boldsymbol{\Lambda}^{1/2}\right) \left(\mathbf{Q} \boldsymbol{\Lambda}^{1/2}\right)^{T} = \mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^{T}.\tag{29}$$

Let us see these two scenarios in detail:

#### **4.2.1 General scenario**

If any *λ<sup>i</sup>* in **Λ** is negative, then we can approximate the matrix **Λ** substituting every negative eigenvalue for 0. The resulting matrix may be written as

$$\mathbf{R} = \mathbf{A} + \mathbf{Y}\_{\prime} \tag{30}$$

Finally, if the original auto-covariance matrix is not symmetric, then the Schur decomposition

<sup>439</sup> Upper Burst Error Bound for Atmospheric Correlated

where **<sup>U</sup>** ∈ �*nxn* is upper triangular as was shown in (Golub & van Loan, 1996). In these circumstances, we can form the matrix **Λ** taking into account solely the diagonal elements of

**K** = **Q** (**Λ** + **Υ**)

= **Q** (**Λ** + **Υ**)

*T*

where **Im** is the 'm'-element identity matrix. By calculating **w**[*n*] = **Kz**[*n*], being

the desired auto-covariance matrix of the log-amplitude fluctuations is obtained since

Furthermore, every *<sup>w</sup>*(*i*)[*n*], <sup>∀</sup>*<sup>i</sup>* <sup>=</sup> 1...*<sup>m</sup>* in Eq. (40) remains statistically Gaussian so it is possible to filter them in the way proposed in (Jurado-Navas et al., 2007) in order to obtain *<sup>α</sup>*(*i*)

each of the *m* receivers. The filter employed in each branch is the Gaussian filter shown in (Jurado-Navas et al., 2007). Therefore, the output signal of these filters is the log-amplitude fluctuations, *<sup>χ</sup>*(*i*)(*t*) <sup>∀</sup>*<sup>i</sup>* <sup>=</sup> 1...*m*, for each receiver. These outputs exhibit a Gaussian probability density function and have the desired covariance matrix, **C**χ, among sequences but without its inherent time-dependencies that would be enforced by the frozen-in hypothesis, according to the space-time separable statistics approximation imposed in this chapter. Next, its probability density function is converted from Gaussian to lognormally distributed, generally accepted for the irradiance fluctuations under weak turbulence conditions. Figure 1 shows

In this respect, it is possible to employ the same procedure explained in this paper but knowing that, instead of the log-normal employed here, a Beckman probability density much

in order to generate the correlated log-normal scintillation samples. Then, it follows that

as was shown in Eq. (28), where **C**<sup>χ</sup> = **C**<sup>χ</sup> if there does not exist any negative element in diagonal matrix **Λ** written in (26). Equation (37) shows why we need that **C**<sup>χ</sup> were approximated by a positive semi-definite matrix if any eigenvalue of **C**<sup>χ</sup> is negative even when a Schur decomposition could be always accomplished if the starting matrix is squared.

*w*(1)[*n*], *w*(2)[*n*], ... *w*(*m*)[*n*]

1 2 (**Λ** + **Υ**) 1 2 *T*

**RZZ** <sup>=</sup> *<sup>E</sup>*{**z**[*n*]**z***T*[*n*]} <sup>=</sup> **Im**, (39)

*T*

**U** before applying any of the two approximations proposed in this paper.

1 2 *T*

**4.3 Generation of turbulent spatially-correlated channels**

Optical Communications Using an Alternative Matrix Decomposition

**Q** (**Λ** + **Υ**)

*z*(1)[*n*], *z*(2)[*n*], ... *z*(*m*)[*n*]

**<sup>w</sup>**[*n*] =

*<sup>E</sup>*{**w**[*n*] · **<sup>w</sup>***T*[*n*]} <sup>=</sup> *<sup>E</sup>*{**<sup>K</sup>** · **<sup>z</sup>**[*n*]·**z**[*n*]

Now we can form the coloring matrix, **<sup>K</sup>** ∈ �*mxm*

(Jurado-Navas & Puerta-Notario, 2009)

Its correlation matrix, *RZZ*, is given by

1 2

**Q** (**Λ** + **Υ**)

**C**<sup>χ</sup> = **QUQ***T*, (36)

1/2 , (37)

as a set of uncorrelated white Gaussian signals.

**<sup>Q</sup>***<sup>T</sup>* <sup>=</sup> **<sup>Q</sup>**(**<sup>Λ</sup>** <sup>+</sup> **<sup>Υ</sup>**)**Q***<sup>T</sup>* <sup>=</sup> **<sup>C</sup>**<sup>χ</sup>,

, (40)

*sc* (*t*) for

*<sup>T</sup>* · **<sup>K</sup>***T*} <sup>=</sup> **KK***<sup>T</sup>* <sup>=</sup> **<sup>C</sup>**<sup>χ</sup>. (41)

(38)

is written as

**KK***<sup>T</sup>* =

Next, define **Z** =

clearly the overall process.

where **<sup>Υ</sup>** is a diagonal perturbation matrix with diagonal elements, {*υi*}*<sup>m</sup> <sup>i</sup>*=1, of

$$\upsilon\_{i} = \begin{cases} -\lambda\_{i\prime} \text{ if } \lambda\_{i} < 0; \\ 0, \text{ if } \lambda\_{i} \ge 0; \end{cases} \qquad \forall \mathbf{i} = 1, 2...\text{m}. \tag{31}$$

In this case, the approximation **<sup>C</sup>**�<sup>χ</sup> <sup>=</sup> **QRQ***<sup>T</sup>* is proven to optimize the Frobenius norm if (*vi* + *λi*) → 0, and with the constraint of having a positive semi-definite solution matrix, **C**�χ; where the Frobenius norm is detailed in (Golub & van Loan, 1996) and represented by

$$\|\mathbf{A}\|\_{F} = \sqrt{\sum\_{i=1}^{m} \sum\_{j=1}^{m} |\mathbf{c}\_{\mathcal{X}}(\mathbf{i}, \mathbf{j}) - \widetilde{\mathbf{c}}\_{\mathcal{X}}(\mathbf{i}, \mathbf{j})|^{2}}. \tag{32}$$

#### **4.2.2 Specific scenario**

If instead of employing the Frobenius norm, the expression (27) is optimized for the relevant 2−norm matrix defined as:

$$\|\mathbf{C}\_{\mathbf{X}} - \tilde{\mathbf{C}}\_{\mathbf{X}}\|\_{2} = \sqrt{\rho \left\{ \left( \mathbf{C}\_{\mathbf{X}} - \tilde{\mathbf{C}}\_{\mathbf{X}} \right)^{H} \left( \mathbf{C}\_{\mathbf{X}} - \tilde{\mathbf{C}}\_{\mathbf{X}} \right) \right\}}\,\tag{33}$$

then we can obtain a new and different solution for the same problem. In equation (33), *ρ* �� **C**<sup>χ</sup> − **C**�<sup>χ</sup> �*H*� **C**<sup>χ</sup> − **C**�<sup>χ</sup> �� is the spectral radius of the matrix � **C**<sup>χ</sup> − **C**�<sup>χ</sup> �*H*� **C**<sup>χ</sup> − **C**�<sup>χ</sup> � , defined as the supremum among the absolute values of its eigenvalues. Obviously, the hermitian operator can be substituted by a transpose operator since the samples are all real. As a characteristic feature, optimizing the 2−norm entails, in general, multiple solutions for **C**�χ, as a main difference with respect to the previous situation.

Again, we define **R** = **Λ** + **Υ**, being **Υ** the perturbation matrix. One of the possible solutions for **Υ** that let **R** be a semi-definite positive matrix is obtained by constructing **Υ** as a diagonal matrix whose elements are defined as {*υ<sup>i</sup>* = −*λ* ∀ i = 1...*m*}, being *λ* the highest among the absolute values of the negative eigenvalues of the matrix **C**χ. Obviously, as in the above scenario, we obtain greater accuracies when smaller magnitudes were summed to **Λ**. The limit, again, is the necessary magnitude that let the negative eigenvalue be canceled and so, **C**<sup>χ</sup> were approximated to the closest positive semi-definite matrix. In this way, this second scenario is proven to minimize the *p*-norms (with *p* = 1, 2, ∞) since **Υ** is a diagonal matrix in which all its entries have the same value. These norms can be consulted in (Golub & van Loan, 1996).

As final comment to this second scenario, we can remark that the great resemblance obtained between **C**<sup>χ</sup> and **C**�<sup>χ</sup> is based on the orthogonality of the rows of the matrix **Q**. For example, suppose that the covariance matrix is non-positive definite expressed as

$$\mathbf{C}\_{\mathbf{X}} = \begin{pmatrix} 0.100 \ 0.079 \ 0.010 \\ 0.079 \ 0.100 \ 0.079 \\ 0.010 \ 0.079 \ 0.100 \end{pmatrix}\_{3 \times 3} \tag{34}$$

Hence, applying the decomposition shown in scenario 2, we obtain that **C**�<sup>χ</sup> = **C**χ, except for the principal diagonal, that now is �*cχ*(*i*, *<sup>i</sup>*) = 0.106 <sup>∀</sup>*<sup>i</sup>* <sup>=</sup> 1...3, i.e.:

$$\mathbf{C\_{X}} = \begin{pmatrix} 0.106 \ 0.079 \ 0.010 \\ 0.079 \ 0.106 \ 0.079 \\ 0.010 \ 0.079 \ 0.106 \end{pmatrix}\_{3 \times 3} \tag{35}$$

Finally, if the original auto-covariance matrix is not symmetric, then the Schur decomposition is written as

$$\mathbf{C}\_{\mathbf{X}} = \mathbf{Q} \mathbf{U} \mathbf{Q}^{T} \,. \tag{36}$$

where **<sup>U</sup>** ∈ �*nxn* is upper triangular as was shown in (Golub & van Loan, 1996). In these circumstances, we can form the matrix **Λ** taking into account solely the diagonal elements of **U** before applying any of the two approximations proposed in this paper.

#### **4.3 Generation of turbulent spatially-correlated channels**

Now we can form the coloring matrix, **<sup>K</sup>** ∈ �*mxm*

10 Numerical Simulations / Book 1

In this case, the approximation **<sup>C</sup>**�<sup>χ</sup> <sup>=</sup> **QRQ***<sup>T</sup>* is proven to optimize the Frobenius norm if (*vi* + *λi*) → 0, and with the constraint of having a positive semi-definite solution matrix, **C**�χ; where the Frobenius norm is detailed in (Golub & van Loan, 1996) and represented by

> *m* ∑ *j*=1

If instead of employing the Frobenius norm, the expression (27) is optimized for the relevant

then we can obtain a new and different solution for the same problem. In equation (33),

defined as the supremum among the absolute values of its eigenvalues. Obviously, the hermitian operator can be substituted by a transpose operator since the samples are all real. As a characteristic feature, optimizing the 2−norm entails, in general, multiple solutions for

Again, we define **R** = **Λ** + **Υ**, being **Υ** the perturbation matrix. One of the possible solutions for **Υ** that let **R** be a semi-definite positive matrix is obtained by constructing **Υ** as a diagonal matrix whose elements are defined as {*υ<sup>i</sup>* = −*λ* ∀ i = 1...*m*}, being *λ* the highest among the absolute values of the negative eigenvalues of the matrix **C**χ. Obviously, as in the above scenario, we obtain greater accuracies when smaller magnitudes were summed to **Λ**. The limit, again, is the necessary magnitude that let the negative eigenvalue be canceled and so, **C**<sup>χ</sup> were approximated to the closest positive semi-definite matrix. In this way, this second scenario is proven to minimize the *p*-norms (with *p* = 1, 2, ∞) since **Υ** is a diagonal matrix in which all its entries have the same value. These norms can be consulted in (Golub & van Loan,

As final comment to this second scenario, we can remark that the great resemblance obtained between **C**<sup>χ</sup> and **C**�<sup>χ</sup> is based on the orthogonality of the rows of the matrix **Q**. For example,

> 0.100 0.079 0.010 0.079 0.100 0.079 0.010 0.079 0.100

Hence, applying the decomposition shown in scenario 2, we obtain that **C**�<sup>χ</sup> = **C**χ, except for

0.106 0.079 0.010 0.079 0.106 0.079 0.010 0.079 0.106

⎞ ⎠ 3×3

⎞ ⎠ 3×3

**C**<sup>χ</sup> − **C**�<sup>χ</sup>

is the spectral radius of the matrix �

�*H*�

**C**<sup>χ</sup> − **C**�<sup>χ</sup>

*<sup>i</sup>*=1, of

<sup>|</sup>c*χ*(i,j) <sup>−</sup>�c*χ*(i,j)|2. (32)

��

**C**<sup>χ</sup> − **C**�<sup>χ</sup>

, (33)

�*H*�

. (34)

. (35)

**C**<sup>χ</sup> − **C**�<sup>χ</sup>

� ,

0, if *<sup>λ</sup><sup>i</sup>* <sup>≥</sup> 0; <sup>∀</sup><sup>i</sup> <sup>=</sup> 1, 2...m. (31)

where **<sup>Υ</sup>** is a diagonal perturbation matrix with diagonal elements, {*υi*}*<sup>m</sup>*

�**A**�*<sup>F</sup>* =

�**C**<sup>χ</sup> − **C**�χ�<sup>2</sup> =

**C**�χ, as a main difference with respect to the previous situation.

suppose that the covariance matrix is non-positive definite expressed as

⎛ ⎝

⎛ ⎝

**C**<sup>χ</sup> =

the principal diagonal, that now is �*cχ*(*i*, *<sup>i</sup>*) = 0.106 <sup>∀</sup>*<sup>i</sup>* <sup>=</sup> 1...3, i.e.:

**C**<sup>χ</sup> =

��

−*λi*, if *λ<sup>i</sup>* < 0;

���� *m* ∑ *i*=1

> � *ρ* ��

*υ<sup>i</sup>* = �

**4.2.2 Specific scenario**

**C**<sup>χ</sup> − **C**�<sup>χ</sup>

*ρ* ��

1996).

2−norm matrix defined as:

�*H*�

**C**<sup>χ</sup> − **C**�<sup>χ</sup>

$$\mathbf{K} = \mathbf{Q} \left(\mathbf{A} + \mathbf{Y}\right)^{1/2},\tag{37}$$

in order to generate the correlated log-normal scintillation samples. Then, it follows that (Jurado-Navas & Puerta-Notario, 2009)

$$\mathbf{K}\mathbf{K}^T = \left(\mathbf{Q}\left(\boldsymbol{\Lambda} + \mathbf{Y}\right)^{\frac{1}{2}}\right)\left(\mathbf{Q}\left(\boldsymbol{\Lambda} + \mathbf{Y}\right)^{\frac{1}{2}}\right)^T = \mathbf{Q}\left(\boldsymbol{\Lambda} + \mathbf{Y}\right)^{\frac{1}{2}}\left(\left(\boldsymbol{\Lambda} + \mathbf{Y}\right)^{\frac{1}{2}}\right)^T \mathbf{Q}^T = \mathbf{Q}\left(\boldsymbol{\Lambda} + \mathbf{Y}\right)\mathbf{Q}^T = \tilde{\mathbf{C}}\_{\mathbf{X}'} \tag{38}$$

as was shown in Eq. (28), where **C**<sup>χ</sup> = **C**<sup>χ</sup> if there does not exist any negative element in diagonal matrix **Λ** written in (26). Equation (37) shows why we need that **C**<sup>χ</sup> were approximated by a positive semi-definite matrix if any eigenvalue of **C**<sup>χ</sup> is negative even when a Schur decomposition could be always accomplished if the starting matrix is squared.

Next, define **Z** = *z*(1)[*n*], *z*(2)[*n*], ... *z*(*m*)[*n*] *T* as a set of uncorrelated white Gaussian signals. Its correlation matrix, *RZZ*, is given by

$$\mathbf{R}\_{\mathbf{ZZ}} = \mathbb{E}\left\{ \mathbf{z}[n] \mathbf{z}^{T}[n] \right\} = \mathbf{I}\_{\mathbf{m}\prime} \tag{39}$$

where **Im** is the 'm'-element identity matrix. By calculating **w**[*n*] = **Kz**[*n*], being

$$\mathbf{w}[n] = \begin{bmatrix} w^{(1)}[n] \ \_ \boldsymbol{w}^{(2)}[n] \ \_ \times \ \_ \boldsymbol{w}^{(m)}[n] \end{bmatrix}^T \tag{40}$$

the desired auto-covariance matrix of the log-amplitude fluctuations is obtained since

$$E\{\mathbf{w}[n]\cdot\mathbf{w}^T[n]\} = E\{\mathbf{K}\cdot\mathbf{z}[n]\cdot\mathbf{z}[n]^T\cdot\mathbf{K}^T\} = \mathbf{K}\mathbf{K}^T = \tilde{\mathbf{C}}\_{\mathbf{X}}.\tag{41}$$

Furthermore, every *<sup>w</sup>*(*i*)[*n*], <sup>∀</sup>*<sup>i</sup>* <sup>=</sup> 1...*<sup>m</sup>* in Eq. (40) remains statistically Gaussian so it is possible to filter them in the way proposed in (Jurado-Navas et al., 2007) in order to obtain *<sup>α</sup>*(*i*) *sc* (*t*) for each of the *m* receivers. The filter employed in each branch is the Gaussian filter shown in (Jurado-Navas et al., 2007). Therefore, the output signal of these filters is the log-amplitude fluctuations, *<sup>χ</sup>*(*i*)(*t*) <sup>∀</sup>*<sup>i</sup>* <sup>=</sup> 1...*m*, for each receiver. These outputs exhibit a Gaussian probability density function and have the desired covariance matrix, **C**χ, among sequences but without its inherent time-dependencies that would be enforced by the frozen-in hypothesis, according to the space-time separable statistics approximation imposed in this chapter. Next, its probability density function is converted from Gaussian to lognormally distributed, generally accepted for the irradiance fluctuations under weak turbulence conditions. Figure 1 shows clearly the overall process.

In this respect, it is possible to employ the same procedure explained in this paper but knowing that, instead of the log-normal employed here, a Beckman probability density much

<sup>0</sup> 0.01 0.02 0.03 0.04 0.05 0.06 <sup>0</sup>

<sup>0</sup> 0.01 0.02 0.03 0.04 0.05 0.06 <sup>0</sup>

Combining scintillation sequence, space−time separable statistics

*<sup>χ</sup>* = 0.25

Envelope, Rx 1 Envelope, Rx 2 Envelope, Rx 3

**Time (s.)**

(b)

1

Fig. 2. Scintillation sequences generated from (a) an AR model; (b) a space-time separable statistics model, for a 3-receivers system. The EGC sequence is displayed in thicker solid line.

(Deutsch & Miller, 1981) definition of a burst error with lengths of 192, and 64 bits, not containing more than 4 consecutive correct bits (*Lb* = 5, as explained in (Deutsch & Miller,

Then, pulses with on-off keying format (OOK) and Gaussian shape (OOK-GS) with a duty cycle (d.c.) of 100% are adopted, where identical average optical power is transmitted for every simulated case representing the reference condition to establish the comparative analysis (Jurado-Navas et al., 2010). All these features are included in the system model, where its remarkable elements are: first, the channel model presented in (Jurado-Navas et al., 2007) corresponding to a turbulent atmospheric environment, but included in a *m*-branch reception as in (Jurado-Navas & Puerta-Notario, 2009), which represent the *m* different correlated turbulence-induced fadings at each of the optical receivers; secondly, a 500 kHz three-pole Bessel highpass filter for natural and artificial light adverse effects suppression; and, thirdly, a five-pole, Bessel, lowpass filter employed as a matched filter. As said above, a conventional EGC is implemented whereas the detection procedure considered is a maximum likelihood (ML) detection. The receivers employed in this paper are point receivers whereas the weather-induced attenuation is neglected so that we concentrate our attention on turbulence effects. Furthermore, the atmospheric-induced beam spreading that causes a power reduction at the receiver is also neglected because, in our specific case, we are considering a terrestrial link where beam divergence is typically on the order of 10 *μ*Rad. This temporal spreading may be considered at high data rate, as in (Jurado-Navas et al., 2009), particularly when operating in special scenarios where dust particles are likely present.

The first set of results are displayed in Figure 3 for a three-receivers system and for *σ*<sup>2</sup>

*<sup>χ</sup>* = 0.1, where the CC coefficients between scintillation sequences are:

2

3

**Scintillation envelope**

4

5

6

<sup>441</sup> Upper Burst Error Bound for Atmospheric Correlated

Combining scintillation sequence, AR model Envelope, Rx 1 Envelope, Rx 2 Envelope, Rx 3

Optical Communications Using an Alternative Matrix Decomposition

**Time (s.)**

(a)

1981)) any sequence of burst error.

**6. Numerical results**

• *ρ*12=*ρ*21=*ρ*23=*ρ*32=0.30 and *ρ*13=*ρ*31=0.008 • *ρ*12=*ρ*21=*ρ*23=*ρ*32=0.56 and *ρ*13=*ρ*31=0.1 • *ρ*12=*ρ*21=*ρ*23=*ρ*32=0.79 and *ρ*13=*ρ*31=0.1

and *σ*<sup>2</sup>

1

2

3

**Scintillation envelope**

4

5

6

Fig. 1. Block diagram representing the generation of 'm' equal power log-normal scintillation sequences with a space-time separable statistics approximation.

more accurately reflects the statistics of the intensity scintillations (Hill & Frehlich, 1997) if Rytov variance increases even beyond the limits of the weak turbulence regime; or even the recently discovered Málaga probability density function (Jurado-Navas et al., 2011b), also very accurate with the statistics of the intensity scintillations.

#### **5. System model**

To study the performance of different values of spatial CC coefficients, *ρij*, between points *i* and *j* in the receiver plane, intensity modulation and direct detection (IM/DD) links are assumed operating at a laser wavelength of 830 nm through a 250 m horizontal path at a bit rate of 50 Mbps. Assume three receivers in the system, where a conventional equal-gain combining (EGC) of diversity branches is implemented. Figure 2 shows the temporal behavior of correlated scintillation sequences for these three different receivers with a CC coefficient of *ρ*<sup>12</sup> = 0.56 and *ρ*<sup>13</sup> = 0.1 respectively, for this two extreme scenarios: an AR model, as in (Jurado-Navas & Puerta-Notario, 2009), and the space-time separable statistics approximation employed in this chapter. The component of the wind velocity, *<sup>u</sup>*⊥, transverse to the propagation direction is fixed to *<sup>u</sup>*<sup>⊥</sup> = 20 m/s, where *<sup>ρ</sup>ij* is the normalized CC coefficient between points *i* and *j* in the receiver plane. Clearly, the election of a 20 m/s wind speed in this chapter will not be by far the typical operational scenario unless one of the terminals is in motion or the optical link is settled in particular geographical locations specifically affected by strong wind (Campins et al., 2007); but, however, such value is selected for numerical convenience so that we can work with a lower computational cost, because the total amount of samples needed to obtain a scintillation sequence with the desired lognormal statistics is smaller than using more conventional values for *<sup>u</sup>*⊥, but letting us extract the same conclusions than employing a standard velocity.

Hence, a different temporal variability is shown in Figure 2 (a) and (b) for the EGC-combined scintillation sequence, being faster for the obtained sequence from a space-time separable statistics model. This different variability must have a different repercussion in the associated burst error rate curves, as shown below, where we followed Deutsch and Miller's 12 Numerical Simulations / Book 1

sc

sc

sc

Fig. 1. Block diagram representing the generation of 'm' equal power log-normal scintillation

more accurately reflects the statistics of the intensity scintillations (Hill & Frehlich, 1997) if Rytov variance increases even beyond the limits of the weak turbulence regime; or even the recently discovered Málaga probability density function (Jurado-Navas et al., 2011b), also

To study the performance of different values of spatial CC coefficients, *ρij*, between points *i* and *j* in the receiver plane, intensity modulation and direct detection (IM/DD) links are assumed operating at a laser wavelength of 830 nm through a 250 m horizontal path at a bit rate of 50 Mbps. Assume three receivers in the system, where a conventional equal-gain combining (EGC) of diversity branches is implemented. Figure 2 shows the temporal behavior of correlated scintillation sequences for these three different receivers with a CC coefficient of *ρ*<sup>12</sup> = 0.56 and *ρ*<sup>13</sup> = 0.1 respectively, for this two extreme scenarios: an AR model, as in (Jurado-Navas & Puerta-Notario, 2009), and the space-time separable statistics approximation employed in this chapter. The component of the wind velocity, *<sup>u</sup>*⊥, transverse to the propagation direction is fixed to *<sup>u</sup>*<sup>⊥</sup> = 20 m/s, where *<sup>ρ</sup>ij* is the normalized CC coefficient between points *i* and *j* in the receiver plane. Clearly, the election of a 20 m/s wind speed in this chapter will not be by far the typical operational scenario unless one of the terminals is in motion or the optical link is settled in particular geographical locations specifically affected by strong wind (Campins et al., 2007); but, however, such value is selected for numerical convenience so that we can work with a lower computational cost, because the total amount of samples needed to obtain a scintillation sequence with the desired lognormal statistics is smaller than using more conventional values for *<sup>u</sup>*⊥, but letting us extract the same

Hence, a different temporal variability is shown in Figure 2 (a) and (b) for the EGC-combined scintillation sequence, being faster for the obtained sequence from a space-time separable statistics model. This different variability must have a different repercussion in the associated burst error rate curves, as shown below, where we followed Deutsch and Miller's

sequences with a space-time separable statistics approximation.

very accurate with the statistics of the intensity scintillations.

conclusions than employing a standard velocity.

**5. System model**

Fig. 2. Scintillation sequences generated from (a) an AR model; (b) a space-time separable statistics model, for a 3-receivers system. The EGC sequence is displayed in thicker solid line.

(Deutsch & Miller, 1981) definition of a burst error with lengths of 192, and 64 bits, not containing more than 4 consecutive correct bits (*Lb* = 5, as explained in (Deutsch & Miller, 1981)) any sequence of burst error.

Then, pulses with on-off keying format (OOK) and Gaussian shape (OOK-GS) with a duty cycle (d.c.) of 100% are adopted, where identical average optical power is transmitted for every simulated case representing the reference condition to establish the comparative analysis (Jurado-Navas et al., 2010). All these features are included in the system model, where its remarkable elements are: first, the channel model presented in (Jurado-Navas et al., 2007) corresponding to a turbulent atmospheric environment, but included in a *m*-branch reception as in (Jurado-Navas & Puerta-Notario, 2009), which represent the *m* different correlated turbulence-induced fadings at each of the optical receivers; secondly, a 500 kHz three-pole Bessel highpass filter for natural and artificial light adverse effects suppression; and, thirdly, a five-pole, Bessel, lowpass filter employed as a matched filter. As said above, a conventional EGC is implemented whereas the detection procedure considered is a maximum likelihood (ML) detection. The receivers employed in this paper are point receivers whereas the weather-induced attenuation is neglected so that we concentrate our attention on turbulence effects. Furthermore, the atmospheric-induced beam spreading that causes a power reduction at the receiver is also neglected because, in our specific case, we are considering a terrestrial link where beam divergence is typically on the order of 10 *μ*Rad. This temporal spreading may be considered at high data rate, as in (Jurado-Navas et al., 2009), particularly when operating in special scenarios where dust particles are likely present.

#### **6. Numerical results**

The first set of results are displayed in Figure 3 for a three-receivers system and for *σ*<sup>2</sup> *<sup>χ</sup>* = 0.25 and *σ*<sup>2</sup> *<sup>χ</sup>* = 0.1, where the CC coefficients between scintillation sequences are:


−15 −10 −5 <sup>0</sup> <sup>5</sup> <sup>10</sup> 10−6

(a)

formats have been displayed in (b) for *σ*<sup>2</sup>

σχ 2=0.1

 **Normalized average optical power, dB**

σχ 2=0.5

10−5

10−4

 OOK−GScc (d.c. 100%)

OOK−GSc (d.c. 100%)

OOK−GS (d.c. 100%)

Space time separable statistics AR model

 **Burst Error Rate (average)**

Fig. 4. Burst error rate of a system with 3 receivers versus normalized average optical power using OOK format, EGC and ML detection, assuming *ρ*<sup>12</sup> = 0.56 and *ρ*<sup>13</sup> = 0.1, being *<sup>u</sup>*<sup>⊥</sup> = 20 m/s, whereas the burst error length is established to 192 bits. An OOK-GS format

*<sup>χ</sup>* = 0.15.

the inherent cuadratic response of the optical detector and the average power constraints previously mentioned (Jurado-Navas et al., 2010). For this reason, we use the increase in PAOPR as a figure of merit to compare the different pulses performance, taking into account that a burst error rate analysis must be performed to take into account the temporal variability of different combined sequences and verify that the increase in distortion of shorter pulses does not counteract the PAOPR benefits. In this sense, pulses have been modified by varying their statistics of the amplitude sequence with the purpose of increasing this PAOPR, for example with a reduced d.c. of 25%; or even using OOK-GS formats with memory to avoid the appearance of more than one pulse in sets of two (OOK-GSc) and three (OOK-GScc) consecutive symbol periods (Jurado-Navas et al., 2010) but maintaining the average optical power at the same constant level in all cases. The inclusion of memory on OOK formats allows to increase the separation in burst error rate between the AR model and space-time separable statistics one. The same conclusion may be deduced when different combining techniques are employed (an EGC technique offers more distant results than a select best scheme); or when a different number of receivers make part of the system, or even when the definition of a burst error is modified allowing to contain, for instance, up to 9 consecutive correct bits any sequence of burst error, as displayed in Figure 4 (b). Thus, when we get more complicated the transmission format with the aim of obtaining a better performance, then the difference in burst error rate between considering (AR model) or not (separable statistics model) the frozen-in hypothesis is getting increased in a more meaningful way, but offering this separable statistics model a reasonably accurate upper error limit in terms of burst error

When technical specifications may not permit sufficient receiver spacing, scintillation sequences may be spatially correlated. For these cases, an efficient method for generating an

10−3

10−2

10−1

<sup>443</sup> Upper Burst Error Bound for Atmospheric Correlated

−15 −10 −5 <sup>0</sup> <sup>5</sup> 10−6

(b)

 **Normalized average optical power, dB**

*<sup>χ</sup>*. Different transmission

OOK−GS y SELECT BEST (d.c. 100%)

> OOK−GS (d.c. 25%)

OOK−GS 2 Rx (d.c. 100%)

OOK−GS (d.c. 100%) Lb = 10

σχ 2=0.25

with d.c. of 100% is represented in (a), for different values of *σ*<sup>2</sup>

rate in all simulated cases with a reduced computational load.

**7. Concluding remarks**

σχ 2=0.15

Optical Communications Using an Alternative Matrix Decomposition

10−5

Space time separable statistics AR model

σχ 2=0.01

10−4

10−3

 **Burst Error Rate (average)**

10−2

10−1

Fig. 3. Burst error rate of a system with 3 receivers versus normalized average optical power using OOK and ML detection, for different values of *ρij* and *σ*<sup>2</sup> *<sup>χ</sup>*, with a wind velocity of *<sup>u</sup>*<sup>⊥</sup> = 20 m/s. The burst error length is established to (a) 192 bits, (b) 64 bits.

where, for the last case, we need to calculate the positive semidefinite approximation to the original autocovariance matrix of the log-amplitude, as was explained in this chapter for a space-time separable statistics model and in (Jurado-Navas & Puerta-Notario, 2009) for an AR model, because the original autocovariance matrix is not positive definite and can not be factorized by a classical Cholesky decomposition method. Anyway, the obtainment of the covariance matrix involving the process in both scenarios must be factored to finally build the set of Gaussian log-amplitude sequences measured at the different receivers in the system. As shown in Figure 3, there only exists a slight difference (approx. 1–2 optical dB at a burst error rate of 10−<sup>6</sup> for a length of burst of 192 bits) between the obtained performance from the AR model (*ρ<sup>l</sup>* → 0) and the one obtained from the space-time separable statistics model (*ρ<sup>l</sup>* → 1) proposed in this paper. In fact, this latter model offers an accurate upper error bound for the link performance in terms of burst error rate. Similar conclusions may be deduced from the curves simulated with a length of burst of 64 bits. Nevertheless, from this Figure 3, we can observe that diversity reception can improve the performance of the link, but a gain penalty is shown when the CC coefficient between two receivers is high, especially if the log-amplitude variance, *σ*<sup>2</sup> *<sup>χ</sup>*, is larger. This fact can be a plausible reason to include the consideration of the spatial channel coherence as a key factor to fully evaluate the performance of atmospheric optical communication systems.

Furthermore, Figure 4 shows some obtained results when the CC coefficient has been established to *ρ*12=*ρ*21=*ρ*23=*ρ*32=0.56 and *ρ*13=*ρ*31=0.1, for a length of burst of 192 bits. Figure 4 (a) represents the behavior of an OOK-GS format with a 100% d.c. for different intensities of turbulence. As in Figure 3, the difference in performance is not significant for the two extreme scenarios studied in this chapter. Only when the log-amplitude variance, *σ*<sup>2</sup> *<sup>χ</sup>*, is getting stronger, such difference between scenarios gets higher, for instance, from 0.17 to 1.2 optical dB for *σ*<sup>2</sup> *<sup>χ</sup>* = 0.01 and *σ*<sup>2</sup> *<sup>χ</sup>* = 0.5 respectively at a burst error rate of 10−6.

Next, the pair of curves displayed for *σ*<sup>2</sup> *<sup>χ</sup>* = 0.15 are taken as a reference in Figure 4 (b). To show the superiority of the different pulse, we analyze the performance of Gaussian pulse shapes in terms of the peak-to-average optical power ratio (PAOPR) and burst error rate. Note that the PAOPR is a favorable characteristic in IM/DD infrared links due to 14 Numerical Simulations / Book 1

10−2

Space time separable statistics AR model

 **Burst Error Rate (average)**

Fig. 3. Burst error rate of a system with 3 receivers versus normalized average optical power

where, for the last case, we need to calculate the positive semidefinite approximation to the original autocovariance matrix of the log-amplitude, as was explained in this chapter for a space-time separable statistics model and in (Jurado-Navas & Puerta-Notario, 2009) for an AR model, because the original autocovariance matrix is not positive definite and can not be factorized by a classical Cholesky decomposition method. Anyway, the obtainment of the covariance matrix involving the process in both scenarios must be factored to finally build the set of Gaussian log-amplitude sequences measured at the different receivers in the system. As shown in Figure 3, there only exists a slight difference (approx. 1–2 optical dB at a burst error rate of 10−<sup>6</sup> for a length of burst of 192 bits) between the obtained performance from the AR model (*ρ<sup>l</sup>* → 0) and the one obtained from the space-time separable statistics model (*ρ<sup>l</sup>* → 1) proposed in this paper. In fact, this latter model offers an accurate upper error bound for the link performance in terms of burst error rate. Similar conclusions may be deduced from the curves simulated with a length of burst of 64 bits. Nevertheless, from this Figure 3, we can observe that diversity reception can improve the performance of the link, but a gain penalty is shown when the CC coefficient between two receivers is high, especially if the log-amplitude

*<sup>χ</sup>*, is larger. This fact can be a plausible reason to include the consideration of the

spatial channel coherence as a key factor to fully evaluate the performance of atmospheric

Furthermore, Figure 4 shows some obtained results when the CC coefficient has been established to *ρ*12=*ρ*21=*ρ*23=*ρ*32=0.56 and *ρ*13=*ρ*31=0.1, for a length of burst of 192 bits. Figure 4 (a) represents the behavior of an OOK-GS format with a 100% d.c. for different intensities of turbulence. As in Figure 3, the difference in performance is not significant for the two

stronger, such difference between scenarios gets higher, for instance, from 0.17 to 1.2 optical

*<sup>χ</sup>* = 0.5 respectively at a burst error rate of 10−6.

To show the superiority of the different pulse, we analyze the performance of Gaussian pulse shapes in terms of the peak-to-average optical power ratio (PAOPR) and burst error rate. Note that the PAOPR is a favorable characteristic in IM/DD infrared links due to

extreme scenarios studied in this chapter. Only when the log-amplitude variance, *σ*<sup>2</sup>

10−1

100

−12 −10 −8 −6 −4 −2 <sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>6</sup> 10−3

(b)

*<sup>χ</sup>* = 0.15 are taken as a reference in Figure 4 (b).

2=0.1

 **Normalized average optical power, dB**

*<sup>χ</sup>*, with a wind velocity of

<sup>2</sup> <sup>σ</sup> =0.25 <sup>χ</sup>

<sup>ρ</sup>12=0.79 y ρ13=0.1 <sup>ρ</sup>12=0.56 y ρ13=0.1 <sup>ρ</sup>12=0.3 y ρ13=0.008

σχ

*<sup>χ</sup>*, is getting

<sup>ρ</sup>12=0.79 y ρ13=0.1 <sup>ρ</sup>12=0.56 y ρ13=0.1 <sup>ρ</sup>12=0.3 y ρ13=0.008

2=0.25

−12 −10 −8 −6 −4 −2 <sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>6</sup> 10−6

(a)

 **Normalized average optical power, dB**

2=0.1 σχ

using OOK and ML detection, for different values of *ρij* and *σ*<sup>2</sup>

*<sup>u</sup>*<sup>⊥</sup> = 20 m/s. The burst error length is established to (a) 192 bits, (b) 64 bits.

10−5

variance, *σ*<sup>2</sup>

dB for *σ*<sup>2</sup>

optical communication systems.

*<sup>χ</sup>* = 0.01 and *σ*<sup>2</sup>

Next, the pair of curves displayed for *σ*<sup>2</sup>

Space time separable statistics AR model

σχ

10−4

 **Burst Error Rate (average)**

10−3

10−2

10−1

Fig. 4. Burst error rate of a system with 3 receivers versus normalized average optical power using OOK format, EGC and ML detection, assuming *ρ*<sup>12</sup> = 0.56 and *ρ*<sup>13</sup> = 0.1, being *<sup>u</sup>*<sup>⊥</sup> = 20 m/s, whereas the burst error length is established to 192 bits. An OOK-GS format with d.c. of 100% is represented in (a), for different values of *σ*<sup>2</sup> *<sup>χ</sup>*. Different transmission formats have been displayed in (b) for *σ*<sup>2</sup> *<sup>χ</sup>* = 0.15.

the inherent cuadratic response of the optical detector and the average power constraints previously mentioned (Jurado-Navas et al., 2010). For this reason, we use the increase in PAOPR as a figure of merit to compare the different pulses performance, taking into account that a burst error rate analysis must be performed to take into account the temporal variability of different combined sequences and verify that the increase in distortion of shorter pulses does not counteract the PAOPR benefits. In this sense, pulses have been modified by varying their statistics of the amplitude sequence with the purpose of increasing this PAOPR, for example with a reduced d.c. of 25%; or even using OOK-GS formats with memory to avoid the appearance of more than one pulse in sets of two (OOK-GSc) and three (OOK-GScc) consecutive symbol periods (Jurado-Navas et al., 2010) but maintaining the average optical power at the same constant level in all cases. The inclusion of memory on OOK formats allows to increase the separation in burst error rate between the AR model and space-time separable statistics one. The same conclusion may be deduced when different combining techniques are employed (an EGC technique offers more distant results than a select best scheme); or when a different number of receivers make part of the system, or even when the definition of a burst error is modified allowing to contain, for instance, up to 9 consecutive correct bits any sequence of burst error, as displayed in Figure 4 (b). Thus, when we get more complicated the transmission format with the aim of obtaining a better performance, then the difference in burst error rate between considering (AR model) or not (separable statistics model) the frozen-in hypothesis is getting increased in a more meaningful way, but offering this separable statistics model a reasonably accurate upper error limit in terms of burst error rate in all simulated cases with a reduced computational load.

#### **7. Concluding remarks**

When technical specifications may not permit sufficient receiver spacing, scintillation sequences may be spatially correlated. For these cases, an efficient method for generating an

**z**[*n*] Vector of independent zero mean Gaussian variates with unit variance.

<sup>445</sup> Upper Burst Error Bound for Atmospheric Correlated

*ρij* Normalized cross-correlation coefficient between points *i* and *j* in the receiver

*ρ<sup>l</sup>* Degree of randomness as effect of the dynamic evolution of the turbulence

Abramovich, Y. I.; Spencer, N.K. & Gorokhov, A. Y. (2001). Detection-estimation of more

Alefeld, G. & Mayer, G. (1993). The Cholesky method for interval data. *Linear Algebra and its Applications,* Vol. 194, No.22, (November 1993), pp. 161-182, ISSN 0024-3795. Al Naboulsi, M. & Sizun, H. (2004). Fog attenuation prediction for optical and infrared waves.

Andrews, L. C. & Phillips, R. L. (1998). *Laser Beam Propagation Through Random Media*, SPIE

Andrews, L. C.; Phillips, R. L. & Hopen, C. Y. (2000). Aperture Averaging of Optical

Anguita, J. A.; Neifeld, M. A. & Vasic, B. V. (2007). Spatial correlation and irradiance statistics

Baddour, K.E.. & Beaulieu, N.C. (2002). Accurate simulation of multiple cross-correlated

267–271, ISBN 0-7803-7400-2, New York, NY, USA, April 28 - May 2, 2002. Beaulieu, N. C. (1999). Generation of correlated Rayleigh fading envelopes. *IEEE Communication Letters*, Vol. 3, No. 6 (June 1999), pp. 172–174, ISSN 1089-7798. Beaulieu, N. C. & Merani, M. L. (2000). Efficient simulation of correlated diversity channels.

*Optics*, Vol. 46, No. 26 (September 2007), pp. 6561–6571, ISSN 1559-128X. Arnon, S. (2003). Effects of atmospheric turbulence and building sway on optical

*Media,* Vol. 10, No. 1 (2000), pp. 53 – 70, ISSN 1745-5049.

207–210, ISBN 0-7803-6596-8, Chicago, IL September 2000.

uncorrelated Gaussian sources than sensors in nonuniform linear antenna arrays part I: fully augmentable arrays. *IEEE Transactions on Signal Processing,* Vol. 49, No. 5

*SPIE Optical Engineering,* Vol. 43, No. 2 (February 2004), pp. 319–329, ISSN 091-3286.


Scintillations: Power Fluctuations and the Temporal Spectrum. *Waves in Random*

in a multiple-beam terrestrial free-space optical communication link. *OSA Applied*

wireless-communication systems. *Optics Letters,* Vol. 28, No. 2 (January 2003), pp.

fading Channels. *IEEE International Conference on Communications (ICC) 2002*, pp.

*IEEE Wireless Communications and Networking Conference (WCNC 2000)*, Vol. 1, pp.

*αsc*(*t*) Time-varying atmospheric scintillation sequence. *χ*(*t*) Log-amplitude fluctuation of scintillation.

Optical Communications Using an Alternative Matrix Decomposition

**Λ** Diagonal matrix containing the eigenvalues of **C***χ*. *�* Turbulent kinetic energy dissipation rate (*m*2/*s*3).

*<sup>I</sup>* Scintillation index (normalized irradiance variance).

*ψ*(**r**, *L*) Phase perturbations of Rytov approximation.

ω[*n*] Driving noise vector process of an AR model. **Υ** Perturbation matrix employed to obtain **C**<sup>χ</sup>.

(May 2001), pp. 959–971, ISSN 1053-587X.

*λ* Wavelength.

*σ*2

*σ*2

**10. References**

plane (m).

(=*ρ<sup>l</sup>* = *τ*0/*τe*).

*<sup>χ</sup>* Log-amplitude variance. *τ*<sup>0</sup> Turbulence correlation time. *τe* Lifetime of turbulent eddies.

Washington, USA.

129–131, ISSN 0146-9592.

accurate approximation of 'm' equal power log-normal scintillation sequences with any CC coefficient is proposed in this paper, overcoming the restrictions of a Cholesky decomposition. Hence, the AR-model proposed in (Jurado-Navas & Puerta-Notario, 2009) and its inherent numerically ill-conditioned covariance matrix (Baddour & Beaulieu, 2002) may be avoided in many cases when calculating burst error rate curves due to the difference between the two extreme scenarios studied in this chapter is usually limited to approximately 2 dB at a burst rate of 10−6. In this sense, the space-time separable statistics model proposed here can be used to consider spatial correlations among scintillation sequences without fear of making big mistakes and with the advantage of a reduced computational time. Thus, such separable statistics model may be seen as a highly accurate upper error bound of the whole model detailed in (Jurado-Navas & Puerta-Notario, 2009).
