Markov speckle for efficient random bit generation

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 Optical Society of America OCIS codes: (030.6140) Speckle; (110.6150) Speckle Imaging. 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 Signal Process. 36(10), 1628–1640 (1988). 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). 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). 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). 6. D. P. Kelly, J. E. Ward, U. Gopinathan, and J. T. Sheridan, ‘‘Controlling speckle using lenses and free space,’’ Opt. Lett. 32, 23 3394–3396 (2007). 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,’’ Nature Photon. 6, 312–315 (2012). 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). 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, Third Edition (Oxford University Press, 2001). 14. H. Derin and P. A. Kelly, ‘‘Discrete-index Markov-type random processes,’’ Proc. IEEE 77, 1485–1510 (1989). 15. H. Rue and L. Held, Gaussian Markov Random Fields: Theory and Applications (Chapman and Hall, 2005). 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). 17. Y. Ait-Sahalia, ‘‘Do interest rates really follow continuous-time Markov diffusions?,’’ Tech Rp., 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). 19. S. Park and V. S. Pande, ‘‘Validation of Markov state models using Shannon’s entropy,’’ J. Chem. Phys 124, 054118 (2006). 20. B. Chen and Y. Hong, ‘‘Testing for the Markov property in time series,’’ Econ. Theory 28, 130–178 (2012). 21. T. W. Anderson and L. A. Goodman, ‘‘Statistical inference about Markov chains,’’ Ann. Math. Statist. 28(1), 89–110 (1957). 22. I. Yamaguchi and T. Zhang, ‘‘Phase-shifting digital holography,’’ Opt. Lett. 22(16), 1268–1270 (1997). 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). 24. T. M. Cover and J. A. Thomas, Elements of Information Theory (John Wiley and Sons, Inc., 1991), Ch. 11. 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). 26. A. W. Marshall and I. Olkin, ‘‘A multivariate exponential distribution’’ J. Amer. Statist. Assoc. 62, 30–44 (1967).


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.

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.

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) with σ 0 = lim , where |φ | is a random phasor [11].Since we are concerned with correlations between multiple pixels, Eq. ( 1) can be transformed to a complex random vector A A A of the speckle field at a set of N neighboring pixels along x.As each element of A A A remains circular symmetric in the presence of correlations, the pdf of the circular symmetric complex Gaussian speckle vector A A A is zero-mean Gaussian: where J J J is the covariance matrix of the complex field A A A representing the correlations between all N pixels.By definition, J J J is a positive definite symmetric matrix.J J J −1 is often referred to as a precision matrix.Within the context of optics, J J 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 J 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 J J(x 1 , x 2 ) is a weighted integral of a function M(η) that describes the illuminated area's shape at the scattering surface: with wavelength λ 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: where F is the Fourier transform operator.Finally, it is useful to define the average speckle size l c as the width of the autocorrelation J J 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 J J(Δx).Also, from Eq. (4) it is clear the mask function M(η) uniquely determines J J J(Δx).Thus, our goal is to first establish a sufficient condition on J J J(Δx) that guarantees Markov speckle, and then to manipulate M(η) to ensure J J 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.

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 A = (A 1 , . . ., A N ).The 2D case is examined in Section 4. The discretized complex values each pixel detects define the finite state space S of the random process.A detected speckle field A A A with values in S is a first-order unilateral Markov process, or a Markov chain, if the conditional probability satisfies The one-sided conditioning in Eq. ( 5) must hold for all 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: Also, a Markov process corresponding to a spatially stationary speckle field is homogeneous, that is, ) for all n.For completeness, assumed Markov processes are also aperiodic, irreducible and reversible.Finally, a transition matrix P P P is useful when representing the evolution of a homogeneous unilateral Markov process.P P P tabulates the conditional probabilities of transitioning from any state α to any state β at an immediately neighboring pixel: ).An important equation that the matrix P P P must satisfy is the Chapman-Kolmogorov equation [13], which we use in Section 6: P P P mn = P P P mv • P P P vn .(7) Here, v is an intermediate pixel between pixels m and n, and P P P mv and P P P vn are the transition matrices from pixel 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 for 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 {ρ k } are coefficients and {U n } are independent Gaussian random variables such that U n and A m are independent for 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 J J, fully describe its dependence relationships.Specifically, if J −1 nm = J J J −1 (n, m) = 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 n at pixel n, given field values A m at all other pixels: A similar expression is in [1] and [15].To satisfy the first-order Markov property, Eq. ( 10)'s conditional probability of A n must only depend on 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 −1 nm are zero except those between the matrix J J J −1 's super-diagonal and sub-diagonal.Strategies to ensure that J J J −1 is tridiagonal are examined next.

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 J J −1 .While small speckle satisfies this condition trivially, for large speckle we modify J J 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 J J J and precision matrix J J J −1 : where I is the identity matrix.A diagonal J J J −1 indicates no conditioned variables appear on the left side of Eq. (10).A discrete speckle field A A A with diagonal J J J is a purely random process.

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 J J −1 is tridiagonal.A tridiagonal J J J −1 transforms Eq. (10) to where This conditional probability follows the bilateral Markov property Eq. ( 6).We show in Appendix B that a spatially stationary (A 1 , . . ., A N ) also satisfies the unilateral Markov property Eq. ( 5).Thus, any stationary speckle field with a tridiagonal precision matrix J J J −1 is a first order Markov process of either type.We note the precision matrix of g th order Markov processes must have non-zero entries only between the positive and negative g th diagonals.A covariance matrix J J J e of the following exponential form generates a tridiagonal precision matrix: First Order Neighborhood Second Order Neighborhood Conditionally Independent of (i,j) Fig. 3. Speckle as a second-order Markov process in 2D with a neighborhood defined over 8 pixels (here only speckle amplitude is displayed).Independent of average speckle size, the conditional probability of each pixel in this Markov speckle field only depends on these 8 neighbors.
The corresponding inverse J J J −1 e is Here, σ = √ 2σ 0 and ρ are constants.For shift-invariant speckle, the exponential covariance matrix in Eq. ( 13) can be expressed as a one-dimensional autocorrelation function: where γ o = − ln(ρ).As expressed in Eq. ( 4), J J 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 J J J −1 e : where γ = − 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 J 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.

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 Ψ = {Ψ i j ⊂ L : (i, j) ∈ L}, where Ψ i j is a set of neighbors of (i, j).Informally, the neighborhood of pixel (i, j) is a set Ψ i j of nearby pixels that does not include (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: Here, Ω is a finite subset of pixels in the lattice that contains Ψ i, j but not (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, Here, {ρ k,l } is a set of correlation coefficients and {U i j } are independent random processes such that U i, j and A k,l are independent for (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 This separable correlation function is optically created by a separable apodizing mask, where again γ = − 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 nonseparable and higher-order Markov processes is possible following further investigation into GMRF theory [14,15].

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.

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.

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, for all values 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 P P P n,n−2 = P P P n,n−1 • P P P n−1,n−2 (22) 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 J 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 P P(α, β ) from a set of statistics is P P 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 P P n,n−1 , P P P n−1,n−2 ) or at alternate pixels (P P 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): where |S| represents the size of the state space.When r is zero, the sequence obeys the Markov property exactly.

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  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 A 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 P P P and r in Eq. ( 23).

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 J J(Δx)/J J J(0)|.The correlation area a and size of main lobe l c are proportional to propagation distance z for a given mask (see Eq. ( 4)).Likewise, a and l c are both inversely proportional to mask parameters ρ and w.We use two different Markov state space definitions to populate a transition matrix P P 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 P P n,n−1 , P P P n,n−1 P P P n−1,n−2 , P P P n,n−2 and (P P P n,n−2 − P P P n,n−1 P P 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 P 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 R R R 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 = 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.

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.

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  Fig. 9. Given a fixed desired speckle size, the entropy of a detected speckle field can be maximized using a Cauchy mask with an easily determined autocorrelation parameter ρ.
Table 1 offers calculated field entropies h(A) for several different autocorrelation functions assuming N = 100 and ρ = .9. *Note since the sinc autocorrelation's covariance matrix J J J is singular, its entropy is calculated using only J J J's nonzero eigenvalues.
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 J J is given by, h(N (0, J J J)) = 1  2 log (2πe) N • det(J J J) , 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 Here, we apply the assumption that the covariance matrix J J J is real.A real J J J is created by a symmetric apodizing mask function M(η), already assumed in Section 4. Since this modified covariance matrix indicates X X X and Y Y Y are uncorrelated, the entropy of a correlated speckle field We use Eq. ( 25) to calculate the entropy of speckle fields generated by several different mask functions, shown in Table 1.Parameter selection for these masks is discussed below.Burg's maximum entropy theorem [24] states that the maximum entropy rate stochastic process satisfying the autocorrelation constraints J J J(Δx) = ρ Δx for Δx = 0, 1...g is the g th order Gauss-Markov process in the form of Eq. ( 9), where the U n are i.i.d.∼ N (0, σ 2 0 ).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 M(η).To select the appropriate ρ Δx to satisfy this theorem's conditions, we note that speckle is often approximately assigned the 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 c , related to the aperture width 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, Given a speckle field generated through an open aperture of fixed width 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.

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 higherorder 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 δ : where 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 where the hatted functions represent a Fourier transform to spatial frequency coordinates ν x .A discrete field is represented as a sampled version of 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 F of the sampled autocorrelation function J J J 0 (Δx) in terms of A shows effects of pixelization: = F A(x) Π(x) A * (x) Π * (x) = Ĵ J J(ν x ) • sinc 2 (δ ν x ).
Here, Ĵ J J is the Fourier transform of the un-sampled autocorrelation J J J and denotes convolution.The squared 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 c spanning many pixels (i.e., large average speckle size) have a band-limited power spectrum with support narrower than 1/δ .In this limit the effect of discretization is negligible: lim l c →∞ Ĵ J J 0 (ν x ) ≈ Ĵ J J(ν x ).As long as the average speckle size extends across several pixels (l c > δ ), approximating J J J(Δx) with a discretely sampled J J 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 J 0 (ν x ) = sinc 2 (δ ν x ).The modified autocorrelation width in this limit becomes l c = δ , the pixel width.The discretized covariance matrix in the small speckle limit thus becomes J J J = σ 2 0 • 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 A A A will not change (the sum of correlated or uncorrelated Gaussian random variables remains Gaussian).Thus, for very small speckle with l c ≤ δ , the uncorrelated multivariate distribution will transform Eq. ( 2

) into p(A
where each A i follows the circular symmetric complex univariate Gaussian distribution Eq. ( 1).The factorization of Eq. ( 31) indicates small speckle with l c ≤ δ follows an i.i.d.complex Gaussian process, which obeys the Markov property.

Fig. 1 .
Fig. 1. Outline of an optical speckle field as a Markov process.(a) A 2D complex speckle field (phase in color) is examined along one dimension.If the speckle field obeys the firstorder Markov property, the conditional dependence of the field at pixel n (blue dot) will only depend on its immediate neighbors (red and green dots), and no other pixels.(b) This relationship can be visualized as a transition process between field values in complex space, or through transitions over an undirected graph.

Fig. 2 .
Fig. 2. Speckle is designed to follow a first-order Markov process using a Cauchydistributed apodizing mask M placed either (a) directly at the scatterer surface (z = 0) or (b) in the aperture plane of an imaging lens with focal length f .

Fig. 4 .
Fig. 4. Aperture windowing of the Cauchy mask function M(η) slightly shifts the desired autocorrelation function J J J(Δx) from an ideal exponential curve.

Fig. 5 .
Fig. 5. Diagram of the experimental setup used to generate correlated speckle field measurements.(a) K random speckle fields are generated via phase-shifting holography by imaging K random SLM patterns onto a volumetric scatterer.(b) Estimates of the three transition matrices in Eq. (22) are formed by vectorizing and processing the detected speckle fields.

Fig. 6 .
Fig. 6.An example set of transition matrices for the case of speckle modulated by a square aperture mask (unmodified speckle).(a) Transition matrices generated through simulation of square-masked speckle.(b) Transition matrices found experimentally.(c) Speckle field autocorrelation functions for simulated and experimental data.(d) Example square-masked speckle field data used for these plots.

Fig. 7 .
Fig.7.An example set of transition matrices for the case of speckle modulated by a Cauchy function mask in the same layout as Fig.6.Note the R 32 R 21 and R 31 matrices match more closely in width and slant than those generated by the square mask in Fig.6, leading to difference matrices ΔP 31 and ΔR 31 that are closer to 0.

TVFig. 8 .
Fig. 8. Plots comparing Markov TV error vs. speckle size for a Cauchy, Gaussian and square apodizing mask.Speckle generated via the Cauchy mask exhibits a lower TV error and thus is in closer agreement with a Markov process.