Some Statistical Definitions
Ensemble averages (also called expected values) of a random process are defined as $$ \begin{equation} \left\langle f\left(u\right) \right\rangle = \int\limits_{-\infty}^{\infty} f\left(u\right) p\left(u\right) \, d u, \end{equation} $$ where \(u\) is the random variable, \(p\left(u\right)\) is its probability density function (PDF), and \(f\left(u\right)\) is an arbitrary function of \(u\). For example, the mean and mean-square value of \(u\) are given by $$ \begin{equation} \mu_u = \left\langle u \right\rangle = \int\limits_{-\infty}^{\infty} u \, p\left(u\right) \, d u \end{equation} $$ and $$ \begin{equation} \left\langle u^2 \right\rangle = \int\limits_{-\infty}^{\infty} u^2 \, p\left(u\right) \, d u, \end{equation} $$ respectively. We can combine these to compute variance \(\sigma_u^2\), given by $$ \begin{equation} \sigma_u^2 = \left\langle \left(u - \left\langle u \right\rangle \right)^2 \right\rangle. \end{equation} $$ These are called point statistics.
If we have two random variables \(u\) and \(v\) described by a joint PDF \(p\left(u,v\right)\), we can describe how they are related on average by their correlation, defined as $$ \begin{equation} \Gamma_{uv} = \left\langle u v \right\rangle = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} u v \, p\left(u,v\right) \, d u \, dv. \end{equation} $$ A similar measure of this relationship between \(u\) and \(v\) is the covariance \(C_{uv} = \Gamma_{uv} - \left\langle u \right\rangle \left\langle v \right\rangle \). Finally, we define the correlation coefficient \( \rho_{uv} = \Gamma_{uv} / \left( \left\langle u \right\rangle \left\langle v \right\rangle \right) \). When \(v = u \), these are called the autocorrelation and autocovariance, respectively.
Temporal Random Processes
Random processes are random variables that evolve over time, like \(u\left(t\right)\) and \(v\left(t\right)\). In general, the PDF and expected values are functions of time, e.g., $$ \begin{equation} \mu_u\left(t\right) = \left\langle u\left(t\right) \right\rangle = \int\limits_{-\infty}^{\infty} u\left(t\right) \, p\left(u, t\right) \, du, \end{equation}$$ $$ \begin{equation} \left\langle u^2\left(t\right) \right\rangle = \int\limits_{-\infty}^{\infty} u^2\left(t\right) \, p\left(u, t\right) \, du, \end{equation} $$ and $$ \begin{equation} \sigma_u^2\left(t\right) = \left\langle \left[u\left(t\right) - \left\langle u\left(t\right) \right\rangle \right]^2 \right\rangle. \end{equation} $$ Statistical correlation of a random process involves averages at two different moments in time, \(t_1\) and \(t_2\). Temporal correlation is a good example: $$ \begin{equation} \Gamma_{u}\left(t_1,t_2\right) = \left\langle u\left(t_1\right) u\left(t_2\right) \right\rangle = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} u\left(t_1\right) u\left(t_2\right) \, p\left[u\left(t_1\right),u\left(t_2\right)\right] \, d u\left(t_1\right) \, d u\left(t_2\right). \end{equation} $$ Also, we have the closely-related temporal covariance, which is the correlation of the process \(u\left(t\right)\) with its mean subtracted according to $$ \begin{align} C_{u}\left(t_1,t_2\right) &= \left\langle \left[u\left(t_1\right)-\mu_u\left(t_1\right)\right] \left[u\left(t_2\right)-\mu_u\left(t_2\right)\right] \right\rangle \\ &= \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} \left[u\left(t_1\right)-\mu_u\left(t_1\right)\right] \left[u\left(t_2\right)-\mu_u\left(t_2\right)\right] \, p\left[u\left(t_1\right),u\left(t_2\right)\right] \, d u\left(t_1\right) \, d u\left(t_2\right). \end{align} $$ Covariance is related to correlation according to $$ \begin{equation} C_{u}\left(t_1,t_2\right) = \Gamma_{u}\left(t_1,t_2\right) - \left\langle u\left(t_1\right) \right\rangle \left\langle u\left(t_2\right)\right\rangle. \end{equation} $$ Finally, there is another statistic that appears in optical turbulence often, called the structure function \(D_u\left(t_1,t_2\right).\) Generally, it involves only one process, and it is defined as $$ \begin{equation} D_u\left(t_1,t_2\right) = \left\langle \left[ u\left(t_1\right) - u\left(t_2\right) \right]^2 \right\rangle. \end{equation} $$ The structure function is related to the correlation according to $$ \begin{equation} D_{u}\left(t_1,t_2\right) = \left\langle u^2\left(t_1\right) \right\rangle + \left\langle u^2\left(t_2\right)\right\rangle - 2\Gamma_{u}\left(t_1,t_2\right). \end{equation} $$
Ergodicity and Stationarity
There are a few important categories of random processes. These are WSS, SSS, ergodic, and stationary increments. Processes with these restrictions apply to many physical phenomena and are useful for their mathematical tractibility.
Ergodicity
Ergodicity is the most restrictive property and has the most useful of mathematical properties. For ergodic processes, ensemble averages of ergodic processes are equal to their time averages. We write the time average as $$ \begin{equation} \overline{f\left[u\left(t\right)\right]} = \lim\limits_{T\to\infty} \frac{1}{T} \int\limits_{-T/2}^{T/2} f\left[u\left(t\right)\right] \, dt. \end{equation} $$ Ergodicity is expressed mathematically as $$ \begin{equation} \left\langle f\left[u\left(t\right)\right] \right\rangle = \overline{f\left[u\left(t\right)\right]}. \end{equation} $$
Strict-Sense Stationary
SSS processes satisfy a slightly less restrictive criterion. In this case, the \(N^{th}\)-order PDF is independent of time origin for all \(N\). Thus, we write this as $$ \begin{equation} p\left(u_1, u_2, \ldots, u_N; t_1, t_2, \ldots, t_N\right) = p\left(u_1, u_2, \ldots, u_N; t_1-T, t_2-T, \ldots, t_N-T\right) \end{equation} $$ for all \(T\). As simple examples, a first-order PDF is independent of time altogether, i.e., \( p\left(u\right) \), and a second-order PDF depends only on a time difference, i.e., \( p\left(u_1,u_2;\tau\right) \), where \( \tau=t_2-t_1 \).
Wide-Sense Stationary
WSS processes are slightly less restrictive than SSS processes and have two key properties. First, the mean is independent of time such that $$ \begin{equation} \mu_u\left(t\right) = \mu_u. \end{equation} $$ Second, the correlation depends only on the time difference \(\tau = t_2 - t_1\) such that $$ \begin{equation} \Gamma_{u}\left(t_1,t_2\right) = \Gamma_{u}\left(\tau\right) \end{equation} $$ and $$ \begin{equation} C_{u}\left(t_1,t_2\right) = C_u\left(\tau\right) = \Gamma_{u}\left(\tau\right) - \mu_u^2. \end{equation} $$ Note that the variance is \(\sigma_u^2 = C_u\left(0\right)\), indicating that is constant. Finally, the structure function has a similar relationship to the correlation and covariance given by $$ \begin{align} D_u\left(\tau\right) &= 2\Gamma_{u}\left(0\right) - 2\Gamma_{u}\left(\tau\right)\\ &= 2C_{u}\left(0\right) - 2C_{u}\left(\tau\right) = 2\sigma_u^2 - 2C_{u}\left(\tau\right). \end{align} $$
Practical Forms of Time Averages for Ergodic Processes
To make the time averages more evident, we write the first- and second-order statistics for an ergodic random process \(u\left(t\right)\) as $$ \begin{align} \mu_u &= \lim\limits_{T\to\infty} \frac{1}{T} \int\limits_{-T/2}^{T/2} u\left(t\right) \, dt, \\ \sigma_u^2 &= \lim\limits_{T\to\infty} \frac{1}{T} \int\limits_{-T/2}^{T/2} \left[ u\left(t\right) - \mu_u\right]^2 \, dt, \\ \Gamma_u\left(\tau\right) &= \lim\limits_{T\to\infty} \frac{1}{T} \int\limits_{-T/2}^{T/2} u\left(t\right) u\left(t-\tau\right) \, dt, \\ C_u\left(\tau\right) &= \lim\limits_{T\to\infty} \frac{1}{T} \int\limits_{-T/2}^{T/2} \left[u\left(t\right) - \mu_u\right] \left[u\left(t-\tau\right) - \mu_u\right] \, dt. \end{align} $$ These are very useful forms. Often, we have one or more discretely sampled time series of \(u\left(t\right)\), and so we use discrete versions of these equations in practice. In Matlab, the mean and var functions can be used to compute the statistical mean and variance of an array of values, whether or not the values constitute a time series. The xcorr and xcov can be used to compute the statistical correlation and covariance of an array of values. For these two functions, the array must represent a time series. Also, xcorr and xcov can compute the cross-correlation and cross-covariance between two different random processes. Although cross-correlation and cross-covariance are not specifically covered in this article, the equations generalize to cross-correlation and cross-covariance in a straighforward way.
Of course, mean, var, xcorr and xcov operate on a finite number of samples. Thus, \(T\) cannot go to infinity, and the continuous integrals with infinite limits are approximated with discrete sums. More samples produces a better approximation. For computing correlation and covariance, each separation \(\tau\) has a different number of samples in the discrete sum. This produces a biased estimate of the correlation and covariance. This is explained in Matlab's documentation for xcorr and xcov. There is an optional argument to specify whether the calculation should produce a biased or unbiased estimate. I have another article on this supplemental information website to illustrate the differnces because it is important to understand the differnces and do the correct calculation.
Stationary Increments
Finally, stationary increments is the least restrictive category discussed here. In this case, we can define a new process from \( u\left(t\right) \) such that $$ \begin{equation} v\left(t_1, t_2\right) = u\left(t_2\right) - u\left(t_1\right). \end{equation} $$ If \( v\left(t_1, t_2\right) \) is SSS, we say that \( u\left(t\right) \) has stationary increments. When we encounter a random process that is not WSS, the structure function can be a useful and more appropriate way of characterizing a process than the correlation.
Random Fields
Random fields are random variables that are a function of the position in space, \( \mathbf{r} \) and sometimes time. In general, the ensemble average values are a function of position and time, and their definitions are very similar to the corresponding definitions for temporal random processes, e.g., $$ \begin{align} \mu_u\left(\mathbf{r},t\right) &= \left\langle u\left(\mathbf{r},t\right) \right\rangle = \int\limits_{-\infty}^{\infty} u\left(\mathbf{r},t\right) \, p\left(u,\mathbf{r},t\right) \, d u,\\ &\left\langle u^2\left(\mathbf{r},t\right) \right\rangle = \int\limits_{-\infty}^{\infty} u^2\left(\mathbf{r},t\right) \, p\left(u,\mathbf{r},t\right) \, d u,\\ \sigma_u^2\left(\mathbf{r},t\right) &= \left\langle \left[u\left(\mathbf{r},t\right) - \left\langle u\left(\mathbf{r},t\right) \right\rangle \right]^2 \right\rangle,\\ \Gamma_{u}\left(\mathbf{r}_1,t_1,\mathbf{r}_2,t_2\right) &= \left\langle u\left(\mathbf{r}_1,t_1\right) u\left(\mathbf{r}_2,t_2\right) \right\rangle \\ &\qquad = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} u\left(\mathbf{r}_1,t_1\right) u\left(\mathbf{r}_2,t_2\right) \, p\left[u\left(\mathbf{r}_1,t_1\right), u\left(\mathbf{r}_2,t_2\right)\right] \, d u\left(\mathbf{r}_1,t_1\right) \, d u\left(\mathbf{r}_2,t_2\right),\\ C_{u}\left(\mathbf{r}_1,t_1,\mathbf{r}_2,t_2\right) &= \Gamma_{u}\left(\mathbf{r}_1,t_1,\mathbf{r}_2,t_2\right) - \left\langle u\left(\mathbf{r}_1,t_1\right) \right\rangle \left\langle u\left(\mathbf{r}_2,t_2\right)\right\rangle \\ D_u\left(\mathbf{r}_1,t_1,\mathbf{r}_2,t_2\right) &= \left\langle \left[ u\left(\mathbf{r}_1,t_1\right) - u\left(\mathbf{r}_2,t_2\right) \right]^2 \right\rangle. \end{align} $$
Homogeneity and Isotropy
The principles of ergodicity for temporal random processes generalize to spatial processes. When a random field is stationary in 2-D or 3-D space, it is called homogeneous. In the case of homogeneous random fields, spatial statistics simplify to $$ \begin{align} \mu_u\left(\mathbf{r},t\right) &= \mu_u,\\ \sigma_u^2\left(\mathbf{r},t\right) &= \sigma_u^2,\\ \Gamma_{u}\left(\mathbf{r}_1,t_1,\mathbf{r}_2,t_2\right) &= \Gamma_{u}\left(\mathbf{r}_1-\mathbf{r}_2,t_2-t_1\right),\\ C_{u}\left(\mathbf{r}_1,t_1,\mathbf{r}_2,t_2\right) &= C_{u}\left(\mathbf{r}_1-\mathbf{r}_2,t_2-t_1\right) \\ &= \Gamma_{u}\left(\mathbf{r}_1-\mathbf{r}_2,t_2-t_1\right) - \mu_u^2\\ D_u\left(\mathbf{r}_1,t_1,\mathbf{r}_2,t_2\right) &= D_u\left(\mathbf{r}_1-\mathbf{r}_2,t_2-t_1\right)\\ &= 2\Gamma_{u}\left(\mathbf{0},0\right) - 2\Gamma_{u}\left(\mathbf{r}_1-\mathbf{r}_2,t_2-t_1\right)\\ &= 2C_{u}\left(\mathbf{0},0\right) - 2C_{u}\left(\mathbf{r}_1-\mathbf{r}_2,t_2-t_1\right)\\ &= 2\sigma_u^2 - 2C_{u}\left(\mathbf{r}_1-\mathbf{r}_2,t_2-t_1\right). \end{align} $$ Furthermore with spatially ergodic random fields, ensemble averages are equal to their spatial averages. A spatial average can be written as $$ \begin{equation} \overline{f\left[u\left(\mathbf{r}\right)\right]} = \lim\limits_{A\to\infty} \frac{1}{A} \int\limits_{\mathcal{A}} f\left[u\left(\mathbf{r}\right)\right] \, d\mathbf{r}, \end{equation} $$ where \(\mathcal{A}\) is the spatial domain of the integration, an \(A\) is its area. Spatial ergodicity is expressed mathematically as $$ \begin{equation} \left\langle f\left[u\left(\mathbf{r}\right)\right] \right\rangle = \overline{f\left[u\left(\mathbf{r}\right)\right]}. \end{equation} $$
Because random fields depend on a vector argument, there is one more limiting case, namely isotropy. An isotropic random field's moments do not depend on direction, just the magnitude. In this case, second-order moments that depend on separation simplify to these forms: $$ \begin{align} \Gamma_{u}\left(\mathbf{r}_1,t_1,\mathbf{r}_2,t_2\right) &= \Gamma_{u}\left(\left\vert\mathbf{r}_1-\mathbf{r}_2\right\vert,t_2-t_1\right),\\ C_{u}\left(\mathbf{r}_1,t_1,\mathbf{r}_2,t_2\right) &= C_{u}\left(\left\vert\mathbf{r}_1-\mathbf{r}_2\right\vert,t_2-t_1\right) \\ D_u\left(\mathbf{r}_1,t_1,\mathbf{r}_2,t_2\right) &= D_u\left(\left\vert\mathbf{r}_1-\mathbf{r}_2\right\vert,t_2-t_1\right) \end{align} $$ Notice that the relationships between the correlation, covariance, and structure function extend to homogeneous random fields. Often, optical turbulence (refractive index, phase, etc.) is treated as a homogeneous and isotropic random process. This simplifies calculations, although the accuracy of these approximations depends on the time of day, location in the atmosphere, etc.
Just like with temporal random processes, we can write useful forms of first- and second-order ergodic fields by being explicit about the spatial averages. These are given by $$ \begin{align} \mu_u &= \lim\limits_{A\to\infty} \frac{1}{A} \int\limits_{\mathcal{A}} u\left(\mathbf{r}\right) \, d\mathbf{r}, \\ \sigma_u^2 &= \lim\limits_{A\to\infty} \frac{1}{A} \int\limits_{\mathcal{A}} \left[ u\left(\mathbf{r}\right) - \mu_u\right]^2 \, d\mathbf{r}, \\ \Gamma_u\left(\Delta r\right) &= \lim\limits_{A\to\infty} \frac{1}{A} \int\limits_{\mathcal{A}} u\left(\mathbf{r}_1\right) u\left(\mathbf{r}_1-\Delta\mathbf{r}\right) \, d\mathbf{r}_1, \\ C_u\left(\Delta\mathbf{r}\right) &= \lim\limits_{A\to\infty} \frac{1}{A} \int\limits_{\mathcal{A}} \left[u\left(\mathbf{r}_1\right) - \mu_u\right] \left[u\left(\mathbf{r}_1-\Delta\mathbf{r}\right) - \mu_u\right] \, d\mathbf{r}_1, \end{align} $$
Again, these are generally the most useful forms for computing these statistics from discretely sampled data. For 2-D random fields, the mean and var functions can still be used. If we have a 2-D array u, we can comput its sample mean using mean(mean(u)) or mean(u(:)). To compute the sample variance, the simplest way is var(u(:)).
For computing correlation and covariance of 2-D fields, there is no Matlab function like xcorr or xcov. That is why Section 3.2 in my book explains how this calculation works in 2-D. Again, there is a different number of samples for each spatial separation \(\Delta r\). This is more complicated in 2-D than 1-D because in 2-D, and in optics, we often work with data that has a "mask" associated with it. This is due to the combination of circular lenses and mirrors with rectangular sensors. The corr2_ft function from my book is intended to handle this, but the code is a bit flawed. In another article on this website, I provide improved code with a full explanation.
References
- Joseph W. Goodman, Statistical Optics Ch. 3, Wiley, New York, NY (1985)
- Larry C. Andrews and Ronald L. Phillips, Laser Beam Propagation through Random Media, 2nd Ed., SPIE Press, Bellingham, WA (2005)