## Abstract

Optical speckle is commonly observed in measurements using coherent radiation. While lacking experimental validation, previous work has often assumed that speckle’s random spatial pattern follows a Markov process. Here, we present a derivation and experimental confirmation of conditions under which this assumption holds true. We demonstrate that a detected speckle field can be designed to obey the first-order Markov property by using a Cauchy attenuation mask to modulate scattered light. Creating Markov speckle enables the development of more accurate and efficient image post-processing algorithms, with applications including improved de-noising, segmentation and super-resolution. To show its versatility, we use the Cauchy mask to maximize the entropy of a detected speckle field with fixed average speckle size, allowing cryptographic applications to extract a maximum number of useful random bits from speckle images.

© 2012 OSA

## 1. Introduction

The simplest probabilistic description of an optical speckle field assumes it as an uncorrelated random process across space. This description is valid when the field is sampled at a rate less than or equal to the average speckle size, and is useful in many scenarios including denoising [1], entropy analysis [2] and laser speckle imaging [3], to name a few. Often, however, this sampling condition is not satisfied. Setups that benefit from detecting speckle enlarged across multiple sensor pixels include those for optical encryption [4], optical phase conjugation [5], speckle shape analysis [6], and structured illumination [7], among others. The correlations that arise between neighboring pixels are rarely modeled exactly, complicating attempts to calculate a detected field’s entropy, determine specifics about a scattering source or remove unwanted speckle noise, for example. An accurate Markov model representation of these inter-pixel dependencies can increase the accuracy and efficiency of such computational procedures.

Plainly put, a Markov process is a random sequence of states, such that the probability of observing a particular state only depends on the properties of directly neighboring states. In the context of speckle, the states we will be concerned with are the particular value a pixel takes on when an optical field is detected. Neighboring states are the values detected by adjacent pixels. Considering a speckle field in 1D, the Markov condition requires the probability of detecting a certain value at one pixel depends only on the value detected at its two neighboring pixels, and no others (Fig. 1(a)).

There has been some previous interest in attempting to model the speckle “noise” in synthetic aperture radar (SAR) images as a Markov process [1, 8 – 10]. Doing so enables both removal of speckle noise and SAR image segmentation. While the above references offer some mathematical support for assuming Markov speckle (most notably [1]), their derivations are approximate for two reasons. First, this prior work is only concerned with processing images post-capture, and does not physically modify the imaging system to produce exact Markov speckle. Second, this work aims to fit a Markov model to the detected speckle’s intensity, while we find an exact solution only for the speckle’s complex field. Here we show how a simple modification to any speckle detection setup, in the form of an apodizing mask, leads to a detected speckle field that obeys a Markov process. Moreover, the Markov property holds as the average speckle size is varied to cover an arbitrarily large number of pixels.

We first review the Gaussian distribution of a 1D speckle field, explain the Markov property related to Gaussian distributions, and provide a sufficient condition for a detected speckle field to be Markov in Section 2. In Section 3, we show one approach to realize Markov speckle by adding a designed apodizing mask. Section 4 extends this model to 2D, and Section 5 offers a discussion of practical limitations. The accuracy of using an apodizing mask to create Markov speckle is then experimentally investigated in Section 6.

Finally, Section 7 explores one possible application of Markov speckle: maximizing the number of extractable random bits from a digitally detected speckle field. Optically probing a volumetric scatterer produces sets of speckle patterns that can be turned into unclonable encryption keys [2, 4]. Generating and detecting large speckle spanning several pixels ensures these random keys are reproducible. In Section 7 we demonstrate that Markov speckle of arbitrary size exhibits an entropy-maximizing property. We then argue that such maximum-entropy Markov speckle may lead to larger random keys for a fixed average speckle size, potentially increasing the efficiency of current optical encryption setups.

## 2. Mathematical background

The following analysis is based on derivations in [11]. A coherent, monochromatic, polarized field with many de-phased contributions leaving a scattering region is assumed as the initial speckle field. Furthermore, attention will be restricted to a discretized representation of the field at a set of pixels of finite size *δ*, assuming that the average speckle size is equal to or greater than *δ*. The effects of discretization, the case of speckle size being smaller than *δ* and their connection to a Markov process are discussed in Appendix A. Finally, it is assumed that the joint probability distribution of the speckle field does not change across the detector area of interest, leading to a homogeneous Markov model and stationary statistics.

#### 2.1. Speckle field covariance and the attenuation mask

A complex speckle field *A*(*x,y,z*) measured at a discrete pixel location (*x,y*) on a distant plane *z* perpendicular to the direction of propagation is the random sum of many independent phasor components. Evolution of the field over space is a random walk on the complex plane, as in Fig. 1(b). We first consider a 1D detector along *x* at fixed *z*. The field *A* at one pixel is a circular symmetric complex Gaussian random variable with probability density function (pdf)

*ϕ*| is a random phasor [11]. Since we are concerned with correlations between multiple pixels, Eq. (1) can be transformed to a complex random vector

**of the speckle field at a set of**

*A**N*neighboring pixels along

*x*. As each element of

**remains circular symmetric in the presence of correlations, the pdf of the circular symmetric complex Gaussian speckle vector**

*A***is zero-mean Gaussian:**

*A***is the covariance matrix of the complex field**

*J***representing the correlations between all**

*A**N*pixels. By definition,

**is a positive definite symmetric matrix.**

*J*

*J*^{−1}is often referred to as a precision matrix.

Within the context of optics, ** J** is also referred to as the complex mutual intensity function, or degree of coherence function, of the speckle field at plane

*z*. This coherence function is expressed as

**(**

*J**x*

_{1},

*x*

_{2}) = 〈

*A*(

*x*

_{1})

*A**(

*x*

_{2})〉, where angular brackets denote an ensemble average of detected fields scattered from

*z*= 0 (Fig. 2). As shown in [12], the coherence function

**(**

*J**x*

_{1},

*x*

_{2}) is a weighted integral of a function

*M*(

*η*) that describes the illuminated area’s shape at the scattering surface:

*λ*and wavenumber

*k*. Equation (3) neglects a constant phase factor and assumes an unresolvable microstructure surface from plane

*z*. The illuminated area shape at the scatterer is determined by an amplitude-attenuating mask

*M*(

*η*) placed at the scattering material’s surface. A nearly identical relationship is found if an imaging setup is used instead of free-space propagation and the mask

*M*(

*η*) is inserted at the aperture plane (Fig. 2(b)). Substituting Δ

*x*=

*x*

_{1}−

*x*

_{2}into Eq. (3) leads to a scaled Fourier transform relationship between the mask

*M*(

*η*) and the speckle field autocorrelation at plane

*z*:

*l*as the width of the autocorrelation

_{c}**(Δ**

*J**x*)’s main lobe for an open (unmodified) aperture of width

*w*. Notice that speckle’s multivariate distribution Eq. (2) is fully characterized by the autocorrelation function

**(Δ**

*J**x*). Also, from Eq. (4) it is clear the mask function

*M*(

*η*) uniquely determines

**(Δ**

*J**x*). Thus, our goal is to first establish a sufficient condition on

**(Δ**

*J**x*) that guarantees Markov speckle, and then to manipulate

*M*(

*η*) to ensure

**(Δ**

*J**x*) satisfies this condition. We derive the correct mask function

*M*(

*η*) for Markov speckle in Section 3, and experimentally place this mask at a scatterer’s surface to test its performance in Section 6.

#### 2.2. The Markov property

A sequence of random variables obeying the homogeneous Markov property is described by a fixed transition probability that dictates transitions between immediately neighboring events. A good introduction to Markov random processes can be found in [13]. We continue to limit our attention to a 1D optical field *A* across an adjacent set of *N* pixels ** A** = (

*A*

_{1},...,

*A*). The 2D case is examined in Section 4. The discretized complex values each pixel detects define the finite state space

_{N}*S*of the random process. A detected speckle field

**with values in**

*A**S*is a first-order unilateral Markov process, or a Markov chain, if the conditional probability satisfies

*n*≥ 1. We are primarily interested in first-order conditioning since it offers the most compact description of a correlated random process and lends useful properties to entropy maximization. Unless otherwise stated, future references to Markovity will imply first-order conditioning. Since speckle is a spatial process, the first-order bilateral Markov property for values on

*either side*of each pixel must be satisfied:

*A*|

_{n}*A*

_{n}_{−1},

*A*

_{n}_{+1}) = (

*A*

_{2}|

*A*

_{1},

*A*

_{3}) for all

*n*. For completeness, assumed Markov processes are also aperiodic, irreducible and reversible. Finally, a transition matrix

**is useful when representing the evolution of a homogeneous unilateral Markov process.**

*P***tabulates the conditional probabilities of transitioning from any state**

*P**α*to any state

*β*at an immediately neighboring pixel:

**(**

*P**α*,

*β*) = (

*A*=

_{n}*β*|

*A*

_{n}_{−1}=

*α*). An important equation that the matrix

**must satisfy is the Chapman-Kolmogorov equation [13], which we use in Section 6: Here,**

*P**v*is an intermediate pixel between pixels

*m*and

*n*, and

*P**and*

_{mv}

*P**are the transition matrices from pixel*

_{vn}*m*to

*v*and pixel

*v*to

*n*, respectively.

**Three characterizations of Gaussian Markov processes:** A sequence of random variables that satisfies Eq. (6) is called a strict-sense Markov (SSM) process. A weaker sense of Markovity, termed wide-sense Markov (WSM), relies upon a neighbor-dependent conditioning defined through the conditional mean. A random process is first-order bilateral WSM if and only if the linear minimum mean squared error (MSE) estimate satisfies

*p*,

*q*> 1, where 𝔼̂ is the minimum MSE estimator of the random process, given the conditioning. We can equivalently characterize a first-order bilateral WSM process in terms of an autoregressive representation: where the {

*ρ*} are coefficients and {

_{k}*U*} are independent Gaussian random variables such that

_{n}*U*and

_{n}*A*are independent for

_{m}*n*≠

*m*.

Equation (2) shows that a speckle field is a multivariate Gaussian random process. As shown in detail in [14], a Gaussian process is SSM if and only if it is WSM. Thus, Eq. (6), Eq. (8) and Eq. (9) offer equivalent definitions of a Gaussian speckle field’s Markovity. We will use Eq. (8) to assist us in experimental confirmation of Markov speckle, while Eq. (9) will lead us to a definition of maximum speckle entropy in Section 7.

**Conditional probability of Gaussian processes:** A multivariate Gaussian process has the unique property that its second-order statistics, tabulated by the covariance matrix ** J**, fully describe its dependence relationships. Specifically, if
${J}_{nm}^{-1}={\mathit{J}}^{-1}\left(n,m\right)=0$ then pixel

*n*and

*m*are independent, conditioned on all other pixels [15]. This special property is exhibited by the conditional pdf of speckle field value

*A*at pixel

_{n}*n*, given field values

*A*at all other pixels:

_{m}*A*must only depend on

_{n}*A*

_{n}_{−1}and

*A*

_{n}_{+1}. Thus, the summands corresponding to

*m*with |

*m*−

*n*| > 1 on the right hand side of Eq. (10) should disappear. This is achieved when all entries ${J}_{nm}^{-1}$ are zero except those between the matrix

*J*^{−1}’s super-diagonal and sub-diagonal. Strategies to ensure that

*J*^{−1}is tridiagonal are examined next.

## 3. Markov speckle in 1D

This section discusses two scenarios under which a detected speckle field obeys the Markov property: a limiting case where the average speckle size does not exceed one pixel, and a designed case that is independent of speckle size. As noted above, Markov speckle must have a tridiagonal precision matrix *J*^{−1}. While small speckle satisfies this condition trivially, for large speckle we modify *J*^{−1} by shaping the optical field with a designed attenuation mask *M*(*η*).

**A. Limit of small speckle:** Speckle with an average size *l _{c}* less than the detector pixel size

*δ*obeys the Markov property (Appendix A). The value detected by each pixel will be uncorrelated, leading to a diagonal covariance matrix

**and precision matrix**

*J*

*J*^{−1}:

**I**is the identity matrix. A diagonal

*J*^{−1}indicates no conditioned variables appear on the left side of Eq. (10). A discrete speckle field

**with diagonal**

*A***is a purely random process.**

*J***B. Speckle with designed correlation:** A speckle field with average size *l _{c}* larger than one pixel width

*δ*can obey the first-order Markov property only when its precision matrix

*J*^{−1}is tridiagonal. A tridiagonal

*J*^{−1}transforms Eq. (10) to

*A*

_{1},...,

*A*) also satisfies the unilateral Markov property Eq. (5). Thus, any stationary speckle field with a tridiagonal precision matrix

_{N}

*J*^{−1}is a first order Markov process of either type. We note the precision matrix of

*g*order Markov processes must have non-zero entries only between the positive and negative

^{th}*g*diagonals. A covariance matrix

^{th}

*J**of the following exponential form generates a tridiagonal precision matrix:*

_{e}*ρ*are constants. For shift-invariant speckle, the exponential covariance matrix in Eq. (13) can be expressed as a one-dimensional autocorrelation function:

*γ*= − ln(

_{o}*ρ*). As expressed in Eq. (4),

*J**(Δ*

_{e}*x*) is fully characterized by the attenuation function

*M*(

*η*) of an apodizing mask. Substituting Eq. (15) into Eq. (4) and taking the inverse Fourier transform, we obtain the following sufficient condition on

*M*(

*η*) to guarantee a tridiagonal precision matrix ${\mathit{J}}_{e}^{-1}$:

*γ*= −ln(

*ρ*)/

*λz*. Apart from the constant pre-factor, Eq. (16) describes a Cauchy distribution with respect to

*η*. Placing an amplitude apodizing mask

*M*(

*η*) following Eq. (16) at the scatterer surface in a free-space speckle detection setup leads to an exponential autocorrelation

*J**(Δ*

_{e}*x*) of the speckle field a sufficient distance

*z*away, which obeys the Markov property. The same holds for speckle in the imaging setup in Fig. 2(b) with a Cauchy mask placed at the aperture plane. This amplitude-only mask simply shapes the speckle pattern in such a way that its autocorrelation function across many pixels can be recursively described using one parameter (

*ρ*), leading to first-order Markov conditional probability relationships.

## 4. Markov speckle in 2D

The above analysis also extends to form 2D speckle patterns into Markov processes. To explain how, we introduce the concept of a Markov Random Field (MRF). We define an MRF over a rectangular lattice *L*, which here is the 2D pixel array that detects the speckle field (indexed by (*i*, *j*)). Also, a neighborhood set Ψ associated with *L* is defined (formally) as Ψ = {Ψ* _{ij}* ⊂

*L*: (

*i*,

*j*) ∈

*L*}, where Ψ

*is a set of neighbors of (*

_{ij}*i*,

*j*). Informally, the neighborhood of pixel (

*i*,

*j*) is a set Ψ

*of nearby pixels that does not include (*

_{ij}*i*,

*j*). For example, the first-order neighborhood of any pixel (

*i*,

*j*) contains the 4 pixels it shares an immediate border with (Fig. (3)). We can extend our 1D Markov process definition to describe a 2D MRF by limiting the conditional dependence of pixel (

*i*,

*j*) to its neighborhood:

*but not (*

_{i,j}*i*,

*j*). Equation (17) describes a special type of SSM field. As with the 1D case, we can also characterize a WSM process in 2D in terms of a bilateral autoregressive form,

*ρ*} is a set of correlation coefficients and {

_{k,l}*U*} are independent random processes such that

_{ij}*U*and

_{i,j}*A*are independent for (

_{k,l}*i*,

*j*) ≠ (

*k,l*).

Since speckle is Gaussian, a 2D Markov speckle pattern on *L* that satisfies Eq. (17) belongs to a subclass of MRF, known as a Gaussian Markov Random Field (GMRF). As with the SSM and WSM equivalence in 1D, Eq. (17) and Eq. (18) also offer equivalent characterizations of 2D speckle as a GMRF [15]. For a second-order 2D Markov process, the right-hand side of Eq. (18) sums over the 8 immediately surrounding neighbors of pixel (*i*, *j*). Assuming generation from a spatially symmetric *x*–*y* separable correlation function, the sum simplifies to

*γ*= −ln(

*ρ*)/

*λz*. From Eq. (19), a 2D speckle Markov process will depend on its 4 immediate neighbors (first-order neighborhood) and to a lesser extent its 4 diagonal neighbors (second-order neighborhood) assuming

*ρ*< 1. The separable form Eq. (20) of the apodizing mask is the most direct method of achieving Markovity for 2D speckle. Designing non-separable and higher-order Markov processes is possible following further investigation into GMRF theory [14, 15].

## 5. Practical considerations for Markov speckle

Creating a speckle field that exactly follows the Markov property is limited in practice by several experimental conditions beyond the introduction of noise. First, the derivation of Eq. (16) and Eq. (20) assume an optical geometry that extends to infinity in both directions *η* and *ξ*. A limited mask baseline *w* (i.e., a finite aperture width) cuts off the tail of the Cauchy function equation, generating an exponential correlation function that deviates from an ideal curve (Fig. 4). Although deviations are small, a modified Cauchy function may be solved for via an optimization procedure to better approximate this desired exponential autocorrelation.

Second, spatial discretization effects prevent realization of an exact Markov relationship. Effects caused by a discrete pixel size are discussed in detail in Appendix A. For large speckle, these effects are minimized with an increase in average speckle size for a fixed pixel size. On the other hand, small speckle (less than one pixel) obeys the Markov property regardless of the shape of its correlation function. Digitization of the optical field into discrete *values* does not fundamentally limit the creation of an exact Markov sequence. Transition probabilities can be determined for a discrete state space of any size, often set by the bit depth of the sensor.

Third, we note that since the speckle field in the above derivations is complex, it has a complex Markov state space. This does not present any fundamental limitations in establishing Markovity. A transition matrix can be created by labeling each complex state with a particular value to jointly describe the field amplitude and phase, or the state space may be defined as the real and/or imaginary field value. Note that the mask in Eq. (20) will not cause speckle *intensity* to behave as Markov. The observation that speckle’s complex field will generally follow Markov assumptions more closely than its intensity was first suggested in [16]. We explore in detail how the intensity connects to an unobservable Markov process through a squaring operation in Appendix C.

## 6. Experimental verification

This section first introduces an intuitive test that measures the degree to which a large dataset follows the Markov property. Then, both simulated and experimental data demonstrate how Cauchy-masked speckle performs better on this test than speckle fields attenuated by other mask functions.

#### 6.1. Chapman-Kolmogorov validation equation

In general, it is difficult to offer an exhaustive proof that a large sequence of data obeys the Markov property. A full test of bilateral Markovity for a 1D data sequence must consider the validity of,

*m*and

*n*. Considering all possible lags

*m*is computationally infeasible for large random sequences. Instead of an exhaustive proof, a test previously explored with large datasets [17 – 19] uses the Chapman-Kolmogorov equation to check for properties consistent with a Markov process. Referring to Eq. (7), we will test if the equality holds. Equation (22) is a necessary condition for a homogeneous Markov process. This test’s main benefit is it only requires computing and storing first-order conditional relationships, reaching statistical significance with less data than other second-order tests [20]. Second, it offers the ability to visualize performance errors in the three transition matrices.

A detailed discussion of the validity of using the Chapman-Komolgorov (CK) test to verify Markov speckle is presented in Appendix B. The first assumption required to derive the CK test Eq. (22) from Eq. (21) is a monotonically decreasing correlation function ** J**(Δ

*x*), valid for all mask functions we test (although counter-examples can be constructed [20]). Second, we assume a spatially stationary field to show speckle obeys both the unilateral (Eq. (5)) and bilateral (Eq. (6)) Markov property.

The three transition matrices that comprise the CK test must be estimated from recorded speckle field data. The maximum likelihood (ML) estimator for a stationary transition matrix ** P**̂(

*α*,

*β*) from a set of statistics is

*P*^{̂}(

*α*,

*β*) =

*T*(

*α*,

*β*)/∑

*′*

_{β}*T*(

*α*,

*β*′), where

*T*(

*α*,

*β*) represents the number of observed transitions from state

*α*to state

*β*(i.e., pixel value

*α*to pixel value

*β*) [21]. This linear estimate is simply an average transition rate between different pixel values. It is direct to show the CK test equation based on the ML estimator is equivalent to our WSM definition in Eq. (8) considering three states.

ML estimates for each of the three transition matrices in Eq. (22) can be constructed by sweeping through detected speckle and counting the number of transitions *T*(*α*, *β*) from speckle field value *α* to *β*, either at adjacent pixels (*P*_{n,n}_{−1}, *P*_{n}_{−1}_{,n}_{−2}) or at alternate pixels (*P*_{n,n}_{−2}). The ergodic theorem guarantees this sweeping process is equivalent to constructing a conditional mean estimator over many independent speckle field realizations. A robust expectation measurement is created by repeating the sweep process over a large set of independent images following identical statistics (Fig. 5).

Given a set of transition matrix estimates, we define the following error metric *r* based on total variation (TV error) to measure how well the data sequence satisfies the CK test Eq. (22):

*S*| represents the size of the state space. When

*r*is zero, the sequence obeys the Markov property exactly.

#### 6.2. Experimental setup and procedure

The experimental setup used to verify a masked speckle field follows a first-order Markov process via Eq. (22) is diagrammed in Fig. 5(a). A solid state 532nm CW laser is split into an object and reference arm (both spatially filtered and collimated). The object arm is first incident onto an amplitude-modulating spatial light modulator (SLM) displaying a random binary pattern *b _{t}* (3.3cm LCD, 1920x1080 pixels). After removal of higher diffraction orders with a filtered 4

*f*setup, the randomized object field is imaged onto the back surface of a volumetric scattering material (25mm

^{2}opal diffusing glass). The field on the opposite side of the scatterer (front surface) is assumed to be delta-correlated, serving as the speckle field source.

A patterned grayscale apodizing mask *M*(*η*, *ξ*) is positioned directly adjacent to the scatterer front surface (Kodak LVT-exposed on film at 2032dpi, Bowhaus Printing). After passing through the mask, the speckle field propagates a distance *z* to a CMOS detector (1936 x 1456 pixels, 4.54*μ*m width). An electro-optic phase modulator (EOM) is used to phase-shift the reference plane wave by *π*/2 four times before recombination in a phase-shifting digital holography setup. An estimate of the speckle field phase is generated from four phase-shifted intensity images via the phase recovery equation [22]. The amplitude of the object wave is solved for with a similar equation (pixels where a division by 0 occurs are ignored).

Displaying a different binary pattern *b _{t}*

_{+1}on the amplitude SLM leads to the measurement of an independent speckle field

*A*

_{t}_{+1}. Displaying 100 different random amplitude SLM images builds a 1936 x 1456 x 100 pixel dataset of complex field measurements, which are unwrapped into a vector of approximately 10

^{8}elements after windowing out nonuniform image areas. The scanning procedure discussed in the previous subsection is then applied to this vector (ignoring transitions at image edges and between images) to generate the three transition matrix estimates

**̂ and**

*P**r*in Eq. (23).

#### 6.3. Results

Here we demonstrate that Cauchy-masked speckle follows Markov statistics more closely than speckle that passes through other mask functions. Since the Markov TV error *r* depends both on speckle size and correlation function shape, it is measured as a function of the speckle’s normalized correlation area, *a* = |∑ ** J**(Δ

*x*)/

**(0)|. The correlation area**

*J**a*and size of main lobe

*l*are proportional to propagation distance

_{c}*z*for a given mask (see Eq. (4)). Likewise,

*a*and

*l*are both inversely proportional to mask parameters

_{c}*ρ*and

*w*.

We use two different Markov state space definitions to populate a transition matrix ** P** with complex field measurements. First, the complex field’s amplitude and phase is concatenated into a single state for a complete description of the Markov process (complex state space). Second, only the real value of the field is used to create a state space half as large (real state space). Transition matrices for the real state space are easier to visualize but do not offer a complete description of the Markov process.

Figure 6 displays complex and real transition matrices used by the CK test Eq. (22), including (from left to right) *P*_{n,n}_{−1}, *P*_{n,n}_{−1}*P*_{n}_{−1}_{,n}_{−2}, *P*_{n,n}_{−2} and (*P*_{n,n}_{−2} − *P*_{n,n}_{−1}*P*_{n}_{−1}_{,n}_{−2}). The right-most matrices offer a visualization of error under the Markov assumption. Transition matrices in Fig. 6(b) are created experimentally using a square open aperture of width *w* = 0.8cm with the detector at *z* = 22cm, making *a* = 5.54. The ** P** matrices contain transition data for the complex state space after the field is discretized into 16 amplitude (A) bins and 8 phase (

*ϕ*) bins (|

*S*| = 128). The

**matrices display the same data for the real state space, discretized into 64 states. The matrices in Fig. 6(a) are created via simulation of speckle through a square mask with similar parameters (details below). The speckle correlation functions obtained through experiment and simulation at this distance are in Fig. 6(c). TV error for the square-masked speckle is**

*R**r*= 0.090 in experiment and

*r*= 0.088 in simulation.

Figure 7 contains an identical set of plots except with the square mask replaced by a Cauchy apodizing mask following Eq. (20) (*w* = 1cm, *γ* = *w*/8). The detector distance was slightly varied to *z* = 20cm to achieve roughly the same correlation area as with the square mask (*a* = 5.56). The difference matrices on the right, proportional to *r*, are plotted on the same scale as those for the square mask in Fig. 6 and demonstrate reduced variation. TV error for the Cauchy-masked speckle is *r* = 0.061 in experiment and *r* = 0.058 in simulation.

The simulation model used above is based on the transmission matrix formalism of scattering [23]. Scattering is modeled by multiplication of a complex random Gaussian matrix with many independent random incident field vectors. The effects of the apodizing mask at the scatterer surface, including the windowed aperture effects discussed in Section 5, are added by convolution along the scattering matrix columns. The correlation area *a* is set by the area of the convolution kernel. Simulations were run over approximately 10^{8} random variables.

To demonstrate that Cauchy-masked speckle of arbitrary average size remains a Markov process, we measure the error metric *r* for speckle obtained with 5 different correlation areas *a*. This is experimentally achieved by varying *z* to 5 distances per mask. The results of this experiment are shown in Fig. 8. A *w* = 1cm Gaussian apodizing mask is also compared to the square mask and Cauchy mask described above (all with similar total transmission) to demonstrate that Markovity is highly dependent upon mask shape. In both simulation and experiment, *r* slowly decreases with an increase in *a*, which is a result of the transition matrices becoming increasingly diagonal with larger speckle. However, the derived Cauchy mask creates speckle that more closely follows Markov statistics, independent of speckle size.

#### 6.4. Discussion

In both simulation and experiment, Cauchy-masked speckle has a significantly lower TV error *r* than un-masked or Gaussian-masked speckle. This demonstrates that speckle can be optically designed to follow the Markov property with increased accuracy. This trend is independent of average speckle size. However, we observe that unmodified speckle and speckle apodized by radially decreasing masks are also roughly Markov. Dominant sources of error include the finite extent of the apodizing mask and pixel discretization effects (included in simulation). Differences between simulation and experiment can be mostly attributed to the unstable nature of the phase-shifting digital holography setup, especially when measuring speckle with a small correlation area *a*. Further issues include a possible global phase variation across the sensor, loss of dynamic range from using a reference beam, imperfect absorption by mask elements, and inaccuracies in modeling the volumetric scatterer under the transmission matrix formalism.

## 7. Application: Entropy maximization

Markov speckle can be immediately applied to improve speckle removal algorithms (where it is a common assumption [1, 8, 10]), or offer an additional constraint to enhance speckle-based super-resolution reconstruction [7], for example. In this section, we focus on the application of Markov speckle to random bit generation. Recent work demonstrates that one can use speckle to create highly random yet reproducible sequences of bits [2, 4]. These random speckle keys are used in cryptographic applications including secure identification, authentication and communication key establishment. Detecting speckle with an average size greater than several pixels is required by these setups to ensure the field pattern is experimentally reproducible in the presence of noise. The number of useful random bits that can be extracted from a detected speckle field is commonly assumed to be an increasing function of its entropy, as discussed in [2]. We demonstrate that for a fixed average speckle size *l _{c}*, the entropy of a detected speckle field is maximized when a Cauchy mask is used to create Markov speckle. This, in turn, can allow speckle encryption setups to increase their efficiency of random bit generation.

The entropy of a real joint Gaussian distribution of *N* variables in 1D with covariance matrix ** J** is given by,
$h\left(\mathcal{N}\left(0,\mathit{J}\right)\right)=\frac{1}{2}\text{log}\left[{\left(2\pi \text{e}\right)}^{N}\cdot \text{det}\left(\mathit{J}\right)\right]$, where det denotes determinant [24]. Similar to the real Gaussian case, the covariance matrix of a circular symmetric

*complex*Gaussian process also determines its entropy. For a complex speckle process

**=**

*A***+ j**

*X***, the concatenated vector [**

*Y***] follows a multivariate Gaussian distribution with covariance matrix**

*X;Y***is real. A real**

*J***is created by a symmetric apodizing mask function**

*J**M*(

*η*), already assumed in Section 4. Since this modified covariance matrix indicates

**and**

*X***are uncorrelated, the entropy of a correlated speckle field**

*Y***over**

*A**N*pixels is

Burg’s maximum entropy theorem [24] states that the maximum entropy rate stochastic process satisfying the autocorrelation constraints ** J**(Δ

*x*) =

*ρ*

_{Δ}

*for Δ*

_{x}*x*= 0, 1...

*g*is the

*g*order Gauss-Markov process in the form of Eq. (9), where the

^{th}*U*are i.i.d. $~\mathcal{N}\left(0,{\sigma}_{0}^{2}\right)$. Since a circular symmetric complex Gaussian process is separable into two uncorrelated real Gaussian processes given a real covariance matrix, this theorem directly applies to speckle fields generated from a symmetric mask function

_{n}*M*(

*η*). To select the appropriate

*ρ*

_{Δ}

*to satisfy this theorem’s conditions, we note that speckle is often approximately assigned the*

_{x}*single*parameter of “average size” in many applications (e.g., [2,8]), indicating just

*one*parameter

*ρ*

_{1}can be fixed. Many setups to this point generate speckle through an unmodified aperture, which in 1D corresponds to a rect function. Thus, “average size” is often assumed to imply the main lobe width of the resulting speckle’s sinc autocorrelation function

*l*, related to the aperture width

_{c}*w*by

*l*=

_{c}*λz*/

*w*(Fig. 9). Assuming this relationship, a single

*ρ*

_{1}parameter tied both to average speckle size and aperture width is determined as,

*w*, we conclude that first-order Markov speckle designed by inserting the Cauchy apodizing mask in Eq. (16) using

*ρ*from Eq. (26) will maximize the speckle’s entropy. Parameters for other correlation functions in Table 1 are determined through a similar requirement that the first neighboring pixel’s correlation equals

*ρ*

_{1}, and all exhibit lower entropy. Since the CK test in Section 6 suggests but does not prove Markovity, we may interpret the trends in Table 1 as additional evidence that the Cauchy mask is able to produce optimally Markov speckle. Using the similarity of Eq. (9) and Eq. (18), the above analysis also extends to support a maximum entropy theorem in 2D.

## 8. Conclusions and future work

The addition of an amplitude-modulating mask following a Cauchy distribution at a scattering surface creates speckle satisfying first-order Markov conditions. An experimental test offers support to this claim. The entropy of a speckle field with a fixed average speckle size is maximized when this mask is used, leading to more efficient speckle-based random bit generation. Other applications that may benefit from Markov-designed speckle include those in which speckle is viewed as a source of noise (e.g., SAR imagery, ultrasound, and digital holography). Optically modifying the speckle to follow Markov statistics at the detector will assist in its digital removal, as discussed in [1, 9, 10]. Furthermore, general modification of speckle’s autocorrelation function to a desired curve, as with modifications to the point-spread function of a camera, may enable improved depth estimation or superresolution by image post-processing. Finally, several fundamental properties of Markov speckle, including masks that generate higher-order Markov processes and methods of modeling speckle intensity as a hidden Markov process, have yet to be fully explored and may lead to interesting insights.

## Appendix A: Pixel sampling effects on speckle’s autocorrelation

The detection of a continuous function *A*(*x*) by a discrete pixel lattice causes the second-order statistics of *A*(*x*) to deviate slightly. We represent the detection process as a convolution operation with a pixel transfer function Π(*x*), assumed to be a rect function of width *δ*:

*A*(

*x*) and

*A*

_{0}(

*x*) represent continuous and detected optical intensity, respectively. Since our experiment uses four discrete intensity measurements to compute a complex field, we will assume this convolution relationship also holds for a computed complex

*A*and

*A*

_{0}. Pixel modulation in the spatial frequency domain is

*ν*. A discrete field is represented as a sampled version of

_{x}*A*

_{0}and

*Â*

_{0}following Shannon’s sampling theorem. The mean of the detected complex field is unchanged by the above sampling process (Eq. (28) indicates 〈

*A*

_{0}(

*x*′)〉 = 0 given 〈

*A*(

*x*)〉 = 0). Expressing the Fourier transform ℱ of the sampled autocorrelation function

*J*_{0}(Δ

*x*) in terms of

*A*shows effects of pixelization:

**̂ is the Fourier transform of the un-sampled autocorrelation**

*J***and ★ denotes convolution. The squared**

*J**sinc*function represents pixelation effects with approximate main lobe width 1/

*δ*. Sampling effects become clear in two limiting cases. First, speckle fields with an autocorrelation width

*l*spanning many pixels (i.e., large average speckle size) have a band-limited power spectrum with support narrower than 1/

_{c}*δ*. In this limit the effect of discretization is negligible: lim

_{lc→∞}

**̂**

*J*_{0}(

*ν*) ≈

_{x}**̂(**

*J**ν*). As long as the average speckle size extends across several pixels (

_{x}*l*>

_{c}*δ*), approximating

**(Δ**

*J**x*) with a discretely sampled

*J*_{0}(Δ

*x*) remains accurate.

Second, in the limit of a small average speckle size, the pixel power spectrum cuts off the speckle field’s autocorrelation: lim*l _{c}*

_{→0}

**̂**

*J*_{0}(

*ν*) = sinc

_{x}^{2}(

*δν*). The modified autocorrelation width in this limit becomes

_{x}*l*′

*=*

_{c}*δ*, the pixel width. The discretized covariance matrix in the small speckle limit thus becomes $\mathit{J}={\sigma}_{0}^{2}\cdot \mathbf{I}$, the diagonal covariance matrix used to justify Markovity in the small speckle limit in Section 3.

Note that although each pixel integrates over multiple speckles when *l _{c}* <

*δ*at the detection plane, the first-order statistics of the random process

**will not change (the sum of correlated or uncorrelated Gaussian random variables remains Gaussian). Thus, for very small speckle with**

*A**l*≤

_{c}*δ*, the uncorrelated multivariate distribution will transform Eq. (2) into

*A*follows the circular symmetric complex univariate Gaussian distribution Eq. (1). The factorization of Eq. (31) indicates small speckle with

_{i}*l*≤

_{c}*δ*follows an i.i.d. complex Gaussian process, which obeys the Markov property.

## Appendix B: The Chapman-Kolmogorov test

To avoid checking the equality of Eq. (21) for all values of *m*, we rely on a first assumption that the pixel correlation function ** J**(Δ

*x*) is monotonically decreasing, valid for all speckle correlation functions experimentally tested. This common assumption is used in other “necessary condition" Markov tests for large sets of data [19, 20, 25]. Checking Eq. (21) with

*m*= 2,

*m*= 2 but not for some

*m*> 2.

For a random speckle field, we can further simplify the bilateral Markovity of Eq. (32) to the unilateral Markovity of Eq. (5). The following argument shows that bilateral Gaussian Markovity implies unilateral Gaussian Markovity. Consider the precision matrix Eq. (14) of an *N*-variate first-order bilateral process (extension to *g ^{th}* order is direct). The only non-zero off-diagonal entry in
${\mathit{J}}_{e}^{-1}$’s last row is at (

*N*,

*N*− 1), which by the conditional probability Eq. (10) indicates the variable

*A*also satisfies the unilateral Markov property. One can check that the unilateral Markov property holds for any

_{N}*A*(1 ≤

_{p}*p*<

*N*) as follows. Since the unilateral Markov property of variable

*A*does not concern any

_{p}*A*with

_{q}*q*>

*p*, we only need to look at the joint distribution of (

*A*

_{1},...,

*A*). Because the distribution is Gaussian, the

_{p}*p*×

*p*principal minor

*J**of*

_{p}

*J**is the covariance matrix of (*

_{e}*A*

_{1},...,

*A*). A direct computation shows the corresponding precision matrix ${\mathit{J}}_{p}^{-1}$ has the the same bandwidth as ${\mathit{J}}_{e}^{-1}$. The last row of ${\mathit{J}}_{p}^{-1}$ also has one non-zero off diagonal entry at (

_{p}*p*,

*p*−1). Thus,

*A*also satisfies the unilateral Markov property for all

_{p}*p*. In general, the precision matrices

*J*^{−1}and ${\mathit{J}}_{p}^{-1}$ have equal bandwidths under the assumption of a spatially stationary random process. We stated this assumption in Section 2.

This simplification allows us to replace Eq. (32) with the unilateral condition

While Eq. (33) is a useful necessary condition test, its estimation requires construction of a three-dimensional conditional probability matrix, which scales poorly with a large state space. The CK test in Eq. (22) is derived from Eq. (33) by multiplying either side by (*A*

_{n}_{−2})(

*A*

_{n}_{−1}|

*A*

_{n}_{−2}) and simplifying to

*A*

_{n}_{−1}and dividing by

*P*(

*A*

_{n}_{−2}) in Eq. (34) then leads to

*m*> 2) by replacing

*A*

_{n}_{−2}with

*A*

_{n}_{−}

*and selecting any intermediate pixel to integrate over.*

_{m}Finally, since the above derivations are only valid for a 1D Markov process, the 2D speckle fields we capture in experiment are only tested along one dimension. Successfully applying a 1D Markov test to a 2D image follows from our assumption of an *x*–*y* separable autocorrelation function. While a similar derivation may be used to create a complete unilateral test in 2D, the result is a significantly more complex conditional probability relationship between 4 variables.

## Appendix C: Speckle intensity and the Markov property

The derivation that Cauchy-masked speckle fields follow a Markov process does not directly extend to speckle intensity. Intensity’s integration over phase changes Eq. (2)’s multivariate complex Gaussian to a multivariate exponential density (e.g. see the derivation of speckle intensity’s bivariate exponential density in Chapter 4 of [11]). While multivariate Gaussians are conveniently described in full by their mean vector and covariance matrix (i.e., the covariance matrix ** J** incorporates all conditional probability relationships), multivariate exponentials are not. The general form of a multivariate exponential density function is [26]

*ω*is a different correlation parameter. As with the speckle field, if

*n ≤*2 then Eq. (36) describes uncorrelated (

*n*= 1) or nearly uncorrelated (

*n*= 2) speckle intensity and simplifies to fulfill the first-order Markov condition. For any

*n*> 2, Eq. (36) relies on

*n*-order correlation relationships that cannot be controlled optically. These

^{th}*n*-order correlations prevent determination of Markov speckle intensity with average size extending over

^{th}*n*pixels.

Instead, modeling speckle intensity as a hidden Markov process (HMP) offers a more direct Markovity relationship. Underlying this HMP is a complex first-order Markov speckle field generated using the proposed Cauchy mask. The |*S*| detectable complex field values form the HMP’s discrete state space. The HMP’s observation space is the discrete set of |*V*| detectable speckle intensity values. The |*S*| x |*V*| emission matrix contains the conditional probability of observing intensity value *I _{v}* from any complex field value

*A*and takes a very simple form: most matrix entries will be zero except for those following the deterministic relationship

_{s}*I*= |

_{v}*A*|

_{s}^{2}.

## Acknowledgments

RH was supported in part by the National Defense Science and Engineering Graduate Fellowship Program. RYC was supported by Joel A. Tropp under AFOSR award FA9550-09-1-0643.

## References and links

**1. **P. A. Kelly, H. Derin, and K. D. Hartt, “Adaptive segmentation of speckle images using a hierarchical random field model,” IEEE Trans. Acoust., Speech Sig. Process. **36**(10), 1628–1640 (1988). [CrossRef]

**2. **B. Skoric, “On the entropy of keys derived from laser speckle: statistical properties of Gabor-transformed speckle,” J. Opt. A: Pure Appl. Opt **10**, 055304 (2008). [CrossRef]

**3. **H. J. Rabal and R. A. Braga, *Dynamic Laser Speckle and Applications* (CRC Press, 2009).

**4. **R. Pappu, B. Recht, J. Taylor, and N. Gershenfeld, “Physical one-way functions,” Science **297**, 1074376 (2002). [CrossRef]

**5. **Y. M. Wang, B. Judkewitz, C. DiMarzio, and C. Yang,“Deep-tissue focal fluorescence imaging with digitally time-reversed ultrasound-encoded light,” Nature Commun. **3**, 928 (2012). [CrossRef]

**6. **D. P. Kelly, J. E. Ward, U. Gopinathan, and J. T. Sheridan, “Controlling speckle using lenses and free
space,” Opt. Lett. **32**, 3394–3396 (2007). [CrossRef]

**7. **E. Mundry, K. Belkebir, J. Girard, J. Savatier, E. Moal, C. Nocoletti, M. Allain, and A. Sentenac, “Structured illumination microscopy using unknown speckle patterns,” Nat. Photonics **6**, 312–315 (2012). [CrossRef]

**8. **O. Lankoande, M. M. Hayat, and B. Santhanam, “Scene estimation from speckled synthetic aperture radar imagery: Markov random-field approach,” J. Opt. Soc. Am. A **23**, 1269–1272 (2006). [CrossRef]

**9. **R. T. Frankot and R. Chellappa, “Lognormal random-field models and their applications to radar image synthesis,” IEEE Trans. Geosci. Remote Sens. **25**, 2196–2212 (2002).

**10. **H. Xie, L. E. Pierce, and F. T. Ulaby, “SAR speckle reduction using wavelet denoising and Markov random field modeling,” IEEE Trans. Geosci. Remote Sens. **40**, 195–208 (1987).

**11. **J. Goodman, *Speckle Phenomena in Optics* (Ben Roberts and Company, 2007).

**12. **J. C. Dainty, *Topics in Applied Physics: Laser Speckle and Related Phenomena* (Springer-Verlag, 1984).

**13. **J. Grimmett and D. Stirzaker, *Probability and Random Processes*, 3rd ed. (Oxford University Press, 2001).

**14. **H. Derin and P. A. Kelly, “Discrete-index Markov-type random processes,” Proc. IEEE **77**, 1485–1510 (1989). [CrossRef]

**15. **H. Rue and L. Held, *Gaussian Markov Random Fields: Theory and Applications* (Chapman and Hall, 2005). [CrossRef]

**16. **H. Derin, P. A. Kelly, G. Veniza, and S. G. Labitt, “Modeling and segmentation of speckle images using complex data,” IEEE Trans. Geosci. Remote Sens. **40**(1), 76–87 (1990). [CrossRef]

**17. **Y. Ait-Sahalia, “Do interest rates really follow continuous-time Markov diffusions?,” *Tech Rep.*, University of Chicago (1997).

**18. **A. de Matos and M. Fernandes, “Testing the Markov property with high frequency data,” J. Econometrics **141**, 44–64 (2007). [CrossRef]

**19. **S. Park and V. S. Pande, “Validation of Markov state models using Shannon’s entropy,” J. Chem. Phys **124**, 054118 (2006). [CrossRef]

**20. **B. Chen and Y. Hong, “Testing for the Markov property in time series,” Econ. Theory **28**, 130–178 (2012). [CrossRef]

**21. **T. W. Anderson and L. A. Goodman, “Statistical inference about Markov chains,” Ann. Math. Statist. **28**(1), 89–110 (1957). [CrossRef]

**22. **I. Yamaguchi and T. Zhang, “Phase-shifting digital holography,” Opt. Lett. **22**(16), 1268–1270 (1997). [CrossRef]

**23. **M. C. W. van Rossum and T. M. Nieuwenhuizen, “Multiple scattering of classical waves: microscopy, mesoscopy and diffusion,” Rev. Mod. Phys. **71**, 313–369 (1999). [CrossRef]

**24. **T. M. Cover and J. A. Thomas, *Elements of Information Theory* (John Wiley and Sons, Inc., 1991), chap. 11. [CrossRef]

**25. **W. C. Swope, J. W. Pitera, and F. Suits, “Describing protein folding kinetics by molecular dynamics simulations 1. theory,” J. Phys. Chem. B **108**, 6571–6581 (2004). [CrossRef]

**26. **A. W. Marshall and I. Olkin, “A multivariate exponential distribution” J. Amer. Statist. Assoc. **62**, 30–44 (1967). [CrossRef]