Speckle is the spatial fluctuation of irradiance seen when coherent light is reflected from a rough surface. It is due to light reflected from the surface’s many nooks and crannies accumulating vastly discrepant time delays, spanning much more than an optical period, en route to an observation point. Although speckle with continuous-wave (cw) illumination is well understood, the emerging interest in non-line-of-sight (NLoS) imaging using coherent light has created the need to understand the higher-order speckle that results from multiple rough-surface reflections, viz., speckled speckle and speckled speckled speckle. Moreover, the recent introduction of phasor-field (${{\mathcal {P}}}$-field) NLoS imaging—which relies on amplitude-modulated coherent illumination—requires pushing beyond cw scenarios for speckle and higher-order speckle. In this paper, we take first steps in addressing the foregoing needs using a three-diffuser transmissive geometry that is a proxy for three-bounce NLoS imaging. In the small-diffusers limit, we show that the irradiance variances of cw and modulated nth-order speckle coincide and are (2n − 1)-times those of ordinary (first-order) speckle. The more important case for NLoS imaging, however, involves extended diffuse reflectors. For our transmissive geometry with extended diffusers, we treat third-order cw speckle and first-order modulated speckle. Our results there imply that speckle is unlikely to impede successful operation of coherent-illumination cw imagers, and they suggest that the same might be true for ${{\mathcal {P}}}$-field imagers.
1. Introduction
When continuous-wave (cw) laser light that has been diffusely reflected by a rough surface illuminates an observation plane some distance away, a speckle pattern is visible in the measured irradiance. This speckle is due to wave-optical interference between reflections from independent wavelength-scale surface patches that vary in height by many optical wavelengths. This phenomenon, which we call ordinary (first-order) speckle, is well studied [1]. The irradiance produced by diffuse reflection is exponentially distributed at any observation point, and thus has a variance equal to its squared mean. The irradiance’s covariance function for diffuse reflection from an extended surface is also well understood, as are limited properties of the second-order speckle seen when the speckle pattern from a first diffuse reflection is reflected from a second rough surface and observed, as speckled speckle, on a new observation plane. However, with the growing interest in non-line-of-sight (NLoS) imaging, colloquially referred to as “seeing around corners”, there is now an unfulfilled need to understand higher-order speckle effects, e.g., speckled speckled speckle, such as results from laser light sequentially reflecting off three rough surfaces. Moreover, the advent of phasor-field (${{\mathcal {P}}}$-field) NLoS imagers [2–4]—which rely on amplitude-modulated coherent illumination—dictates that the preceding speckle questions be addressed for amplitude-modulated as well as for cw illumination.
In this paper, we take first steps in addressing the preceding issues using a paraxial, scalar-wave, transmissive geometry—which is a proxy for three-bounce NLoS imaging—that we have employed in our earlier treatments of ${{\mathcal {P}}}$-field imaging with quasimonochromatic coherent illumination [3,5]. After some preliminaries, which allow us to obtain the complete statistics of cw and modulated speckled speckled speckle in the small-diffusers limit, we begin our analyses in earnest with third-order speckle for cw illumination of extended diffusers. There, although small-diffuser speckled speckled speckle is seven times stronger than ordinary (first-order) speckle—i.e., its irradiance variance is seven times its squared mean—we find that cw speckled speckled speckle is highly mitigated by the geometry of the problem. In particular, our closed-form expression for the irradiance covariance of cw third-order speckle proves that the geometry of typical NLoS imaging scenarios reduces that third-order speckle to the ordinary cw speckle produced by the final diffuser. Furthermore, the speckle fluctuations that remain will be suppressed in power collection over any reasonable detector area. To quantify the impact of those residual power fluctuations, we evaluate the signal-to-noise ratio (SNR) for direct detection and show that power fluctuations act as an excess noise—above the fundamental shot-noise limit—that sets a maximum attainable SNR. For typical parameter values, we find this saturation SNR to be quite generous and thus conclude that the impact of third-order speckle is unlikely to be significant in NLoS imaging with cw coherent illumination.
Next, we move on to the speckle produced by the modulated illumination of extended diffusers, as used in ${{\mathcal {P}}}$-field imaging. Here, our analysis is limited to first-order speckle. We establish an upper bound on the zero-frequency component of the speckle when the initial illumination is space-time factorable finding that, at worst, such speckle is of ordinary strength. This result is seemingly at odds with Teichman’s analysis [4] for factorable, single-frequency modulation, which finds the modulation-frequency-component speckle to be stronger than ordinary speckle. To resolve the apparent discrepancy, we analyze a single-frequency-modulation limiting case in our ${{\mathcal {P}}}$-field framework and show that it recovers both Teichman’s result for the modulation-frequency speckle and our upper bound for the zero-frequency speckle. Then, using realistic parameter values for NLoS imaging scenarios, we conclude that the speckle enhancement effect reported by Teichman is likely to be minimal. Moving further, we analyze the first-order-speckle size for the modulated case and find that it is close to that of the cw case’s first-order speckle. If these similarities—for both speckle strength and speckle size—between modulated and cw speckle from extended diffusers persist in their second-order and third-order cases, then the adverse effects of speckled speckled speckle on ${{\mathcal {P}}}$-field NLoS imaging may be inconsequential.
2. Preliminaries
As developed in our earlier analysis of ${{\mathcal {P}}}$-field imaging [3,5,6], we use paraxial, scalar-wave optics in a transmissive geometry that serves as a proxy for a typical reflective, three-bounce NLoS geometry. The light at each plane is characterized by its baseband complex-field envelope $E_z({\boldsymbol \rho }_z,t)$, which modulates an optical carrier of frequency ${\omega _0}$ to produce a $\textrm {W}^{1/2}/\textrm {m}$-units optical field $\textrm {Re}[E_z({\boldsymbol \rho }_z,t)e^{-i{\omega _0} t}]$, where ${\boldsymbol \rho }_z$ is the 2D transverse spatial coordinate in the plane indicated by $z$. For cw speckle $E_z({\boldsymbol \rho }_z,t) = E_z({\boldsymbol \rho }_z)$ will have no time dependence, whereas for modulated speckle $E_z({\boldsymbol \rho }_z,t)$ will have bandwidth $\Delta \omega \ll \omega _0$. In both cases, $I_z({\boldsymbol \rho }_z,t) \equiv |E_z({\boldsymbol \rho }_z,t)|^{2}$ will be the short-time-average (STA) irradiance at the $z$-plane [7].
The geometry for our third-order speckle analysis is depicted in Fig. 1. Coherent space-time factorable illumination, $E_0({\boldsymbol \rho }_0,t) = E_0({\boldsymbol \rho }_0)S(t)$ with
(1)$$E_0({\boldsymbol \rho}_0) = \sqrt{I_0}\,e^{-4\left| {{\boldsymbol \rho}_0} \right|^{2}/d_0^{2}},$$
and $S(t) = 1$ for cw illumination, is incident at plane 0, which contains a diffuser with thickness profile $h_0({\boldsymbol \rho }_0)$. In the standard NLoS imaging configuration, this diffuser represents the visible wall at which reflection into the hidden space occurs. Plane 1 contains a second diffuser, whose thickness profile is $h_1({\boldsymbol \rho }_1)$ and whose size is modeled by a Gaussian pupil with $e^{-1}$-field-attenuation-diameter $d_1$. This diffuser represents a finite-sized, planar, diffuse target in the hidden scene, where, for simplicity, we have ignored any albedo variations across the target. Plane 2 contains a final diffuser, whose thickness profile is $h_2({\boldsymbol \rho }_2)$ and whose finite size is modeled by a Gaussian pupil with $e^{-1}$-field-attenuation-diameter $d_2$. In the NLoS scenario, it represents the visible wall where light reflects back to the imager. That imager’s entrance pupil lies in plane 3. Note that plane 2’s finite pupil enables us to obtain convergent paraxial-regime results for the variance of third-order cw speckle. A finite pupil at plane 0 is not needed for that purpose, because the initial illumination is self-limited to within that wall’s boundaries. The distances between planes 0, 1, 2, and 3 are all $L$, a choice made for convenience rather than necessity.2.1 Basic principles
All the analysis to follow rests on four basic principles: Fresnel diffraction for monochromatic light; the van Cittert–Zernike theorem for propagating the mutual coherence function (MCF) of spatially-incoherent light; the central limit theorem for sums of large numbers of independent random variables; and the law of iterated expectation.
To see how these principles come into play in the cw case, we start with how the initial illumination, $E_0({\boldsymbol \rho }_0)$ from Eq. (1), first passes through the plane-0 diffuser to become $E_0'({\boldsymbol \rho }_0)$ and then diffracts over a length $L$ free-space path to become the illumination, $E_1({\boldsymbol \rho }_1)$, at plane 1. We have that
(2)$$E'_0({\boldsymbol \rho}_0) = E_0({\boldsymbol \rho}_0)e^{i\omega_0h_0({\boldsymbol \rho}_0)/c},$$
where $c$ is light speed and we have ignored the diffuser’s refractive index. The essence of ${{\mathcal {P}}}$-field imaging [2] is that the diffusers are rough at the optical wavelength, $\lambda _0 = 2\pi c/\omega _0$, but smooth at the modulation wavelength, $\Delta \lambda = 2\pi c/\Delta \omega$ [8]. Our previous work [3,5,6] enforced this behavior by taking the $\{h_n({\boldsymbol \rho }_n): n = 0,1,2\}$ to be statistically-independent [9], identically-distributed, zero-mean Gaussian random processes whose covariance function (3)$$K_{hh}({\boldsymbol \rho}_n,{\boldsymbol{\tilde{\rho}}}_n) \equiv \langle h_n({\boldsymbol \rho}_n)h_n({\boldsymbol{\tilde{\rho}}}_n)\rangle = \sigma_h^{2} e^{-|{\boldsymbol \rho}_n-{\boldsymbol{\tilde{\rho}}}_n|^{2}/\rho_h^{2}},$$
has a standard deviation $\sigma _h$ satisfying $\lambda _0 \ll \sigma _h \ll \Delta \lambda$ and a correlation length $\rho _h$ satisfying $\rho _h \sim \lambda _0$. For the cw case, with $\langle \cdot \rangle _0$ denoting ensemble averaging over $h_0({\boldsymbol \rho }_0)$, these statistics imply that (4)$$\langle E'_0({\boldsymbol \rho}_0)\rangle_0 = E_0({\boldsymbol \rho}_0) e^{-{\omega_0}^{2}\sigma_h^{2}/2c^{2}} \approx 0,$$
(5)$$\langle E'_0({\boldsymbol \rho}_0)E'_0(\tilde{{\boldsymbol \rho}}_0)\rangle_0 = E_0({\boldsymbol \rho}_0)E_0({\boldsymbol{\tilde{\rho}}}_0) \exp\!\left[-{\omega_0}^{2}\sigma_h^{2}\left( {1+e^{-|{\boldsymbol \rho}_0-{\boldsymbol{\tilde{\rho}}}_0|^{2}/\rho_h^{2}}} \right)/c^{2}\right] \approx 0,$$
and (6)$$\langle E'_0({\boldsymbol \rho}_0) E^{\prime *}_0(\tilde{{\boldsymbol \rho}}_0)\rangle_0 = E_0({\boldsymbol \rho}_0)E^{*}_0(\tilde{{\boldsymbol \rho}}_0)\exp\!\left[-{\omega_0}^{2}\sigma_h^{2}\left( {1-e^{-|{\boldsymbol \rho}_0-{\boldsymbol{\tilde{\rho}}}_0|^{2}/\rho_h^{2}}} \right)/c^{2}\right] $$
(7)$$ \approx E_0({\boldsymbol \rho}_0)E^{*}_0(\tilde{{\boldsymbol \rho}}_0)\lambda_0^{2}\delta({\boldsymbol \rho}_0-\tilde{{\boldsymbol \rho}}_0), $$
where $\delta (\cdot )$ is the unit impulse and the approximation uses the fact that the smallest correlation area for a wavelength-$\lambda _0$ propagating wave is $\sim \lambda _0^{2}$ [10].Fresnel diffraction at frequency $\omega _0$ now gives
(8)$$E_1({\boldsymbol \rho}_1) = \frac{e^{i\omega_0L/c}}{i\lambda_0L}\int\!\textrm{d}^{2}{\boldsymbol \rho}_0\,E'_0({\boldsymbol \rho}_0)e^{i\omega_0 |{\boldsymbol \rho}_1-{\boldsymbol \rho}_0|^{2}/2cL},$$
and we can use this result, in conjunction with the approximation in (7), to obtain (9)$$\langle E_1({\boldsymbol \rho}_1)E_1^{*}(\tilde{{\boldsymbol \rho}}_1)\rangle_0 = \frac{e^{i\omega_0(|{\boldsymbol \rho}_1|^{2}-|\tilde{{\boldsymbol \rho}}_1|^{2})/2cL}}{L^{2}}\int\!\textrm{d}^{2}{\boldsymbol \rho}_0\,|E_0({\boldsymbol \rho}_0)|^{2}e^{i\omega_0(\tilde{{\boldsymbol \rho}}_1-{\boldsymbol \rho}_1)\cdot {\boldsymbol \rho}_0/cL},$$
via the van Cittert–Zernike theorem [1]. This result, together with Eq. (1), immediately proves that the cw case’s diffuser-averaged STA irradiance at plane 1 is independent of ${\boldsymbol \rho }_1$ and given by $\langle I_1\rangle \equiv \langle I_1({\boldsymbol \rho }_1)\rangle _0 = \pi d_0^{2}I_0/8L^{2}$ [11]. Now, by using Eq. (2) in Eq. (8), the central limit theorem tells us that $E_1({\boldsymbol \rho }_1)$ will be a zero-mean, complex-valued, Gaussian random process that is completely characterized by the MCF from Eq. (9). This result then implies that $I_1({\boldsymbol \rho }_1)$ is exponentially distributed and so has $\textrm {Var}[I_1({\boldsymbol \rho }_1)] = \langle I_1\rangle ^{2}$.The final principle our analysis will need—iterated expectation—comes in at this point. Conditioned on knowledge of $|E_1({\boldsymbol \rho }_1)|^{2}$, the procedure we have just employed can be used to show that $E_2({\boldsymbol \rho }_2)$ is a zero-mean, complex-valued, Gaussian random process that is completely characterized by
(10)$$\langle E_2({\boldsymbol \rho}_2)E_2^{*}(\tilde{{\boldsymbol \rho}}_2)\rangle_1 = \frac{e^{i\omega_0(|{\boldsymbol \rho}_2|^{2}-|\tilde{{\boldsymbol \rho}}_2|^{2})/2cL}}{L^{2}}\int\!\textrm{d}^{2}{\boldsymbol \rho}_1\,|E_1({\boldsymbol \rho}_1)|^{2}e^{-8|{\boldsymbol \rho}_1|^{2}/d_1^{2}}e^{i\omega_0(\tilde{{\boldsymbol \rho}}_2-{\boldsymbol \rho}_2)\cdot {\boldsymbol \rho}_1/cL},$$
where we have used $E_1'({\boldsymbol \rho }_1) = E_1({\boldsymbol \rho }_1)e^{-4|{\boldsymbol \rho }_1|^{2}/d_1^{2}}e^{i\omega _0h_1({\boldsymbol \rho }_1)/c}$ for the field emerging from the diffuser at plane 1 and $\langle \cdot \rangle _1$ denotes ensemble averaging over that diffuser. Iterated expectation now gives us (11)$$\langle E_2({\boldsymbol \rho}_2)E_2^{*}(\tilde{{\boldsymbol \rho}}_2)\rangle_{0,1} = \frac{e^{i\omega_0(|{\boldsymbol \rho}_2|^{2}-|\tilde{{\boldsymbol \rho}}_2|^{2})/2cL}}{L^{2}}\int\!\textrm{d}^{2}{\boldsymbol \rho}_1\,\langle |E_1({\boldsymbol \rho}_1)|^{2}\rangle_0\, e^{-8|{\boldsymbol \rho}_1|^{2}/d_1^{2}}e^{i\omega_0(\tilde{{\boldsymbol \rho}}_2-{\boldsymbol \rho}_2)\cdot {\boldsymbol \rho}_1/cL},$$
for the unconditional (ensemble averaged over the plane-0 and plane-1 diffusers) MCF of $E_2({\boldsymbol \rho }_2)$. Hence we get $\langle I_2\rangle \equiv \langle I_2({\boldsymbol \rho }_2) \rangle _{0,1} = \pi d_1^{2} \langle I_1\rangle /8L^{2}$ for the unconditional, diffuser-averaged STA irradiance at plane 2. By now it should be clear that we can pursue a similar argument to that just completed and show that $\langle I_3\rangle \equiv \langle I_3({\boldsymbol \rho }_3)\rangle _{0,1,2} = \pi d_2^{2}\langle I_2\rangle /8L^{2}$.We can obtain a further result for cw speckled speckle in the small-diffuser limit, wherein ($d_0d_1/4\lambda _0 L)^{2} \ll 1$. Small diffusers are especially interesting because they give rise to large speckles that prevent speckle-fluctuation suppression via spatial averaging, i.e., they imply worst-case signal-to-noise ratio behavior, as we will see later. In particular, within this regime we have that
(12)$$|E_1({\boldsymbol \rho}_1)|^{2} e^{-8|{\boldsymbol \rho}_1|^{2}/d_1^{2}} \approx |E_1(\textbf{0})|^{2} e^{-8|{\boldsymbol \rho}_1|^{2}/d_1^{2}}.$$
Consequently, the unconditional probability density function (pdf) for $I_2({\boldsymbol \rho }_2) = |E_2({\boldsymbol \rho }_2)|^{2}$ in the small-diffuser limit, i.e., the pdf for cw speckled speckle in that regime, is [1] (13)$$p_{I_2}(\mathcal{I}_2) = \int_0^{\infty}\!\textrm{d}\mathcal{I}_1\, \frac{\exp(-\mathcal{I}_1/\langle I_1\rangle)}{\langle I_1\rangle} \frac{\exp[-\mathcal{I}_2/(\pi d_1^{2}\mathcal{I}_1/8L^{2})]}{(\pi d_1^{2}\mathcal{I}_1/8L^{2})}u(\mathcal{I}_2) = \left[2K_0\!\left( {2\sqrt{\mathcal{I}_2/\langle I_2\rangle}} \right)/\langle I_2\rangle\right]\!u(\mathcal{I}_2),$$
where $K_0(\cdot )$ is the zeroth-order modified Bessel function of the second kind and $u(\cdot )$ is the unit-step function.The results in this subsection can be found in Goodman’s monograph [1]. We have reviewed them for two reasons. First, in Sec. 2.2, they will let us analyze both cw and modulated speckled speckle speckle in the small-diffusers limit. Second, in Secs. 3 and 4, respectively, they will be generalized to treat cw third-order speckle from extended diffusers and modulated first-order speckle from an extended diffuser.
2.2 Third-order speckle in the small-diffusers regime
Guided by the small-diffuser result for cw second-order speckle, and motivated by the desire to quantify worst-case speckle behavior for three-bounce NLoS imaging, let us consider cw third-order speckle when $(d_0d_1/4\lambda _0L)^{2} \ll 1$ and $(d_1d_2/4\lambda _0L)^{2} \ll 1$, i.e., when the Fresnel-number products for propagation between planes 0 and 1 and between planes 1 and 2 are both very small. The work from the previous subsection immediately shows us that
(14)$$|E_2({\boldsymbol \rho}_2)|^{2} e^{-8|{\boldsymbol \rho}_2|^{2}/d_2^{2}} \approx |E_2(\textbf{0})|^{2} e^{-8|{\boldsymbol \rho}_2|^{2}/d_2^{2}},$$
so that, conditioned on knowledge of $|E_2(\textbf {0})|^{2}$, we have that $E_3({\boldsymbol \rho }_3)$ is a zero-mean, complex-valued, Gaussian random process that is completely characterized by its conditional MCF, (15)$$\langle E_3({\boldsymbol \rho}_3)E_3^{*}(\tilde{{\boldsymbol \rho}}_3)\rangle_2 = \frac{e^{i\omega_0(|{\boldsymbol \rho}_3|^{2}-|\tilde{{\boldsymbol \rho}}_3|^{2})/2cL}}{L^{2}}|E_2(\textbf{0})|^{2}\int\!\textrm{d}^{2}{\boldsymbol \rho}_2\,e^{-8|{\boldsymbol \rho}_2|^{2}/d_2^{2}}e^{i\omega_0(\tilde{{\boldsymbol \rho}}_3-{\boldsymbol \rho}_3)\cdot {\boldsymbol \rho}_2/cL}.$$
The unconditional pdf for $I_3({\boldsymbol \rho }_3) = |E_3({\boldsymbol \rho }_3)|^{2}$ in the small-diffusers limit, i.e., the pdf for cw speckled speckled speckle in that regime, is therefore (16)$$p_{I_3}(\mathcal{I}_3) = \int_0^{\infty}\!\textrm{d}\mathcal{I}_2\, \frac{\displaystyle 2K_0\!\left( {2\sqrt{\mathcal{I}_2/\langle I_2\rangle}} \right)}{\displaystyle \langle I_2\rangle}\frac{\exp[-\mathcal{I}_3/(\pi d_2^{2}\mathcal{I}_2/8L^{2})]}{(\pi d_2^{2}\mathcal{I}_2/8L^{2})}u(\mathcal{I}_3).$$
Recourse to integral tables yields (17)$$p_{I_3}(\mathcal{I}_3) = [G_{0,3}^{3,0}(\mathcal{I}_3/\langle I_3\rangle \vert 0,0,0)/\langle I_3\rangle]u(\mathcal{I}_3),$$
where $G_{0,3}^{3,0}(\cdot \vert 0,0,0)$ is a Meijer $G$-function [12]. In Fig. 2 we have plotted the pdfs of the normalized irradiances $\tilde {I}_n \equiv I_n(\textbf {0})/\langle I_n\rangle$ for $n = 1, 2, 3$, which show the increasing randomness that occurs in progressing from speckle to speckled speckle to speckled speckled speckle. Indeed, more pdf iterations and properties of the Meijer $G$-function can be used to show that the unconditional pdf for $n$th-order speckle in the small-diffusers regime is [6] (18)$$p_{I_n}(\mathcal{I}_n) = [G_{0,n}^{n,0}(\mathcal{I}_n/\langle I_n\rangle \vert 0,\ldots,0)/\langle I_n\rangle]u(\mathcal{I}_n),$$
where $\langle I_n\rangle = (\Pi _{m=0}^{n-1}\pi d_m^{2}/8L^{2})I_0$. The normalized variance then turns out to be [6] (19)$$\textrm{NVar}_{I_n} \equiv \textrm{Var}[I_n(\textbf{0})]/\langle I_n\rangle^{2} = (2^{n}-1).$$
Interestingly, this subsection’s results for high-order cw speckle in the small-diffusers limit have immediate translations into corresponding results for high-order modulated speckle in that regime. In particular, with our assumption of quasimonochromatic, space-time factorable, modulated illumination, the small-diffuser assumption $(d_0d_1/4\lambda _0L)^{2} \ll 1$ leads to
(20)$$E_1({\boldsymbol \rho}_1,t)e^{-4|{\boldsymbol \rho}_1|^{2}/d_1^{2}} = e^{i\omega_0L/c}\left(\int\!\textrm{d}^{2}{\boldsymbol \rho}_0\,\frac{E'_0({\boldsymbol \rho}_0)}{i\lambda_0L}\right)S(t-L/c)e^{-4|{\boldsymbol \rho}_1|^{2}/d_1^{2}} $$
(21)$$ = E_1(\textbf{0})S(t-L/c)e^{-4|{\boldsymbol \rho}_1|^{2}/d_1^{2}}, $$
where $E_1(\textbf {0})$ is the cw-illumination complex envelope whose squared magnitude appears in Eq. (12). Similarly, because we have just shown that $E_1({\boldsymbol \rho }_1,t)$ for quasimonochromatic, space-time factorable, initial illumination is itself quasimonochromatic and space-time factorable when $(d_0d_1/4\lambda _0L)^{2} \ll 1$, we have that the additional small-diffuser assumption $(d_1d_2/4\lambda _0L)^{2} \ll 1$ leads to (22)$$E_2({\boldsymbol \rho}_2,t)e^{-4|{\boldsymbol \rho}_2|^{2}/d_2^{2}} = e^{i\omega_0L/c}\left(\int\!\textrm{d}^{2}{\boldsymbol \rho}_1\,\frac{E'_1({\boldsymbol \rho}_1)}{i\lambda_0L}\right)S(t-2L/c)e^{-4|{\boldsymbol \rho}_2|^{2}/d_2^{2}} $$
(23)$$ = E_2(\textbf{0})S(t-2L/c)e^{-4|{\boldsymbol \rho}_2|^{2}/d_2^{2}}, $$
where $E_2(\textbf {0})$ is the cw-illumination complex envelope whose squared magnitude appears in Eq. (14). So, except for their having time-delayed temporal modulations $|S(t-2L/c)|^{2}$ and $|S(t-3L/c)|^{2}$, the behaviors of second-order and third-order modulated speckle in the small-diffusers regime are identical to what we found for their cw counterparts.3. Third-order cw speckle from extended diffusers
First-order speckle has long been an issue for line-of-sight laser radars. For a rough-surfaced target, the single-pulse, single-pixel SNR of a heterodyne-detection laser radar asymptotes to a saturation SNR of 1—set by first-order speckle—with increasing target-return strength [13]. Direct-detection laser radars are largely immune to first-order speckle because each pixel is configured to contain sufficient speckles to average out their individual fluctuations without unduly compromising spatial resolution. NLoS laser imagers could potentially suffer third-order speckle’s seven-fold increased fluctuation strength that prevails with small diffuse reflectors. If unabated, this increase would result in a saturation SNR of 1/7. Whether or not such will be the case requires understanding the statistics of third-order speckle from extended targets, which is the case of interest for NLoS laser imagers. That task, for the cw case, is this section’s mission. We begin by relating direct detection’s saturation SNR to third-order speckle’s irradiance statistics.
Suppose plane 0 in Fig. 1 is illuminated with cw light from Eq. (1) and that a direct-detection system integrates the optical power transmitted through a diameter-$D$ circular pupil in plane 3 over the time interval $0\le t\le T$. We will neglect technical noises, e.g., thermal noise, and normalize the detector’s output to represent the number of detected photons, $N$, in that time interval. By the conditional Poissonian nature of photon-counting statistics for randomized laser light [14], we have that the resulting SNR is
(24)$$\textrm{SNR} \equiv \frac{\langle N\rangle^{2}}{\textrm{Var}(N)} = \frac{(\eta \langle P_3\rangle T/\hbar\omega_0)^{2}}{\eta \langle P_3\rangle T/\hbar\omega_0 + \eta^{2}\textrm{Var}(P_3)T^{2}/(\hbar\omega_0)^{2}},$$
where $\eta$ is the detector’s quantum efficiency, $\hbar \omega _0$ is the photon energy, and (25)$$P_3 = \int_{|{\boldsymbol \rho}_3|\le D/2}\!\textrm{d}^{2}{\boldsymbol \rho}_3\,I_3({\boldsymbol \rho}_3)$$
is the detected power. The first term in the SNR’s denominator is due to shot noise—the fundamental noise of semiclassical photodetection [14], which is always present—and the second term in that denominator is the excess noise associated with randomness in the detector’s illumination. We can rewrite Eq. (24) as (26)$$\textrm{SNR} = \frac{\textrm{SNR}_\textrm{sat}}{\textrm{SNR}_\textrm{sat}/\langle N\rangle + 1},$$
where the saturation SNR, (27)$$\textrm{SNR}_\textrm{sat} \equiv \frac{\langle P_3\rangle^{2}}{\textrm{Var}(P_3)} = \frac{(\pi D^{2}\langle I_3\rangle/4)^{2}}{\int_{|{\boldsymbol \rho}_3|\le D/2}\!\textrm{d}^{2}{\boldsymbol \rho}_3\int_{|\tilde{{\boldsymbol \rho}}_3|\le D/2}\!\textrm{d}^{2}\tilde{{\boldsymbol \rho}}_3\,\textrm{Covar}[I_3({\boldsymbol \rho}_3),I_3(\tilde{{\boldsymbol \rho}}_3)]},$$
is the maximum achievable SNR, and it is only approached (from below) as $\langle N\rangle \rightarrow \infty$. From Sec. 2.1 we have that $\langle I_3\rangle = \pi ^{3}d_0^{2}d_1^{2}d_2^{2} I_0/512 L^{6}$. All that remains, before we can evaluate $\textrm {SNR}_\textrm {sat}$, is to find plane 3’s normalized irradiance covariance, (28)$$\textrm{NCovar}_{I_3}({\boldsymbol \rho}_3-\tilde{{\boldsymbol \rho}}_3) \equiv \frac{\textrm{Covar}[I_3({\boldsymbol \rho}_3),I_3(\tilde{{\boldsymbol \rho}}_3)]}{\langle I_3\rangle^{2}},$$
where we are anticipating its being spatially homogeneous, as indeed will turn out to be the case. We will find this normalized irradiance covariance in the next subsection, using Gaussian moment factoring and iterated expectation. Along the way we will get the normalized covariances for $I_1({\boldsymbol \rho }_1)$ and $I_2({\boldsymbol \rho }_2)$, whose behaviors aid our understanding of how extended diffusers mitigate high-order speckle.3.1 Irradiance covariance of cw third-order speckle
The normalized covariance of plane 1’s irradiance is easily obtained. We know that $E_1({\boldsymbol \rho }_1)$ is a zero-mean, complex-valued, Gaussian random process that is completely characterized by the MCF from Eq. (9). Gaussian moment factoring gives us
(29)$$\textrm{Covar}[I_1({\boldsymbol \rho}_1),I_1(\tilde{{\boldsymbol \rho}}_1)] = \langle |E_1({\boldsymbol \rho}_1)|^{2}|E_1(\tilde{{\boldsymbol \rho}}_1)|^{2}\rangle_0-\langle I_1\rangle^{2} $$
(30)$$= |\langle E_1({\boldsymbol \rho}_1)E_1^{*}(\tilde{{\boldsymbol \rho}}_1)\rangle_0|^{2} $$
(31)$$=\left|\frac{e^{i\omega_0(|{\boldsymbol \rho}_1|^{2}-|\tilde{{\boldsymbol \rho}}_1|^{2})/2cL}}{L^{2}}\int\!\textrm{d}^{2}{\boldsymbol \rho}_0\,|E_0({\boldsymbol \rho}_0)|^{2}e^{i\omega_0(\tilde{{\boldsymbol \rho}}_1-{\boldsymbol \rho}_1)\cdot {\boldsymbol \rho}_0/cL}\right|^{2}, $$
which is spatially homogeneous, as presumed earlier. Plane 1’s normalized irradiance covariance is then found to be (32)$$\textrm{NCovar}_{I_1}({\boldsymbol \rho}_1-{\boldsymbol{\tilde{\rho}}}_1) = e^{-4 \Omega_{01} \left| {{\boldsymbol \rho}_1-{\boldsymbol{\tilde{\rho}}}_1} \right|^{2} / d_1^{2}},$$
where $\Omega _{01} \equiv (\pi d_0d_1/4\lambda _0L)^{2}$.Proceeding now toward obtaining plane 2’s irradiance covariance, we start from $E_2({\boldsymbol \rho }_2)$’s being—conditioned on knowledge of $E_1({\boldsymbol \rho }_1)$—a zero-mean, complex-valued, Gaussian random process that is completely characterized by its conditional MCF from Eq. (10). Gaussian moment factoring now gives us $I_2({\boldsymbol \rho }_2)$’s conditional correlation function,
(33)$$\langle {I_2({\boldsymbol \rho}_2)I_2({\boldsymbol{\tilde{\rho}}}_2)} \rangle_{1} = \langle {E_2({\boldsymbol \rho}_2)E_2^{*}({\boldsymbol \rho}_2)E_2({\boldsymbol{\tilde{\rho}}}_2)E_2^{*}({\boldsymbol{\tilde{\rho}}}_2)} \rangle_1 $$
(34)$$= \langle {I_2({\boldsymbol \rho}_2)} \rangle_1 \langle {I_2({\boldsymbol{\tilde{\rho}}}_2)} \rangle_1 + |\langle {E_2({\boldsymbol \rho}_2)E_2^{*}({\boldsymbol{\tilde{\rho}}}_2)} \rangle_1|^{2}. $$
Using Fresnel propagation, these terms expand to give (35)$$\begin{aligned}\langle {I_2({\boldsymbol \rho}_2)I_2({\boldsymbol{\tilde{\rho}}}_2)} \rangle_{1} =& \frac{1}{L^{4}} \Bigg[ \int\!\textrm{d}^{2}{{\boldsymbol \rho}_1}\,\int\!\textrm{d}^{2}{{\boldsymbol{\tilde{\rho}}}_1}\, I_1({\boldsymbol \rho}_1)I_1({\boldsymbol{\tilde{\rho}}}_1) e^{-8 (\left| {{\boldsymbol \rho}_1} \right|^{2} + \left| {{\boldsymbol{\tilde{\rho}}}_1} \right|^{2}) / d_1^{2}}\\ &\times \left( {1 + e^{-i\omega_0({\boldsymbol \rho}_1-{\boldsymbol{\tilde{\rho}}}_1)\cdot({\boldsymbol \rho}_2-{\boldsymbol{\tilde{\rho}}}_2) / c L}} \right) \Bigg]. \end{aligned}$$
Now, using the law of iterated expectation and taking advantage of the linearity of expectation, averaging over the first-diffuser’s statistics yields (36)$$\begin{aligned} \langle {I_2({\boldsymbol \rho}_2)I_2({\boldsymbol{\tilde{\rho}}}_2)} \rangle_{0,1} =& \frac{1}{L^{4}} \Bigg[ \int\!\textrm{d}^{2}{{\boldsymbol \rho}_1}\,\int\!\textrm{d}^{2}{{\boldsymbol{\tilde{\rho}}}_1}\, \langle {I_1({\boldsymbol \rho}_1)I_1({\boldsymbol{\tilde{\rho}}}_1)} \rangle_0 e^{-8 (\left| {{\boldsymbol \rho}_1} \right|^{2} + \left| {{\boldsymbol{\tilde{\rho}}}_1} \right|^{2}) / d_1^{2}}\\ &\times \left( {1 + e^{-i\omega_0({\boldsymbol \rho}_1-{\boldsymbol{\tilde{\rho}}}_1)\cdot({\boldsymbol \rho}_2-{\boldsymbol{\tilde{\rho}}}_2) / c L}} \right) \Bigg], \end{aligned}$$
for $I_2({\boldsymbol \rho }_2)$’s unconditional correlation function. Using this result we get (37)$$\textrm{NCovar}_{I_2}({\boldsymbol \rho}_2-{\boldsymbol{\tilde{\rho}}}_2) = e^{-4 \Omega_{01} \left| {{\boldsymbol \rho}_2-{\boldsymbol{\tilde{\rho}}}_2} \right|^{2} / d_0^{2}} + \frac{1}{1+\Omega_{01}} \left( {1 + e^{-4\Omega_{01} \left| {{\boldsymbol \rho}_2-{\boldsymbol{\tilde{\rho}}}_2} \right|^{2} / d_0^{2}({1+\Omega_{01})}}} \right),$$
which again is spatially homogeneous. This normalized covariance has interesting behavior with an interesting interpretation. When $\Omega _{01} \gg 1$, it reduces to (38)$$\textrm{NCovar}_{I_2}({\boldsymbol \rho}_2-{\boldsymbol{\tilde{\rho}}}_2) = e^{-(\pi d_1/2\lambda_0L)^{2} \left| {{\boldsymbol \rho}_2-{\boldsymbol{\tilde{\rho}}}_2} \right|^{2} }$$
which resembles the normalized covariance for first-order speckle, cf. (39)$$\textrm{NCovar}_{I_1}({\boldsymbol \rho}_1-{\boldsymbol{\tilde{\rho}}}_1) = e^{-(\pi d_0/2\lambda_0L)^{2} \left| {{\boldsymbol \rho}_1-{\boldsymbol{\tilde{\rho}}}_1} \right|^{2}},$$
which follows from Eq. (32). This is not an accidental coincidence. When $\Omega _{01} \gg 1$, the speckle size in $I_1({\boldsymbol \rho }_1)$ is much smaller than $d_1$. Moreover, $I_2({\boldsymbol \rho }_2)$ is conditionally exponential, given $I_1({\boldsymbol \rho }_1)$, with conditional mean $\langle I_2\rangle _1 = \int \!\textrm {d}^{2}{\boldsymbol \rho }_1\,I_1({\boldsymbol \rho }_1)e^{-8|{\boldsymbol \rho }_1|^{2}/d_1^{2}}/L^{2}$. By the law of large numbers, we have that $\int \!\textrm {d}^{2}{\boldsymbol \rho }_1\,I_1({\boldsymbol \rho }_1)e^{-8|{\boldsymbol \rho }_1|^{2}/d_1^{2}}/L^{2} \approx \pi d^{2}_1\langle I_1\rangle /8L^{2}$, equivalently $\langle I_2({\boldsymbol \rho }_2)\rangle _1 \approx \langle I_2\rangle$, because of speckle averaging over the plane-1 pupil. Furthermore, this means we can take $E_2'({\boldsymbol \rho }_2) \equiv E_2({\boldsymbol \rho }_2)e^{i\omega _0h_2({\boldsymbol \rho }_2)/c}$ to be a zero-mean, complex-valued, Gaussian random process insofar as calculating the statistics of $E_3({\boldsymbol \rho }_3)$ is concerned, i.e., $\Omega _{01} \gg 1$ has totally suppressed the speckle generated in propagation from plane 0 to plane 1 insofar as evaluating the speckle incurred in propagating from plane 1 to plane 2.Putting aside, for now, the $\Omega _{01} \gg 1$ condition and its consequences, it should be clear that even without that condition, we can proceed with an iterated-expectation procedure to find $\textrm {NCovar}_{I_3}({\boldsymbol \rho }_3-{\boldsymbol {\tilde {\rho }}}_3)$. For the sake of brevity, we will omit the details and just give the final answer:
(40)$$\begin{aligned} \textrm{NCovar}_{I_3}({\boldsymbol \rho}_3-{\boldsymbol{\tilde{\rho}}}_3) =& e^{-4 \Omega_{12} \left| {{\boldsymbol \rho}_3-{\boldsymbol{\tilde{\rho}}}_3} \right|^{2} / d_1^{2}} + \frac{1}{1+\Omega_{01}}\left( {1 + e^{-4 \Omega_{12} \left| {{\boldsymbol \rho}_3-{\boldsymbol{\tilde{\rho}}}_3} \right|^{2} / d_1^{2}}} \right)\\ &+ \frac{1}{1+\Omega_{12}} \left( { 1 + e^{-4\Omega_{12} \left| {{\boldsymbol \rho}_3-{\boldsymbol{\tilde{\rho}}}_3} \right|^{2} / d_1^{2}({1+\Omega_{12})}}} \right)\\ &+ \frac{1}{1+\Omega_{01}+\Omega_{12}} \left( {1 + e^{-4\Omega_{12}({1+\Omega_{01})} \left| {{\boldsymbol \rho}_3-{\boldsymbol{\tilde{\rho}}}_3} \right|^{2} / d_1^{2}({1+\Omega_{01}+\Omega_{12})}}} \right), \end{aligned}$$
where $\Omega _{12} \equiv (\pi d_1d_2/4\lambda _0L)^{2}$.When $\Omega _{01} \gg 1$ and $\Omega _{12} \gg 1$, the preceding normalized covariance becomes
(41)$$\textrm{NCovar}_{I_3}({\boldsymbol \rho}_3-{\boldsymbol{\tilde{\rho}}}_3) = e^{-(\pi d_2/2\lambda_0L)^{2} \left| {{\boldsymbol \rho}_3-{\boldsymbol{\tilde{\rho}}}_3} \right|^{2} },$$
which is the normalized covariance for first-order speckle produced by propagation from plane 2 to plane 3. Indeed, when $\Omega _{01} \gg 1$ and $\Omega _{12} \gg 1$, the law of large numbers implies $\langle I_2({\boldsymbol \rho }_2)\rangle _1 \approx \langle I_2\rangle$ and $\langle I_3({\boldsymbol \rho }_3)\rangle _2 \approx \langle I_3\rangle$, so what would have been speckled speckled speckle in $I_3({\boldsymbol \rho }_3)$ reduces to the first-order speckle for propagation from plane 2 to plane 3 [15].Taking values close to what we might expect in practice—$\lambda _0 = 532{\,\textrm {nm}}$ optical wavelength, $L = 1{\,\textrm {m}}$ to 10 m scene depth and standoff, $d_0=$1 mm to 1 cm spot size, $d_1=$3 cm to 2 m target size, and $d_2=$1 m to 10 m wall size—we find that $20 \le \Omega _{01} \le 10^{9}$ and $2\times 10^{9} \le \Omega _{12} \le 10^{15}$. The least favorable attenuation factor in Eq. (40), $1/(1+\Omega _{01}) \approx 0.05$, is already small enough to make cw third-order speckle reduce to first-order speckle.
3.2 Saturation signal-to-noise ratio
The saturation SNR is related to the normalized covariance $I_3({\boldsymbol \rho }_3)$ as follows,
(42)$$\textrm{SNR}_\textrm{sat} = \frac{(\pi D^{2}/4)^{2}}{\int_{|{\boldsymbol \rho}_3|\le D/2}\!\textrm{d}^{2}{\boldsymbol \rho}_3\int_{|\tilde{{\boldsymbol \rho}}_3|\le D/2}\!\textrm{d}^{2}\tilde{{\boldsymbol \rho}}_3\,\textrm{NCovar}_{I_3}({\boldsymbol \rho}_3-\tilde{{\boldsymbol \rho}}_3)}.$$
Switching to sum and difference coordinates, ${\boldsymbol \rho }_+=\left ( {{\boldsymbol \rho }_3+{\boldsymbol {\tilde {\rho }}}_3} \right )/2$ and ${\boldsymbol \rho }_-={\boldsymbol \rho }_3-{\boldsymbol {\tilde {\rho }}}_3$, the ${\boldsymbol \rho }_+$ integration yields (43)$$\textrm{SNR}_\textrm{sat} = \frac{(\pi D^{2}/4)^{2}}{\int_{|{\boldsymbol{\rho}_-}| \le D}\textrm{d}^{2}{\boldsymbol \rho}_-\, \textrm{NCovar}_{I_3}({\boldsymbol{\rho}_-}) O({\boldsymbol{\rho}_-},D)}$$
where (44)$$O({\boldsymbol \rho}_-,D) =\frac{\displaystyle D^{2}}{\displaystyle 2}\left[ {\cos^{-1}\!\left( {\frac{\displaystyle |{\boldsymbol \rho}_-|}{\displaystyle D}} \right)-\frac{\displaystyle |{\boldsymbol \rho}_-|}{\displaystyle D}\sqrt{1-\frac{\displaystyle |{\boldsymbol \rho}_-|^{2}}{\displaystyle D^{2}}}} \right], \textrm{ for}\, 0 \leq |{\boldsymbol \rho}_-| \leq D,$$
is the two-circle overlap function. An exact evaluation of Eq. (43) is tedious but results in (45)$$\begin{aligned}\textrm{SNR}_\textrm{sat} = &\left[ \frac{1}{1+\Omega_{01}}+ \frac{1}{1+\Omega_{12}} + \frac{1}{1+\Omega_{01}+\Omega_{12}} + \frac{d_1^{2}}{D^{2}\Omega_{12}}\left( {1+\sqrt{\Omega_{12}}} \right)\left( {1+\frac{1}{1+\Omega_{01}}} \right)\right.\\ & - \frac{d_1^{2} (2+\Omega_{01})}{D^{2}\sqrt{\Omega_{12}} (1+\Omega_{01})} B\!\left( {\frac{2D^{2} \sqrt{\Omega_{12}}}{d_1^{2}}} \right) - \frac{d_1^{2}}{D^{2}\Omega_{12}} B\!\left( {\frac{2D^{2} \Omega_{12}}{d_1^{2} (1+\Omega_{12})}} \right)\\ & \left. - \frac{d_1^{2}}{D^{2}\Omega_{12} (1+\Omega_{01})} B\!\left( {\frac{2D^{2} \Omega_{12} (1+\Omega_{01})}{d_1^{2} (1+\Omega_{01} + \Omega_{12})}} \right) \right]^{-1}, \end{aligned}$$
where (46)$$B(x) \equiv e^{-x}[{\textrm{BI}}_0(x) + {\textrm{BI}}_1(x)],$$
with ${\textrm {BI}}_n$ being the $n$th-order modified Bessel function of the first kind.In Fig. 3 we have plotted $\textrm {SNR}_\textrm {sat}$ in dB versus $D/d_1$ for four representative $(\Omega _{01},\Omega _{12})$ pairs: both being at their maximum values given at the end of Sec. 3.1, both at their minimum values given there, both at intermediate values between those maxima and minima, and both having unit values. As a specific example, let $\lambda _0 = 532$ nm, $L = 1$ m, $d_0 = 1$ mm, $d_1 = 10$ cm, and $d_2 = 1$ m. We then have $\Omega _{01} = 2.18 \times 10^{4}$ and $\Omega _{12}= 2.18 \times 10^{10}$, implying that $\textrm {SNR}_\textrm {sat} = 31.4$ dB for $D=1$ cm. In short, modest detector diameters will reduce cw speckled speckled speckle’s strength to a miniscule amount within the range of typical Fresnel-number products. In contrast, the $(\Omega _{01},\Omega _{12}) = (1,1)$ curve in Fig. 3 shows the approach to the small-diffusers, $\textrm {SNR}_\textrm {sat} = 1/7 \approx 0.143$, limit, viz., for unit Fresnel-number products we have $\textrm {SNR}_\textrm {sat} \approx 0.27$ for $D$ between 1 mm and 1 cm.
It should be noted, however, that speckle averaging at the detector is not always desirable. This is especially true for Willomitzer et al.’s [16] synthetic-wavelength-holography approach to NLoS imaging. It is a variant of ${{\mathcal {P}}}$-field imaging that: (1) uses sequential, cw illumination at two optical frequencies; (2) heterodyne detects each frequency’s $E_3({\boldsymbol \rho }_3)$ at high spatial resolution, using a detector array, to obtain its speckle pattern; and (3) forms a ${{\mathcal {P}}}$-field image of the hidden-space’s target plane by processing the two speckle patterns. Spatial integration over multiple speckles at each detector element degrades the speckle-pattern measurements and is thus undesirable [5].
4. Modulated first-order speckle from an extended diffuser
Our reason for studying modulated speckle is its potential relevance to ${{\mathcal {P}}}$-field NLoS imaging. With the exception of Teichman’s work [4], prior theoretical treatments of ${{\mathcal {P}}}$-field imaging have ignored the possible ill-effects of high-order speckle on such systems. The ${{\mathcal {P}}}$ field, as defined in [3], is the temporal Fourier transform of the diffuser-averaged STA irradiance,
(47)$${{\mathcal{P}}}_z({\boldsymbol \rho}_z,{\omega_-}) \equiv \int\!\textrm{d}t\, \langle I_z({\boldsymbol \rho}_z,t)\rangle e^{i{\omega_-} t} = \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\, \langle {{\mathcal{E}}}_z({\boldsymbol \rho}_z,\omega_+ +{\omega_-}/2) {{\mathcal{E}}}_z^{*}({\boldsymbol \rho}_z,\omega_+-{\omega_-}/2)\rangle,$$
where ${{\mathcal {E}}}_z({\boldsymbol \rho }_z,\omega )\equiv \int {\textrm {d}t}\, E_z({\boldsymbol \rho }_z,t) e^{i\omega t}$ and $\langle \cdot \rangle$ denotes ensemble averaging over all relevant diffusers. To analyze high-order speckle’s impact on estimating ${{\mathcal {P}}}_z({\boldsymbol \rho }_z,{\omega _-})$ from experimental data, we introduce (48)$${\tilde{\mathcal{P}}}_z({\boldsymbol \rho}_z,{\omega_-}) \equiv \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\, {{\mathcal{E}}}_z({\boldsymbol \rho}_z,\omega_++{\omega_-}/2) {{\mathcal{E}}}_z^{*}({\boldsymbol \rho}_z,\omega_+-{\omega_-}/2),$$
whose ensemble average equals ${{\mathcal {P}}}_z({\boldsymbol \rho }_z,{\omega _-})$. Unfortunately, despite our having found quantitative—indeed favorable—results for the cw third-order speckle from extended diffusers, analysis of the modulated speckle for extended diffusers is far more challenging. Accordingly, we will limit ourselves to the first-order case, i.e., characterizing the ${\tilde {\mathcal {P}}}$-field fluctuations at plane 1, where the hidden target would be located. Our goals will be the same as those we set for cw speckle: determining ${\tilde {\mathcal {P}}}_1({\boldsymbol \rho }_1,{\omega _-})$’s speckle strength and speckle size.4.1 Speckle strength
The spatial incoherence created by plane 0’s diffuser plus the central limit theorem imply that $E_1({\boldsymbol \rho }_1,t)$ is a zero-mean, complex-valued, Gaussian random process that is completely characterized by its space-time MCF $\langle E_1({\boldsymbol \rho }_1,t)E_1^{*}(\tilde {{\boldsymbol \rho }}_1,\tilde {t})\rangle _0$. It follows that $\mathcal {E}_1({\boldsymbol \rho }_1,{\omega})$ is also a zero-mean, complex-valued, Gaussian random process, and its complete characterization is given by its space-frequency MCF $\langle \mathcal {E}_1({\boldsymbol \rho }_1,\omega )\mathcal {E}_1^{*}(\tilde {{\boldsymbol \rho }}_1,\tilde {\omega })\rangle _0$. Gaussian moment factoring now gives us
(49)$$\begin{aligned}&\langle {|{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,{\omega_-})|^{2}} \rangle_0\\ &= \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\, \int\!\frac{\textrm{d}{{\tilde{\omega}_+}}}{2\pi}\,\\ & \langle { {{\mathcal{E}}}_1({\boldsymbol \rho}_1,\omega_++{\omega_-}/2) {{\mathcal{E}}}_1^{*}({\boldsymbol \rho}_1,\omega_+-{\omega_-}/2) {{\mathcal{E}}}_1^{*}({\boldsymbol \rho}_1,{\tilde{\omega}_+}+{\omega_-}/2) {{\mathcal{E}}}_1({\boldsymbol \rho}_1,{\tilde{\omega}_+}-{\omega_-}/2)} \rangle_0\\ &= \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\, \int\!\frac{\textrm{d}{{\tilde{\omega}_+}}}{2\pi}\,\\ & \Bigg [ \langle {{{\mathcal{E}}}_1({\boldsymbol \rho}_1,\omega_++{\omega_-}/2) {{\mathcal{E}}}_1^{*}({\boldsymbol \rho}_1,\omega_+-{\omega_-}/2)} \rangle_0 \langle {{{\mathcal{E}}}_1^{*}({\boldsymbol \rho}_1,{\tilde{\omega}_+}+{\omega_-}/2) {{\mathcal{E}}}_1({\boldsymbol \rho}_1,{\tilde{\omega}_+}-{\omega_-}/2)} \rangle_0\\ &+ \langle {{{\mathcal{E}}}_1({\boldsymbol \rho}_1,\omega_++{\omega_-}/2) {{\mathcal{E}}}_1^{*}({\boldsymbol \rho}_1,{\tilde{\omega}_+}+{\omega_-}/2)} \rangle_0 \langle {{{\mathcal{E}}}_1^{*}({\boldsymbol \rho}_1,\omega_+-{\omega_-}/2) {{\mathcal{E}}}_1({\boldsymbol \rho}_1,{\tilde{\omega}_+}-{\omega_-}/2)} \rangle_0 \Bigg]. \end{aligned}$$
The first term, after integration, is $|{{\mathcal {P}}}_1({\boldsymbol \rho }_1,{\omega _-})|^{2}$. Hence the remaining term, which can be expanded via the Fresnel-diffraction formula, is the ${\tilde {\mathcal {P}}}$ field’s variance, viz., (50)$$\begin{aligned}& \textrm{Var}[{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,{\omega_-})] = \langle {|{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,{\omega_-})|^{2}} \rangle_0 - \left| {{{\mathcal{P}}}_1({\boldsymbol \rho}_1,{\omega_-})} \right|^{2}\\& \quad = \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\, \int\!\frac{\textrm{d}{{\tilde{\omega}_+}}}{2\pi}\,\\ & \quad \quad \langle {{{\mathcal{E}}}_1({\boldsymbol \rho}_1,\omega_++{\omega_-}/2) {{\mathcal{E}}}_1^{*}({\boldsymbol \rho}_1,{\tilde{\omega}_+}+{\omega_-}/2)} \rangle_0 \langle {{{\mathcal{E}}}_1^{*}({\boldsymbol \rho}_1,\omega_+-{\omega_-}/2) {{\mathcal{E}}}_1({\boldsymbol \rho}_1,{\tilde{\omega}_+}-{\omega_-}/2)} \rangle_0 \\& \quad = \frac{1}{L^{4}} \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\, \int\!\frac{\textrm{d}{{\tilde{\omega}_+}}}{2\pi}\, \int\!\textrm{d}^{2}{{\boldsymbol \rho}_0}\, \int\!\textrm{d}^{2}{{\boldsymbol{\tilde{\rho}}}_0}\,\\ & \quad \quad {{\mathcal{E}}}_0({\boldsymbol \rho}_0,\omega_++{\omega_-}/2) {{\mathcal{E}}}_0^{*}({\boldsymbol \rho}_0,{\tilde{\omega}_+}+{\omega_-}/2) {{\mathcal{E}}}_0^{*}({\boldsymbol{\tilde{\rho}}}_0,\omega_+-{\omega_-}/2) {{\mathcal{E}}}_0({\boldsymbol{\tilde{\rho}}}_0,{\tilde{\omega}_+}-{\omega_-}/2)\\ \end{aligned}$$
(51)$$\times \; e^{i \left( {\omega_+-{\tilde{\omega}_+}} \right)\left( {\left| {{\boldsymbol \rho}_1-{\boldsymbol \rho}_0} \right|^{2}-\left| {{\boldsymbol \rho}_1-{\boldsymbol{\tilde{\rho}}}_0} \right|^{2}} \right)/ 2 c L}. $$
Invoking $E_0({\boldsymbol \rho }_0,t)$’s being space-time factorable, Eq. (51) becomes (52)$$\begin{aligned} & \textrm{Var}[{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,{\omega_-})]\\& \quad = \frac{1}{L^{4}} \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\, \int\!\frac{\textrm{d}{{\tilde{\omega}_+}}}{2\pi}\, \int\!\textrm{d}^{2}{{\boldsymbol \rho}_0}\, \int\!\textrm{d}^{2}{{\boldsymbol{\tilde{\rho}}}_0}\, \left| {E_0({\boldsymbol \rho}_0)} \right|^{2} \left| {E_0({\boldsymbol{\tilde{\rho}}}_0)} \right|^{2} e^{i \left( {\omega_+-{\tilde{\omega}_+}} \right)\left( {\left| {{\boldsymbol \rho}_1-{\boldsymbol \rho}_0} \right|^{2}-\left| {{\boldsymbol \rho}_1-{\boldsymbol{\tilde{\rho}}}_0} \right|^{2}} \right)/ 2 c L}\\ & \quad \times{{\mathcal{S}}}(\omega_++{\omega_-}/2) {{\mathcal{S}}}^{*}({\tilde{\omega}_+}+{\omega_-}/2) {{\mathcal{S}}}^{*}(\omega_+-{\omega_-}/2) {{\mathcal{S}}}({\tilde{\omega}_+}-{\omega_-}/2), \end{aligned}$$
where $\mathcal {S}(\omega ) \equiv \int \!\textrm {d}t\,S(t)e^{i\omega t}$.Equation (52) does not lend itself to further evaluation, but it does allow us to establish an upper bound on the variance of ${\tilde {\mathcal {P}}}_1({\boldsymbol \rho }_1,0)$,
(53)$$\begin{aligned} & \textrm{Var}[{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,0)]\\& \quad = \frac{1}{L^{4}} \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\, \int\!\frac{\textrm{d}{{\tilde{\omega}_+}}}{2\pi}\, \int\!\textrm{d}^{2}{{\boldsymbol \rho}_0}\, \int\!\textrm{d}^{2}{{\boldsymbol{\tilde{\rho}}}_0}\, \left| {E_0({\boldsymbol \rho}_0)} \right|^{2} \left| {E_0({\boldsymbol{\tilde{\rho}}}_0)} \right|^{2} \left| {{{\mathcal{S}}}(\omega_+)} \right|^{2} \left| {{{\mathcal{S}}}({\tilde{\omega}_+})} \right|^{2}\\ & \quad \quad \times e^{i \left( {\omega_+-{\tilde{\omega}_+}} \right)\left( {\left| {{\boldsymbol \rho}_1-{\boldsymbol \rho}_0} \right|^{2}-\left| {{\boldsymbol \rho}_1-{\boldsymbol{\tilde{\rho}}}_0} \right|^{2}} \right)/ 2 c L}\\ & \quad \leq \frac{1}{L^{4}} \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\, \int\!\frac{\textrm{d}{{\tilde{\omega}_+}}}{2\pi}\, \int\!\textrm{d}^{2}{{\boldsymbol \rho}_0}\, \int\!\textrm{d}^{2}{{\boldsymbol{\tilde{\rho}}}_0}\, \left| {E_0({\boldsymbol \rho}_0)} \right|^{2} \left| {E_0({\boldsymbol{\tilde{\rho}}}_0)} \right|^{2} \left| {{{\mathcal{S}}}(\omega_+)} \right|^{2} \left| {{{\mathcal{S}}}({\tilde{\omega}_+})} \right|^{2}\\ \end{aligned}$$
(54)$$ =\left| {\frac{1}{L^{2}}\int\!\textrm{d}^{2}{{\boldsymbol \rho}_0}\,{{\mathcal{P}}}_0({\boldsymbol \rho}_0,0)} \right|^{2} = \left| {{{\mathcal{P}}}_1({\boldsymbol \rho}_1,0)} \right|^{2}, $$
which shows that the modulated first-order speckle at zero frequency from an extended diffuser is never stronger than ordinary (cw first-order) speckle.The preceding bound is seemingly at odds with Teichman’s analysis [4] for factorable, single-frequency modulation, which finds the modulation-frequency-component speckle to be stronger than ordinary speckle. It turns out, however, that this apparent discrepancy is because Teichman’s ${{\mathcal {P}}}$-field speckle analysis applies to the modulation-frequency component, whereas our bound applies to the zero-frequency component. We will harmonize the two results by assuming
(55)$$S(t) = e^{-t^{2}/T^{2}}\cos({\Omega} t) $$
(56)$$ {{\mathcal{S}}}(\omega) = \frac{T\sqrt{\pi}}{2} \left( {e^{-T^{2}(\omega-{\Omega})^{2}/4} + e^{-T^{2}(\omega+{\Omega})^{2}/4}} \right), $$
with ${\Omega } T \gg 1$, so that $\mathcal {S}(\omega )$ is confined to narrow bandwidths about $\omega = \pm \Omega$. Then, because $|S(t)|^{2} = e^{-2t^{2}/T^{2}}[1+\cos (2{\Omega } t)]/2$, we will be concerned with the speckle strengths in ${\tilde {\mathcal {P}}}_1({\boldsymbol \rho }_1,2\Omega )$ and ${\tilde {\mathcal {P}}}_1({\boldsymbol \rho }_1,0)$.The narrowband nature of $\mathcal {S}(\omega )$ allows us to use
(57)$${{\mathcal{S}}}({\omega_+}+{\Omega}) {{\mathcal{S}}}^{*}({\tilde{\omega}_+}+{\Omega}) {{\mathcal{S}}}^{*}({\omega_+}-{\Omega}) {{\mathcal{S}}}({\tilde{\omega}_+}-{\Omega}) \approx \frac{\pi^{2} T^{4}}{16} e^{-T^{2}({\omega_+}^{2}+{\tilde{\omega}_+}^{2})/2}$$
in finding $\textrm {Var}[{\tilde {\mathcal {P}}}_1({\boldsymbol \rho }_1,2\Omega )]$ from Eq. (52), where it will function as a low-pass filter. As a result, we have $\left | {\omega_+ -{\tilde {\omega }_+}} \right |/2cL \ll {\Omega }/2cL = \pi / \Lambda L$, where $\Lambda =2\pi c/{\Omega }$ is the modulation wavelength. Then, if $d_0^{2}<\Lambda L$, so that $E_0({\boldsymbol \rho }_0) = \sqrt {I_0}\,e^{-4|{\boldsymbol \rho }_0|^{2}/d_0^{2}} \approx 0$ for $\left | {{\boldsymbol \rho }_0} \right | > \sqrt {\Lambda L}$, Fresnel diffraction reduces to Fraunhofer diffraction, yielding (58)$$\begin{aligned} \textrm{Var}[{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,2{\Omega})]&= \left(\frac{\pi I_0 T^{2}}{4L^{2}}\right)^{2}\! \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\, \int\!\frac{\textrm{d}{{\tilde{\omega}_+}}}{2\pi}\, \int\!\textrm{d}^{2}{{\boldsymbol \rho}_0}\, \int\!\textrm{d}^{2}{{\boldsymbol{\tilde{\rho}}}_0}\, e^{-8(|{\boldsymbol \rho}_0|^{2}+|{\boldsymbol{\tilde{\rho}}}_0|^{2})/d_0^{2}}\\ &\times e^{-T^{2}({\omega_+}^{2}+{\tilde{\omega}_+}^{2})/2} e^{-i \left( {\omega_+-{\tilde{\omega}_+}} \right) {\boldsymbol \rho}_1 \cdot ({\boldsymbol \rho}_0-{\boldsymbol{\tilde{\rho}}}_0) / c L} \end{aligned}$$
(59)$$ = \frac{\pi P_0^{2} T^{2}}{32 L^{4}} \sqrt{\frac{1}{1+\left( {d_0 / 2 c L T} \right)^{2}\left| {{\boldsymbol \rho}_1} \right|^{2}}}, $$
where $P_0 \equiv \pi d_0^{2}I_0/8$ is the power illuminating plane 0. Equation (59) applies to reasonably practical scenarios, e.g., $\sqrt {\Lambda L} = 22\,$cm for $\Lambda = 5\,$cm and $L=1\,$ m, and interestingly is independent of the modulation frequency. This result is maximized on axis, where the square-root term vanishes. Say $\left | {{\boldsymbol \rho }_1} \right | < 2L$, as is certainly necessary for paraxial operation. From $d_0^{2} < \Lambda L$ it then follows that (60)$$(d_0/2 c L T)^{2}\left| {{\boldsymbol \rho}_1} \right|^{2} < \;(2\pi/{\Omega} T)^{2} L/\Lambda \ll 1,$$
provided ${\Omega } T \gg 2\pi \sqrt {L/\Lambda } \approx 28.1$ for $\Lambda =5{\,\textrm {cm}}$ and $L=1\,$m. In this reasonable regime, the square-root term in Eq. (59) can be neglected entirely so that (61)$$\textrm{Var}[{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,2{\Omega})] = \pi P_0^{2} T^{2} / 32 L^{4}.$$
Using the assumptions employed thus far, but strengthening $d_0 < \sqrt {\Lambda L}$ to $d_0 \ll \sqrt {\Lambda L}$, we can show that (62)$$\begin{aligned}\left| {{{\mathcal{P}}}_1({\boldsymbol \rho}_1,2{\Omega})} \right|^{2} =& \left(\frac{\pi I_0T^{2}}{4L^{2}}\right)^{2} \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\, \int\!\frac{\textrm{d}{{\tilde{\omega}_+}}}{2\pi}\, \int\!\textrm{d}^{2}{{\boldsymbol \rho}_0}\, \int\!\textrm{d}^{2}{{\boldsymbol{\tilde{\rho}}}_0}\,\\ & \times e^{-8(|{\boldsymbol \rho}_0|^{2} + |{\boldsymbol{\tilde{\rho}}}_0|^{2})/d_0^{2}} e^{-T^{2}({\omega_+}^{2}+{\tilde{\omega}_+}^{2}) / 2} e^{-i 2{\Omega} {\boldsymbol \rho}_1 \cdot ({\boldsymbol \rho}_0-{\boldsymbol{\tilde{\rho}}}_0) / c L} \end{aligned}$$
(63)$$ = \frac{\pi P_0^{2} T^{2}}{32 L^{4}} e^{-\left( {d_0{\Omega} / 2 c L} \right)^{2}\left| {{\boldsymbol \rho}_1} \right|^{2}}. $$
Equations (61) and (63) indicate that the modulation-frequency speckle is as strong as ordinary speckle on axis and is stronger off axis, in complete agreement with Teichman’s analysis. However, we can see that in the worst case the increase is only a factor of (64)$$e^{(d_0{\Omega}/ 2 c L)^{2}|{\boldsymbol \rho}_1|^{2}} < \; e^{(2\pi d_0 / \Lambda)^{2}} \approx 1.08,$$
for $\left | {{\boldsymbol \rho }_1} \right | < 2L$, $\Lambda = 5\,\textrm {cm}$, and $d_0 = 2.2\,\textrm {mm}$. So, Teichman’s result for first-order modulated speckle from an extended diffuser is qualitatively correct, but in the paraxial regime that speckle has approximately ordinary strength. It must be emphasized, however, that our analysis for modulated speckle from extended diffusers does not extend beyond the first-order case.To illustrate our zero-frequency bound for narrowband modulation, we use
(65)$$\left| {{{\mathcal{S}}}(\omega_+)} \right|^{2} \left| {{{\mathcal{S}}}({\tilde{\omega}_+})} \right|^{2} \approx \frac{\pi^{2} T^{4}}{16} \left( {e^{-T^{2}(\omega_+ - {\Omega})^{2}/2} + e^{-T^{2}(\omega_+ + {\Omega})^{2}/2}} \right) \left( {e^{-T^{2}({\tilde{\omega}_+} - {\Omega})^{2}/2} + e^{-T^{2}({\tilde{\omega}_+} + {\Omega})^{2}/2} } \right),$$
in Eq. (53) and parallel what we just did for the modulation-frequency component. We find that (66)$$\begin{aligned} \textrm{Var}[{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,0)] = & \frac{\pi P_0^{2} T^{2}}{16 L^{4} \sqrt{1+\alpha^{2} |{\boldsymbol \rho}_1|^{2}}}\Bigg( 1 + e^{-{\Omega}^{2} T^{2} \alpha^{2} |{\boldsymbol \rho}_1|^{2} / \left( {1+\alpha^{2} |{\boldsymbol \rho}_1|^{2}} \right)} + 2 e^{-{\Omega}^{2} T^{2}}\\ & + 4 e^{-{\Omega}^{2} T^{2}\left( {2+3\alpha^{2} |{\boldsymbol \rho}_1|^{2}} \right) / 4 \left( {1+\alpha^{2} |{\boldsymbol \rho}_1|^{2}} \right)} \Bigg), \end{aligned}$$
where $\alpha \equiv d_0 / 2 c L T$. Equation (66) has its unique maximum on axis, where it takes the value (67)$$\textrm{Var}[{\tilde{\mathcal{P}}}_1({{\mathbf 0}},0)] = \frac{\pi P_0^{2} T^{2}}{8 L^{4}} \left( {1 + e^{-{\Omega}^{2} T^{2}} + 2 e^{-{\Omega}^{2} T^{2} / 2}} \right).$$
Moreover, we have that (68)$$\left| {{{\mathcal{P}}}_1({\boldsymbol \rho}_1,0)} \right|^{2} = \frac{\pi P_0^{2} T^{2}}{8 L^{4}} \left( {1 + e^{-{\Omega}^{2} T^{2}} + 2 e^{-{\Omega}^{2} T^{2} / 2}} \right),$$
which together with the variance result implies that the zero-frequency speckle from the extended diffuser has ordinary strength on axis and is attenuated off axis, in complete agreement with our bound from Eq. (54).4.2 Speckle size
As seen in the cw case for third-order speckle from extended diffusers, ${\tilde {\mathcal {P}}}_1({\boldsymbol \rho }_1,{\omega _-})$’s variance for the modulated speckle from an extended diffuser is not the sole determinant of whether those speckle fluctuations will severely limit estimating ${{\mathcal {P}}}_1({\boldsymbol \rho }_1,{\omega _-})$ from experimental data. The speckle covariance is equally important, if not more so. Deriving modulated first-order speckle’s irradiance covariance, however, is more difficult than obtaining its cw counterpart, as we shall soon see. To start, by using $E_0({\boldsymbol \rho }_0,t)$’s being space-time factorable, Fresnel-diffraction integrals, and Gaussian moment factoring we can obtain
(69)$$\begin{aligned}& \textrm{Covar}[{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,{\omega_-}),{\tilde{\mathcal{P}}}_1({\boldsymbol{\tilde{\rho}}}_1,{\omega_-})]= \langle {{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,{\omega_-}){\tilde{\mathcal{P}}}^{*}_1({\boldsymbol{\tilde{\rho}}}_1,{\omega_-})} \rangle_0-{{\mathcal{P}}}_1({\boldsymbol \rho}_1,{\omega_-}) {{\mathcal{P}}}^{*}_1({\boldsymbol{\tilde{\rho}}}_1,{\omega_-})\\ & = \frac{I_0^{2}}{L^{4}} \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\,\int\!\frac{\textrm{d}{{\tilde{\omega}_+}}}{2\pi}\,\int\!\textrm{d}^{2}{{\boldsymbol \rho}_0}\,\int\!\textrm{d}^{2}{{\boldsymbol{\tilde{\rho}}}_0}\, e^{-8(|{\boldsymbol \rho}_0|^{2}+|{\boldsymbol{\tilde{\rho}}}_0|^{2})/d_0^{2}}\\ & \times \mathcal{S}(\omega_++{\omega_-}/2)\mathcal{S}^{*}({\tilde{\omega}_+}+{\omega_-}/2)\mathcal{S}^{*}(\omega_+-{\omega_-}/2)\mathcal{S}({\tilde{\omega}_+}-{\omega_-}/2)\\ & \times e^{i [ \left( {{\omega_0}+\omega_++{\omega_-}/2} \right)\left| {{\boldsymbol \rho}_1-{\boldsymbol \rho}_0} \right|^{2} - \left( {{\omega_0}+{\tilde{\omega}_+}+{\omega_-}/2} \right)\left| {{\boldsymbol{\tilde{\rho}}}_1-{\boldsymbol \rho}_0} \right|^{2} - \left( {{\omega_0}+\omega_+-{\omega_-}/2} \right)\left| {{\boldsymbol \rho}_1-{\boldsymbol{\tilde{\rho}}}_0} \right|^{2} + \left( {{\omega_0}+{\tilde{\omega}_+}-{\omega_-}/2} \right)\left| {{\boldsymbol{\tilde{\rho}}}_1-{\boldsymbol{\tilde{\rho}}}_0} \right|^{2} ] / 2 c L}. \end{aligned}$$
The differing spatial-coordinate combinations in the four Fresnel-diffraction exponents makes it impossible to group them in any useful way, regardless of their frequency coefficients. Moreover, expanding these exponents’ squares (not shown) does not provide a route to simplification. Furthermore, deleting the exponentials altogether to obtain an upper bound, as we did for $\textrm {Var}[{\tilde {\mathcal {P}}}_1({\boldsymbol \rho }_1,0)]$, is not an option because it removes all spatial dependence, i.e., it suppresses the very behavior we are seeking.Despite the preceding difficulties, useful insight into $\textrm {Covar}[{\tilde {\mathcal {P}}}_1({\boldsymbol \rho }_1,2\Omega ),{\tilde {\mathcal {P}}}_1({\boldsymbol {\tilde {\rho }}}_1,2\Omega )]$ can be obtained in the narrowband-modulation, $d_0^{2} < \Lambda L$ case considered earlier. Paralleling the work done there we get
(70)$$\begin{aligned} &\textrm{Covar}[{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,2{\Omega}),{\tilde{\mathcal{P}}}_1({\boldsymbol{\tilde{\rho}}}_1,2{\Omega})] = \left(\frac{\pi I_0T^{2}}{4L^{2}}\right)^{2}\! e^{i{\Omega}\left( {\left| {{\boldsymbol \rho}_1} \right|^{2}-\left| {{\boldsymbol{\tilde{\rho}}}_1} \right|^{2}} \right)/ c L}\\ & \quad \times \int\!\frac{\textrm{d}{\omega_+}}{2\pi}\,\int\!\frac{\textrm{d}{{\tilde{\omega}_+}}}{2\pi}\,\int\!\textrm{d}^{2}{{\boldsymbol \rho}_0}\,\int\!\textrm{d}^{2}{{\boldsymbol{\tilde{\rho}}}_0}\, e^{-8(|{\boldsymbol \rho}_0|^{2}+|{\boldsymbol{\tilde{\rho}}}_0|^{2})/d_0^{2}}e^{-T^{2}\left( {{\omega_+}^{2}+{\tilde{\omega}_+}^{2}} \right)/2}\\ & \quad \times e^{i [-{\boldsymbol \rho}_0\cdot{\boldsymbol \rho}_1\left( {{\omega_0}+\omega_++{\Omega}} \right)+{\boldsymbol \rho}_0\cdot{\boldsymbol{\tilde{\rho}}}_1\left( {{\omega_0}+{\tilde{\omega}_+}+{\Omega}} \right)+{\boldsymbol{\tilde{\rho}}}_0\cdot{\boldsymbol \rho}_1\left( {{\omega_0}+\omega_+-{\Omega}} \right)-{\boldsymbol{\tilde{\rho}}}_0\cdot{\boldsymbol{\tilde{\rho}}}_1\left( {{\omega_0}+{\tilde{\omega}_+}-{\Omega}} \right)]/ c L}. \end{aligned}$$
This 6D integral can be evaluated in closed form, resulting in (71)$$\begin{aligned} \textrm{Covar}[{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,2\Omega),&{\tilde{\mathcal{P}}}_1({\boldsymbol{\tilde{\rho}}}_1,2\Omega)] = \frac{\pi P_0^{2}T^{2}}{32 L^{4}}\frac{e^{i\Omega(|{\boldsymbol \rho}_1|^{2}-|{\boldsymbol{\tilde{\rho}}}_1|^{2})/cL}e^{-\alpha^{2}\Omega^{2}T^{2}|{\boldsymbol \rho}_-|^{2}/4}}{\sqrt{(1+ \alpha^{2}|{\boldsymbol \rho}_+|^{2})(1 + \alpha^{2}|{\boldsymbol \rho}_-|^{2}/4) - \alpha^{4}({\boldsymbol \rho}_+\cdot{\boldsymbol \rho}_-)^{2}/4}}\\ &\times \exp\!\left\{-\frac{{\omega_0}^{2}T^{2}}{4}\frac{\alpha^{2}|{\boldsymbol \rho}_-|^{2}+\alpha^{4}\left[ {|{\boldsymbol \rho}_+|^{2}|{\boldsymbol \rho}_-|^{2}-({\boldsymbol \rho}_+\cdot{\boldsymbol \rho}_-)^{2}} \right]}{(1+ \alpha^{2}|{\boldsymbol \rho}_+|^{2})(1 + \alpha^{2}|{\boldsymbol \rho}_-|^{2}/4) - \alpha^{4}({\boldsymbol \rho}_+\cdot{\boldsymbol \rho}_-)^{2}/4}\right\}, \end{aligned}$$
where, as before, ${\boldsymbol \rho }_+ \equiv ({\boldsymbol \rho }_1+{\boldsymbol {\tilde {\rho }}}_1)/2$, ${\boldsymbol \rho }_- \equiv {\boldsymbol \rho }_1-{\boldsymbol {\tilde {\rho }}}_1$, and $\alpha \equiv d_0/2cLT$.Despite its formidable length, Eq. (71) readily yields important insights. First, we have that
(72)$$\lim_{T\rightarrow \infty}\frac{|\textrm{Covar}[{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,2\Omega),{\tilde{\mathcal{P}}}_1({\boldsymbol{\tilde{\rho}}}_1,2\Omega)]|} {\sqrt{\textrm{Var}[{\tilde{\mathcal{P}}}_1({\boldsymbol \rho}_1,2\Omega)]\textrm{Var}[{\tilde{\mathcal{P}}}_1({\boldsymbol{\tilde{\rho}}}_1,2\Omega)]}} = e^{-({\omega_0}^{2}+\Omega^{2})d_0^{2}|{\boldsymbol \rho}_-|^{2}/16c^{2}L^{2}}.$$
Second, comparing with Eq. (32)—and remembering that $\Omega \ll {\omega _0}$—shows that modulated first-order speckle has approximately the same speckle size as cw first-order speckle in the $T\rightarrow \infty$ limit. Moreover, because $\max \!\left ( {|{\boldsymbol \rho }_+|^{2}} \right ) = d_1^{2}/4$, $\max \!\left ( {|{\boldsymbol \rho }_-|^{2}} \right ) = d_1^{2}$, $\max \!\left [ {|{\boldsymbol \rho }_+|^{2}|{\boldsymbol \rho }_-|^{2} - ({\boldsymbol \rho }_+\cdot {\boldsymbol \rho }_-)^{2}} \right ] = d_1^{4}/16$, and $\max \!\left [ {\left ( {{\boldsymbol \rho }_+\cdot {\boldsymbol \rho }_-} \right )^{2}} \right ] = d_1^{4}/64$, this speckle-size near coincidence holds for $T$ satisfying $d_0^{2}d_1^{2}/16c^{2}L^{2}T^{2} \ll 1$. Taking $d_0 = 3\,$mm, $d_1 =10$ cm, and $L = 1$ m we have that $d_0^{2}d_1^{2}/16c^{2}L^{2}T^{2} \le 0.0625$ for $T \ge 1$ ps, which is much smaller than the minimum $T$ value for validating the narrowband ($\Omega T \gg 1$) quasimonochromatic ($\Omega \ll {\omega _0}$) modulation assumed in this covariance analysis. Furthermore, as we have already shown, the first-order speckle strength of the modulated case is close to that of the cw case. Thus, should these behaviors hold for second-order and third-order speckle, the situation would be quite favorable for speckle-suppression in ${{\mathcal {P}}}$-field NLoS imaging, i.e., it would mean that the speckled speckled speckle from extended diffuse reflectors would reduce to the speckle produced by the last reflection and that final speckle could easily be averaged out in power collection over realistically-sized detectors.5. Summary and discussion
Increasing interest in NLoS imaging with coherent illumination is driving the need to understand high-order speckle, principally the speckled speckled speckle that arises in three-bounce NLoS imaging. This need is especially pressing for ${{\mathcal {P}}}$-field NLoS imaging because initial ${{\mathcal {P}}}$-field experiments have afforded some of the best NLoS coherent-illumination imagery obtained to date [17]. Inasmuch as ${{\mathcal {P}}}$-field imaging relies on modulated coherent illumination, the understanding to be sought should encompass modulated speckled speckled speckle.
In this paper, we took first steps in addressing these issues using space-time factorable initial illumination in a three-diffuser transmissive geometry that is a proxy for three-bounce NLoS imaging. In the small-diffusers limit we showed that the irradiance variances of cw and modulated $n$th-order speckle coincide and are $(2^{n}-1)$-times those of ordinary (first-order) speckle. If not mitigated, this result implies that the maximum SNR in three-bounce NLoS imaging, i.e., its $\textrm {SNR}_\textrm {sat}$, would be 1/7. The more important case for NLoS imaging, however, involves extended diffuse reflectors. For our transmissive geometry with extended diffusers, we treated third-order cw speckle and first-order modulated speckle. It turned out that speckle is unlikely to impede successful operation of coherent-illumination cw imagers because typical parameter values for NLoS scenarios reduce cw third-order speckle in our transmissive-geometry proxy to the ordinary (first-order) speckle produced by the last diffuser. In addition, speckle averaging in optical power collection over typical detector sizes suppresses that residual first-order speckle to the point that its impact on the cw imager’s SNR is quite benign, i.e., $\textrm {SNR}_\textrm {sat} \gg 1$ is achieved. More importantly, insofar as ${{\mathcal {P}}}$-field NLoS imaging is concerned, our analysis of modulated first-order speckle revealed that its speckle strength and speckle size were very similar to those of the cw case. So, should the same correspondence apply to modulated third-order speckle, then ${{\mathcal {P}}}$-field NLoS imagers would be largely immune to the adverse effects of speckle.
It remains to work out the behavior—in speckle strength and speckle size—of modulated third-order speckle for extended diffusers. Based on what we have accomplished, that important task appears to be quite formidable. A second remaining task of significance is the design and execution of experiments that establish the extent to which speckle effects are actually discernible in coherent-illumination NLoS imaging.
Funding
Defense Advanced Research Projects Agency (HR0011-16-C-0030).
Acknowledgments
The authors acknowledge fruitful interactions with members of the DARPA REVEAL teams from the University of Wisconsin and Southern Methodist University. They also thank Dr. Ravi Athale for organizing a valuable workshop on phasor-field imaging, and Dr. Jeremy Teichman for sharing an early version of his analysis of speckle effects in ${{\mathcal {P}}}$-field imaging.
Disclosures
The authors declare no conflicts of interest.
References
1. J. W. Goodman, Speckle Phenomena in Optics, (Roberts & Company, 2007).
2. S. A. Reza, M. La Manna, and A. Velten, “A physical light transport model for non-line-of-sight imaging applications,” arXiv:1802.1823 [physics.optics].
3. J. Dove and J. H. Shapiro, “Paraxial theory of phasor-field imaging,” Opt. Express 27(13), 18016–18037 (2019). [CrossRef]
4. J. A. Teichman, “Phasor field waves: a mathematical treatment,” Opt. Express 27(20), 27500–27506 (2019). [CrossRef]
5. J. Dove and J. H. Shapiro, “Paraxial phasor-field physical optics,” Opt. Express 28(14), 21095–21109 (2020). [CrossRef]
6. J. Dove, “Theory of phasor-field imaging,” Ph.D. thesis, Massachusetts Institute of Technology (2020).
7. The STA irradiance $I_z(\boldsymbol {\rho}_z, t)$ is the power density, averaged over the time interval [t − Ta, t] where 2π/ω0 ≪ Ta ≪ 2π/Δω, that illuminates the transverse coordinates $\boldsymbol {\rho}_z $ in the z-plane.
8. For cw illumination Δλ is infinite, hence in that case we need only assume that the diffuser is rough at the optical wavelength.
9. In the NLoS scenario, both $h_0(\boldsymbol {\rho}_0) $ and $h_2(\boldsymbol {\rho}_2) $ represent the visible wall. Nevertheless, they can be taken to be statistically independent in our transmissive-geometry proxy if the NLoS imager’s visible-wall illumination falls on a different portion of that wall than what lies in its camera’s field of view.
10. Interestingly, having λ0 ≪ σh ≪ Δλ and ρh ∼ λ0 results in the first and second moments of ${E_0}^\prime(\boldsymbol {\rho}_0) $ being independent of both σh and ρh.
11. That the diffuser-averaged STA irradiance is independent of the transverse coordinate is an artifact of Fresnel diffraction and total spatial incoherence produced by our delta-function approximation to $\langle e^{i\omega _0[h_0({\boldsymbol \rho }_0)-h_0(\tilde {{\boldsymbol \rho }}_0)]/c}\rangle$. Note that 〈I1〉 < 〈I0〉 holds because d0 ≪ L is required for Fresnel diffraction to be valid.
12. E. W. Weisstein, “Meijer G-function,” http://mathworld.wolfram.com/MeijerG-Function.html.
13. J. H. Shapiro, B. A. Capron, and R. C. Harney, “Imaging and target detection with a heterodyne-reception optical radar,” Appl. Opt. 20(19), 3292–3313 (1981). [CrossRef]
14. J. H. Shapiro, “The quantum theory of optical communications,” IEEE J. Sel. Top. Quantum Electron. 15(6), 1547–1569 (2009). [CrossRef]
15. We used this condition in our analysis of synthetic-wavelength holography [5].
16. F. Willomitzer, P. V. Rangarajan, F. Li, M. M. Balaji, M. P. Christensen, and O. Cossairt, “Synthetic wavelength holography: an extension of Gabor’s holographic principle to imaging with scattered wavefronts,” arXiv:1912.11438 [physics.optics].
17. X. Liu, I. Guillén, M. La Manna, J. H. Nam, S. A. Reza, T. H. Le, A. Jarabo, D. Gutierrez, and A. Velten, “Non-line-of-sight imaging using phasor-field virtual wave optics,” Nature 572(7771), 620–623 (2019). [CrossRef]
Open Access
Add to BibSonomy
