Random vectors, covariance & correlation
Second-order statistics rebuilt from the expectation operator up: covariance, the correlation coefficient and how badly it is estimated from few samples, the covariance matrix of a complex IQ vector, correlation functions in time, and finally the channel covariance matrix across an array. That last object, through the van Cittert–Zernike theorem, gives the spatial coherence curve that sound-speed estimators maximize and that aberration flattens.
You will be able to
- Manipulate E[·], var, cov as linear and bilinear operators, and say when a spatial or temporal average may replace the ensemble average
- Read a 2×2 covariance matrix as an ellipse, and explain why the two regression lines and the principal axis differ
- Quote the sampling spread of \(\hat{\rho }\) (Fisher z) and recognise spurious correlation at small or effective N
- Use the covariance matrix of a complex circular Gaussian vector, its rank when snapshots are scarce, and diagonal loading
- Disambiguate the six things called "correlation" or "coherence" in the ultrasound literature
- Derive the VCZ triangle for a focused uniform aperture and predict how noise, aberration and a wrong sound speed change it
07.1Random variables, vectors and the expectation operator
A random variable x is a function that assigns a number to each outcome of an experiment. In ultrasound the experiment is "draw a tissue realization": the positions and reflectivities of the sub-resolution scatterers. One RF sample at a fixed depth is a random variable; the 64 IQ samples across the channels at one depth form a random vector \(\mathbf{x} \in \mathbb{C}^{64}\).
Everything in this module is built on one operator, the expectation (ensemble average over realizations):
- \(x\)
- random variable (e.g. one RF sample at a fixed depth)
- \(p(x)\)
- probability density of x
- \(\mu\)
- mean of x
- \(g\)
- any function of x
- \(\mathrm{E}[\cdot]\)
- expectation (ensemble average) of the bracketed quantity
- \(\int\)
- integral over all values of x
- \(dx\)
- integration element
- \(x,\ y\)
- random variables
- \(a,\ b,\ c\)
- fixed scalar constants
- \(\mathbf{x}\)
- random vector
- \(\mathbf{A},\ \mathbf{b}\)
- fixed matrix and fixed vector
- \(\pmb{\mu}\)
- mean vector, \(\mathrm{E}[\mathbf{x}]\)
- \(\mathrm{E}[\cdot]\)
- expectation, as in (07.1)
Variance is the expected squared deviation, and it is not linear:
- \(x,\ y\)
- random variables (real or complex)
- \(\operatorname{var}(\cdot)\)
- variance
- \(\sigma_x^2\)
- variance of x
- \(\mu\)
- mean of x, (07.1)
- \(a,\ b\)
- fixed constants
- \(\operatorname{cov}(x,y)\)
- covariance, (07.5)
- \(\operatorname{Re}\)
- real part
- \(\lvert\cdot\rvert\)
- magnitude
- \(\mathrm{E}[\cdot]\)
- expectation, (07.1)
Ensemble versus sample averages; stationarity and ergodicity
You never observe the ensemble. You observe one tissue realization and average over space (pixels in an ROI), time (slow-time frames) or channels. That replacement is legitimate only under two conditions:
- Stationarity (wide-sense, WSS): the mean is constant and the covariance between two samples depends only on their separation, \(C(t_{1}, t_{2}) = C(t_{1} - t_{2})\). For speckle this means a statistically homogeneous region: same scatterer density, same attenuation-compensated level, same PSF.
- Ergodicity (in the mean or in correlation): the sample average along one realization converges to the ensemble average. For the mean, a sufficient condition is that the autocovariance decays fast enough to be summable, \(\sum _{k} |C(k)| < \infty\).
The price of correlated samples is quantified by the variance of the sample mean. For N samples of a WSS sequence with correlation coefficient \(\rho (k)\) between samples k apart:
- \(\bar{x}\)
- sample mean of N samples of a WSS sequence
- \(\sigma^2\)
- variance of one sample
- \(\rho(k)\)
- correlation coefficient between samples k apart
- \(N_{\mathrm{eff}}\)
- effective number of independent samples (\(\approx N/\ell\) for a correlation length of \(\ell\) samples)
- \(\operatorname{var}(\cdot)\)
- variance of the sample mean
- \(N\)
- number of samples averaged
- \(k\)
- lag between samples, in samples
- \(\sum_{k}\)
- sum over all lags from \(-(N-1)\) to \(N-1\)
- \(\equiv\)
- defines \(N_{\mathrm{eff}}\)
Two different "repeats" are often confused. Re-acquiring a static phantom at the same probe position gives the same speckle realization: frame averaging removes electronic noise but not speckle. Independent realizations come from different scatterer draws: moving the probe by more than a resolution cell, a different phantom region, elevational translation, or tissue motion. A spatial average over an ROI stands in for the ensemble only if the ROI is homogeneous (stationary) and spans many speckle cells. This is the assumption behind every QUS estimate (backscatter coefficient, attenuation, envelope statistics) and behind the spatial-coherence estimates in section 07.6.
07.2Covariance and the correlation coefficient
Covariance measures how two random variables move together after their means are removed:
- \(\mu_x,\ \mu_y\)
- means of x and y
- \({}^*\)
- complex conjugate (IQ data are complex)
- \(\operatorname{cov}(x,y)\)
- covariance of x and y (units of x times y)
- \(x,\ y\)
- random variables
- \(\mathrm{E}[\cdot]\)
- expectation, (07.1)
The correlation coefficient removes the scale:
- \(\sigma_x,\ \sigma_y\)
- standard deviations
- \(\rho_{xy}\)
- correlation coefficient, dimensionless
- \(\operatorname{cov}(x,y)\)
- covariance, (07.5)
- \(\lvert\cdot\rvert\)
- magnitude
- \(\le\)
- bound set by Cauchy–Schwarz, (07.7)
The bound is Cauchy–Schwarz
Zero-mean random variables with finite variance form a vector space with inner product \(\langle x, y\rangle = \mathrm{E}[x y^{*}]\). The norm is \(\| x\| = \sqrt{\mathrm{E}[|x|^{2}]} = \sigma _{x}\). Cauchy–Schwarz from Module 05 then gives the bound directly:
- \(\langle\cdot,\cdot\rangle\)
- inner product of zero-mean variables, \(\mathrm{E}[\tilde{x}\tilde{y}^{*}]\)
- \(\lVert\cdot\rVert\)
- norm, \(\sqrt{\mathrm{E}\lvert\tilde{\cdot}\rvert^{2}}\) (the standard deviation)
- \(\tilde{x},\ \tilde{y}\)
- mean-removed variables
- \(\mu_x,\ \mu_y\)
- means, (07.5)
- \(\sigma_x,\ \sigma_y\)
- standard deviations, (07.6)
- \(\rho\)
- correlation coefficient, (07.6)
- \(\angle\)
- angle between two zero-mean random variables viewed as vectors
- \(\iff\)
- if and only if
- \(a\)
- fixed scalar
- \(\operatorname{cov}\)
- covariance, (07.5)
Uncorrelated is not independent
Independence means \(p(x, y) = p(x)p(y)\), which forces \(\mathrm{E}[f(x)g(y)] = \mathrm{E}[f(x)]\mathrm{E}[g(y)]\) for every f, g. Uncorrelated asks this only for f, g linear. Example: \(x \sim N(0,\ 1), y = x^{2}\). Then \(\operatorname{cov} (x, y) = \mathrm{E}[x^{3}] = 0\), yet y is a deterministic function of x. The ultrasound version: the in-phase component I and the envelope |z| of circular speckle are uncorrelated (by symmetry) but obviously dependent.
Gaussian exception. If x and y are jointly Gaussian, uncorrelated implies independent, because the joint density factorizes when the off-diagonal of \(\Sigma\) is zero. Two marginally Gaussian variables are not enough; joint Gaussianity is the requirement.
Mean removal
\(\mathrm{E}[xy^{*}]\) (the "correlation" in the signal-processing sense) and cov(x, y) differ by \(\mu _{x}\mu _{y}^{*}\). RF and IQ samples of speckle are zero-mean, so the two coincide. Envelopes, intensities and B-mode pixels are not zero-mean: correlating them without removing the mean gives a large coefficient for two speckle patterns that have nothing to do with each other: \(\mathrm{E}[a]^{2}/\mathrm{E}[a^{2}] = \pi/4 \approx 0.79\) for two independent Rayleigh envelopes, 0.5 for two independent intensities, and close to 1 for log-compressed (dB) B-mode pixels.
Covariance ellipse, principal axes and the two regression lines
For a 2D Gaussian with covariance \(\Sigma\), the contours of constant density are ellipses \((x - \mu )^{T}\Sigma ^{-1}(x - \mu ) = k^{2}\). Their axes are the eigenvectors of \(\Sigma\) (Module 06), with semi-axis lengths \(k\sqrt{\lambda _{i}}\). Three different lines go through the cloud:
- \(\theta\)
- angle of the major axis from the x axis
- \(\Sigma_{xx},\ \Sigma_{xy},\ \Sigma_{yy}\)
- entries of the 2 × 2 covariance matrix \(\Sigma\)
- \(\rho,\ \sigma_x,\ \sigma_y\)
- correlation coefficient and standard deviations, (07.6)
- \(\mu_x,\ \mu_y\)
- means
- \(\mathrm{E}[y \mid x]\)
- conditional mean of y given x
- \(\tan\)
- tangent
- \(\Sigma\)
- 2 × 2 covariance matrix of (x, y)
2D Gaussian: covariance ellipse, eigenvectors and regression lines
SimulationSeeded Gaussian draws (JS); the ellipses, principal axes and regression lines are the closed forms computed from \(\Sigma\).A "1σ ellipse" in 2D contains \(1 - e^{-1/2}\) = 39.3% of the probability and the \(2\sigma\) ellipse 86.5%, because the squared Mahalanobis distance is \(\chi ^{2}\) with 2 degrees of freedom: \(P(d \le k) = 1 - e^{-k^{2}/2}\). The 68% / 95% rule is one-dimensional. A 95% confidence region in 2D needs \(k = 2.45\).
Fitting a line to a calibration plot (measured vs true attenuation, estimated vs reference sound speed) with ordinary least squares regresses one variable on the other and assumes the regressor is noise-free. If both carry error, the OLS slope is biased toward zero (attenuation bias, or regression dilution); the principal axis (total least squares) treats both symmetrically. Module 05 has the least-squares geometry; here you see it as the difference between the red and teal lines.
07.3Estimating covariance from data
Given N independent pairs \((x_{n}, y_{n})\), the natural estimators are
- \((x_n, y_n)\)
- N independent pairs, n = 1 … N
- \(\hat{C}_{xy}\)
- sample covariance
- \(\hat{\rho}\)
- sample correlation coefficient
- \(\bar{x},\ \bar{y}\)
- sample means
- \(N\)
- number of pairs
- \(n\)
- sample index
- \(\hat{C}_{xx},\ \hat{C}_{yy}\)
- sample variances (the same estimator with y = x or x = y)
- \((\cdot)^{*}\)
- complex conjugate
- \(\sum_n\)
- sum over n = 1 … N
Why 1/(N − 1)
- \(x_n\)
- samples, i.i.d. with mean \(\mu\) and variance \(\sigma^{2}\)
- \(\bar{x}\)
- sample mean, (07.9)
- \(\mu,\ \sigma^{2}\)
- true mean and variance
- \(N\)
- number of samples
- \(\mathrm{E}[\cdot]\)
- expectation; result \((N-1)\sigma^{2}\) is why the estimator divides by \(N-1\)
How noisy the estimates are
- \(K\)
- number of independent snapshots
- \(R_{ij},\ \hat{R}_{ij}\)
- true and sample covariance between channels i and j
- \(\hat{C}_{xx}\)
- sample variance, (07.9)
- \(\sigma^{2}\)
- true variance
- \(N\)
- number of real samples
- \(\operatorname{var}(\cdot)\)
- variance of the estimator
- \(\mathrm{E}\lvert\cdot\rvert^{2}\)
- mean-squared error of the estimate
- \(R_{ii},\ R_{jj}\)
- channel powers (diagonal of R)
The sampling distribution of ρ̂: Fisher's z
\(\hat{\rho }\) is bounded in [−1, 1], so its distribution is skewed when \(|\rho |\) is large, and its spread depends on the unknown \(\rho\). Fisher (1915, 1921) found a variance-stabilizing transform:
- \(\hat{\rho}\)
- sample correlation coefficient from N independent bivariate-normal pairs
- \(\rho\)
- true correlation coefficient
- \(z\)
- Fisher's z, approximately normal
- \(\mathcal{N}(m, v)\)
- normal distribution with mean m and variance v
- \(\operatorname{atanh}\)
- inverse hyperbolic tangent
- \(\ln\)
- natural logarithm
- \(N\)
- number of pairs, (07.9)
- \(\approx\)
- approximately distributed as
- \(z,\ \rho,\ \hat{\rho}\)
- Fisher's z and correlation coefficients, (07.12)
- \(N\)
- number of pairs
- \(\operatorname{std}(\cdot)\)
- standard deviation of the estimate
- \(dz/d\rho\)
- derivative of \(z=\operatorname{atanh}\rho\) with respect to \(\rho\)
- \(\Rightarrow\)
- implies
Spurious correlation
When \(\rho = 0, \operatorname{std} (\hat{\rho }) \approx 1/\sqrt{N - 1}\). At \(N = 10, a |\hat{\rho }| \ge 0.5\) happens about 14% of the time by chance alone (t-test: \(t = \hat{\rho }\sqrt{N - 2}/\sqrt{1 - \hat{\rho }^{2}}\) with \(N - 2\) degrees of freedom). Worse, all of the above assumes independent pairs. With serially correlated samples, replace N by \(N_{\mathrm{eff}}\) from section 07.1. Two smooth random curves (or two speckle images over the same ROI) routinely show \(|\hat{\rho }| > 0.5\) with "N = thousands of pixels" because \(N_{\mathrm{eff}}\) is only the number of speckle cells.
Sampling distribution of ρ̂ and the Fisher-z prediction
SimulationMonte Carlo: 500 experiments of N bivariate-normal pairs (JS). The prediction curve is the Fisher-z normal approximation (07.12), itself an approximation.In speckle tracking (Module 11) and elastography, the correlation coefficient of two kernels is computed from a window of a few wavelengths. A 2 mm axial kernel at 5 MHz, 60% bandwidth contains roughly 6 axial speckle cells even if it holds 130 RF samples at 50 MHz. The exact effective count depends on the statistic (rule 11 in 07.8): for a correlation coefficient between RF kernels it is about 2.7 times the cell count, ≈ 17 here. Either way it is of order ten, not 130, and the jitter of the correlation peak, and hence of displacement, scales with it, not with the sample count. The same applies to the per-pixel coherence estimates in SLSC imaging: the axial kernel length sets their variance.
07.4The covariance matrix of a random vector
For a random vector \(\mathbf{x} \in \mathbb{C}^{N}\) with mean \(\pmb{\mu}\), all pairwise covariances are collected into one matrix:
- \(\mathbf{x} \in \mathbb{C}^N\)
- random vector, e.g. the IQ samples of N channels at one depth
- \(\pmb{\mu}\)
- mean vector
- \({}^H\)
- conjugate transpose
- \(\mathbf{R}\)
- \(N \times N\) covariance matrix of the channels (Hermitian)
- \(R_{ij}\)
- covariance of channels i and j
- \(R_{ii},\ \sigma_i^{2}\)
- power (variance) of channel i
- \(\mathrm{E}[\cdot]\)
- expectation, (07.1)
- \(\operatorname{cov}\)
- covariance, (07.5)
Properties, and why each matters
- Hermitian: \(\mathbf{R}^{H} = \mathbf{R}\). Real eigenvalues, orthonormal eigenvectors (spectral theorem, Module 06).
- Positive semidefinite: for any weight vector \(\mathbf{w}, \mathbf{w}^{H}\mathbf{R}\mathbf{w} = \mathrm{E}|\mathbf{w}^{H}(\mathbf{x} - \pmb{\mu})|^{2} \ge 0\). The quadratic form is the output power of a beamformer with weights \(\mathbf{w}\). DAS with unit weights has output power \(\mathbf{1}^{H}\mathbf{R}\mathbf{1} = \sum _{ij}R_{ij}\): the sum of all entries of the channel covariance matrix.
- Linear maps: \(\operatorname{cov} (\mathbf{A}\mathbf{x}) = \mathbf{A}\mathbf{R}\mathbf{A}^{H}\). Apodization, subarray sums, SVD filters all transform \(\mathbf{R}\) this way.
- Eigendecomposition = PCA: \(\mathbf{R} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^{H}\). The first eigenvector is the weight vector that maximizes output power for unit norm (Rayleigh quotient), \(\lambda _{1}\) that power. In Module 06 the SVD of a Casorati matrix was this PCA done on data instead of on \(\mathbf{R}\).
- \(\mathbf{U},\ \mathbf{\Lambda}\)
- eigenvectors and eigenvalues of R, \(\mathbf{R} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^H\)
- \(\mathbf{R}\)
- covariance matrix, (07.14)
- \(\mathbf{C}\)
- correlation matrix, entries \(C_{ij}=\rho_{ij}\)
- \(\mathbf{D}\)
- diagonal matrix of the channel variances
- \(\operatorname{diag}\)
- keeps only the diagonal
- \(\mathbf{x},\ \pmb{\mu}\)
- random vector and its mean, (07.14)
- \(\mathbf{y}\)
- whitened (decorrelated, unit-variance) vector
- \(\mathbf{I}\)
- identity matrix
- \(\lVert\cdot\rVert\)
- 2-norm; \(\lVert\mathbf{y}\rVert^{2}\) is the Mahalanobis distance
- \((\cdot)^{-1/2}\)
- inverse matrix square root
Complex circular Gaussian: the statistics of IQ data
Write \(z = I\) + jQ. The real pair (I, Q) has a 2×2 real covariance with three free numbers, but \(\mathrm{E}[|z|^{2}]\) is only one number. The missing information is in the pseudo-covariance (complementary covariance):
- \(z = I + jQ\)
- complex IQ sample (in-phase I, quadrature Q)
- \(R\)
- covariance, E|z|²
- \(P\)
- pseudo-covariance (complementary covariance)
- \(I,\ Q\)
- in-phase and quadrature components (real)
- \(z\)
- complex sample, (07.16)
- \(\mathbf{z}\)
- vector of complex samples (one per channel)
- \(\mathbf{R},\ \mathbf{P}\)
- covariance and pseudo-covariance matrices
- \(j\)
- imaginary unit
- \(\mathrm{E}[\cdot]\)
- expectation
- \((\cdot)^{H},\ (\cdot)^{T}\)
- conjugate transpose, transpose
- \(\iff\)
- if and only if
- \(\mathbf{0}\)
- zero matrix
- \(a\)
- envelope amplitude |z|
- \(U(-\pi, \pi]\)
- uniform distribution of the phase
- \(\sigma^{2}\)
- mean power \(\mathrm{E}\lvert z\rvert^{2}\)
- \(z\)
- complex IQ sample
- \(p(\cdot)\)
- probability density
- \(\arg z\)
- phase of z, rad
- \(\sim\)
- is distributed as
- \(\mathrm{E}\lvert z\rvert,\ \operatorname{std}\lvert z\rvert\)
- mean and standard deviation of the envelope
- \(\pi\)
- 3.14159…
Fully developed speckle is circular because each IQ sample is a sum of many phasors with independent, uniformly distributed phases: the central limit theorem makes it Gaussian, the uniform phases make it circular. The speckle-statistics module (Module 26) develops Rayleigh, Rician and K distributions from this starting point. Circularity breaks with I/Q gain or phase imbalance in the demodulator, with a coherent (specular, point-target) component that adds a non-zero mean, and for real RF samples (\(Q = 0\), so \(P = R\)).
IQ samples: circularity, Rayleigh envelope, uniform phase
SimulationMonte Carlo IQ samples: circular complex Gaussian passed through a gain and phase imbalance, plus a constant coherent term (JS); theory curves are the circular Rayleigh and uniform laws.Sample covariance with few snapshots: rank deficiency and diagonal loading
With K snapshots stacked as columns of \(\mathbf{X} \in \mathbb{C}^{N\times K}\), the sample covariance is \(\hat{\mathbf{R}} = (1/K)\mathbf{X}\mathbf{X}^{H}\). Its rank is at most min(K, N). With \(K < N\) it is singular, so \(\hat{\mathbf{R}}^{-1}\), which MVDR and whitening need, does not exist. Even for K moderately above N, the eigenvalues of \(\hat{\mathbf{R}}\) are spread: the large ones biased up, the small ones biased down (Marchenko–Pastur), and the small ones dominate the inverse. Reed, Mallett and Brennan (1974) quantified how many snapshots are enough. For the sample-matrix-inverse (SMI) weights \(\mathbf{w} = \hat{\mathbf{R}}^{-1}\mathbf{s}\), with K i.i.d. complex Gaussian training snapshots that do not contain the desired signal, the SINR relative to the optimum \(\mathbf{R}^{-1}\mathbf{s}\) is a random variable whose distribution depends only on K and N:
- \(N\)
- number of channels (array elements)
- \(K\)
- number of i.i.d. complex Gaussian training snapshots, free of the desired signal
- \(\mathbf{s}\)
- steering vector
- \(\hat{\mathbf{R}}\)
- sample covariance \((1/K)\,\mathbf{X}\mathbf{X}^H\)
- \(\rho\)
- SINR of the sample-matrix-inverse weights relative to the optimum (here ρ is not a correlation coefficient)
- \(\mathrm{SINR}(\mathbf{w})\)
- signal-to-interference-plus-noise ratio of weights \(\mathbf{w}\)
- \(\mathbf{R}\)
- true covariance
- \(\mathrm{Beta}(\cdot,\cdot)\)
- beta distribution
- \(\mathrm{E}[\cdot]\)
- expectation
Diagonal loading is the standard repair:
- \(\delta\)
- loading level, power units
- \(\kappa\)
- condition number \(\lambda_{\max}/\lambda_{\min}\)
- \(\hat{\mathbf{R}}_\delta\)
- diagonally loaded sample covariance
- \(\hat{\mathbf{R}}\)
- sample covariance, (07.18)
- \(\mathbf{I}\)
- identity matrix
- \(\lambda_i(\cdot)\)
- i-th eigenvalue
- \(\lambda_{\max}\)
- largest eigenvalue of \(\hat{\mathbf{R}}\)
Eigenvalues of the sample covariance: K snapshots vs N = 32 channels
SimulationSample covariances of K complex Gaussian snapshots drawn from a true R (VCZ triangle plus white noise at 20 dB), eigenvalues by Jacobi; SMI SINR loss by Monte Carlo (JS).Spatial (sub-aperture) averaging and diagonal loading in practice
Ultrasound MVDR has two snapshot problems at once: a pixel offers only a few axial samples (K of order 1–20, set by the pulse length), and those snapshots contain the desired echo itself, which the RMB analysis excludes and which makes small-K behaviour worse (the adaptive weights partly cancel the signal they should pass). Two fixes are combined in every published implementation:
- \(L\)
- subarray length, elements
- \(2K + 1\)
- axial samples in the range gate (here K is a half-length, not the RMB snapshot count)
- \(\mathbf{a}\)
- steering vector (a = 1 after focusing delays)
- \(\Delta\)
- relative loading, typically \(10^{-2}\)
- \(\tilde{\mathbf{R}}\)
- sub-aperture and range-averaged covariance, \(L \times L\)
- \(N\)
- number of array elements
- \(K\)
- half-length of the axial window (samples)
- \(k\)
- axial sample offset within the window
- \(l\)
- sub-aperture start index
- \(\mathbf{x}_l[k]\)
- snapshot of L consecutive elements starting at element l, axial sample k
- \(\mathbf{w}\)
- beamformer weights (MV, with distortionless constraint)
- \(\operatorname{tr}\)
- trace
- \(\mathbf{I}\)
- identity matrix
- \((\cdot)^{H}\)
- conjugate transpose
In adaptive beamforming the "snapshots" for one pixel are the few axial samples in a short range gate (and possibly a few transmit events), typically far fewer than the 64–128 channels. Hence subarray averaging plus loading in every published ultrasound MVDR implementation (Module 25). The same rank argument appears in the reflection-matrix and distortion-matrix methods for aberration correction: the number of independent speckle realizations (isoplanatic patches times transmit angles) bounds the rank you can trust.
07.5Correlation functions in time and the power spectrum
For a WSS process x(t), covariance between two samples depends only on the lag \(\tau\), so it becomes a function:
- \(\tau\)
- time lag, µs
- \(x(t),\ y(t)\)
- WSS random processes
- \(R_{xx}(\tau)\)
- autocorrelation
- \(C_{xx}(\tau)\)
- autocovariance
- \(R_{xy}(\tau)\)
- cross-correlation
- \(\rho(\tau)\)
- normalized autocovariance (dimensionless)
- \(\mu\)
- mean of x
- \(\mathrm{E}[\cdot]\)
- expectation
- \((\cdot)^{*}\)
- complex conjugate
Wiener–Khinchin
- \(R_{xx}(\tau)\)
- autocorrelation of the WSS process x(t)
- \(S_{xx}(f)\)
- power spectral density, power per Hz
- \(f\)
- frequency, MHz
- \(S_{xx}\)
- power spectral density
- \(j\)
- imaginary unit
- \(d\tau,\ df\)
- integration elements
- \(\int\)
- integral over all lags or all frequencies
- \(h,\ H(f)\)
- impulse response and frequency response
- \(\ast\)
- convolution
- \(N_0\)
- white-noise power spectral density
- \(y\)
- filter output
- \(x\)
- input process
- \(R_{yy},\ R_{xx}\)
- output and input autocorrelations
- \(R_{hh}\)
- deterministic autocorrelation of the impulse response
- \(S_{xx},\ S_{yy}\)
- input and output power spectra
- \(t,\ \tau\)
- time and lag, µs
- \(\lvert\cdot\rvert\)
- magnitude
- \((\cdot)^{*}\)
- complex conjugate
For a Gaussian-envelope pulse with centre frequency \(f_{0}\) and −6 dB (amplitude) bandwidth B, the autocorrelation envelope is Gaussian with full width at half maximum
- \(f_0\)
- centre frequency, MHz
- \(\sigma_t\)
- standard deviation of the Gaussian envelope, µs
- \(B\)
- −6 dB (amplitude) bandwidth, MHz
- \(\tau_{\mathrm{FWHM}}\)
- full width at half maximum of the autocorrelation envelope, µs
- \(c\)
- sound speed, 1540 m/s
- \(p(t)\)
- Gaussian-windowed pulse
- \(t\)
- time, µs
- \(\tau\)
- lag, µs
- \(R_{pp}\)
- pulse autocorrelation
- \(\propto\)
- proportional to
- \(\ln\)
- natural logarithm
- \(\cos\)
- cosine
- \(\text{axial correlation length}\)
- distance over which echoes stay correlated, mm; \(c/2\) accounts for the round trip
Estimating R(τ) and S(f)
- \(\hat{R}(k)\)
- estimated autocorrelation at lag k samples
- \(x[n]\)
- sampled signal
- \(N\)
- number of samples
- \(k\)
- lag, samples
- \(n\)
- sample index
- \(X(f)\)
- discrete-time Fourier transform of the N samples
- \(f\)
- frequency, cycles/sample
- \(\lvert X(f)\rvert^{2}/N\)
- periodogram, a spectrum estimate
- \(\mathrm{DTFT}\)
- discrete-time Fourier transform
- \((\cdot)^{*}\)
- complex conjugate
Speckle RF line: autocorrelation, PSD and correlation length vs bandwidth
SimulationOne simulated RF line: white Gaussian reflectivity (4096 samples at 50 MHz) filtered by the Gaussian pulse spectrum with FFTs (JS); the theory curves are the closed forms (07.23)–(07.24).One word, six objects: a terminology table
The ultrasound literature uses "correlation" and "coherence" for at least six distinct quantities. They differ in what is averaged, over which index, and how it is normalized.
| Name | Definition | Normalization / range | Ultrasound use |
|---|---|---|---|
| Correlation coefficient | \(\rho = \operatorname{cov} (x, y)/(\sigma _{x}\sigma _{y})\), one number | mean removed, scale removed; [−1, 1] (complex: \(|\rho | \le 1\)) | Kernel match quality, decorrelation, NCC peak height in speckle tracking |
| Correlation function | \(R_{xy}(\tau ) = \mathrm{E}[x(t + \tau )y^{*}(t)]\), a function of lag | raw (power units) or normalized \(\mathrm{NCC}(\tau )\) ∈ [−1, 1] | Time-delay estimation (Module 11), Kasai Doppler: \(\bar{f}\) = PRF·arg \(R(1)/(2\pi )\) on slow-time IQ |
| Correlation matrix | \(\mathrm{E}[\mathbf{x}\mathbf{x}^{H}]\) (signal processing) or \(\mathbf{D}^{-1/2}\mathbf{RD}^{-1/2}\) (statistics) | first: no mean removal; second: unit diagonal | Channel or slow-time ensembles; identical to \(\mathbf{R}\) for zero-mean IQ |
| Covariance matrix | \(\mathbf{R} = \mathrm{E}[(\mathbf{x} - \pmb{\mu})(\mathbf{x} - \pmb{\mu})^{H}]\) | Hermitian PSD, power units | MVDR, PCA / SVD clutter filtering, eigen-based Doppler |
| Spatial coherence | \(C(m) = \langle R_{i,i+m}/\sqrt{R_{ii}R_{i+m,i+m}}\rangle _{i}\), function of element lag | normalized per element pair; 1 at \(m = 0\) | VCZ, SLSC imaging, coherence factor CF, sound-speed and aberration metrics (07.6) |
| Spectral coherence (magnitude-squared) | \(\gamma ^{2}_{xy}(f) = |S_{xy}(f)|^{2}/(S_{xx}(f)S_{yy}(f))\) | [0, 1] per frequency; needs averaging (=1 for one segment) | Frequency-resolved similarity, e.g. of RF echoes before and after deformation; system identification |
Magnitude-squared coherence estimated from a single segment is identically 1 for any two signals, because \(|X Y^{*}|^{2} = |X|^{2}|Y|^{2}\). Coherence of any kind is a statement about an average. The same holds for spatial coherence from one snapshot: \(\hat{R}_{ij}/\sqrt{\hat{R}_{ii}\hat{R}_{jj}}\) has unit magnitude for every pair. SLSC uses an axial kernel (several snapshots) for exactly this reason, and takes the real part of the normalized correlation.
07.6Spatial covariance across an array: the van Cittert–Zernike theorem
Fix one focal point. After receive focusing delays, the IQ samples \(s_{i}\) across the N elements form a random vector whose ensemble is "all possible speckle realizations in the focal zone". Its covariance matrix \(\mathbf{R} = \mathrm{E}[\mathbf{s}\mathbf{s}^{H}]\) is the channel covariance. What does it look like for diffuse scatterers?
Derivation (focal zone, narrowband)
Let \(a(\xi )\) be the scatterer reflectivity across the focal plane at depth z, delta-correlated: \(\mathrm{E}[a(\xi )a^{*}(\xi ')] = \sigma _{a}^{2} \delta (\xi - \xi ')\). In the focal plane the transmit field is the Fraunhofer pattern of the transmit apodization t(x), and the focused receive signal on element i is the same Fourier kernel at element position \(x_{i}\):
- \(\xi\)
- lateral position in the focal plane, mm
- \(z\)
- focal depth, mm
- \(x_n,\ x_i\)
- element positions, mm
- \(t_n\)
- transmit apodization weights
- \(a(\xi)\)
- scatterer reflectivity, delta-correlated
- \(k\)
- wavenumber 2π/λ
- \(s_i\)
- focused IQ sample on receive element i
- \(P_t(\xi)\)
- transmit field at lateral position \(\xi\) in the focal plane
- \(j\)
- imaginary unit
- \(\sum_n\)
- sum over transmit elements
- \(\int d\xi\)
- integral over the focal plane
- \(\sigma_a^2\)
- reflectivity power per unit length
- \(\lvert P_t(\xi)\rvert^2\)
- transmit intensity across the focal plane (the incoherent source)
- \(\mathrm{E}[\cdot]\)
- expectation over scatterer realizations
- \(s_i,\ s_j\)
- focused IQ samples on receive elements i and j, (07.26)
- \((\cdot)^{*}\)
- complex conjugate
- \(k,\ x_i,\ x_j,\ z,\ \xi,\ j\)
- as in (07.26)
- \(t \star t\)
- autocorrelation of the transmit apodization
- \(m\)
- element lag (separation in elements)
- \(N_t\)
- number of transmit elements
- \(C(m)\)
- spatial coherence, per-pair normalized covariance averaged along a diagonal of R
- \(\mathrm{FT}\)
- Fourier transform
- \(R_{ij}\)
- channel covariance, (07.14)
- \(\propto\)
- proportional to
- \(x_i - x_j\)
- element separation, mm
- \(i,\ j\)
- element indices
- \(\lvert\cdot\rvert\)
- magnitude
- \(P_t\)
- transmit field, (07.26)
Assumptions: diffuse, delta-correlated scatterers (fully developed speckle); focal zone (away from the focus the insonified region broadens and C(m) narrows); Fresnel approximation; isoplanatic medium; no noise; a known transmit aperture, whose autocorrelation sets the curve. Mallart and Fink (1991) derived and measured this pulse-echo form of VCZ, and in 1994 used it to define a focusing criterion based on the same coherence (the coherence factor below is a closely related, later formulation used as an image weighting, usually credited to Hollman, Rigby and O'Donnell, 1999).
hanning(64) transmit weights. Receive apodization does not change them: it cancels in the per-pair normalization.Three things that pull the curve down
- \(\mathbf{R}_s\)
- speckle (signal) covariance
- \(\sigma_n^2\)
- noise power per element, uncorrelated between elements
- \(\mathrm{SNR}\)
- per-element signal-to-noise power ratio
- \(C_s(m)\)
- noise-free coherence
- \(\mathbf{R}\)
- measured covariance, (07.14)
- \(\mathbf{I}\)
- identity matrix
- \(C(m)\)
- coherence at lag m, (07.28)
- \(m\)
- element lag
- \(\phi_i\)
- phase error on element i, rad
- \(P_t^{\mathrm{ab}}\)
- aberrated transmit field
- \(\langle\cdot\rangle_i\)
- average along the m-th diagonal
- \(D_\phi(m)\)
- phase structure function \(\mathrm{E}[(\phi_i - \phi_{i+m})^2]\), rad²
- \(t_n\)
- transmit weights, (07.26)
- \(s_i\)
- focused sample on element i, (07.26)
- \(\to\)
- replaced by
- \(R_{ij}\)
- channel covariance, (07.14)
- \(\xi,\ k,\ x_i,\ z\)
- as in (07.26)
- \(C_{\mathrm{VCZ}}(m)\)
- coherence without aberration (van Cittert–Zernike), (07.28)
- \(N_t,\ N_r\)
- number of transmit and receive elements
- \(C(m)\)
- coherence at lag m
- \(e^{-D_\phi(m)}\)
- coherence loss from phase errors, dimensionless
- \(\hat{c}\)
- sound speed assumed by the beamformer, m/s
- \(\hat{z}\)
- depth the beamformer assigns to the arrival time, mm
- \(\phi_i(\hat{c})\)
- phase error on element i caused by the wrong sound speed, rad
- \(f_0\)
- centre frequency, MHz
- \(x_i\)
- element position, mm
- \(z\)
- true focal depth, mm
- \(c\)
- true sound speed, m/s
- \(\sqrt{\cdot}\)
- square root: path length to the element, mm
- \(\ne\)
- not equal
Channel covariance and VCZ coherence: N = 64, focus 30 mm, 5 MHz, pitch 0.3 mm
SimulationMonte Carlo of the narrowband focal-plane (Fraunhofer, paraxial) model of 07.6: delta-correlated scatterers on a 256-point lateral grid, weighted by the transmit beam, a thin phase screen and white noise (JS). Not a full time-domain channel-data simulation; the MATLAB script m07_vcz_coherence.m does that with exact geometry.Why maximizing coherence estimates sound speed
For diffuse scatterers in a homogeneous medium, the triangle is the upper bound on C(m) for a given aperture: any residual phase across the aperture (wrong \(\hat{c}\), aberration) or noise can only lower it. With \(\hat{c}\) as the only free parameter, the residual phase vanishes at \(\hat{c} = c\), so coherence-based metrics (SLSC, CF, or the full curve) as functions of \(\hat{c}\) peak at the true speed. The sweep panel above computes the expected SLSC and CF analytically for every \(\hat{c}\). The curvature of the peak sets the estimator's sensitivity: with a 19.2 mm aperture focused at 30 mm (F/1.6), a 40 m/s error costs ≈ 4% of SLSC (lags 1–6) but ≈ 28% of CF, which uses every lag. Estimation noise on SLSC from a short axial kernel is of the same order as that 4%, so practical estimators average many pixels. The paraxial edge defocus phase is ≈ \((\pi D^{2}/(2\lambda z))\cdot (\Delta c/c)\), about 1.6 rad here for \(\Delta c = 40\) m/s (1.45 rad with the exact hyperbolic delays): it grows with the aperture squared over \(\lambda z\), so larger apertures and smaller F-numbers sharpen the peak.
Why deviation from the triangle measures aberration
Because the VCZ triangle is known a priori for a given aperture and apodization, the ratio \(C_{\mathrm{measured}}(m)/C_{\mathrm{VCZ}}(m)\) isolates everything that is not ideal: the noise factor at \(m \ge 1\) (a constant), and the two-way aberration factor (decaying with lag, roughly \(\exp (-D_{\varphi}(m))\), half from receive and half from transmit broadening). The lag at which it decays gives the screen's correlation length; its depth gives the rms phase. Correcting the screen (estimating \(\varphi _{i}\) from neighbouring-element cross-correlation, as in Flax and O'Donnell 1988, or from the distortion matrix) restores the curve toward the triangle, and the restored coherence is the standard figure of merit for an aberration-correction method.
Target B (sound speed): global estimators sweep \(\hat{c}\) and maximize a focus-quality metric (coherence factor, SLSC, image sharpness); local estimators (tomographic or differentiable-beamforming methods) need the same coherence logic per region plus an inversion. Target C (aberration): the measured coherence deficit is both the diagnostic and the cost function; the distortion-matrix approach (Module 06) can be read as finding the phase law that makes the de-scanned field maximally coherent. ULM (Target A) uses coherence indirectly: microbubbles are point-like, not diffuse, so their channel covariance is rank one (no VCZ triangle), and that contrast is what coherence-based bubble detection exploits.
The coherence factor \(\mathrm{CF} = |\sum _{i}s_{i}|^{2}/(N \sum _{i}|s_{i}|^{2})\) for diffuse speckle with a uniform aperture has an ensemble value (ratio of expectations) of \(\sum _{m}(N - |m|)^{2}/N^{3}\) \(= 2/3 + 1/(3N^{2}) \approx 2/3\). For white noise it is exactly 1/N. The average of the per-snapshot ratio is different (\(\approx 0.57\) for \(N = 64\)), because \(\mathrm{E}[X/Y] \neq \mathrm{E}[X]/\mathrm{E}[Y]\). Papers differ on which they report, and on whether CF is computed per sample or over an axial kernel; the readout shows both.
Lediju et al. (2011) compute, for each lag m, the normalized correlation of the delayed RF signals of elements i and \(i + m\) over a short axial kernel, average it over the \(N - m\) element pairs, and sum lags 1 to \(M: \mathrm{SLSC} = \sum _{m=1}^{M} \hat{R}(m)\). With IQ data the real part of the complex normalized correlation is used. Check three things before comparing numbers across papers: (1) whether "lag" counts element separation (adjacent elements = lag 1, as here) or is expressed as a fraction of the aperture; (2) whether M is quoted in elements or as \(Q = M/N\), so the same Q on a 64- and a 128-element array means different element counts; (3) whether the curve is averaged with \(1/(N - m)\) and whether the sum is normalized (by M or by the ideal VCZ value, as in the demo's "% of ideal"). The ideal value depends on the transmit aperture through the VCZ triangle, so an SLSC value is meaningful only together with \(N_{t}, M\) and the apodization.
07.7Visual summary
Twelve pictures for second-order statistics, each with the five-second answer and its number. All curves come from the closed forms in the cards (Fisher z, Rayleigh, Marchenko–Pastur, the VCZ lag sums) or from a fixed-seed draw of 260 Gaussian pairs. The slider sets the correlation coefficient \(\rho\) used by the four scatter and density cards.
07.8How experts think about it
The expert model has three layers. Speckle is a circular complex Gaussian field, so one covariance matrix is the whole story. That covariance is set by the imaging system, not the tissue: in time by the pulse spectrum (Wiener–Khinchin), across the array by the transmit intensity at the focus (van Cittert–Zernike). Every departure from the system-predicted covariance (noise, aberration, a wrong sound speed, a point target, motion) is a measurable signal, and every estimate of it carries a sampling error set by the number of independent speckle cells, not by the number of samples.
Delta-correlated scatterers make every point of the focal plane an independent source, weighted by the transmit intensity \(|P_{t}(\xi )|^{2}\). After focusing, the receive kernel is \(e^{-jkx_{i}\xi /z}\), so
- \(\mathrm{E}[\cdot]\)
- expectation over scatterer realizations
- \(s_i,\ s_j\)
- focused IQ samples on receive elements i and j
- \(P_t(\xi)\)
- transmit field across the focal plane
- \(\xi,\ z\)
- lateral position and focal depth, mm
- \(x_i,\ x_j\)
- element positions, mm
- \(k\)
- wavenumber \(2\pi/\lambda\), rad/mm
- \(j\)
- imaginary unit
- \(\propto\)
- proportional to (constant factor is the reflectivity power \(\sigma_a^2\)), as in (07.27)
Because \(|P_{t}|^{2}\) is itself the transform of the transmit-aperture autocorrelation, C(m) is that autocorrelation in element units, independent of frequency. It is the optical VCZ theorem with the transmit beam as the incoherent source (Mallart and Fink 1991). The phase of \(P_{t}\) drops out: only the insonified intensity matters.
The autocorrelation of a rect is a triangle. The receive array only samples it: receive apodization cancels in the per-pair normalization, and the receive width sets only the largest available lag, \(N_{r} - 1\). A 32-element transmit received on 64 elements gives \(1 - m/32\): the MATLAB lab with Ntx = 32 prints \(C(16) = 0.51\) and \(C(32) = 0.01\). A 64-element Hann transmit gives 0.67, 0.18, 0.01 at lags 16, 32, 48. A single-element transmit insonifies a wide region and the coherence is gone within a lag or two. See (07.28).
Receive focusing still forms the VCZ integral at any depth z, but the source is now the defocused transmit beam. Geometrically it is about \(D|1 - z/z_{f}|\) wide, and the FT of a source of width w has its first null at a separation \(\lambda z/w\). Exact-geometry check (64 elements, 0.3 mm pitch, 5 MHz, \(z_{f} = 30\) mm, F/1.6): at 27 mm the −6 dB beam is 1.5 mm and \(C(16) = 0.18\) instead of 0.73; at 40 mm \(C(4) = 0.50\); at 15 mm \(C(1) = 0.59\). The depth of field \(7\lambda F\#^{2}\) is only ≈ 5 mm here.
Noise independent across channels adds \(\sigma _{n}^{2}\mathbf{I}\): only the diagonal grows, so every normalized off-diagonal entry is divided by \(1 + 1/\mathrm{SNR}\). The step between lag 0 and lag 1 is a channel-SNR meter: See (07.29).
- \(\mathrm{SNR}\)
- per-element signal-to-noise power ratio (dimensionless; dB = 10 log)
- \(C(1)\)
- measured lag-one coherence, as in (07.29)
- \(C_{\mathrm{VCZ}}(1)\)
- noise-free lag-one coherence, as in (07.28)
- \(\rho,\ \rho_{\mathrm{meas}}\)
- true and noise-biased correlation coefficient between two kernels
- \(\mathrm{SNR}_1,\ \mathrm{SNR}_2\)
- SNRs of the two kernels
\(N = 128\), measured \(C(1) = 0.80,\) \(C_{\mathrm{VCZ}}(1) = 0.992: \mathrm{SNR} = 4.2\) (6.2 dB). This is lag-one coherence as an image-quality metric (Long, Bottenus and Trahey 2018). The second formula is the same bias for a correlation coefficient between two noisy kernels.
A near-field screen costs ≈ \(e^{-D_{\varphi}/2}\) on receive (phase differences averaged along each diagonal) and the same on transmit (beam broadening). \(D_{\varphi}(m) = \mathrm{E}[(\varphi _{i} - \varphi _{i+m})^{2}]\) grows with lag and saturates at \(2\sigma ^{2}\). \(\sigma = 0.5\) rad, screen correlation \(e^{-m^{2}/(2\ell ^{2})}, \ell = 6\) elements, \(N = 64: C(1) = 0.977\) against \(0.984, C(8) = 0.65\) against 0.875. For \(\varphi _{i} = \alpha i\) the receive factor \(e^{-j\alpha m}\) and transmit factor \(e^{+j\alpha m}\) cancel: the focus moves sideways and the curve stays on the triangle. See (07.30).
For a given transmit aperture the VCZ curve is an upper bound: any residual phase can only lower \(|\langle e^{j\Delta \varphi}\rangle |\). A wrong \(\hat{c}\) leaves a mostly quadratic residual with edge phase ≈ \((\pi D^{2}/(2\lambda z))(\Delta c/c),\ 1.6\) rad for \(D = 19.2\) mm, \(z = 30\) mm, 5 MHz, \(\Delta c = 40\) m/s. With \(\hat{c}\) the only free parameter the residual vanishes at \(\hat{c} = c\), so CF, SLSC or the full curve peak there. CF weights all lags and is sharp (−28% in the paraxial demo, −30 to −36% with exact geometry); short-lag SLSC is flat (−4%) but robust to noise. The MATLAB lab recovers \(c_{\mathrm{true}} = 1510\) m/s to within 1 m/s.
Mainlobe speckle keeps C(m) near the triangle at small m whatever its amplitude, while thermal noise and much reverberation clutter have almost no coherence beyond lag 0. Summing only short lags keeps the speckle signal and avoids the long lags that aberration and finite kernels make noisy. A hypoechoic lesion filled with clutter turns dark in SLSC because its coherence is low, not because its echo is weak (Lediju et al. 2011).
Each IQ sample sums many phasors with independent phases uniform over \(2\pi\): the central limit theorem gives a Gaussian, the uniform phases make it rotation invariant. Then \(\mathbf{R}\) alone is the complete description, \(p(\mathbf{z}) = \exp (-\mathbf{z}^{H}\mathbf{R}^{-1}\mathbf{z})/(\pi ^{N} \det \mathbf{R})\). Envelope \(\mathrm{SNR} \sqrt{\pi /(4 - \pi )} = 1.91\); the MATLAB lab gets 1.915 from 30 random phasors per cell. VCZ, MVDR and the coherence estimators all assume it.
Correlation tests only linear dependence: for symmetric \(x, \operatorname{cov} (x, x^{2}) = \mathrm{E}[x^{3}] = 0\). A jointly Gaussian density with diagonal covariance factorizes, so there uncorrelated implies independent: I and Q of circular speckle are independent. I and the envelope are uncorrelated by symmetry yet dependent: \(\operatorname{var} (I^{2}) = 2\sigma _{I}^{4}\) and \(\operatorname{var} (|z|^{2}) = 4\sigma _{I}^{4}\), so \(\operatorname{corr} (I^{2}, |z|^{2}) = 2/\sqrt{8} = 1/\sqrt{2}\).
The delta method on z gives \(d\rho /dz = 1 - \rho ^{2}\): the spread collapses as \(|\rho | \to 1\) and \(\hat{\rho }\) is skewed toward 0, so symmetric ± intervals on \(\hat{\rho }\) fail near \(\pm 1\). Fisher's z is close to normal with a variance that does not depend on \(\rho\). \(\hat{\rho } = 0.6\) from 12 patients: \(z = 0.693 \pm 1.96/3\) \(= [0.040,\ 1.347]\), so \(\rho \in [0.04,\ 0.87]\). The MATLAB lab checks it: \(\operatorname{std} z = 0.2415\) against \(1/\sqrt{17} = 0.2425\), coverage 95.0%. See (07.12).
\(\operatorname{var} (\bar{x}) = (\sigma ^{2}/N)\sum _{k}(1 - |k|/N)\rho (k)\), so \(N_{\mathrm{eff}} = T/\int \rho _{q}(\tau )d\tau\), where \(\rho _{q}\) is the correlation of the quantity being averaged. For an intensity mean \(\rho _{q}\) = |envelope of \(r|^{2}, \int = \sqrt{2\pi}\sigma _{t} = 0.94/B\). For \(\hat{\rho }\) between RF kernels (Bartlett) \(\rho _{q} = r^{2}, \int \approx 0.47/B\), since the carrier averages \(\cos ^{2}\) to \(\tfrac{1}{2}\). A 5 mm window (T·B = 19.5) holds 15.6 cells of 1.25/B but \(N_{\mathrm{eff}} = 20.8\) for an intensity mean and 41.5 for \(\hat{\rho }\); Monte Carlo 21.8 and 40.0. See (07.4).
White scatterers convolved with the PSF give a speckle autocorrelation equal to the PSF autocorrelation. Gaussian pulse: the RF autocorrelation envelope has FWHM \(4\sqrt{2} \ln 2/(\pi B) = 1.25/B\) (MATLAB: 1.246). The intensity autocovariance is \(|R|^{2}\) (Isserlis), FWHM \(1/\sqrt{2}\) of that, 0.88/B (0.23 mm here). Laterally, the two-way uniform-aperture PSF \(\operatorname{sinc} ^{2}\) has a −6 dB width of \(0.89 \lambda F\#\) and an intensity-speckle FWHM of \(0.83 \lambda F\#\).
Wiener–Khinchin for WSS processes; equivalently the DFT approximately diagonalizes a Toeplitz covariance, with PSD samples as eigenvalues. White reflectivity through the system gives \(S_{yy} = N_{0}|H|^{2}\), which is why spectral QUS divides by a reference-phantom spectrum to remove the system. For a finite record the periodogram is exactly the DFT of the biased \(\hat{R}\) (MATLAB: agreement to \(5 \times 10^{-16}\)), and each bin stays exponential, std/mean ≈ 1, at 256 or 4096 samples. See (07.22).
\(\hat{\mathbf{R}} = \mathbf{XX}^{H}/K\) is a sum of K rank-one terms. Reed, Mallett and Brennan showed that the SINR of \(\hat{\mathbf{R}}^{-1}\mathbf{s}\) relative to \(\mathbf{R}^{-1}\mathbf{s}\) is \(\mathrm{Beta}(K - N + 2, N - 1)\) distributed for any \(\mathbf{R}\) and \(\mathbf{s}\). Monte Carlo, \(N = 32: K = 61\) gives \(0.501, K = 128\) gives 0.759 (theory 0.760), \(K = 32\) gives \(0.061 = 2/33\). See (07.18).
Loading adds \(\delta\) to every eigenvalue with eigenvectors unchanged, capping \(\kappa\) at \((\lambda _{\max} + \delta )/\delta\) and bounding the white-noise gain, hence the sensitivity to delay errors (Cox, Zeskind and Owen 1987). In Demo 07.4 \((N = 32, K = 64)\) loading at −20 dB raises the mean SINR ratio from 0.52 to 0.71. Sub-aperture averaging of the \(N - L + 1\) overlapping \(L \times L\) covariances multiplies the snapshots and decorrelates coherent echoes (Shan, Wax and Kailath 1985): 64 elements, \(K = 40, L = 32\) gives 1320 snapshots for a 32 × 32 matrix.
\(\mathrm{CF} = |\sum s_{i}|^{2}/(N\sum |s_{i}|^{2})\), by Parseval the fraction of aperture-domain energy at zero spatial frequency. In lag form its ensemble value is \(\mathrm{CF} = \sum _{m}(N - |m|)C(m)/N^{2}\): the uniform triangle gives \(\sum (N - |m|)^{2}/N^{3}\) \(= 2/3 + 1/(3N^{2}) = 0.667\) for \(N = 64\). A point at focus has identical channels, \(\mathrm{CF} = 1\). For white noise only \(m = 0\) survives, \(\mathrm{CF} = N/N^{2} = 1/N\); its per-snapshot value is \(\mathrm{Beta}(1, N - 1)\), mean also 1/N. With noise: \((N\cdot \mathrm{CF}_{s} + \sigma _{n}^{2})/(N(1 + \sigma _{n}^{2})),\) \(0.34\) at 0 dB.
The Cramér–Rao bound for delay estimation between partially correlated speckle signals (Walker and Trahey 1995) scales as \(\sqrt{1/\rho ^{2} - 1}\) at high SNR, times a factor that falls with kernel length, centre frequency and bandwidth. Small decorrelation from strain, out-of-plane motion or noise therefore dominates displacement error long before \(\rho\) "looks low". The same number is used directly as a quality metric: \(1 - \hat{\rho }\) flags out-of-plane motion and tracking failure, lag-one coherence grades channel data, and decorrelation rate measures flow or motion.
For independent circular Gaussian channels, \(|\hat{\rho }|^{2}\) from K snapshots is \(\mathrm{Beta}(1, K - 1)\), so \(\mathrm{E}|\hat{\rho }|^{2} = 1/K\) and \(\mathrm{E}|\hat{\rho }| = \Gamma (K)\Gamma (3/2)/\Gamma (K + \tfrac{1}{2})\) \(\approx \sqrt{\pi /(4K)}\) (Monte Carlo at \(K = 16:\ 0.222\) against 0.223). With \(K = 1\) every pair has unit magnitude, and the same holds for magnitude-squared coherence from one segment. SLSC therefore takes the real part, which has zero mean for noise, over an axial kernel.
07.9Cheat sheet
| Concept | Formula | Ultrasound use |
|---|---|---|
| Expectation | \(\mathrm{E}[a x + b y]\) \(= a\mathrm{E}[x] + b\mathrm{E}[y];\) \(\mathrm{E}[\mathbf{Ax}] = \mathbf{A}\pmb{\mu}\) | Mean of any linear beamformer output |
| Effective N | \(\operatorname{var} (\bar{x}) = \sigma ^{2}/N \cdot \sum (1 - |k|/N)\rho (k)\) \(= \sigma ^{2}/N_{\mathrm{eff}}\) | Number of speckle cells in a QUS ROI or tracking kernel |
| Correlation coefficient | \(\rho = \operatorname{cov} /(\sigma _{x}\sigma _{y})\) \(= \cos \angle (\tilde{x}, \tilde{y}), |\rho | \le 1\) | NCC in speckle tracking, decorrelation |
| Regression vs PCA | \(\mathrm{slope}_{y|x} = \rho \sigma _{y}/\sigma _{x}; \tan 2\theta = 2\rho \sigma _{x}\sigma _{y}/(\sigma _{x}^{2} - \sigma _{y}^{2})\) | Calibration fits with error in both variables |
| Sample covariance | \(1/(N - 1)\) with mean removal; 1/K if zero-mean known | Channel covariance from IQ snapshots |
| Fisher z | \(\operatorname{atanh} \hat{\rho } \sim N(\operatorname{atanh} \rho ,\ 1/(N - 3));\) \(\operatorname{std} \hat{\rho } \approx (1 - \rho ^{2})/\sqrt{N - 1}\) | Confidence on correlation peaks and coherence values |
| Covariance matrix | \(\mathbf{R} = \mathrm{E}[(\mathbf{x} - \pmb{\mu})(\mathbf{x} - \pmb{\mu})^{H}]\), Hermitian PSD; \(\mathbf{w}^{H}\mathbf{Rw}\) = output power | DAS power = \(\sum _{ij}R_{ij}\); MVDR; PCA |
| Circular complex Gaussian | \(\mathrm{E}[\mathbf{zz}^{T}] = 0; |z|\) Rayleigh, SNR 1.91; phase uniform | Fully developed speckle IQ |
| Rank of \(\hat{R}\), loading | rank ≤ \(\min (K, N); \lambda _{i} + \delta ; K \approx 2N\) for 3 dB (RMB) | Adaptive beamforming with few range samples |
| Wiener–Khinchin | \(S(f) = \mathrm{FT}\{R(\tau )\}\); white in, \(R_{yy} = N_{0}R_{hh}\) | Speckle autocorrelation = pulse autocorrelation |
| Correlation length | \(\tau _{\mathrm{FWHM}} \approx 1.25/B\) (Gaussian); axial \(c\tau /2\) | Speckle cell size, \(N_{\mathrm{eff}}\), tracking kernel length |
| VCZ coherence | \(C(m) = (t \star t)(m)/(t \star t)(0)\); uniform: \(1 - |m|/N\) | Ideal spatial coherence of diffuse speckle at focus |
| Noise on coherence | \(C(m \ge 1) \times 1/(1 + 1/\mathrm{SNR}),\) \(C(0) = 1\) | SLSC noise rejection; channel SNR estimate |
| Aberration on coherence | ≈ × \(\exp (-D_{\varphi}(m))\) two-way (receive \(e^{-D_{\varphi}/2}\) × transmit \(e^{-D_{\varphi}/2}\)) | Aberration diagnosis and correction metric |
| SLSC, CF | \(\sum _{m=1}^{M} \hat{C}(m);\) \(\mathrm{CF} = |\sum s|^{2}/(N\sum |s|^{2})\): point 1, speckle \(2/3 + 1/(3N^{2})\), noise 1/N | Coherence imaging; sound-speed optimization |
| Off-focus coherence | insonified width \(w \approx D|1 - z/z_{f}|\); first null at separation ≈ \(\lambda z/w\) | Why coherence metrics need the focal zone or synthetic focusing |
| Noise bias on \(\rho\) | \(\rho _{\mathrm{meas}} = \rho /\sqrt{(1 + 1/\mathrm{SNR}_{1})(1 + 1/\mathrm{SNR}_{2})}\); channel \(\mathrm{SNR} = C(1)/(C_{\mathrm{VCZ}}(1) - C(1))\) | Lag-one coherence; NCC peak height |
| Coherence of pure noise | \(\mathrm{E}|\hat{\rho }|^{2} = 1/K;\) \(\mathrm{E}|\hat{\rho }| \approx \sqrt{\pi /(4K)}\); one snapshot: \(|\hat{\rho }| = 1\) | Bias floor of magnitude coherence; why SLSC uses Re |
| Effective N by statistic | \(N_{\mathrm{eff}} = T/\int \rho _{q}d\tau\): intensity mean T·B/0.94, RF \(\hat{\rho }\) T·B/0.47, cells T·B/1.25 | Variance of QUS, tracking and coherence estimates |
| Decorrelation and jitter | \(\sigma _{\Delta t} \propto \sqrt{1/\rho ^{2} - 1}\) (CRLB, high SNR) | Displacement error budget; decorrelation as quality metric |
07.10Interview and qualifying-exam questions
I1Derive the van Cittert–Zernike result for pulse-echo speckle intuitively. What does the coherence curve of a uniform focused aperture look like? · onsite
Diffuse scatterers are independent sources, so their contributions add in power: the focal plane acts as an incoherent source whose brightness is the transmit intensity \(|P_{t}(\xi )|^{2}\). After receive focusing, a source at lateral position \(\xi\) reaches element i with phase \(-kx_{i}\xi /z\), a plane wave across the aperture. The covariance of elements i and j is the sum over sources of \(e^{-jk(x_{i} - x_{j})\xi /z}\) weighted by \(|P_{t}|^{2}\): a Fourier transform of the source intensity evaluated at the separation. Since \(|P_{t}|^{2}\) is the transform of the transmit-aperture autocorrelation, \(C(m) = (t \star t)(m)/(t \star t)(0)\). Uniform: a triangle \(1 - m/N_{t}\), zero at the transmit width; \(N_{t} = 64\) gives \(C(32) = 0.5\).
Follow-up: 32-element transmit, 64-element receive? Zero at lag 32; receive apodization does not change the normalized curve. Away from the focus? The insonified width grows to ≈ \(D|1 - z/z_{f}|\) and the first null moves to a separation ≈ \(\lambda z/w\): at 27 mm with a 30 mm focus, C(16) drops from 0.73 to 0.18.
I2Why does short-lag spatial coherence suppress noise and clutter while keeping speckle? · screen
Channel noise is uncorrelated between elements, so it adds only to the diagonal of \(\mathbf{R}\): it raises lag 0 and scales every lag ≥ 1 by \(\mathrm{SNR}/(1 + \mathrm{SNR})\). SLSC never uses lag 0, so noise does not add to it; at 0 dB SNR the lag-1 value halves (0.98 → 0.49) instead of being swamped. Much reverberation clutter is also nearly incoherent beyond lag 0. Mainlobe speckle, in contrast, keeps \(C(m) \approx 1 - m/N\) at short lags regardless of its brightness, so a clutter-filled anechoic region turns dark. Short lags are used because long lags are both small (triangle) and noisy, and are the first to be destroyed by aberration.
Follow-up: why the real part? The magnitude of a normalized correlation of pure noise is biased up to ≈ \(\sqrt{\pi /(4K)}\) for K snapshots; the real part has zero mean.
I3You sweep the assumed sound speed and plot coherence (CF or SLSC). What do the position, height and width of the peak tell you? · research
Position: in a homogeneous medium the residual focusing phase vanishes only at \(\hat{c} = c\), so the argmax estimates c (the MATLAB lab recovers 1510 m/s to within 1 m/s). In layered media it is a path-averaged speed. Height: compared with the VCZ ideal (SLSC 5.67 for N = 64, M = 6; CF 0.667), a low peak means noise (constant factor at lags ≥ 1) or aberration (lag-dependent loss). Width: set by the edge defocus phase \((\pi D^{2}/(2\lambda z))(\Delta c/c):\ 1.6\) rad at 40 m/s for F/1.6, giving −4% SLSC but about −30% CF. A broad, flat peak means low sensitivity (large F-number, short lags only) or aberration. A shifted, flattened peak with a screen present (+25 m/s for SLSC in the lab) is the signature of aberration bias; noise alone lowers the peak without moving it.
Follow-up: how to make it local? Estimate per region and invert along paths (tomographic or differentiable beamforming), because the focusing delay integrates slowness along the path.
I4How many snapshots do you need for MVDR on a 64-element array, and what do people actually do in ultrasound? · onsite
\(\hat{R} = XX^{H}/K\) has rank ≤ K, so \(K \ge N = 64\) just to invert it. Reed–Mallett–Brennan: with i.i.d. signal-free snapshots the SINR ratio is \(\mathrm{Beta}(K - N + 2, N - 1)\), mean \((K - N + 2)/(K + 1)\), so \(K = 2N - 3 = 125\) for −3 dB and \(K = N\) gives \(2/(N + 1)\) (−15.1 dB for \(N = 64\)). A pixel offers perhaps 1–20 axial samples, so ultrasound MVDR uses sub-aperture averaging (\(L \le N/2\), multiplying snapshots by \(N - L + 1\)) plus diagonal loading \(\delta = \Delta \cdot \operatorname{tr} (\hat{R})/L\) with \(\Delta \approx 10^{-2}\).
Follow-up: why is it worse than RMB in practice? The snapshots contain the desired echo, so errors in steering or delays make the adaptive weights cancel signal; loading bounds the white-noise gain and restores robustness.
I5Distinguish correlation coefficient, correlation function, covariance matrix and coherence. · screen
The correlation coefficient is one number, \(\operatorname{cov} (x, y)/(\sigma _{x}\sigma _{y})\), mean and scale removed, \(|\rho | \le 1\) by Cauchy–Schwarz (NCC peak in speckle tracking). The correlation function \(R_{xy}(\tau )\) is the same inner product as a function of lag (delay estimation, Kasai's arg R(1)). The covariance matrix collects all pairwise covariances of a vector; it is Hermitian PSD, and \(w^{H}Rw\) is beamformer output power. Spatial coherence C(m) is the per-pair normalized covariance averaged along a diagonal of that matrix, a function of element lag. Spectral coherence \(\gamma ^{2}(f)\) normalizes the cross-spectrum per frequency. All coherences are averages: from one segment or snapshot they are identically 1.
Follow-up: is "correlation matrix" the normalized one? It depends on the field: statistics normalizes to unit diagonal; array processing means \(\mathrm{E}[xx^{H}]\), identical to R for zero-mean IQ.
I6State and sketch the Wiener–Khinchin theorem. Why does it matter for speckle and QUS? · onsite
For a WSS process the PSD is the Fourier transform of the autocorrelation, \(S(f) = \int R(\tau )e^{-j2\pi f\tau}d\tau\). Sketch: \(\mathrm{E}|X_{T}(f)|^{2}/T\) for a record of length T is the transform of the triangularly weighted \(R(\tau )(1 - |\tau |/T)\), which tends to S(f) as \(T \to \infty\). Finite-record version: the periodogram is exactly the DFT of the biased \(\hat{R}(k)\). For speckle, white reflectivity through the pulse gives \(S = N_{0}|P(f)|^{2}\) and \(R = N_{0}R_{pp}\): the speckle autocorrelation is the pulse autocorrelation (FWHM 1.25/B for a Gaussian pulse), which is why spectral QUS divides by a reference-phantom spectrum.
Follow-up: why average periodograms? Each bin of a raw periodogram is exponential with std/mean = 1 whatever the record length; Welch averaging of K segments cuts that by about \(\sqrt{K}\) at the cost of resolution.
I7What is the coherence factor for a point target, for diffuse speckle and for noise, with a uniform 64-element aperture? · onsite
\(\mathrm{CF} = |\sum s_{i}|^{2}/(N\sum |s_{i}|^{2})\). Point at focus: identical channels, \(\mathrm{CF} = 1\). Diffuse speckle (ratio of expectations): \(\sum _{m}(N - |m|)C(m)/N^{2}\) with C the triangle, = \(2/3 + 1/(3N^{2}) = 0.667\). Noise: only lag 0 contributes, \(\mathrm{CF} = 1/N = 0.016\). Speckle at 0 dB channel SNR: (N·0.667 + 1)/(2N) = 0.34.
Follow-up: the per-pixel average of CF in speckle is ≈ 0.57, not 0.667, because \(\mathrm{E}[X/Y] \neq \mathrm{E}[X]/\mathrm{E}[Y]\); check which a paper reports, and whether CF was computed per sample or over an axial kernel.
I8Two speckle kernels have true correlation 0.95, each with 10 dB SNR of independent noise. What correlation do you measure, and how would you get channel SNR from coherence? · screen
Noise adds to the variances but not the covariance: \(\rho _{\mathrm{meas}} = \rho /\sqrt{(1 + 1/\mathrm{SNR}_{1})(1 + 1/\mathrm{SNR}_{2})}\) \(= 0.95 \times 10/11\) = 0.864. The same factor scales spatial coherence at lags ≥ 1, so with the VCZ value known, \(\mathrm{SNR} = C(1)/(C_{\mathrm{VCZ}}(1) - C(1)):\) \(N = 128, C(1) = 0.80\) gives \(\mathrm{SNR} = 4.2\) (6.2 dB). This is lag-one coherence (Long, Bottenus and Trahey 2018).
Follow-up: what inflates C(1) artificially? Crosstalk or common-mode electronic noise, which is correlated between neighbouring channels.
I9A QUS parameter correlates with fat fraction at \(\hat{\rho } = 0.6\) in 12 patients. Give a 95% interval. What changes if each patient contributes 50 correlated ROIs? · screen
Fisher \(z: \operatorname{atanh} 0.6 = 0.693\), se = \(1/\sqrt{N - 3} = 1/3\), interval 0.693 ± 0.653 = [0.040, 1.347] in \(z, \rho\) ∈ [0.04, 0.87]. Barely excludes zero. Adding 50 ROIs per patient does not give \(N = 600\): the ROIs within a patient share the same fat fraction and are strongly correlated, so the effective N stays near 12 (use a mixed model or per-patient averages).
Follow-up: what is \(\operatorname{std} (\hat{\rho })\) at \(\rho = 0.9\) with \(N = 16\) independent pairs? \((1 - 0.81)/\sqrt{15} = 0.049\).
I10Why is speckle IQ circular complex Gaussian, and what are the envelope statistics? Are I and |z| independent? · screen
Each IQ sample sums many phasors with independent, uniform phases: the central limit theorem gives a Gaussian, uniform phases give rotation invariance, so \(\mathrm{E}[z^{2}] = 0, I\) and Q uncorrelated with equal variance, and (jointly Gaussian) independent. Then |z| is Rayleigh with mean/std = \(\sqrt{\pi /(4 - \pi )} = 1.91, |z|^{2}\) is exponential with mean/std = 1, the phase is uniform. I and |z| are uncorrelated (symmetry under I → −I) but dependent \((|z| \ge |I|);\) \(\operatorname{corr} (I^{2}, |z|^{2}) = 1/\sqrt{2}\).
Follow-up: what breaks circularity, and would the envelope show it? I/Q gain or phase imbalance gives \(|P|/R \approx 0.42\) for 20% gain error and 10° skew, but the envelope SNR only moves from 1.91 to 1.82: the envelope is a poor detector.
I11Your measured coherence curve falls below the VCZ triangle. How do you tell noise from aberration? · research
Divide by the known triangle. Noise gives a ratio that is flat for all \(m \ge 1\) (equal to \(\mathrm{SNR}/(1 + \mathrm{SNR})\)) with a discontinuity at lag 0. Aberration gives a ratio that decays with lag, ≈ \(e^{-D_{\varphi}(m)}\) two-way, where \(D_{\varphi}\) is the phase structure function; the lag where it falls gives the screen correlation length and the depth gives the rms phase (one screen realization in the MATLAB lab, \(\sigma = 0.7\) rad, \(\ell = 4\) elements: \(C(4) = 0.57\) against 0.94; the ensemble value \(e^{-D_{\varphi}(4)}\cdot 0.94\) is 0.64). A wrong sound speed gives a similar lag-dependent loss, so sweep \(\hat{c}\) first. A linear phase error gives no loss at all.
Follow-up: what figure of merit would you report for an aberration-correction method? The restored coherence curve (or SLSC, CF) relative to the triangle, at matched SNR, plus point-target resolution.
I12An RF kernel is 5 mm long at 5 MHz, 60% bandwidth, 50 MHz sampling. How many independent samples does it contain? · research
It holds 325 samples but about 16 speckle cells of length c·1.25/(2B) = 0.32 mm. The effective N depends on the statistic: \(N_{\mathrm{eff}} = T/\int \rho _{q}(\tau )d\tau\) with \(\rho _{q}\) the correlation of the averaged quantity. For an intensity mean, \(\int |r_{\mathrm{env}}|^{2} = 0.94/B\) gives 20.8; for a correlation coefficient between two RF kernels, \(\int r^{2} \approx 0.47/B\) gives 41.5 (Monte Carlo 40). So \(\operatorname{std} (\hat{\rho })\) at \(\rho = 0.9\) is ≈ 0.03, between the cell-count estimate 0.05 and the naive sample-count estimate 0.011.
Follow-up: why does the RF correlation coefficient get twice the intensity count? The carrier makes \(r(\tau )^{2}\) oscillate as \(\cos ^{2}\), halving its integral.
07.11MATLAB lab
Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run.
m07_covariance_basics.mCovariance ellipse, principal axes, the two regression lines, Fisher-z Monte Carlo · change Sigma, rho_mc, Npairs · runtime ≈ 1 s
%% m07_covariance_basics.m -- covariance ellipse, principal axes, regression lines, Fisher z
% Purpose : draw a 2D Gaussian with a chosen covariance Sigma, estimate the sample
% covariance, read it as an ellipse (eigenvectors = principal axes), compare
% the three lines through the cloud, then Monte-Carlo the sampling
% distribution of the correlation coefficient and check Fisher's z.
% Change : Sigma (any symmetric positive definite 2x2), Npts, rho_mc, Npairs.
% Expect : rho = 0.40, lambda = 9.966 / 3.034, axis angle 68.1 deg,
% slopes y-on-x 0.600, principal 2.488, x-on-y 3.750,
% 39.3 % inside the 1-sigma ellipse, 86.5 % inside 2-sigma;
% Fisher z std ~ 1/sqrt(N-3) and ~95 % interval coverage.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
rng(1);
%% 1. Draw samples with covariance Sigma
Sigma = [4 2.4; 2.4 9]; % covariance (a.u.^2), same as Self-check Q1
Npts = 2000; % number of samples
Lc = chol(Sigma, 'lower'); % Sigma = Lc*Lc'
X = Lc * randn(2, Npts); % 2 x Npts, zero mean
%% 2. Sample covariance (1/(N-1)) written out, then rho
mu = mean(X, 2);
Xc = X - mu; % implicit expansion
S = (Xc * Xc') / (Npts - 1); % sample covariance
rho_true = Sigma(1,2) / sqrt(Sigma(1,1) * Sigma(2,2));
rho_hat = S(1,2) / sqrt(S(1,1) * S(2,2));
fprintf('true Sigma = [%6.3f %6.3f; %6.3f %6.3f]\n', Sigma');
fprintf('sample S = [%6.3f %6.3f; %6.3f %6.3f]\n', S');
fprintf('rho true %.3f rho_hat %.3f\n', rho_true, rho_hat);
%% 3. Eigen-decomposition = principal axes of the ellipse
[U, Lam] = eig(Sigma);
[lam, order] = sort(diag(Lam), 'descend'); U = U(:, order);
if U(1,1) < 0, U(:,1) = -U(:,1); end % fix sign for a readable angle
theta = atan2(U(2,1), U(1,1)) * 180 / pi; % major-axis angle from x axis (deg)
fprintf('eigenvalues = %.3f, %.3f major axis angle = %.1f deg\n', lam(1), lam(2), theta);
%% 4. The three lines through the cloud
slope_yx = Sigma(1,2) / Sigma(1,1); % E[y|x]: minimises vertical errors
slope_pca = U(2,1) / U(1,1); % principal axis: perpendicular errors (TLS)
slope_xy = Sigma(2,2) / Sigma(1,2); % E[x|y] drawn in the x-y plane
fprintf('slopes: y-on-x %.3f principal %.3f x-on-y %.3f\n', slope_yx, slope_pca, slope_xy);
%% 5. Fraction inside the k-sigma ellipses (chi-square with 2 dof)
d2 = sum(X .* (Sigma \ X), 1); % squared Mahalanobis distance
fprintf('inside 1-sigma: %.1f %% (theory %.1f %%) 2-sigma: %.1f %% (theory %.1f %%)\n', ...
100*mean(d2 <= 1), 100*(1-exp(-1/2)), 100*mean(d2 <= 4), 100*(1-exp(-2)));
figure('Position', [100 100 1000 420]);
subplot(1,2,1);
plot(X(1,1:800), X(2,1:800), '.', 'Color', [0.45 0.6 0.85]); hold on;
t = linspace(0, 2*pi, 200);
for k = [1 2]
E = U * diag(k*sqrt(lam)) * [cos(t); sin(t)];
plot(E(1,:), E(2,:), 'k-', 'LineWidth', 1.5);
end
xr = 8;
plot([-xr xr], slope_yx*[-xr xr], 'r-', 'LineWidth', 1.5);
plot([-xr xr], slope_pca*[-xr xr], 'k--', 'LineWidth', 1.2);
plot([-xr xr], slope_xy*[-xr xr], 'Color', [0 0.6 0.5], 'LineWidth', 1.5);
axis equal; xlim([-8 8]); ylim([-10 10]); grid on;
xlabel('x (a.u.)'); ylabel('y (a.u.)');
title('1\sigma/2\sigma ellipses, principal axis (--), y on x (red), x on y (green)');
%% 6. Fisher z Monte Carlo for the sample correlation coefficient
rho_mc = 0.5; Npairs = 20; Nexp = 5000; % true rho, pairs per experiment, experiments
x = randn(Npairs, Nexp);
y = rho_mc * x + sqrt(1 - rho_mc^2) * randn(Npairs, Nexp);
xc = x - mean(x, 1); yc = y - mean(y, 1);
r = sum(xc .* yc, 1) ./ sqrt(sum(xc.^2, 1) .* sum(yc.^2, 1));
z = atanh(r);
se = 1 / sqrt(Npairs - 3);
covered = abs(z - atanh(rho_mc)) <= 1.96 * se;
fprintf('rho = %.2f, N = %d: mean rho_hat %.4f, std %.4f (delta method (1-rho^2)/sqrt(N-1) = %.4f)\n', ...
rho_mc, Npairs, mean(r), std(r), (1 - rho_mc^2) / sqrt(Npairs - 1));
fprintf('Fisher z: std %.4f (1/sqrt(N-3) = %.4f), 95%% interval coverage %.1f %%\n', std(z), se, 100*mean(covered));
subplot(1,2,2);
edges = linspace(-0.3, 1, 66); centers = (edges(1:end-1) + edges(2:end)) / 2; bw = edges(2) - edges(1);
cnt = hist(r, centers);
bar(centers, cnt / (Nexp * bw), 1, 'FaceColor', [0.6 0.7 0.9], 'EdgeColor', 'none'); hold on;
rr = linspace(-0.29, 0.99, 400);
mz = atanh(rho_mc) + rho_mc / (2 * (Npairs - 1));
pdf_r = exp(-0.5 * ((atanh(rr) - mz) / se).^2) ./ (se * sqrt(2*pi) * (1 - rr.^2));
plot(rr, pdf_r, 'r-', 'LineWidth', 1.8);
plot([rho_mc rho_mc], [0 max(pdf_r)*1.1], 'k--');
xlabel('sample correlation coefficient'); ylabel('density');
title(sprintf('\\rho = %.1f, N = %d: histogram vs Fisher-z density', rho_mc, Npairs));
What to look for: for \(\Sigma = [4\;2.4;\ 2.4\;9]\) it prints \(\rho = 0.400, \lambda = 9.966 / 3.034\), a major axis at 68.1°, slopes 0.600 (y on x), 2.486 (principal) and 3.750 (x on y), and 39.6% / 86.2% of points inside the \(1\sigma / 2\sigma\) ellipses (theory 39.3% / 86.5%). The Fisher block \((\rho = 0.5, N = 20)\) gives \(\operatorname{std} z = 0.2415\) against \(1/\sqrt{17} = 0.2425\) and 95.0% interval coverage.
m07_complex_gaussian_speckle.mRandom-phasor sums are circular Gaussian; Rayleigh, exponential, uniform phase; uncorrelated vs independent; I/Q imbalance · change Nscat, gain_err, skew_deg · runtime ≈ 1 s
%% m07_complex_gaussian_speckle.m -- speckle IQ is circular complex Gaussian
% Purpose : build IQ samples as sums of many random phasors (sub-resolution scatterers),
% then check every property of a circular complex Gaussian: I and Q
% uncorrelated with equal variance, zero pseudo-covariance E[z^2],
% Rayleigh envelope (mean/std = 1.91), exponential intensity (mean/std = 1),
% uniform phase. Then break circularity with I/Q imbalance, and show that
% I and |z| are uncorrelated but dependent.
% Change : Nscat (scatterers per cell), gain_err, skew_deg, A_coh (coherent component).
% Expect : |P|/R ~ 0.01 (1/sqrt(Ns) floor), envelope SNR 1.91, intensity SNR 1.00,
% corr(I,|z|) ~ 0, corr(I^2,|z|^2) ~ 0.707; with imbalance |P|/R ~ 0.42
% while the envelope SNR barely moves.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
rng(1);
%% 1. Random phasor sums (central limit theorem)
Ns = 20000; % number of IQ samples (independent resolution cells)
Nscat = 30; % scatterers per resolution cell
amp = rand(Nscat, Ns); % arbitrary positive amplitudes
phs = 2*pi*rand(Nscat, Ns); % independent uniform phases
z = sum(amp .* exp(1i*phs), 1) / sqrt(Nscat); % 1 x Ns complex samples
z = z / sqrt(mean(abs(z).^2)); % normalise E|z|^2 = 1
I = real(z); Q = imag(z);
%% 2. Second-order statistics
C_IQ = [mean(I.^2) mean(I.*Q); mean(I.*Q) mean(Q.^2)]; % real 2x2 covariance (zero mean)
R = mean(abs(z).^2); % E[z z*]
P = mean(z.^2); % pseudo-covariance E[z z]
fprintf('real covariance of (I,Q) = [%.3f %.3f; %.3f %.3f]\n', C_IQ');
fprintf('R = E|z|^2 = %.3f P = E[z^2] = %.4f %+.4fj |P|/R = %.4f (floor 1/sqrt(Ns) = %.4f)\n', ...
R, real(P), imag(P), abs(P)/R, 1/sqrt(Ns));
%% 3. Envelope, intensity, phase
a = abs(z); Iint = a.^2; ph = angle(z);
fprintf('envelope mean/std = %.3f (Rayleigh sqrt(pi/(4-pi)) = %.3f)\n', mean(a)/std(a), sqrt(pi/(4-pi)));
fprintf('intensity mean/std = %.3f (exponential: 1)\n', mean(Iint)/std(Iint));
fprintf('E|z| = %.3f (theory sqrt(pi)/2 = %.3f)\n', mean(a), sqrt(pi)/2);
%% 4. Uncorrelated is not independent
cc = @(u, v) mean((u-mean(u)).*(v-mean(v))) / (std(u,1)*std(v,1));
fprintf('corr(I, Q) = %.4f corr(I, |z|) = %.4f corr(I^2, |z|^2) = %.3f (theory 1/sqrt(2) = 0.707)\n', ...
cc(I, Q), cc(I, a), cc(I.^2, Iint));
%% 5. Break circularity: I/Q gain and phase imbalance in the demodulator
gain_err = 0.2; skew_deg = 10; A_coh = 0; % fractional gain error, phase skew (deg), coherent mean
Ib = (1 + gain_err) * I + A_coh;
Qb = (1 - gain_err) * (cos(skew_deg*pi/180) * Q + sin(skew_deg*pi/180) * I);
zb = Ib + 1i*Qb; zb0 = zb - mean(zb);
Pb = mean(zb0.^2); Rb = mean(abs(zb0).^2);
fprintf('with imbalance: |P|/R = %.3f envelope mean/std = %.3f phase histogram max/min = %.2f\n', ...
abs(Pb)/Rb, mean(abs(zb))/std(abs(zb)), max(hist(angle(zb), 36)) / min(hist(angle(zb), 36)));
%% 6. Figures
figure('Position', [100 100 1100 360]);
subplot(1,3,1);
plot(I(1:3000), Q(1:3000), '.', 'MarkerSize', 3); hold on;
plot(Ib(1:3000), Qb(1:3000), '.', 'MarkerSize', 3, 'Color', [0.85 0.4 0.2]);
axis equal; axis([-3 3 -3 3]); grid on;
xlabel('I (a.u.)'); ylabel('Q (a.u.)'); title('IQ plane: circular (blue) vs imbalanced (orange)');
subplot(1,3,2);
centers = linspace(0.05, 3.45, 35); bw = centers(2) - centers(1);
cnt = hist(a, centers);
bar(centers, cnt / (Ns*bw), 1, 'FaceColor', [0.6 0.7 0.9], 'EdgeColor', 'none'); hold on;
aa = linspace(0, 3.5, 300);
plot(aa, 2*aa/R .* exp(-aa.^2/R), 'r-', 'LineWidth', 1.8);
xlabel('envelope |z|'); ylabel('density'); title('envelope vs Rayleigh');
subplot(1,3,3);
centers = linspace(-pi + pi/36, pi - pi/36, 36); bw = centers(2) - centers(1);
cnt = hist(ph, centers);
bar(centers, cnt / (Ns*bw), 1, 'FaceColor', [0.6 0.7 0.9], 'EdgeColor', 'none'); hold on;
plot([-pi pi], [1 1]/(2*pi), 'r-', 'LineWidth', 1.8);
xlim([-pi pi]); ylim([0 0.25]);
xlabel('phase arg z (rad)'); ylabel('density'); title('phase vs uniform 1/(2\pi)');
What to look for: \(|P|/R = 0.008\) (noise floor 0.007), envelope mean/std 1.915 (1.913), intensity mean/std \(1.005 (1),\) \(\operatorname{corr} (I, |z|) = -0.003\) but \(\operatorname{corr} (I^{2}, |z|^{2}) = 0.704 (0.707)\). With 20% gain error and 10° skew, |P|/R jumps to 0.419 while the envelope SNR only drops to 1.82.
m07_speckle_autocorr_psd.mWiener–Khinchin on finite records, correlation length 1.25/B, periodogram variance, Neff for different statistics · change f0, frac_bw, win_mm, rho · runtime ≈ 3 s
%% m07_speckle_autocorr_psd.m -- Wiener-Khinchin, speckle correlation length and N_eff
% Purpose : simulate RF speckle lines as white scatterer reflectivity filtered by a
% Gaussian pulse (f0 = 5 MHz, 60 % -6 dB bandwidth, fs = 50 MHz). Check that
% (1) the periodogram is exactly the DFT of the biased autocorrelation
% (Wiener-Khinchin for finite records), (2) the RF autocorrelation envelope has
% FWHM 1.25/B and the intensity autocovariance FWHM 0.88/B, (3) the periodogram
% does not get less noisy with a longer record, and (4) the effective number of
% independent samples in a 5 mm window depends on the statistic: the speckle-cell
% count T*B/1.25, the intensity mean T*B/0.94, the RF correlation coefficient
% T*B/0.47.
% Change : f0, frac_bw, win_mm (analysis window), rho (for the correlation estimate).
% Expect : max |periodogram - DFT(R_biased)| ~ 1e-13 (relative); tau_FWHM*B = 1.25 (RF),
% 0.88 (intensity); axial cell 0.32 mm; periodogram std/mean ~ 1.0 for 256 and
% 4096 samples; N_eff for 5 mm = 15.6 (cells), ~21 (intensity mean, theory 20.8),
% ~40 (rho_hat at rho = 0.9, theory 41.5), out of 325 samples.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
rng(1);
%% 1. Pulse spectrum and speckle lines
c0 = 1540; fs = 50e6; f0 = 5e6; frac_bw = 0.6; % m/s, Hz, Hz, -6 dB fractional bandwidth
B = frac_bw * f0; % -6 dB (amplitude) bandwidth (Hz)
sgt = sqrt(2*log(2)) / (pi*B); % Gaussian envelope std of the pulse (s)
Ns = 8192; Nl = 200; % samples per line, number of lines
f = [0:Ns/2-1, -Ns/2:-1] * fs / Ns; % FFT frequency grid (Hz)
Pf = exp(-2*pi^2*sgt^2*(f - f0).^2) + exp(-2*pi^2*sgt^2*(f + f0).^2); % pulse spectrum (amplitude)
rfl = @(n) real(ifft(fft(randn(Ns, n)) .* Pf.')); % n lines of white reflectivity * pulse (circular)
X = rfl(Nl); X = X / sqrt(mean(X(:).^2)); % Ns x Nl RF, unit power
fprintf('pulse sigma_t %.1f ns, B = %.2f MHz, 1/B = %.0f ns\n', 1e9*sgt, B/1e6, 1e9/B);
%% 2. Wiener-Khinchin for a finite record: periodogram = DFT of the biased autocorrelation
x = X(1:1024, 1); Nx = numel(x);
Rb = zeros(2*Nx-1, 1); % biased R(k), k = -(Nx-1)..(Nx-1)
for k = 0:Nx-1
Rb(Nx+k) = sum(x(1+k:Nx) .* x(1:Nx-k)) / Nx; Rb(Nx-k) = Rb(Nx+k);
end
Lf = 2*Nx - 1; kk = (0:Lf-1).';
Sper = abs(fft(x, Lf)).^2 / Nx; % periodogram on a 2Nx-1 point grid
Swk = real(exp(1i*2*pi*kk*(Nx-1)/Lf) .* fft(Rb)); % DFT of R(k), index shifted to k = -(Nx-1)
fprintf('max |periodogram - DFT(R_biased)| / max = %.1e\n', max(abs(Sper - Swk)) / max(Sper));
%% 3. Correlation lengths: RF envelope (1.25/B) and intensity (0.88/B)
Xf = fft(X); Xf(Ns/2+2:end, :) = 0; Xf(2:Ns/2, :) = 2*Xf(2:Ns/2, :);
Z = ifft(Xf); % analytic signal (IQ at RF), Ns x Nl
Ra = mean(ifft(abs(fft(Z, 2*Ns)).^2), 2); % averaged complex autocorrelation (biased)
env = abs(Ra) / abs(Ra(1)); % RF autocorrelation envelope vs lag
Ia = abs(Z).^2; Ia = Ia - mean(Ia(:)); % intensity fluctuation
RI = mean(ifft(abs(fft(Ia, 2*Ns)).^2), 2); RI = real(RI) / real(RI(1));
fwhm = @(r) 2 * interp1(r(find(r < 0.5, 1) + (-1:0)), (find(r < 0.5, 1) + (-1:0)) - 1, 0.5) / fs; % s
tR = fwhm(env); tI = fwhm(RI);
fprintf('RF autocorrelation FWHM %.0f ns, x B = %.3f (theory 4*sqrt(2)*ln2/pi = %.3f)\n', 1e9*tR, tR*B, 4*sqrt(2)*log(2)/pi);
fprintf('intensity autocov. FWHM %.0f ns, x B = %.3f (theory /sqrt(2) = %.3f)\n', 1e9*tI, tI*B, 4*log(2)/pi);
fprintf('axial speckle cell c*tau/2 = %.3f mm (RF), %.3f mm (intensity)\n', 1e3*c0*tR/2, 1e3*c0*tI/2);
%% 4. The periodogram does not converge: each bin is ~exponential, std/mean ~ 1
for Nw = [256 4096]
P = abs(fft(X(1:Nw, :))).^2 / Nw; % Nw x Nl periodograms
b = round([4.5 5 5.5] * 1e6 / fs * Nw) + 1;
fprintf('periodogram, %4d samples: std/mean across lines at 4.5/5/5.5 MHz = %.2f %.2f %.2f\n', ...
Nw, std(P(b, :), 0, 2) ./ mean(P(b, :), 2));
end
%% 5. Effective number of independent samples in a 5 mm window
win_mm = 5; W = round(2*win_mm*1e-3/c0*fs); T = W/fs; % window: samples, duration (s)
nw = floor(Ns / W);
Iw = reshape(abs(Z(1:nw*W, :)).^2, W, []); % intensity, one column per window
Neff_I = var(abs(Z(:)).^2) / var(mean(Iw, 1)); % sigma_I^2 / var(window mean)
rho = 0.9; Nn = rfl(Nl); Nn = Nn / sqrt(mean(Nn(:).^2)); % independent speckle, unit power
Y = rho*X + sqrt(1 - rho^2) * Nn; % second signal, correlation rho with X
Xw = reshape(X(1:nw*W, :), W, []); Yw = reshape(Y(1:nw*W, :), W, []);
Xw = Xw - mean(Xw, 1); Yw = Yw - mean(Yw, 1);
r = sum(Xw .* Yw, 1) ./ sqrt(sum(Xw.^2, 1) .* sum(Yw.^2, 1));
Neff_r = ((1 - rho^2) / std(r))^2; % from std(rho_hat) = (1 - rho^2)/sqrt(N_eff)
fprintf('5 mm window = %d samples, T*B = %.2f\n', W, T*B);
fprintf('N_eff: speckle cells T/tau_FWHM %.1f | intensity mean %.1f (theory T/(sqrt(2pi) sigma_t) = %.1f) | rho_hat %.1f (theory 2T/(sqrt(2pi) sigma_t) = %.1f)\n', ...
T*B/1.248, Neff_I, T/(sqrt(2*pi)*sgt), Neff_r, 2*T/(sqrt(2*pi)*sgt));
fprintf('std(rho_hat) at rho = %.1f: measured %.4f | from cells %.4f | from samples %.4f\n', ...
rho, std(r), (1 - rho^2)/sqrt(T*B/1.248 - 1), (1 - rho^2)/sqrt(W - 1));
%% 6. Figures
figure('Position', [100 100 1100 380]);
subplot(1,3,1);
t = (0:399) / fs * 1e6; plot(t, X(1:400, 1), 'b-'); grid on;
xlabel('time (\mus)'); ylabel('RF (\sigma units)'); title('speckle RF line (first 8 \mus)');
subplot(1,3,2);
lagk = -80:80; lt = lagk / fs * 1e6;
Rr = real(Ra); Rr = Rr / Rr(1);
idx = mod(lagk, 2*Ns) + 1;
plot(lt, Rr(idx), 'b-', 'LineWidth', 1.4); hold on;
plot(lt, exp(-(lagk/fs).^2 / (4*sgt^2)) .* cos(2*pi*f0*lagk/fs), 'r--', 'LineWidth', 1.2);
plot(lt, env(idx), 'k-'); plot(lt, RI(idx), '-', 'Color', [0 0.6 0.4], 'LineWidth', 1.4);
grid on; xlim([-1.6 1.6]); ylim([-1.05 1.1]);
xlabel('lag \tau (\mus)'); ylabel('normalized correlation');
title('RF autocorr. (blue), theory (red), envelope (black), intensity (green)');
subplot(1,3,3);
Nw = 256;
Xs = reshape(X(1:Nw*floor(Ns/Nw), :), Nw, []); % non-overlapping 256-sample segments
hw = 0.5 - 0.5*cos(2*pi*(0:Nw-1).'/Nw); % Hann window (written out)
Pw = mean(abs(fft(Xs .* hw)).^2, 2); fw = (0:Nw/2) * fs / Nw;
Pth = Pf(1:Ns/2+1).^2; fth = f(1:Ns/2+1);
plot(fw/1e6, 10*log10(Pw(1:Nw/2+1) / max(Pw)), 'b-', 'LineWidth', 1.4); hold on;
plot(fth/1e6, 10*log10(Pth / max(Pth)), 'r--', 'LineWidth', 1.2);
grid on; xlim([0 15]); ylim([-40 3]);
xlabel('frequency (MHz)'); ylabel('PSD (dB re peak)'); title('averaged periodogram vs |P(f)|^2');
What to look for: periodogram and DFT of the biased autocorrelation agree to \(5 \times 10^{-16}; \tau _{\mathrm{FWHM}}\cdot B = 1.246\) for RF (theory 1.248) and 0.879 for intensity (0.883), an axial cell of 0.320 mm; periodogram std/mean ≈ 1 at both 256 and 4096 samples. For a 5 mm window (325 samples) \(N_{\mathrm{eff}}\) is 15.6 speckle cells, 21.8 for an intensity mean (theory 20.8) and 40.0 for \(\hat{\rho }\) (theory 41.5); \(\operatorname{std} (\hat{\rho })\) at \(\rho = 0.9\) is 0.030, against 0.050 from the cell count and 0.011 from the sample count.
m07_vcz_coherence.mExact-geometry channel data from diffuse scatterers: VCZ triangle, noise, phase screen, SLSC, CF, RMB snapshot rule · change Ntx, SNR_dB, sig_phi, ell_phi · runtime ≈ 3 s
%% m07_vcz_coherence.m -- van Cittert-Zernike: channel covariance and coherence vs lag
% Purpose : narrowband channel data from diffuse scatterers in the focal plane of a
% 64-element array (exact spherical geometry, focused transmit and receive).
% Estimate the sample spatial covariance, show it as an image, and compare the
% coherence-vs-lag curve with the VCZ triangle (autocorrelation of the transmit
% aperture). Then add incoherent noise (lag >= 1 scaled by 1/(1+1/SNR)), a
% near-field phase screen (two-way loss ~ exp(-D_phi(m))), compute SLSC and
% the coherence factor, and check the Reed-Mallett-Brennan snapshot rule.
% Change : Ntx (transmit aperture, elements), SNR_dB, sig_phi, ell_phi, M (snapshots).
% Expect : clean C(1), C(16), C(32), C(48) = 0.99, 0.73, 0.48, 0.24 vs the paraxial
% triangle 0.984, 0.750, 0.500, 0.250 (the +-6 mm scatterer window, a stand-in
% for range gating, rounds the kink at lag 0; exact geometry at F/1.6 costs ~0.02);
% 0 dB noise halves lags >= 1; SLSC (lags 1-6) ~ 5.67 clean;
% CF 0.649 (ratio of means; the paraxial model gives 2/3 + 1/(3N^2) = 0.667, the
% exact geometry and scatterer window lower it) and 0.556 (mean of per-snapshot
% ratios; 0.567 in the paraxial model): E[X/Y] differs from E[X]/E[Y];
% Ntx = 32: C(16) = 0.51, C(32) = 0.01 (triangle 0.50, 0), the receive aperture
% only samples the transmit autocorrelation;
% RMB mean SINR loss ~ (K-N+2)/(K+1): 0.50 at K = 2N-3.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
rng(1);
%% 1. Geometry
c0 = 1540; f0 = 5e6; lam = c0/f0; k = 2*pi/lam; % m/s, Hz, m, rad/m
N = 64; pitch = 0.3e-3; z = 30e-3; % elements, m, focal depth (m)
xe = ((0:N-1) - (N-1)/2) * pitch; % element x positions (m), 1 x N
Ntx = 64; % active transmit elements (centred); try 32
tx = abs(((0:N-1) - (N-1)/2)) < Ntx/2; % transmit apodization (uniform), 1 x N
xi = (-6e-3 : 20e-6 : 6e-3); % scatterer grid across the focal plane (m)
Nxi = numel(xi);
M = 2000; % independent speckle realizations
% one-way propagation phasors element -> scatterer, minus the focusing delay to (0,z)
r = sqrt((xe.' - xi).^2 + z^2); % N x Nxi distances (m)
r0 = sqrt(xe.'.^2 + z^2); % N x 1 distances to the focus
G = exp(-1i*k*(r - r0)); % focused one-way kernel, N x Nxi
Pt = (tx * G); % transmit field at each scatterer, 1 x Nxi
B = G .* Pt; % receive kernel weighted by transmit field
in6 = xi(abs(Pt) >= 0.5*max(abs(Pt))); % -6 dB (amplitude) extent of the transmit beam
fprintf('transmit -6 dB beam width %.2f mm (uniform aperture theory 1.21 lambda F# = %.2f mm)\n', ...
1e3*(max(in6) - min(in6)), 1e3*1.207*lam*z/(Ntx*pitch));
%% 2. Channel data, sample covariance, coherence vs lag
A = (randn(Nxi, M) + 1i*randn(Nxi, M)) / sqrt(2); % delta-correlated complex reflectivity
S0 = B * A; % N x M focused channel samples
Ps = mean(abs(S0(:)).^2); % signal power per channel
cohfun = @(S) (S * S') / size(S, 2); % sample covariance R = (1/M) S S^H
nrm = @(R) R ./ sqrt(real(diag(R)) * real(diag(R)).');
lagavg = @(C) arrayfun(@(m) mean(real(diag(C, m))), 0:size(C,1)-1);
tri = max(0, 1 - (0:N-1)/Ntx); % VCZ prediction: transmit autocorrelation
R0 = cohfun(S0); C0 = nrm(R0); coh0 = lagavg(C0);
fprintf('clean: C(1) %.3f C(16) %.3f C(32) %.3f C(48) %.3f (triangle %.3f %.3f %.3f %.3f)\n', ...
coh0([2 17 33 49]), tri([2 17 33 49]));
disp('top-left 4x4 corner of the normalized sample covariance:'); disp(round(100*C0(1:4,1:4))/100);
%% 3. Incoherent noise: lag 0 stays 1, lags >= 1 scale by 1/(1+1/SNR)
SNR_dB = 0; snr = 10^(SNR_dB/10);
Wn = sqrt(Ps/snr) * (randn(N, M) + 1i*randn(N, M)) / sqrt(2);
cohN = lagavg(nrm(cohfun(S0 + Wn)));
fprintf('SNR %g dB: C(1) %.3f, ratio to clean %.3f (theory 1/(1+1/SNR) = %.3f)\n', ...
SNR_dB, cohN(2), cohN(2)/coh0(2), 1/(1 + 1/snr));
%% 4. Near-field phase screen on transmit and receive
sig_phi = 0.7; ell_phi = 4; % rms phase (rad), correlation length (elements)
kk = -15:15; g = exp(-kk.^2 / ell_phi^2); % filter -> corr(m) = exp(-m^2/(2 ell^2))
w = conv(randn(1, N + 30), g, 'same'); phi = w(16:15+N);
phi = phi - mean(phi); phi = sig_phi * phi / sqrt(mean(phi.^2));
Gab = G .* exp(1i*phi.'); % screen at the aperture, both ways
Bab = Gab .* (tx * Gab);
cohA = lagavg(nrm(cohfun(Bab * A)));
avgph = arrayfun(@(m) abs(mean(exp(1i*(phi(1:N-m) - phi(1+m:N)))))^2, 0:N-1);
Dphi = arrayfun(@(m) mean((phi(1:N-m) - phi(1+m:N)).^2), 0:N-1);
fprintf('screen %.1f rad, ell %d el.: lag measured tri*|<e^jdphi>|^2 tri*exp(-D)\n', sig_phi, ell_phi);
for m = [1 2 4 8 16]
fprintf(' %2d %.3f %.3f %.3f\n', ...
m, cohA(m+1), tri(m+1)*avgph(m+1), tri(m+1)*exp(-Dphi(m+1)));
end
%% 5. SLSC (lags 1..6) and coherence factor
Mlag = 6;
slsc = @(coh) sum(coh(2:Mlag+1));
cf_ratio = @(S) mean(abs(sum(S, 1)).^2) / mean(N * sum(abs(S).^2, 1)); % ratio of means
cf_mean = @(S) mean(abs(sum(S, 1)).^2 ./ (N * sum(abs(S).^2, 1))); % mean of ratios
fprintf('SLSC: clean %.2f (ideal %.2f), 0 dB noise %.2f, screen %.2f\n', ...
slsc(coh0), sum(tri(2:Mlag+1)), slsc(cohN), slsc(cohA));
fprintf('CF clean: ratio of means %.3f (paraxial theory %.3f), mean of ratios %.3f\n', ...
cf_ratio(S0), sum((N - abs(-(N-1):(N-1))) .* tri(abs(-(N-1):(N-1)) + 1)) / N^2, cf_mean(S0));
fprintf('CF 0 dB noise %.3f, CF screen %.3f\n', cf_ratio(S0 + Wn), cf_ratio(Bab * A));
%% 6. Reed-Mallett-Brennan: sample-matrix MVDR needs K ~ 2N snapshots for 3 dB
Nr = 32; Rt = toeplitz(1 - (0:Nr-1)/Nr) + 0.01*eye(Nr); % speckle (VCZ) + 20 dB noise
Lr = chol(Rt, 'lower'); s = ones(Nr, 1); opt = real(s' * (Rt \ s));
for K = [Nr, 2*Nr-3, 4*Nr]
loss = zeros(1, 500);
for t = 1:500
X = Lr * (randn(Nr, K) + 1i*randn(Nr, K)) / sqrt(2);
wv = (X*X'/K) \ s; % sample-matrix-inverse weights
loss(t) = abs(wv'*s)^2 / real(wv'*Rt*wv) / opt; % SINR relative to optimum
end
fprintf('RMB N = %d, K = %3d: mean SINR loss %.3f (theory (K-N+2)/(K+1) = %.3f, %.2f dB)\n', ...
Nr, K, mean(loss), (K-Nr+2)/(K+1), 10*log10((K-Nr+2)/(K+1)));
end
%% 7. Figures
figure('Position', [100 100 1050 420]);
subplot(1,2,1);
imagesc(1:N, 1:N, abs(C0)); axis image; colormap(gray); colorbar;
xlabel('element j'); ylabel('element i'); title('|R_{ij}| / sqrt(R_{ii} R_{jj}), clean speckle');
subplot(1,2,2);
plot(0:N-1, tri, 'k--', 'LineWidth', 1.5); hold on;
plot(0:N-1, coh0, 'b-', 'LineWidth', 1.5);
plot(0:N-1, cohN, 'r-', 'LineWidth', 1.2);
plot(0:N-1, cohA, 'Color', [0 0.6 0.4], 'LineWidth', 1.2);
plot(0:N-1, tri .* avgph, ':', 'Color', [0 0.6 0.4], 'LineWidth', 1.5);
grid on; xlim([0 N]); ylim([-0.05 1.05]);
xlabel('element lag m'); ylabel('coherence C(m)');
legend('VCZ triangle', 'clean', sprintf('SNR %g dB', SNR_dB), 'phase screen', 'screen theory', 'Location', 'northeast');
title('spatial coherence vs lag');
What to look for: clean \(C(1), C(16),\) \(C(32), C(48)\) \(= 0.993,\ 0.726,\) \(0.478,\ 0.235\) against the paraxial triangle 0.984, 0.750, 0.500, 0.250 (the finite scatterer window rounds the kink at lag 0, exact geometry costs ≈ 0.02). At 0 dB SNR lag 1 drops to 0.498, ratio 0.501 (theory 0.500). The 0.7 rad screen follows \(\mathrm{tri}\cdot |\langle e^{j\Delta \varphi}\rangle |^{2}\) (lag 4: 0.566 measured, 0.559 predicted). SLSC 5.66 (ideal 5.67), 2.83 at 0 dB, 3.92 with the screen. CF 0.649 as a ratio of means, 0.556 as a mean of per-snapshot ratios. RMB: \(K = 32,\ 61,\ 128\) give 0.061, 0.501, 0.759. Set Ntx = \(32: C(16) = 0.51\) and \(C(32) = 0.01\).
m07_coherence_vs_sound_speed.mToy coherence-based sound-speed estimator: CF and SLSC vs assumed c, clean, with a phase screen, with noise · change c_true, z, sig_phi, SNR_dB · runtime ≈ 2 s
%% m07_coherence_vs_sound_speed.m -- toy coherence-based sound-speed estimator
% Purpose : simulate narrowband focused channel data from diffuse scatterers in a medium
% of sound speed c_true, then re-focus (transmit and receive delays) with an
% assumed speed c_hat swept over 1440-1640 m/s. The coherence factor and the
% short-lag spatial coherence (SLSC) both peak where c_hat = c_true, because any
% residual (quadratic) focusing phase can only lower spatial coherence.
% Repeat with a near-field phase screen and with noise to see the peak drop,
% flatten and (for the screen) shift.
% Change : c_true, z (focal depth), N (aperture), sig_phi, SNR_dB, M (realizations).
% Expect : clean argmax ~ c_true for both metrics; SLSC loses ~4-5 % and CF ~30-36 % at
% |c_hat - c_true| = 40 m/s (F/1.6); the screen lowers and flattens the peak and
% can move it by tens of m/s (here SLSC +25 m/s): aberration biases the estimate.
% Model : narrowband phasors (single frequency f0), scatterers on a fine grid in the
% plane at the true depth z; the beamformer believes the echo comes from
% z_hat = z * c_hat / c_true and focuses there with speed c_hat.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
rng(2);
%% 1. Geometry and data
c_true = 1510; f0 = 5e6; k = 2*pi*f0/c_true; % true speed (m/s), Hz, wavenumber (rad/m)
N = 64; pitch = 0.3e-3; z = 30e-3; % elements, m, true focal depth (m)
xe = ((0:N-1) - (N-1)/2) * pitch; % element positions (m)
xi = -6e-3 : 20e-6 : 6e-3; Nxi = numel(xi); % scatterer grid in the focal plane (m)
M = 300; % speckle realizations (snapshots)
H = exp(-1i*k*sqrt((xe.' - xi).^2 + z^2)); % true one-way propagation, N x Nxi
A = (randn(Nxi, M) + 1i*randn(Nxi, M)) / sqrt(2); % reflectivity, same draws for every c_hat
sig_phi = 0.6; ell = 5; % phase screen: rms (rad), corr. length (el.)
w = conv(randn(1, N + 40), exp(-(-20:20).^2 / ell^2), 'same'); phi = w(21:20+N);
phi = phi - mean(phi); phi = sig_phi * phi / sqrt(mean(phi.^2));
SNR_dB = 10; % channel SNR for the noisy case
Mlag = 6;
lagavg = @(C, mmax) arrayfun(@(m) mean(real(diag(C, m))), 1:mmax);
cvec = 1440:5:1640; nc = numel(cvec);
CF = zeros(3, nc); SL = zeros(3, nc); % rows: clean, screen, noise
Wn = (randn(N, M) + 1i*randn(N, M)) / sqrt(2); % fixed noise draw
%% 2. Sweep the assumed sound speed
for ic = 1:nc
ch = cvec(ic);
zh = z * ch / c_true; % depth the beamformer assigns to the echo
foc = exp(1i*2*pi*f0*sqrt(xe.'.^2 + zh^2)/ch); % focusing phase (delay) per element, N x 1
for cs = 1:3
if cs == 2, Hc = H .* exp(1i*phi.'); else, Hc = H; end
Pt = (foc.' * Hc); % transmit field at the scatterers, 1 x Nxi
S = (foc .* Hc) .* Pt; % receive focusing x propagation x transmit
S = S * A; % N x M focused channel samples
if cs == 3
Ps = mean(abs(S(:)).^2); S = S + sqrt(Ps / 10^(SNR_dB/10)) * Wn;
end
R = (S * S') / M; d = sqrt(real(diag(R))); C = R ./ (d * d.');
coh = lagavg(C, Mlag);
SL(cs, ic) = sum(coh);
CF(cs, ic) = mean(abs(sum(S, 1)).^2) / mean(N * sum(abs(S).^2, 1));
end
end
%% 3. Peak picking with a parabolic fit around the maximum
pfit = @(ym, y0, yp) 0.5 * (ym - yp) / (ym - 2*y0 + yp); % sub-step offset of the vertex
names = {'clean', sprintf('screen %.1f rad', sig_phi), sprintf('SNR %d dB', SNR_dB)};
ideal = sum(1 - (1:Mlag)/N);
i0 = find(cvec == c_true); ip = find(cvec == c_true + 40); im = find(cvec == c_true - 40);
for cs = 1:3
pk = zeros(1, 2); Y = {CF(cs,:), SL(cs,:)};
for q = 1:2
y = Y{q}; ib = min(max(find(y == max(y), 1), 2), nc-1);
pk(q) = cvec(ib) + 5 * pfit(y(ib-1), y(ib), y(ib+1));
end
fprintf('%-16s argmax CF %6.1f m/s SLSC %6.1f m/s | peak CF %.3f SLSC/ideal %.3f | at +-40 m/s: CF %.0f%%/%.0f%%, SLSC %.1f%%/%.1f%% of peak\n', ...
names{cs}, pk(1), pk(2), max(CF(cs,:)), max(SL(cs,:))/ideal, ...
100*CF(cs,im)/max(CF(cs,:)), 100*CF(cs,ip)/max(CF(cs,:)), 100*SL(cs,im)/max(SL(cs,:)), 100*SL(cs,ip)/max(SL(cs,:)));
end
lam = c_true / f0; D = N * pitch;
fprintf('edge defocus phase for 40 m/s, paraxial estimate pi D^2/(2 lambda z) * dc/c = %.2f rad\n', ...
pi * D^2 / (2 * lam * z) * 40 / c_true);
%% 4. Figures
figure('Position', [100 100 1050 400]);
subplot(1,2,1);
plot(cvec, CF(1,:), 'b-', 'LineWidth', 1.6); hold on;
plot(cvec, CF(2,:), 'Color', [0 0.6 0.4], 'LineWidth', 1.6);
plot(cvec, CF(3,:), 'r-', 'LineWidth', 1.6);
plot([c_true c_true], [0 1], 'k--');
grid on; ylim([0 0.75]); xlim([1440 1640]); xlabel('assumed sound speed (m/s)'); ylabel('coherence factor (ratio of means)');
legend(names{1}, names{2}, names{3}, 'true c', 'Location', 'south'); title('CF vs assumed sound speed');
subplot(1,2,2);
plot(cvec, SL(1,:)/ideal, 'b-', 'LineWidth', 1.6); hold on;
plot(cvec, SL(2,:)/ideal, 'Color', [0 0.6 0.4], 'LineWidth', 1.6);
plot(cvec, SL(3,:)/ideal, 'r-', 'LineWidth', 1.6);
plot([c_true c_true], [0 1.05], 'k--');
grid on; ylim([0 1.05]); xlim([1440 1640]); xlabel('assumed sound speed (m/s)'); ylabel('SLSC (lags 1-6) / ideal');
title('SLSC vs assumed sound speed');
What to look for: with \(c_{\mathrm{true}} = 1510\) m/s the clean argmax is 1510.6 m/s (CF) and 1511.0 m/s (SLSC); at \(\pm 40\) m/s CF falls to 64% / 70% of its peak and SLSC only to 94.8% / 96.1%. 10 dB noise lowers both peaks \((\mathrm{CF}\ 0.646 \to 0.584)\) without moving them. A 0.6 rad screen lowers SLSC to 84% of ideal and shifts its argmax to 1535.5 m/s (CF to 1518.8 m/s): aberration biases coherence-based speed estimates. The paraxial edge defocus phase at 40 m/s is 1.69 rad.
07.12Self-check
Q1\(\Sigma = [4\;2.4;\ 2.4\;9]\). Find \(\rho\), the principal-axis angle, and the slopes of the two regression lines.
\(\rho = 2.4/(\sqrt{4}\cdot \sqrt{9}) = 2.4/6\) = 0.40. Eigenvalues: \(\operatorname{tr} = 13, \det\) \(= 36 - 5.76\) \(= 30.24, \lambda\) \(= 6.5 \pm \sqrt{42.25 - 30.24}\) \(= 6.5 \pm 3.47\), so \(\lambda _{1} = 9.97, \lambda _{2} = 3.03\). Principal axis: tan \(2\theta = 4.8/(4 - 9) = -0.96\) with the axis in the quadrant where \(\Sigma _{xy} > 0, \theta\) = 68.1° from the x axis (slope 2.49). Regression y on x: slope \(2.4/4 = 0.60\) (31°). Regression x on y as a line in the plane: slope \(9/2.4 = 3.75\) (75°). The principal axis lies between the two, as it must.
Q2You correlate a QUS parameter with fat fraction in 10 patients and get \(\hat{\rho } = 0.5\). Significant? What is the 95% interval? What if N = 20?
\(t = 0.5\cdot \sqrt{8}/\sqrt{0.75} = 1.63\) with 8 degrees of freedom, two-sided \(p \approx 0.14\): not significant. Fisher: \(z = \operatorname{atanh} 0.5 = 0.549\), se = \(1/\sqrt{7} = 0.378, z\) ± 1.96·se = [−0.192, 1.290], so \(\rho\) ∈ [−0.19, 0.86]. With \(N = 20\): se = \(1/\sqrt{17} = 0.243\), interval [0.074, 1.024] in \(z, \rho \in [0.07,\ 0.77], p \approx 0.025\). The interval is still wide; this is why pilot correlations should be reported with intervals.
Q3For a uniform 64-element aperture focused in diffuse speckle, what is the expected coherence at lag 32, at lag 64, and the ideal SLSC with \(M = 6\) lags?
\(C(m) = 1 - m/N\). Lag 32: 0.50. Lag 64 (= N, beyond the largest pair separation 63): 0; lag 63 gives 1/64. SLSC ideal = \(\sum _{m=1}^{6}(1 - m/64) = 6 - 21/64\) = 5.67. With Hann apodization on transmit the curve is the Hann autocorrelation and lag 32 falls to 0.18.
Q4Channel SNR drops from 30 dB to 0 dB (incoherent electronic noise). What happens to lag-0 and lag-1 coherence? To CF?
Lag 0 is 1 by normalization. Lag \(1: (1 - 1/64)/(1 + 1/\mathrm{SNR})\). At 30 dB, \(0.984/1.001 = 0.983\); at 0 dB, 0.984/2 = 0.492. CF (ratio of expectations) = \((N^{2}\cdot \mathrm{CF}_{s} + N\cdot \sigma _{n}^{2})/(N^{2}(1 + \sigma _{n}^{2}))\) with \(\mathrm{CF}_{s} \approx 0.667\): at 0 dB, \((0.667 + 1/64)/2 = 0.34\). The demo reproduces both within the Monte-Carlo error.
Q5In circular speckle, are I and Q independent? Are I and the envelope |z| independent? Correlated?
I and Q: jointly Gaussian with \(\mathrm{E}[IQ] = 0\) (circularity), therefore independent. I and \(|z|: \operatorname{cov} (I, |z|) = \mathrm{E}[I|z|] = 0\) because the joint distribution is symmetric under I → −I, so they are uncorrelated; but \(|z| \ge |I|\) always, so they are dependent. |z| is not Gaussian, so the Gaussian exception does not apply.
Q6MVDR on a 64-element array with \(K = 40\) axial snapshots per pixel. What is the rank of \(\hat{R}?\) With subarrays of \(L = 32\), how many effective snapshots, and is it enough by the RMB rule?
\(\operatorname{rank} (\hat{R}) \le 40\), so 24 eigenvalues are exactly zero and \(\hat{R}^{-1}\) does not exist. With \(L = 32\) there are \(N - L + 1 = 33\) subarrays, so 40 × 33 = 1320 (correlated) snapshots for a 32×32 matrix, well above \(2L = 64\). Subarray snapshots are not independent, so add diagonal loading (e.g. \(\Delta = 0.01\)) for robustness to steering errors.
Q75 MHz probe, 60% bandwidth, fs = 50 MHz. How many independent samples are in a 5 mm axial window, and what spread do you expect for a \(\hat{\rho }\) computed over it when \(\rho\) = 0.9?
\(B = 3\) MHz, \(\tau _{\mathrm{FWHM}} \approx 1.25/B = 0.42\) µs, axial cell \(c\cdot \tau /2 = 0.32\) mm. 5 mm / 0.32 mm ≈ 16 speckle cells, versus 2·5 mm/1540 m/s·50 MHz ≈ 325 samples. Treating the cells as independent pairs gives \(\operatorname{std} (\hat{\rho }) \approx (1 - 0.81)/\sqrt{15}\) \(\approx 0.05\), a conservative estimate. The exact effective N for a correlation coefficient between RF kernels is \(T/\int r^{2}(\tau )d\tau\) ≈ T·B/0.47 ≈ 41, so \(\operatorname{std} (\hat{\rho }) \approx 0.19/\sqrt{41}\) ≈ 0.03 (Monte Carlo in m07_speckle_autocorr_psd.m: 0.030). Using the sample count, \((1 - 0.81)/\sqrt{324} \approx 0.011\), is wrong by a factor of 3.
Q8A linear phase ramp across the aperture leaves the coherence curve unchanged; a quadratic one lowers it. Why?
A linear phase \(\varphi _{i} = \alpha i\) applied two-way steers the transmit beam and the receive focus by the same angle: \(R_{ij} = e^{j\alpha (i-j)}\cdot g(i - j)\) with g unchanged except for a shift of \(|P_{t}|^{2}\), which only multiplies g by \(e^{-j\alpha (i-j)}\). The two phases cancel, so C(m) is the same triangle; the image shifts laterally instead. A quadratic phase defocuses: the transmit intensity broadens (its FT, the coherence, narrows) and the receive phase difference \(\varphi _{i} - \varphi _{i+m}\) varies along the diagonal, so averaging it reduces |C(m)|. This is why a wrong sound speed, which produces a quadratic residual, is detectable by coherence while a pure steering error is not.
07.13Figures from the literature
Measured coherence curves and coherence images to hold this module's model against. Images are not reproduced here; follow the DOI to the figure. Each entry says what to compare with which equation.
Coherence versus lag: point target, speckle and in vivo thyroid
M. A. Lediju, G. E. Trahey, B. C. Byram, J. J. Dahl, "Short-lag spatial coherence of backscattered echoes: imaging characteristics," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 58(7):1377–1388, 2011. doi:10.1109/TUFFC.2011.1957 (author manuscript: PMC3172134). Fig. 1.
Fig. 1 plots the inter-element correlation of focused RF echoes against receive-element lag for a point target, for a speckle background and for in vivo thyroid tissue. Compare the speckle curve with the triangle of (07.28), and the point-target curve with the rank-one case (identical channels, coherence near 1 out to long lags). Then read the in vivo curve against the speckle curve at lag 1: by (07.29) a constant ratio at lags ≥ 1 is incoherent noise, while a loss that grows with lag points to aberration, (07.30). The lags where speckle and clutter separate most are the ones SLSC sums.
In vivo SLSC image of a thyroid cyst
Same paper, Fig. 10.
A thyroid cyst at about 1.5 cm depth shown as a B-mode image, as SLSC images with short-lag cutoffs Q = 10.4 %, 20.8 % and 31.2 % of the aperture, and as a spatially compounded image. Look inside the cyst: echoes that fill an anechoic region in B-mode carry little coherence beyond lag 0, so the region darkens in SLSC without any change in the tissue brightness that sets B-mode contrast. Note that Q is a fraction of the receive aperture; convert it to elements before comparing with the M = 6 lags of Demo 07.6 (pitfall box in 07.6).
In vivo SLSC of the carotid sheath, and channel-noise robustness
J. J. Dahl, D. Hyun, M. Lediju, G. E. Trahey, "Lesion detectability in diagnostic ultrasound with short-lag spatial coherence imaging," Ultrasonic Imaging 33(2):119–133, 2011. doi:10.1177/016173461103300203. Fig. 8 (in vivo) and Fig. 6 (simulation).
Fig. 8 shows transverse B-mode and SLSC images of the carotid sheath (vagus nerve, jugular vein, arterioles), with SLSC overlaid on B-mode. Look at the vessel lumens: hypoechoic structures whose B-mode interior carries clutter. Fig. 6 is a simulation study: lesion CNR and speckle SNR of B-mode and SLSC as channel noise is added. It is the image-level consequence of (07.29): noise scales every lag ≥ 1 by the same factor, and SLSC never uses lag 0.
Lag-one coherence: the noise model and measured curves in 11 subjects
W. Long, N. Bottenus, G. E. Trahey, "Lag-one coherence as a metric for ultrasonic image quality," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 65(10):1768–1780, 2018. doi:10.1109/TUFFC.2018.2855653 (author manuscript: PMC6378881). Fig. 1 (model) and Fig. 9 (measured).
Fig. 1 is the model of rule 4: theoretical spatial coherence of uniform speckle at channel SNRs of ∞, 6, 0 and −6 dB, with lag-one coherence marked; noise is a delta at lag 0. Fig. 9 shows harmonic coherence curves versus lag measured in a phantom and in 11 subjects at increasing mechanical index. Look at how the lag-1 value climbs with output while the curve shape at longer lags changes less: that is the SNR factor of (07.29), and inverting it gives the channel SNR, \(\mathrm{SNR} = C(1)/(C_{\mathrm{VCZ}}(1) - C(1))\).
07.14Go deeper
- A. Papoulis, S. U. Pillai, Probability, Random Variables and Stochastic Processes, 4th ed., McGraw-Hill, 2002: chapters on two random variables and sequences of random variables (covariance, correlation, Gaussian vectors), and on stationary processes and spectra (Wiener–Khinchin, ergodicity).
- S. M. Kay, Intuitive Probability and Random Processes Using MATLAB, Springer, 2006: covariance matrices, WSS processes and PSD at a gentler pace; and Fundamentals of Statistical Signal Processing, Vol. I: Estimation Theory, Prentice Hall, 1993, for the estimation side (Module 12).
- P. J. Schreier, L. L. Scharf, Statistical Signal Processing of Complex-Valued Data: The Theory of Improper and Noncircular Signals, Cambridge University Press, 2010: pseudo-covariance and circularity.
- J. W. Goodman, Statistical Optics, Wiley, 1985 (2nd ed. 2015): chapter on coherence, including the van Cittert–Zernike theorem. J. W. Goodman, Speckle Phenomena in Optics, Roberts & Co., 2007: first-order speckle statistics.
- H. L. Van Trees, Optimum Array Processing, Wiley, 2002: sample matrix inversion, diagonal loading, spatial smoothing.
- R. A. Fisher, "Frequency distribution of the values of the correlation coefficient in samples from an indefinitely large population", Biometrika 10(4), 1915.
- I. S. Reed, J. D. Mallett, L. E. Brennan, "Rapid convergence rate in adaptive arrays", IEEE Trans. Aerospace and Electronic Systems, 1974.
- R. Mallart, M. Fink, "The van Cittert–Zernike theorem in pulse echo measurements", J. Acoust. Soc. Am. 90(5):2718, 1991. doi:10.1121/1.401867
- R. Mallart, M. Fink, "Adaptive focusing in scattering media through sound-speed inhomogeneities: the van Cittert Zernike approach and focusing criterion", J. Acoust. Soc. Am. 96, 1994.
- M. A. Lediju, G. E. Trahey, B. C. Byram, J. J. Dahl, "Short-lag spatial coherence of backscattered echoes: imaging characteristics", IEEE Trans. Ultrason. Ferroelectr. Freq. Control 58(7):1377–1388, 2011. doi:10.1109/TUFFC.2011.1957
- S. W. Flax, M. O'Donnell, "Phase-aberration correction using signals from point reflectors and diffuse scatterers: basic principles", IEEE Trans. Ultrason. Ferroelectr. Freq. Control 35, 1988.
- C. Kasai, K. Namekawa, A. Koyano, R. Omoto, "Real-time two-dimensional blood flow imaging using an autocorrelation technique", IEEE Trans. Sonics Ultrason. 32, 1985.
- W. Long, N. Bottenus, G. E. Trahey, "Lag-one coherence as a metric for ultrasonic image quality", IEEE Trans. Ultrason. Ferroelectr. Freq. Control 65(10):1768–1780, 2018. doi:10.1109/TUFFC.2018.2855653
- W. F. Walker, G. E. Trahey, "A fundamental limit on delay estimation using partially correlated speckle signals", IEEE Trans. Ultrason. Ferroelectr. Freq. Control 42, 1995.
- H. Cox, R. M. Zeskind, M. M. Owen, "Robust adaptive beamforming", IEEE Trans. Acoust. Speech Signal Process. 35, 1987.
- T.-J. Shan, M. Wax, T. Kailath, "On spatial smoothing for direction-of-arrival estimation of coherent signals", IEEE Trans. Acoust. Speech Signal Process. 33, 1985.
- J.-F. Synnevåg, A. Austeng, S. Holm, "Adaptive beamforming applied to medical ultrasound imaging", IEEE Trans. Ultrason. Ferroelectr. Freq. Control 54, 2007.
- W. Lambert, L. A. Cobus, T. Frappart, M. Fink, A. Aubry, "Distortion matrix approach for ultrasound imaging of random scattering media", Proc. Natl. Acad. Sci. USA 117, 2020.