Distribution of the two-point product of complex amplitudes in the fully saturated scattering regime Open Access

This Letter considers probability distributions for products of complex amplitudes (phasors) of harmonic signals observed at two points in the fully saturated scattering regime. Although the two points most typically are different spatial locations, they may also involve separations in time or frequency. The main significance of the complex-amplitude product is that it corresponds to cross-spectral estimates from time-windowed Fourier transforms (e.g., by Welch's method as discussed in Bendat and Piersol, 1986). As such, this is important for assessing decorrelation of signals across an array of sensors, and the resulting degradations of acoustic beamformer performance (Collier and Wilson, 2003). Multipoint statistical models also enable calculation of acoustic probabilities of detection and classification involving fusion of data across multiple sensors, and the gains achievable by coherent (rather than incoherent) processing between the sensors (Svensson and Lundberg, 2005; Wilson et al., 2020). Finally, the distributions derived here extend our understanding of the fully saturated scattering regime, and for recognizing when the modeling assumptions of this regime are satisfied for sound scattering in the ocean (e.g., Dashen, 1979; Colosi, 2016) and atmosphere (Ostashev and Wilson, 2015).

By complex amplitude, we mean $z=x+iy=Aeiϕ$⁠, where $x$ is the real (observed) part of the signal, $y$ is the quadrature component, $A$ is the amplitude, and $ϕ$ is the phase, such that the actual observed signal is $Re[zeiωt]=A cos (ωt+ϕ)$⁠, where $ω=2πf$⁠, $f$ is the frequency, and $t$ is time. In the quasi-static approximation, $x$⁠, $y$⁠, $A$⁠, and $ϕ$ may slowly vary in time, relative to the period of the wave $T=1/f$⁠. The complex amplitude may represent a time-windowed Fourier transform at a particular frequency, or a signal after application of a narrowband filter and demodulation. By fully saturated regime, we mean that the deterministic component of the signal intensity is negligible compared to the scattered intensity, so that scintillation index attains its saturation value of one, as occurs for relatively strong turbulent fluctuations and long distances from the source (Rytov et al., 1989; Andrews and Phillips, 2005).

Consider a pair of complex amplitudes observed at two different points, designated as $z1=x1+iy1$ and $z2=x2+iy2$⁠. In the fully saturated regime, $x1$ and $y1$ are independent, normally distributed random variables (rvs) with zero mean and variance $σ12$⁠. Likewise, $x2$ and $y2$ are independent and normally distributed with zero mean and variance $σ22$⁠. We allow for the possibility that $σ12≠σ22$ since the observation points may differ. Given these assumptions, it is well known that random samples of the signal power (or time-windowed autospectral estimates), defined as $ζi=|zi|2=xi2+yi2$⁠, have an exponential probability density function (pdf) $p(ζi)=exp (−ζi/mi)/mi$⁠, where $mi=2σi2$ is the mean (Dashen, 1979; Burdic, 1991; Colosi, 2016). However, pdfs involving products of the complex signals at the two points are less well known.

In this Letter, we first derive the pdfs of the real and imaginary parts of the complex amplitude product $ζ12=z1z2∗$⁠, namely $p(Re ζ12)$ and $p(Im ζ12)$⁠, which were apparently unknown previously. The derivation leverages a recent result from Nadarajah and Pogány (2016) for the pdf of the product of zero-mean normal rvs. Second, we derive representative joint pdfs involving $ζ1$⁠, $ζ2$⁠, $Re ζ12$⁠, and $Im ζ12$⁠, namely, $p(ζ1,Re ζ12,ζ2)$⁠, $p(ζ1,Re ζ12,Im ζ12)$⁠, $p(ζ1,ζ2)$⁠, and $p(Re ζ12,Im ζ12)$⁠. The pdf $p(Re ζ12,Im ζ12)$ relates to coherence of the two signals; a larger real part indicates greater signal coherence. The pdf $p(ζ1,ζ2)$ is important for incoherent processing, and was previously derived by Churnside (1989). Here we show that Churnside's equation can be derived efficiently based on the same general formulation used to derive the other pdfs.

In this section, we derive the pdfs for the real and imaginary parts of $ζ12=z1z2∗$ for fully saturated signals. By definition,

Hence, $Re ζ12=x1x2+y1y2$ and $Im ζ12=y1x2−x1y2$⁠. As mentioned already, the $xi$ and $yi$ are distributed normally with zero mean, and $⟨xi2⟩=⟨yi2⟩=σi2$⁠. The cross-correlations of the real parts between the sensors (⁠$x1$ and $x2$⁠), and of the imaginary parts between the sensors (⁠$y1$ and $y2$⁠), are generally non-zero, since the phase variations at the two sensors may be coupled. However, all cross-correlations between real and imaginary parts, at each sensor, and between the sensors, must be zero, as a consequence of the phase being uniformly distributed in the full saturation regime. Mathematically, $⟨x1x2⟩=⟨y1y2⟩=ρσ1σ2$⁠, and $⟨xiyi⟩=0$ for all $i$ and $j$⁠, where $ρ$ is a real-valued correlation coefficient. The value of $ρ$ will depend on the spatial separation between the observation locations (being one when the separation is zero, and decreasing with increasing separation), the spectrum of the turbulence, and the propagation distance.

In the following, we will consider distributions for the sum of $d$ independent, identically distributed (i.i.d.) samples of $ζ12$⁠. Let us first examine just the product $x1x2$⁠. The distribution of this product can be calculated using a result from Nadarajah and Pogány (2016), who consider the mean $ζ¯$ of $d$ i.i.d. samples of the product $n1n2$⁠, where $n1$ and $n2$ are zero-mean normal rvs with unit variance. Using the method of characteristic functions, they show that $ζ¯$ has the following pdf (their Theorem 2.2):

where $β′=d/(1−ρ)$⁠, $γ′=d/(1+ρ)$⁠, $Γ(⋅)$ is the gamma function, and $Kv(⋅)$ is the modified Bessel function of the second kind of order $v$⁠. For present purposes, we wish to find the sum (rather than the mean) of $d$ i.i.d. samples of $x1x2$⁠, for the case when the variances do not necessarily equal one. This sum will be described by the same distribution as Eq. (2), after transforming the random variable (rv) to $ζ=dσ1σ2ζ¯$⁠. With some algebra, we arrive at the following form:

where $α=1/[σ1σ2(1−ρ2)]$⁠, $β=ρ/[σ1σ2(1−ρ2)]$⁠, and $λ=d/2$⁠. This is the customary form for the central variance-gamma distribution (Madan and Seneta 1990). We indicate a rv drawn from this distribution using the notation $ζ∼VG(α,β,λ)$⁠.

The sum of $d$ i.i.d. samples of the product $y1y2$⁠, likewise, has a variance-gamma pdf with the same parameters. Furthermore, since $x1x2$ and $y1y2$ are independent, the sum of these products also has a variance gamma pdf, but with twice the degrees of freedom. Hence the real part of $ζ12$ has the following pdf:

Calculation of the distribution of the imaginary part is somewhat more complicated because the terms $y1x2$ and $x1y2$ are correlated. To address this problem, let us write $x1=σ1n1$ and $y1=σ1n3$⁠, where $ni$ is a unit-variance, zero-mean normal rv. Hence, $⟨x12⟩=⟨y12⟩=σ12$ as required. Furthermore, suppose that $x2$ is a linear combination of $x1$ plus a term proportional to another unit-variance, zero-mean rv, $n2$⁠. That is, $x2=ax1+bn2$⁠, where $a$ and $b$ are constants to be determined. Based on the constraint $⟨x1x2⟩=ρσ1σ2$⁠, we readily find $a=ρσ2/σ1$⁠. Then, since $⟨x22⟩=σ22$⁠, we find $b=σ21−ρ2$⁠. Hence, $x2=σ2(ρn1+1−ρ2n2)$ and, similarly, $y2=σ2(ρn3+1−ρ2n4)$⁠. With some algebra, we find

Note that the pdf of the term $−n1n4$ is not affected by changing the minus sign to a plus, since the pdf of $n4$ is symmetric about the origin. The imaginary part thus reduces to the sum of two independent contributions. Each of these contributions is a product of uncorrelated normal rvs with variance $σ1σ21−ρ2$⁠. This amounts to setting $ρ=0$ in Eq. (4), and then replacing $σ1σ2$ with $σ1σ21−ρ2$⁠. Hence,

Equations (4) and (6) provide the desired distributions for $p(Re ζ12)$ and $p(Im ζ12)$⁠. They are plotted in Fig. 1 for $σ12=σ22=1$⁠, and $d=1$⁠. Four values of $ρ$ are shown: $0.1$⁠, $0.5$⁠, $0.9$⁠, and $1$⁠. For small $ρ$⁠, $Re ζ12$ and $Im ζ12$ have nearly the same pdf, which is approximately exponential but symmetric about the origin (known as a Laplace pdf). This occurs because $x1$⁠, $x2$⁠, $y1$⁠, and $y2$ are negative or positive with equal probability. For $ρ=1$⁠, $x1=x2$ and $y1=y2$⁠, so that the pdf of $Re ζ12$ is exponential as for the power at a single point, whereas the imaginary part is zero.

This section considers joint pdfs involving $ζ1$⁠, $ζ2$⁠, $Re ζ12$⁠, and $Im ζ12$⁠. We begin by recalling that $x1=σ1n1$⁠, $y1=σ1n3$⁠, $x2=σ2(ρn1+1−ρ2n2)$⁠, and $y2=σ2(ρn3+1−ρ2n4)$⁠, where the $ni$ are independent normal rvs with zero mean and unit variance. Changing to cylindrical coordinates $r$⁠, $θ$⁠, $s$⁠, and $ϕ$⁠, such that $n1=r cos θ$⁠, $n2=s cos ϕ$⁠, $n3=r sin θ$⁠, $n4=s sin ϕ$⁠, we have the following:

and

For brevity, let us define the normalized variables $t=ζ1/σ12$⁠, $u=Re ζ12/(σ1σ2)$⁠, $v=Im ζ12/(σ1σ2)$⁠, and $w=ζ2/σ22$⁠. Note from Eqs. (7)–(10) that $t$⁠, $u$⁠, $v$⁠, and $w$ depend only on three independent latent rvs, namely $r$⁠, $s$⁠, and the difference $θ−ϕ$⁠. Since $r2=n12+n32$ and $s2=n22+n42$⁠, these are both exponentially distributed with mean equal to $2$⁠, whereas $θ−ϕ$ is uniformly distributed between $−π$ and $π$⁠. Hence the joint pdf of $η=r2$⁠, $μ=s2$⁠, and $τ=θ−ϕ$ is given by

Joint pdfs involving $t$⁠, $u$⁠, $v$⁠, and $w$ follow by transforming Eq. (11). However, non-singular pdfs can be formulated for only three of these four quantities at a time, due to the aforementioned dependence on three underlying rvs. Let us first consider $p(t,u,v)$⁠, which is of particular interest since it enables calculation $p(u,v)$ by marginalizing over $t$⁠. By the usual formulation for transforming multivariate pdfs [e.g., Stark and Woods (1986), Eqs. (2.11–13a)], we have $p(t,u,v)=[p(η,μ,τ)/|Jt,u,v(η,μ,τ)|]|t,u,v$⁠, where the Jacobian is

Evaluation of the derivatives leads to $|Jt,u,v(η,μ,τ)|=(η/2)(1−ρ2)=(t/2)(1−ρ2)$⁠. Furthermore, it can be shown that $η+μ=(t2+u2+v2−2ρtu)/[t(1−ρ2)]$⁠. Thus, Eq. (11) transforms to

Next, we can integrate (marginalize) over $t$ to find the joint pdf of the real and imaginary parts of $ζ12$⁠. Based on integral (3.471.9) from Gradshteyn and Ryzhik (1980), we find

This equation makes clear that $u$ and $v$ are not independent, since $p(u,v)$ cannot be factored into $p(u)p(v)$⁠. The unnormalized form of the distribution is found by setting $Re ζ12=σ1σ2u$ and $Im ζ12=σ1σ2v$⁠, and dividing by the determinant of the Jacobian (in this case equal to $σ12σ22$⁠), yielding

It can be verified that Eq. (4) with d = 1 follows by integrating over $Im ζ12$⁠, whereas Eq. (6) with d = 1 follows by integrating over $Re ζ12$⁠. Thus, we confirm that the approach in this section leads to results consistent with Nadarajah and Pogány's (2016) Theorem 2.2.

Figure 2 shows $p(Re ζ12,Im ζ12)$ for the case $σ12=σ22=1$ [i.e., $p(u,v)$], with $ρ=0.25$ and $ρ=0.75$⁠. The pdf peaks at $ζ12=0$⁠; as $ρ$ increases, the pdf skews towards positive values of $Re ζ12$⁠. This behavior indicates that the signals at the two observation points are coming into phase.

Let us next consider the joint pdf $p(t,u,w)$⁠. The approach to transforming the pdf $p(η,μ,τ)$ to $p(t,u,w)$ is very similar to the previous derivation of $p(t,u,v)$⁠. However, we must account for the fact that $p(t,u,w)$ is a multi-valued function of $τ$⁠; specifically, since it only depends on $cos τ$ [and not also on $sin τ$ as was the case for $p(t,u,v)$], $−τ$ and $+τ$ both map to the same values of $t$⁠, $u$⁠, and $w$⁠. Thus, we set $p(t,u,w)=2[p(η,μ,τ)/|Jt,u,w(η,μ,τ)|]|t,u,w$ [e.g., Stark and Woods (1986), Eqs. (2.11–15)]. By calculating the derivatives, after some algebra, we arrive at $|Jt,u,w(η,μ,τ)|=(1−ρ2)3/2ημ | sin τ|=(1−ρ2)tw−u2$⁠. It can also be shown that $η+μ=(t+w−2ρu)/(1−ρ2)$⁠. Thus, we have

It can be readily shown from Eqs. (7), (8), and (10) that $tw−u2=(1−ρ2)ημ sin2τ$⁠. Since this quantity is always non-negative, $tw≥u2$⁠, and $u/tw$ is bound between $−1$ and $+1$⁠. This observation enables us to perform the marginalization over $u$ by making the substitution $x=u/tw$⁠, and applying integral (3.387.1) from Gradshteyn and Ryzhik (1980), with result

where $I0(⋅)$ is the modified Bessel function of the first kind of order $0$⁠. Hence, in unnormalized form,

This result agrees with Eq. (19) in Churnside (1989) with the identifications $I1=ζ1$⁠, $I2=ζ2$⁠, $z1=2σ12$⁠, and $z2=2σ22$⁠.

Figure 3 shows $p(ζ1,ζ2)$ for the case $σ12=σ22=1$ [i.e., $p(t,w)$], with $ρ=0.25$ and $ρ=0.75$⁠. The pdf is maximized when $ζ1=ζ2=0$⁠, and decreases gradually as either $ζ1$ or $ζ2$ is increased.

In this Letter, pdfs were derived involving products of complex amplitudes observed at two points, under the assumption that the signal scattering is fully saturated. The real and imaginary parts of the product of the complex amplitude at one point with the conjugate of the complex amplitude at another point were each shown to have variance gamma pdfs, although with somewhat different parameters. Several joint pdfs involving the complex-amplitude product and the signal powers at the two observation points were also derived. The joint pdf of the powers was shown to coincide with a result previously given by Churnside (1989). Taken together, the results here extend our understanding of two-point signal statistics in the fully saturated scattering regime.

This work was supported by the U.S. Army Engineer Research and Development Center, Geospatial Research Engineering basic research program. Permission to publish was granted by the Director, Cold Regions Research and Engineering Laboratory. Any opinions expressed in this paper are those of the authors, and are not to be construed as official positions of the funding agency or the Department of the Army unless so designated by other authorized documents.