Report an errorErrata
Part IV · Module 26

Speckle statistics and image-quality metrics

Speckle is not noise to be wished away; it is a random field whose statistics are fixed by the scatterers and the PSF. Know its distributions, its size and what compounding does to it, and you can measure image quality without fooling yourself or a reviewer.

You will be able to

  • Derive the Rayleigh numbers (\(\mathrm{SNR}_0\) = 1.91, intensity contrast 1, log std 5.57 dB, log-mean offset −2.51 dB) and use them as calibration checks
  • Read envelope statistics: K, Rician, Nakagami \(m\), homodyned K, and how many scatterers make speckle fully developed
  • Predict speckle size from the PSF (Siegert relation) and the gain of incoherent vs coherent compounding from k-space overlap
  • Report contrast, CNR (with a stated definition) and gCNR, and know which ones a display map can inflate
  • Set detection thresholds in Rayleigh and heavier-tailed backgrounds (ULM)

Assumes

  • Module 04 · Sampling, Fourier & LTI systems: beam = FT of aperture, two-way PSF
  • Module 07 · Random vectors, covariance & correlation: circular complex Gaussian speckle, correlation functions, van Cittert–Zernike
  • Module 08 · Analytic signal & IQ demodulation: envelope and log compression
  • Module 23 · Beam patterns, k-space and compounding: DAS and plane-wave compounding
Exact formulaclosed form plotted as is (including numerical quadrature of a closed-form pdf)SimulationMonte Carlo or image simulation run in this page (JS) or in the Octave scriptsApprox. modelsimplified physics; assumptions under the badgeSchematicdrawing for intuition, not computed, not to scaleSpeckle images here are complex white scatterer fields filtered in k-space by a separable pulse-echo transfer function (focal plane, narrowband laterally), computed with a 2D FFT in the page. No Field II or k-Wave output is shown.

26.1Fully developed speckle: the numbers that follow from one assumption

Module 07 showed that the IQ sample \(z\) of a resolution cell that holds many randomly placed scatterers is circular complex Gaussian. Write \(\langle I\rangle = \mathrm{E}\lvert z\rvert^2\) for the mean intensity. Everything in this section follows from that one assumption, and every number is a calibration check you can run on a phantom before you trust a pipeline.

\[ \begin{aligned} p_R(R) &= \frac{2R}{\langle I\rangle}\, e^{-R^2/\langle I\rangle}\ \ (\text{Rayleigh}), \qquad \mathrm{E}R = \frac{\sqrt{\pi}}{2}\sqrt{\langle I\rangle} = 0.886\sqrt{\langle I\rangle}, \qquad \mathrm{E}R^2 = \langle I\rangle \\ \mathrm{var}\,R &= \left(1 - \frac{\pi}{4}\right)\langle I\rangle = (0.463)^2\langle I\rangle \;\;\Rightarrow\;\; \mathrm{SNR}_0 = \frac{\mathrm{E}R}{\mathrm{std}\,R} = \sqrt{\frac{\pi}{4 - \pi}} = 1.913 \end{aligned} \]
(26.1)
\(z\)
beamformed IQ sample (circular complex Gaussian)
\(R = \lvert z\rvert\)
envelope, linear units
\(\langle I\rangle\)
mean intensity \(\mathrm{E}\lvert z\rvert^2\)
\(\mathrm{SNR}_0\)
envelope mean over envelope std (the "speckle SNR")
\(p_R(R)\)
probability density of the envelope
\(\mathrm{E}\)
expected value
\(\mathrm{var},\ \mathrm{std}\)
variance and standard deviation
\(\pi\)
3.14159
\(e\)
exponential function
The ratio does not depend on \(\langle I\rangle\), so gain, TGC level and transmit power cancel. Wagner, Smith, Sandrik and Lopez (1983) introduced \(\mathrm{SNR}_0\) = 1.91 as the benchmark for fully developed speckle.
\[ I = R^2:\quad p_I(I) = \frac{1}{\langle I\rangle}\, e^{-I/\langle I\rangle} \;\;\Rightarrow\;\; \mathrm{std}\,I = \langle I\rangle, \qquad C_I = \frac{\mathrm{std}\,I}{\langle I\rangle} = 1 \]
(26.2)
\(I\)
intensity (squared envelope)
\(C_I\)
intensity speckle contrast, std over mean
\(p_I(I)\)
probability density of the intensity (exponential)
\(\mathrm{std}\,I\)
standard deviation of the intensity
\(\langle I\rangle\)
mean intensity, as in (26.1)
Exponential intensity: the brightest speckle in a region of \(n\) independent cells is about \(\ln n\) times the mean (13.8 for \(n = 10^6\)), which matters for thresholds (Section 26.7).

Log compression. Let \(u = \ln(I/\langle I\rangle)\). For an exponential variable, \(\mathrm{E}u = \Gamma'(1) = -\gamma = -0.5772\) and \(\mathrm{var}\,u = \psi'(1) = \pi^2/6\) (\(\gamma\) is Euler's constant, \(\psi'\) the trigamma function). The displayed value \(D = 10\log_{10}(I/\langle I\rangle) = 20\log_{10}(R/\sqrt{\langle I\rangle})\) is \(u\) times \(10/\ln 10 = 4.343\):

\[ \begin{aligned} \mathrm{E}D &= -4.343\,\gamma = -2.51\ \mathrm{dB}, \qquad \mathrm{std}\,D = 4.343\cdot\frac{\pi}{\sqrt{6}} = 5.57\ \mathrm{dB} \\ p_D(D) &= \frac{\ln 10}{10}\, 10^{D/10} \exp\!\left(-10^{D/10}\right) \end{aligned} \]
(26.3)
\(D\)
log-compressed pixel, dB re \(10\log_{10}\langle I\rangle\)
\(\gamma\)
Euler's constant, 0.5772
4.343
\(10/\ln 10\), dB per neper of intensity
\(\mathrm{E}D,\ \mathrm{std}\,D\)
mean and standard deviation of the log-compressed pixel, dB
\(p_D(D)\)
probability density of \(D\)
\(\ln\)
natural logarithm
\(\exp\)
exponential function
\(\pi\)
3.14159
The log-intensity follows a (reversed) Gumbel / Fisher–Tippett law: skewed, with a long tail toward dark values and the mode at 0 dB. Gain shifts \(D\), so the 5.57 dB spread is gain independent. Module 08 quotes the 5.57 dB; the −2.51 dB offset is the one that biases contrast measurements.
Pitfall

All of these numbers hold for the linear envelope of unprocessed, fully developed speckle. Scanner B-mode pixels are log compressed, clipped, spatially filtered and scan converted: measuring \(\mathrm{SNR}_0\) on them gives values far from 1.91 even in a perfect phantom. Compute speckle statistics on beamformed IQ (or its envelope) before display processing.

26.2Beyond Rayleigh: K, Rician, Nakagami and homodyned K

Two departures from fully developed speckle are common in tissue, and they push the statistics in opposite directions.

  • Too few scatterers per cell, or a fluctuating local mean (pre-Rayleigh): heavier tail, more dark and more very bright pixels, \(C_I > 1\), \(\mathrm{SNR}_0 < 1.91\).
  • A coherent component (post-Rayleigh): a specular interface, a regular scatterer arrangement or a strong point target adds a steady phasor \(A\); the spread narrows, \(C_I < 1\), \(\mathrm{SNR}_0 > 1.91\).

A simple generative model contains both. In one resolution cell there are a Poisson number of scatterers (mean \(\bar N\)), each with a random phase and an amplitude \(a\) (\(\mathrm{E}a^2 = 1\), \(\kappa = \mathrm{E}a^4/(\mathrm{E}a^2)^2\): 1 for equal amplitudes, 2 for Rayleigh-distributed amplitudes), plus a coherent phasor \(A\). With diffuse power \(P = \bar N\) and \(r = A^2/P\):

\[ \begin{aligned} z &= A + \sum_{i=1}^{N} a_i\, e^{j\varphi_i}, \qquad N \sim \mathrm{Poisson}(\bar N) \\ \langle I\rangle &= A^2 + P, \qquad \mathrm{var}\,I = P^2 + 2A^2P + \kappa\bar N \\ m &\equiv \frac{\langle I\rangle^2}{\mathrm{var}\,I} = \frac{(1 + r)^2}{1 + 2r + \kappa/\bar N} \end{aligned} \]
(26.4)
\(\bar N\)
mean number of scatterers per cell
\(a_i, \varphi_i\)
scatterer amplitude and uniform random phase
\(\kappa\)
amplitude kurtosis \(\mathrm{E}a^4/(\mathrm{E}a^2)^2\)
\(A\)
coherent (steady) phasor amplitude
\(P\)
diffuse power, \(\bar N\) here
\(r\)
coherent-to-diffuse power ratio \(A^2/P\)
\(m\)
Nakagami parameter, \(1/C_I^2\)
\(z\)
complex envelope sample, as in (26.1)
\(N\)
number of scatterers in the cell (random)
\(i\)
scatterer index
\(j\)
imaginary unit
\(\sim\ \mathrm{Poisson}(\bar N)\)
N is Poisson distributed with mean \(\bar N\)
\(\langle I\rangle,\ I\)
mean intensity and intensity, as in (26.1)
\(\mathrm{var}\)
variance
\(\equiv\)
defined as
Derivation: \(\lvert A + S\rvert^2 = A^2 + \lvert S\rvert^2 + 2A\,\mathrm{Re}\,S\), odd phase moments vanish, \(\mathrm{E}(\mathrm{Re}\,S)^2 = P/2\), and for a Poisson count \(\mathrm{E}\lvert S\rvert^4 = 2P^2 + \kappa\bar N\). Verified by Monte Carlo in the demo below and in the MATLAB lab.
ModelEnvelope pdf (\(\langle I\rangle = 1\) unless noted)\(m = 1/C_I^2\)Physical picture
Rayleigh\(2R\,e^{-R^2}\)1Many scatterers per cell, no coherent part
Rician\((2R/P)\, e^{-(R^2 + A^2)/P}\, I_0(2AR/P)\)\((1+r)^2/(1+2r)\)Diffuse speckle plus one steady phasor
K (Jakeman and Pusey)\((2b/\Gamma(\alpha))\, (bR/2)^{\alpha}\, K_{\alpha-1}(bR)\), \(b = 2\sqrt{\alpha}\)\(\alpha/(\alpha+2)\)Exponential speckle whose local mean fluctuates as Gamma(\(\alpha\)); exact for a negative-binomial number of steps
Nakagami (Shankar 2000)\((2m^m/\Gamma(m))\, R^{2m-1} e^{-mR^2}\)\(m \ge \tfrac12\)Two-parameter fit to the first two intensity moments; covers both sides of Rayleigh
Homodyned K (Dutt and Greenleaf 1994)coherent \(A\) + K-distributed diffuse part (integral form)\((1+r)^2/(1+2r+2/\alpha)\)General model: both departures at once

For equal-amplitude Poisson scatterers the second moment matches K with \(\alpha = 2\bar N\); for Rayleigh-amplitude scatterers \(\alpha = \bar N\). The moment estimator \(\hat m = (\mathrm{E}R^2)^2/\mathrm{var}(R^2)\) is what most Nakagami imaging papers compute in a sliding window.

Demo 26.1

Random phasor sums: scatterer density, coherent component, envelope statistics

SimulationMonte Carlo phasor sums (Poisson count, uniform phases, equal or Rayleigh amplitudes, plus a coherent phasor), 4000 cells per realization, JS. Rayleigh and Nakagami curves are closed forms.
histogram / measured pointRayleighNakagami fit (m̂) · K branchcoherent phasor · Rician branch
Try: \(\bar N = 1\): 37% (\(e^{-1}\)) of the cells are empty and the histogram piles up at zero (K-like, \(m \approx 0.5\)). Raise \(\bar N\) to 10 and 100: \(\hat m\) climbs toward 1 as \(\bar N/(\bar N + 1)\). Now set \(r\) = 10 dB: \(\hat m \approx 5.8\), \(\mathrm{SNR}_0 \approx 4.6\), and the walk is dominated by the coherent arrow. With \(r\) = 0 dB and \(\bar N = 1\) (equal amplitudes) \(m = 4/(1 + 2 + 1) = 1\) exactly, yet \(\mathrm{SNR}_0 \approx 2.09\) and the histogram is not Rayleigh: one moment ratio cannot tell the two departures apart. \(M\) = 4000 cells per realization, so \(\hat m\) scatters by about \(\pm 2/\sqrt{4000} = \pm 0.03\) near \(m = 1\).
Ultrasound link

Envelope-statistics QUS (Nakagami \(m\), homodyned-K parameters) is used to characterize scatterer concentration and organization, for example in liver and breast tissue characterization; Destrempes and Cloutier (Ultrasound Med. Biol., 2010) review the models. These parameters are fitted to the envelope of beamformed RF after depth (TGC, attenuation) normalization: a mean that drifts across the window itself looks like a K-distribution.

26.3How many scatterers make speckle fully developed?

A resolution cell is not a box: each scatterer at \(x_i\) contributes with the PSF weight \(h(x_i - x_0)\). For a Poisson process of density \(n\) (scatterers per mm² in 2D), Campbell's theorem gives the first two intensity moments directly:

\[ \begin{aligned} \langle I\rangle &= n \int \lvert h\rvert^2, \qquad \mathrm{E}I^2 - 2\langle I\rangle^2 = \kappa\, n \int \lvert h\rvert^4 \\ \Rightarrow\;\; C_I^2 &= 1 + \frac{\kappa}{N_c}, \qquad N_c = \frac{n \left(\int \lvert h\rvert^2\right)^2}{\int \lvert h\rvert^4} \end{aligned} \]
(26.5)
\(n\)
scatterer density, mm⁻² (2D)
\(h\)
PSF (amplitude)
\(N_c\)
mean scatterers per correlation cell
\(V_c\)
correlation cell area \((\int\lvert h\rvert^2)^2/\int\lvert h\rvert^4\), mm²
\(\int\)
integral over the PSF support (area, mm²)
\(\mathrm{E}I^{2}\)
second moment of the intensity
\(\langle I\rangle,\ C_I,\ \kappa\)
as in (26.1), (26.2) and (26.4)
\(\lvert\cdot\rvert\)
magnitude
\(N_c\) is the mean number of scatterers per correlation cell of area \(V_c\). So \(m = N_c/(N_c + \kappa)\): the box model of Section 26.2 with \(\bar N \to N_c\).
\[ \begin{aligned} &\text{Gaussian PSF } h = e^{-x^2/2\sigma_x^2 - z^2/2\sigma_z^2}:\quad V_c = 2\pi\sigma_x\sigma_z, \qquad A_6 = 2\ln 2\cdot\pi\sigma_x\sigma_z \\ &\Rightarrow\;\; N_c = \frac{N_{-6\,\mathrm{dB}}}{\ln 2} = 1.44\, N_{-6\,\mathrm{dB}} \end{aligned} \]
(26.6)
\(\sigma_x, \sigma_z\)
lateral and axial PSF standard deviations, mm
\(A_6\)
area of the −6 dB ellipse, mm²
\(N_{-6\,\mathrm{dB}}\)
mean scatterers inside the −6 dB ellipse
\(h\)
Gaussian amplitude PSF, as in (26.5)
\(x,\ z\)
lateral and axial position, mm
\(\ln\)
natural logarithm
\(N_c\)
scatterers per correlation cell, as in (26.5)
\(V_c\)
correlation cell area, as in (26.5), mm²
At 10 scatterers per −6 dB resolution cell, \(N_c\) = 14.4: \(m\) = 0.935 for equal amplitudes, 0.878 for Rayleigh amplitudes; \(\mathrm{SNR}_0\) = 1.88 and 1.84 (Monte Carlo, 4×10⁵ pixels), within 2–4% of 1.91. The "about 10 scatterers per resolution cell" rule (Oosterveld, Thijssen and Verhoef, Ultrasonic Imaging, 1985, by 3D simulation) is consistent with this, provided you say which cell you mean.

Whether you can see the remaining departure is a separate question. The moment estimator \(\hat m\) from \(n\) independent cells has (delta method on \(\hat\mu^2/\hat v\) with exponential moments \(\mu_2 = 1\), \(\mu_3 = 2\), \(\mu_4 = 9\)):

\[ \mathrm{var}\,\hat m \approx \frac{2^2\cdot 1 + 1^2\cdot(\mu_4 - 1) - 2\cdot 2\cdot 1\cdot\mu_3}{n} = \frac{4 + 8 - 8}{n} = \frac{4}{n} \;\;\Rightarrow\;\; \mathrm{std}\,\hat m \approx \frac{2}{\sqrt{n}} \]
(26.7)
\(\hat m\)
moment estimate \(\hat\mu^2/\hat v\) of the Nakagami parameter
\(n\)
independent speckle cells in the window (not pixels)
\(\mu_k\)
central moments of the unit-mean exponential intensity
\(\mathrm{var}\)
variance
\(\mathrm{std}\)
standard deviation
To tell \(m\) = 0.935 from 1 at two standard deviations you need \(2\cdot 2/\sqrt{n} < 0.065\), i.e. \(n > 3800\) independent cells. "Fully developed" is an operational statement about your estimator as much as about the tissue.
Demo 26.2

PSF-weighted Poisson scatterers: m and SNR0 vs scatterers per resolution cell

SimulationMonte Carlo: Poisson scatterers in a ±3σ box, 2D Gaussian PSF weights, random phases, 2000 pixels per density, JS. Solid curve: closed form \(m = N_c/(N_c + \kappa)\).
theory m = Nc/(Nc+κ)Monte Carlo (2D Gaussian PSF)Rayleigh limit · Nakagami SNR(m)
computing…
Model: scatterers uniformly placed (Poisson) in a ±3σ box around the pixel, Gaussian PSF with \(\sigma_z = 0.6\sigma_x\), random phases; 2000 independent pixels per density. The x axis counts scatterers inside the −6 dB ellipse. Dots scatter by about \(\pm 2/\sqrt{2000} = \pm 0.045\) in \(m\).

26.4Speckle size is set by the PSF, through the Siegert relation

Model the IQ image as \(z = h * s\) with delta-correlated scatterer reflectivity \(s\). Its field autocorrelation is \(\mathrm{E}[z(x)\, z^*(x + \Delta)] = \langle I\rangle\, \mu(\Delta)\), where \(\mu\) is the normalized autocorrelation of the PSF, the inverse Fourier transform of \(\lvert H\rvert^2\) (Module 07, Wiener–Khinchin). For circular Gaussian fields the fourth moment factors (Isserlis), which gives the Siegert relation for the intensity:

\[ \mathrm{E}[I(x)\, I(x + \Delta)] = \langle I\rangle^2 \left(1 + \lvert\mu(\Delta)\rvert^2\right) \;\;\Rightarrow\;\; \mathrm{cov}_I(\Delta) = \langle I\rangle^2 \lvert\mu(\Delta)\rvert^2, \qquad \mu = \frac{\mathcal{F}^{-1}\{\lvert H\rvert^2\}}{\int \lvert H\rvert^2} \]
(26.8)
\(\Delta\)
spatial (or time) lag
\(\mu(\Delta)\)
normalized field autocorrelation of the PSF
\(\mathrm{cov}_I(\Delta)\)
intensity autocovariance
\(H\)
pulse-echo transfer function (k-space support of the PSF)
\(\mathrm{E}[\cdot]\)
expected value over the speckle ensemble
\(I(x)\)
intensity at position x
\(x\)
position, mm
\(\langle I\rangle\)
mean intensity, as in (26.1)
\(\lvert\mu\rvert^2\)
squared magnitude of the normalized field autocorrelation
\(\mathcal{F}^{-1}\)
inverse Fourier transform
\(\int\)
integral over the transfer-function support
Speckle size here = FWHM of the normalized intensity autocovariance \(\lvert\mu\rvert^2\). Other papers use the correlation length \(\int\lvert\mu\rvert^2\) or the envelope ACF; state yours.
\[ \begin{aligned} &\text{axial (Gaussian pulse): } \lvert H(f)\rvert \propto e^{-(f - f_0)^2/2s_f^2}, \qquad B = 2\sqrt{2\ln 2}\; s_f \ \ (-6\ \mathrm{dB}\ \text{two-way bandwidth}) \\ &\lvert\mu(\tau)\rvert^2 = e^{-2\pi^2 s_f^2 \tau^2} \;\;\Rightarrow\;\; \mathrm{FWHM}_\tau = \frac{4\ln 2}{\pi B} = \frac{0.883}{B} \;\;\Rightarrow\;\; \Delta z_{\mathrm{sp}} = 0.441\,\frac{c}{B} \end{aligned} \]
(26.9)
\(s_f\)
std parameter of the Gaussian amplitude spectrum, MHz
\(B\)
−6 dB two-way bandwidth, MHz
\(\tau\)
time lag, µs
\(\Delta z_{\mathrm{sp}}\)
axial speckle size (FWHM of \(\lvert\mu\rvert^2\)), mm
\(f,\ f_0\)
frequency and centre frequency, MHz
\(\lvert H(f)\rvert\)
pulse-echo amplitude spectrum
\(\mathrm{FWHM}_\tau\)
full width at half maximum of \(\lvert\mu\rvert^2\) in time, µs
\(c\)
speed of sound, mm/µs
\(\ln\)
natural logarithm
\(\propto\)
proportional to
This equals the FWHM of the pulse envelope itself. 5 MHz, 60% bandwidth (\(B\) = 3 MHz): 0.226 mm.

Lateral (focal plane, narrowband, same rect aperture on transmit and receive): \(\lvert H(u)\rvert\) is a triangle (the coarray) of half-width \(1/(\lambda F_\#)\), so

\[ \begin{aligned} \mathrm{FWHM}\ \text{of}\ \lvert\mu\rvert^2 &= 0.835\,\lambda F_\#\ \ (\text{rect}), \qquad 1.40\,\lambda F_\#\ \ (\text{Hann}) \\ \text{two-way PSF } {-6}\ \mathrm{dB}\ \text{width} &= 0.886\,\lambda F_\#\ \ (\text{rect}), \qquad 1.44\,\lambda F_\#\ \ (\text{Hann}) \end{aligned} \]
(26.10)
\(\lambda\)
wavelength at \(f_0\), mm
\(F_\#\)
F-number (focal depth over aperture width)
\(u\)
lateral spatial frequency, units of \(1/(\lambda F_\#)\)
\(\mathrm{FWHM}\)
full width at half maximum, mm
\(\lvert\mu\rvert^2\)
as in (26.8), here along the lateral axis
\(\text{rect},\ \text{Hann}\)
uniform and Hann receive apodization
\(\text{two-way PSF}\)
pulse-echo point-spread function
Speckle is slightly smaller than the −6 dB PSF. 5 MHz, F#2 (\(\lambda F_\#\) = 0.616 mm), rect: 0.514 mm. Speckle correlation area \(\iint\lvert\mu\rvert^2\) = 0.241 mm × 0.554 mm, i.e. 7.5 independent speckle cells per mm².
Demo 26.3

Speckle image and its autocovariance vs f0, bandwidth, F# and apodization

Simulation256 × 256 complex white scatterer field filtered in k-space by a separable pulse-echo transfer function (Gaussian axial band, apodization coarray laterally), 2D FFT, JS; theory curves are closed forms.
measured intensity autocovariancetheory |μ|²two-way PSF amplitude
Model: 256 × 256 complex white scatterer field (5.1 mm depth × 10.2 mm lateral, 20 × 40 µm pixels) multiplied in k-space by a separable pulse-echo transfer function: Gaussian axial band, coarray of the apodization laterally (focal plane, narrowband). One realization, so measured widths scatter by a few percent. Try: halve the bandwidth and the axial width doubles; F# 1 → 4 scales the lateral width 4× while the axial width does not move; Hann widens the lateral speckle by 1.68×.
Pitfall

Measure speckle size on intensity (or IQ). The log image has a narrower autocovariance (about 0.77 of the intensity width axially in the MATLAB lab) and the envelope about 0.96 of it, because nonlinear maps add higher harmonics of \(\mu\). Scan-converted images also carry the interpolation kernel's correlation.

26.5Compounding: incoherent averaging, coherent summation and decorrelation

Static tissue gives one speckle realization per imaging geometry, so speckle can only be averaged across geometries: sub-apertures (spatial compounding), steering angles, or frequency sub-bands. The question is how many independent looks you really get.

\[ \begin{aligned} S &= \sum_{k=1}^{N} I_k\ \ (\text{equal } \langle I\rangle):\qquad C_S^2 = \frac{\sum_{k,l} \rho_{I,kl}}{N^2} \equiv \frac{1}{N_{\mathrm{eff}}}, \qquad \rho_{I,kl} = \lvert r_{kl}\rvert^2\ \ (\text{Siegert}) \\ &\text{equicorrelated looks: } N_{\mathrm{eff}} = \frac{N}{1 + (N - 1)\rho_I}, \qquad \text{independent looks: } C = \frac{1}{\sqrt{N}} \end{aligned} \]
(26.11)
\(I_k\)
intensity image of look \(k\)
\(N\)
number of looks
\(r_{kl}\)
complex field correlation between looks \(k\) and \(l\)
\(\rho_{I,kl}\)
intensity correlation, \(\lvert r_{kl}\rvert^2\)
\(C_S\)
speckle contrast of the compound
\(N_{\mathrm{eff}}\)
effective number of independent looks
\(S\)
sum (compound) intensity
\(k,\ l\)
look indices
\(\langle I\rangle\)
mean intensity of each look, as in (26.1)
\(\rho_I\)
intensity correlation between looks, as in (26.13)
\(C\)
speckle contrast of independent looks
\(\equiv\)
defined as
Example: \(N\) = 9 looks with \(\rho_I\) = 0.25 give \(N_{\mathrm{eff}}\) = 3, not 9.
\[ \begin{aligned} &\text{averaging envelopes instead: } \mathrm{SNR}_0 = \frac{1.913\, N}{\sqrt{\sum_{k,l} \rho_{A,kl}}} \\ &\rho_A(r) = \frac{\pi}{4 - \pi}\left[ {}_2F_1\!\left(-\tfrac12, -\tfrac12; 1; \lvert r\rvert^2\right) - 1 \right] \approx 0.915\, \lvert r\rvert^2 \end{aligned} \]
(26.12)
\(\rho_{A}\)
envelope correlation of two Rayleigh variables
\({}_2F_1\)
Gauss hypergeometric function
\(\mathrm{SNR}_0\)
speckle SNR of the averaged envelope
\(1.913\)
single-look envelope speckle SNR, as in (26.1)
\(\rho_A,\ r\)
envelope correlation as a function of the field correlation \(r\)
\(N,\ k,\ l\)
number of looks and look indices, as in (26.11)
\(\pi\)
3.14159
Exact, because means and covariances add; the envelope correlation of two Rayleigh variables with field correlation \(r\) is the hypergeometric expression (verified by Monte Carlo: \(r\) = 0.5 gives 0.2326). Independent looks: \(\mathrm{SNR}_0 = 1.91\sqrt{N}\) exactly.

Coherent compounding (summing IQ, as in plane-wave compounding) is different in kind. Any linear combination of jointly circular Gaussian fields is circular Gaussian, so the compounded envelope is Rayleigh again: \(C_I = 1\) and \(\mathrm{SNR}_0\) = 1.91 for every \(N\). What coherent summation buys is a wider k-space support (finer speckle, lower sidelobes and clutter) and \(10\log_{10}N\) of gain over incoherent noise (Module 23).

Decorrelation is k-space overlap. Two looks of the same delta-correlated scatterers through transfer functions \(H_1\) and \(H_2\) have field correlation

\[ \begin{aligned} &r_{12} = \frac{\int H_1 H_2^*}{\sqrt{\int\lvert H_1\rvert^2 \int\lvert H_2\rvert^2}}, \qquad \rho_I = \lvert r_{12}\rvert^2 \\ &\text{sub-aperture } w \text{ translated by } d:\quad r = 1 - 1.5s^2 + 0.75s^3 \ \ \text{for } s = \frac{2d}{w} \le 1 \\ &\phantom{\text{sub-aperture } w \text{ translated by } d:\quad} r = \frac{(2 - s)^3}{4}\ \ \text{for } 1 < s < 2 \\ &\text{plane waves at } \theta_1, \theta_2,\ \text{receive } F_\#:\quad r = \max\!\left(0,\ 1 - \lvert\sin\theta_1 - \sin\theta_2\rvert\, F_\#\right) \end{aligned} \]
(26.13)
\(H_1, H_2\)
pulse-echo transfer functions of the two looks
\(w, d\)
sub-aperture width and translation (tx and rx), mm
\(s\)
coarray shift \(2d/w\)
\(\theta_{1,2}\)
plane-wave steering angles
\(r_{12},\ r\)
complex field correlation between two looks
\(\rho_I\)
intensity correlation, as in (26.11)
\(H_2^*\)
complex conjugate of \(H_2\)
\(\int\)
integral over the transfer-function support
\(F_\#\)
receive F-number, dimensionless
\(\max\)
larger of the two arguments
Sub-aperture: the coarray triangle shifts by \(2d\). Plane waves: the lateral support is a rect of width \(1/(\lambda F_\#)\) centred at \(\sin\theta/\lambda\). \(d = w/4\): \(\rho_I\) = 0.52; \(d = w/2\): 0.0625 (Burckhardt 1978 and Trahey, Smith and von Ramm 1986 discuss translation decorrelation). Plane waves, F# 1.5: \(\Delta\theta\) = 5° gives \(\rho_I\) = 0.76, 10° gives 0.55; complete decorrelation needs \(\Delta\theta \approx 1/F_\#\) rad. Frequency compounding: disjoint sub-bands have zero overlap, hence independent, at \(N\)× the axial speckle size.
Demo 26.4

Incoherent vs coherent compounding: N looks, inter-look correlation, SNR0 vs √N

Simulation256 × 256 IQ images from a k-space filtered scatterer field (2D FFT, JS): equicorrelated looks by mixing, or plane-wave looks through shifted receive rects; theory curves from (26.11)–(26.13).
incoherent (envelope average): theory · measuredindependent-look limit 1.91√Ncoherent (IQ sum)
Model: 256 × 256 IQ images (6.4 mm × 12.8 mm), 5 MHz, 60% bandwidth. Equicorrelated looks are \(z_k = \sqrt{r}\,z_0 + \sqrt{1 - r}\,w_k\). Plane-wave looks image the same scatterers through receive-aperture rects shifted by \(\sin\theta_k/\lambda\) in lateral k-space. Images show \(20\log_{10}(A/\langle A\rangle)\) on a −25…+8 dB scale. Try: \(r = 0\): incoherent \(\mathrm{SNR}_0\) follows \(1.91\sqrt{N}\) while the coherent sum stays at 1.91. Plane waves at F# 1.5: 15 looks at \(\Delta\theta\) = 1° (14° span) give only \(N_{\mathrm{eff}} \approx 1.3\), and even \(\Delta\theta\) = 4° (56° span) gives 2.7. The coherent image gets finer, never smoother. The angle and F# sliders recompute 15 images when released.
Ultrasound link

Plane-wave compounding papers sum angles coherently for resolution and clutter, then sometimes add incoherent compounding of sub-groups to reduce speckle. The k-space overlap rule says the small angular steps used for coherent compounding (Module 23) give strongly correlated looks: an incoherent average of ±10° at F# 1.5 reduces speckle contrast by much less than \(\sqrt{N}\).

26.6Image-quality metrics: contrast, CNR and gCNR

Definitions used in this module (state yours in every paper; reviewers check):

\[ \begin{aligned} C &= 10\log_{10}\frac{\langle I\rangle_{\mathrm{in}}}{\langle I\rangle_{\mathrm{out}}}\ \ \mathrm{dB} \quad (\text{ratio of mean intensities}) \\ \mathrm{CNR} &= \frac{\lvert\mu_{\mathrm{in}} - \mu_{\mathrm{out}}\rvert}{\sqrt{\sigma_{\mathrm{in}}^2 + \sigma_{\mathrm{out}}^2}} \quad \text{on the linear envelope } A = \lvert z\rvert \\ \mathrm{gCNR} &= 1 - \mathrm{OVL}, \qquad \mathrm{OVL} = \int \min\!\left(p_{\mathrm{in}}(v),\, p_{\mathrm{out}}(v)\right) dv \end{aligned} \]
(26.14)
\(\langle I\rangle_{\mathrm{in}}, \langle I\rangle_{\mathrm{out}}\)
mean intensity inside the lesion ROI and in the background ROI
\(\mu, \sigma\)
mean and std of the envelope in each ROI
\(p_{\mathrm{in}}, p_{\mathrm{out}}\)
pdfs (histograms) of the pixel value \(v\) in each ROI
\(\mathrm{OVL}\)
overlap of the two pdfs
\(C\)
lesion contrast, dB
\(\mathrm{CNR}\)
contrast-to-noise ratio, dimensionless
\(\mu_{\mathrm{in}},\ \mu_{\mathrm{out}},\ \sigma_{\mathrm{in}},\ \sigma_{\mathrm{out}}\)
envelope mean and std inside and outside the lesion, linear units
\(A\)
linear envelope \(\lvert z\rvert\), as in (26.1)
\(\mathrm{gCNR}\)
generalized CNR, 0 (identical pdfs) to 1 (disjoint)
\(v\)
pixel value (integration variable)
\(\min\)
smaller of the two pdfs at each \(v\)
Common variants of CNR: the same ratio on intensity, on the dB image, or with \(\sigma_{\mathrm{out}}\) alone in the denominator. Contrast is also often computed as \(20\log_{10}(\mu_{\mathrm{in}}/\mu_{\mathrm{out}})\) of envelope means, which equals \(C\) only when both regions have the same envelope distribution shape.

For two regions of fully developed speckle (any intensity ratio \(r = \langle I\rangle_{\mathrm{out}}/\langle I\rangle_{\mathrm{in}}\); noise inside a cyst is circular Gaussian too) every metric is a closed-form function of \(r\). With \(q = 1/\sqrt{r}\):

\[ \begin{aligned} \mathrm{CNR}_{\mathrm{env}} &= \frac{\frac{\sqrt{\pi}}{2}(1 - q)}{\sqrt{\left(1 - \frac{\pi}{4}\right)(1 + q^2)}}, \qquad \mathrm{CNR}_{\mathrm{int}} = \frac{1 - 1/r}{\sqrt{1 + 1/r^2}}, \qquad \mathrm{CNR}_{\mathrm{dB}} = \frac{\lvert C\rvert}{\sqrt{2}\cdot 5.57} \\ \mathrm{gCNR} &= r^{-1/(r-1)}\left(1 - \frac{1}{r}\right) \qquad \left(\text{exponential pdfs cross at } t = \frac{r\ln r}{r - 1}\right) \end{aligned} \]
(26.15)
\(r\)
intensity ratio \(\langle I\rangle_{\mathrm{out}}/\langle I\rangle_{\mathrm{in}}\) of two speckle regions, \(\ge 1\)
\(q\)
envelope ratio \(1/\sqrt{r}\)
\(\mathrm{CNR}_{\mathrm{env}},\ \mathrm{CNR}_{\mathrm{int}}\)
CNR on the linear envelope and on the intensity
\(\mathrm{CNR}_{\mathrm{dB}}\)
CNR on the log-compressed (dB) image
\(C\)
contrast in dB, as in (26.14)
\(5.57\)
std of a log-compressed speckle pixel, dB, as in (26.3)
\(\mathrm{gCNR}\)
generalized CNR, as in (26.14)
\(t\)
intensity at which the two exponential pdfs cross
\(\ln\)
natural logarithm
\(\pi\)
3.14159
The dB-image CNR grows without bound as the contrast grows; the envelope CNR saturates at 1.913 and the intensity CNR at 1. gCNR is the total-variation distance between the two distributions; an ideal per-pixel classifier with equal priors succeeds with probability \((1 + \mathrm{gCNR})/2\).
Contrast CCNR envCNR intensityCNR dB (pooled σ)CNR dB (\(\sigma_{\mathrm{out}}\))gCNR
−3 dB0.4560.4460.3810.5390.249
−6 dB0.8530.7260.7621.0770.471
−10 dB1.2470.8961.2691.7950.697
−20 dB1.7130.9902.5393.5910.945
anechoic, no noise1.9131.000set by DR: (DR − 2.51)/5.57same1

Why gCNR is immune to the display map: if \(v = T(A)\) with \(T\) strictly monotone, both densities transform with the same Jacobian, \(p_v(v) = p_A(T^{-1}(v))\, \lvert dT^{-1}/dv\rvert\), so the overlap integral is unchanged. Means and variances are not: CNR on a log image depends on the compression, and on the clipping floor (dynamic range). Clipping maps many values to one, so it can lower gCNR (information is lost) but never raise it.

Schematic Cyst ROI definitions used in Demo 26.5 and m26_cyst_cnr_gcnr.mInner disk, excluded ring and background annulus drawn to scale relative to the cyst radius; cell counts use the speckle correlation area 0.24 × 0.55 mm (5 MHz, 60 %, F#2).
R = 2.5 mm 0.8 R 1.3 R 1.9 R inner ROI background annulus ROI rules inner disk r < 0.8 R = 2.0 mm: away from the edge blur excluded ring 0.8 R–1.3 R: PSF leakage and sidelobe ring background annulus 1.3 R–1.9 R (same depth range) Independent speckle cells (7.5 per mm²) inner disk: π·2.0² ≈ 12.6 mm² → ≈ 94 cells annulus: π(4.75² − 3.25²) ≈ 37.7 mm² → ≈ 280 cells count cells, not pixels: ≈ 100 cells → CNR and gCNR scatter two same-statistics regions of this size read gCNR ≈ 0.1 contrast: C = 10 log₁₀(mean Iᵢₙ / mean Iₒᵤₜ) from these two ROIs
Demo 26.5

Cyst phantom: contrast, four CNRs and gCNR as dynamic range and compression change

Simulationk-space filtered speckle image with an attenuated disk (2D FFT, JS), additive circular Gaussian noise, clipping at −DR and the selected display map; metrics from the ROIs drawn above.
inside ROI (r < 0.8 R)background annulus (1.3–1.9 R)overlap = 1 − gCNR
Model: 5 MHz, 60% bandwidth, F#2 rect apodization; 10.2 mm × 12.8 mm; cyst radius \(R\) = 2.5 mm; displayed value \(v = T(\max(A/A_{\mathrm{rms,out}},\ 10^{-\mathrm{DR}/20}))\). gCNR uses 100 bins at quantiles of the pooled data, so it is exactly invariant to strictly monotone maps. Try: anechoic cyst (slider at −50), log map: drag DR from 80 to 20 dB and the displayed CNR swings between about 7 (near 50 dB) and 3 (at 20 dB) while gCNR barely moves (clipping at 20 dB merges the darkest 1% of background pixels with the cyst). Switch to the linear map: CNR drops to about 1.9. Raise the noise to −10 dB: contrast, CNR and gCNR all fall together, because noise really fills the cyst.
Pitfall

ROI placement and size decide the answer. Keep the inner ROI away from the edge (PSF leakage), the outer annulus out of the cyst's sidelobe ring, and count independent speckle cells, not pixels: a 2.5 mm cyst at 7.5 cells/mm² holds about 95 cells inside the 0.8 \(R\) disk. With that few cells, gCNR between two regions of identical statistics reads about 0.1 (MATLAB lab), and \(\mathrm{CNR}_{\mathrm{env}}\) can exceed the 1.91 ceiling by chance.

26.7Detection in speckle and noise backgrounds

Two detection problems recur: is there a bubble (ULM) or a point target at this pixel, and can an observer see this lesion? Both reduce to the tails and moments derived above.

\[ \begin{aligned} &\text{complex Gaussian background (noise or speckle):} \\ &P_{fa} = P(I > T) = e^{-T/\langle I\rangle} \;\;\Rightarrow\;\; T = -\ln P_{fa}\cdot\langle I\rangle \\ &P_{fa} = 10^{-3}:\ T = 6.9\langle I\rangle\ (8.4\ \mathrm{dB}), \qquad P_{fa} = 10^{-6}:\ T = 13.8\langle I\rangle\ (11.4\ \mathrm{dB}) \\ &\text{target of power } S \text{ (Rician):}\quad P_d = Q_1\!\left(\sqrt{2S/\langle I\rangle},\ \sqrt{2T/\langle I\rangle}\right) \end{aligned} \]
(26.16)
\(T\)
intensity threshold
\(P_{fa}\)
per-pixel false-alarm probability
\(S\)
target (bubble) power
\(P_d\)
detection probability
\(Q_1\)
Marcum Q function of order 1
\(\langle I\rangle\)
mean background intensity, as in (26.1)
\(P(I > T)\)
probability that a background pixel exceeds the threshold
\(\ln\)
natural logarithm
\(Q_1(a, b)\)
Marcum Q function; a = \(\sqrt{2S/\langle I\rangle}\) sets the target strength, b = \(\sqrt{2T/\langle I\rangle}\) the threshold
At \(P_{fa} = 10^{-6}\): \(S/\langle I\rangle\) = 10 dB gives \(P_d\) = 0.25, 13 dB gives 0.87, 15 dB gives 0.997.
\[ \begin{aligned} &\text{K-distributed background } (\langle I\rangle = 1): \\ &P(I > T) = \int \mathrm{Gamma}(x; \alpha)\, e^{-T/x}\, dx = \frac{2}{\Gamma(\alpha)}\, (\alpha T)^{\alpha/2}\, K_\alpha\!\left(2\sqrt{\alpha T}\right) \end{aligned} \]
(26.17)
\(\alpha\)
K shape parameter (Gamma-distributed local mean \(x\), unit mean)
\(K_\alpha\)
modified Bessel function of the second kind
\(I,\ T\)
intensity and threshold, in units of \(\langle I\rangle\)
\(\mathrm{Gamma}(x; \alpha)\)
Gamma density of the local mean \(x\), shape \(\alpha\), unit mean
\(x\)
local mean intensity (integration variable)
\(\Gamma(\alpha)\)
gamma function
\(\langle I\rangle\)
global mean intensity, set to 1
At the Rayleigh \(10^{-6}\) threshold, \(\alpha\) = 10, 4 and 1 give false-alarm rates of 3.4×10⁻⁵, 1.9×10⁻⁴ and 2.1×10⁻³: heavy tails from residual clutter or sparse bright scatterers multiply false detections by orders of magnitude.
\[ d' \approx C\sqrt{N_{\mathrm{cells}}}, \qquad C = 1 - \frac{\langle I\rangle_{\mathrm{in}}}{\langle I\rangle_{\mathrm{out}}} \qquad (\text{low contrast, ideal observer}) \]
(26.18)
\(d'\)
detectability index
\(C\)
fractional lesion contrast (here a linear ratio, not dB)
\(N_{\mathrm{cells}}\)
independent speckle cells covered by the lesion
\(\langle I\rangle_{\mathrm{in}},\ \langle I\rangle_{\mathrm{out}}\)
mean intensity inside the lesion and in the background, as in (26.14)
\(\approx\)
valid for \(\lvert C\rvert \ll 1\) (low contrast)
The mean intensity over \(N\) cells has relative std \(1/\sqrt{N}\) (\(C_I\) = 1 per cell). 5 MHz, 60%, F#2 (7.5 cells/mm²): a 3 mm lesion holds 53 cells; 1 dB contrast (\(C\) = 0.206) gives \(d' \approx 1.5\), 3 dB gives 3.6. This is the speckle-limited contrast-detail trade-off analysed by Smith, Wagner, Sandrik and Lopez (1983).
Demo 26.6

Thresholding a speckle or noise background: Pfa, Pd and heavy tails

Exact formulaExponential, Rician and K pdfs and tail probabilities (26.16)–(26.17), evaluated by numerical quadrature of their closed forms.
Rayleigh backgroundK background (selected α)target + Gaussian background (Rician)threshold
Try: set \(T\) = 11.4 dB (\(P_{fa} = 10^{-6}\) per pixel on a Rayleigh background) and switch \(\alpha\): the same threshold lets through 30 to 2000 times more false detections. A ULM frame of 10⁵ pixels and 10⁴ frames contains 10⁹ pixel tests, so even 10⁻⁶ means about a thousand false bubbles before tracking removes them.
Ultrasound link

After SVD clutter filtering (Module 06) the ULM background is thermal noise plus residual tissue and blood signal. Noise power varies with depth (TGC) and across the aperture, so a fixed global threshold is not a fixed \(P_{fa}\): normalize by a local noise estimate first (a constant-false-alarm-rate detector), and check the background histogram against the exponential law before trusting a threshold.

26.8Visual summary

Fifteen questions an expert answers in five seconds, each with the picture to keep in mind. Every curve is computed from the formula in the answer. The slider rescales the two transform-pair cards: a longer pulse or a larger \(\lambda F_\#\) gives bigger speckle and a narrower spectrum.

26.9How experts think about it

The expert model has three layers. First, fully developed speckle is a circular Gaussian field, so its first-order statistics are fixed numbers (1.91, 1, 5.57 dB, −2.51 dB) and any departure is information: about scatterer density, coherent structure, or your processing. Second, its second-order statistics are the PSF: speckle size, decorrelation and compounding gain all follow from overlaps in k-space. Third, an image-quality metric is a statistic of two distributions, and only metrics invariant to the display map (gCNR) or computed on a stated, linear representation can be compared between papers.

The mean of a dB speckle image sits 2.51 dB below the dB of the mean intensity, at any gain.
\(\mathrm{E}D = -10\gamma/\ln 10\) = −2.51 dB\(\mathrm{std}\,D\) = 5.57 dB

\(u = \ln(I/\langle I\rangle)\) of an exponential variable has mean \(\Gamma'(1) = -\gamma = -0.5772\) and variance \(\pi^2/6\). Multiply by 4.343: −2.507 dB and 5.570 dB (26.3); gain adds a constant to both regions. So a contrast taken as a difference of mean dB values equals 10 log of the intensity ratio only when both regions have the same distribution shape. A Rician target at \(r\) = 10 dB has a log mean of −0.41 dB re its own mean intensity: a mean-dB contrast against Rayleigh background is off by 2.1 dB.

Exact formulaRayleigh and Rician log-intensity pdfs, moments by numerical quadrature.
Fully developed speckle on both sides, no clipping at the display floor or ceiling, statistics on unfiltered data.
Envelope SNR above 1.91 means a coherent component; below 1.91 means too few scatterers or a fluctuating mean.
Rayleigh 1.913Rician 10 dB: 4.65K \(\alpha = 1\): 1.27

Rayleigh sits between two families. A steady phasor narrows the relative spread: \(\mathrm{SNR}_0\) = 2.15, 3.13, 4.65 at \(r = A^2/P\) = 0, 6, 10 dB. Sparse scatterers or a varying local mean fatten the tail: K with \(\alpha\) = 10, 4, 1 gives 1.81, 1.68, 1.27. Check the pipeline on a phantom first: with \(n\) independent cells the estimate scatters by about \(1.36/\sqrt{n}\) (±0.04 for 1000 cells), so 1.91 ± 0.04 is a pass.

Exact formulaRician SNR₀ by quadrature of the Rician pdf; K SNR₀ from Gamma-function moments.
Only on the linear envelope before log compression, smoothing or interpolation. Noise is also Rayleigh, and TGC gradients inside the window lower the number. One moment ratio cannot identify the model.
Nakagami m is one over the squared intensity contrast: m = 1 is Rayleigh, m < 1 pre-Rayleigh, m > 1 post-Rayleigh.
\(\hat m = (\mathrm{E}R^2)^2/\mathrm{var}\,R^2\)\(m \ge \tfrac12\)\(m = (1+r)^2/(1+2r+\kappa/\bar N)\)

Nakagami matches the first two intensity moments with a gamma law for \(I\), so \(m = \langle I\rangle^2/\mathrm{var}\,I = 1/C_I^2\). In the Poisson-plus-coherent model one formula covers both sides (26.4): \(\bar N\) = 10 equal scatterers give 0.91; \(r\) = 10 dB with many scatterers gives 5.76. That is why \(m\) maps are popular: one number, monotone in both departures (Shankar, IEEE TUFFC, 2000).

Exact formula
It does not identify the mechanism: homodyned K with (\(r = 0\), \(\alpha = 18\)) and (\(r\) = 3 dB, \(\alpha = 0.4\)) both give \(m = 0.90\), but \(\mathrm{SNR}_0\) 1.85 vs 2.45 and log std 5.67 vs 3.50 dB. Use higher moments or a full homodyned-K fit when the mechanism matters.
A coherent component 10 dB above the diffuse power takes m from 1 to 5.8 and the log-speckle spread from 5.57 to 2.0 dB.
\(m = (1+r)^2/(1+2r)\)0 dB: 1.336 dB: 2.77

\(I = \lvert A + z\rvert^2\) gives \(\langle I\rangle = A^2 + P\) and \(\mathrm{var}\,I = P^2 + 2A^2P\), hence \(C_I^2 = (1 + 2r)/(1 + r)^2\) and \(m \approx r/2\) for large \(r\). Inverting, \(r = m - 1 + \sqrt{m^2 - m}\). Log-image std (MATLAB lab): 5.09, 3.29, 2.00 dB at \(r\) = 0, 6, 10 dB. Specular interfaces and point targets read as post-Rayleigh, and so does an echo from a strong reflector leaking through a sidelobe.

Exact formulaLog-intensity std of the Rician law by numerical quadrature.
Assumes exactly one deterministic phasor plus circular Gaussian diffuse scattering; a handful of strong scatterers is not one coherent phasor.
K-distributed speckle has \(C_I^2 = 1 + 2/\alpha\), and \(\alpha\) is about twice the number of equal scatterers per correlation cell.
\(m = \alpha/(\alpha+2)\)\(\alpha = 2N_c\) (equal)\(\alpha = N_c\) (Rayleigh amplitudes)

Compound model: exponential speckle whose local mean \(x\) fluctuates as Gamma(\(\alpha\), mean 1), \(I = xE\). Then \(\mathrm{E}I^2 = 2\mathrm{E}x^2 = 2(1 + 1/\alpha)\). Poisson counts give \(C_I^2 = 1 + \kappa/N_c\) (26.5), which matches \(\alpha = 2N_c/\kappa\). The tail is what hurts detection: at the Rayleigh \(10^{-6}\) threshold, \(\alpha = 1\) gives a false-alarm rate of 2.1×10⁻³ (26.17).

Exact formulaTail probabilities (26.16) and (26.17), Bessel K by quadrature.
K is exact for a negative-binomial number of scatterers (Jakeman and Pusey); for Poisson counts only the second moment matches. Mean drift from structure or uncorrected TGC also looks K-like.
Ten scatterers per resolution cell is fully developed speckle to within a few percent.
\(N_c = N_{-6\,\mathrm{dB}}/\ln 2\)\(N_{-6\,\mathrm{dB}} = 10\): \(m\) = 0.935\(\mathrm{SNR}_0\) ≈ 1.84–1.88

Campbell's theorem for Poisson scatterers under a PSF gives \(C_I^2 = 1 + \kappa/N_c\), with \(N_c = n(\int\lvert h\rvert^2)^2/\int\lvert h\rvert^4\). For a Gaussian PSF the correlation area \(2\pi\sigma_x\sigma_z\) is \(1/\ln 2\) = 1.44 times the −6 dB ellipse (26.6). Ten per −6 dB cell: \(N_c\) = 14.4, \(m\) = 0.935 (equal amplitudes) or 0.878 (Rayleigh amplitudes), \(\mathrm{SNR}_0\) 1.88 or 1.84, within 2–4% of 1.91 (Demo 26.2).

Exact formula
The count depends on the cell definition (a factor \(\ln 2\)) and on the amplitude spread \(\kappa\); in 3D the cell includes the elevation slice; clustered scatterers act as fewer, stronger ones.
A Nakagami estimate from n independent speckle cells has a standard deviation of about 2/√n.
\(\mathrm{std}\,\hat m \approx 2/\sqrt{n}\)±0.1 needs \(n\) ≈ 400≈ 53 mm² at 7.5 cells/mm²

Delta method on \(\hat m = \hat\mu^2/\hat v\) for exponential intensity: gradient (2, −1), \(\mathrm{var}\,\hat\mu = 1/n\), \(\mathrm{var}\,\hat v = 8/n\), cov = \(2/n\), so \(\mathrm{var}\,\hat m = (4 + 8 - 8)/n = 4/n\) (26.7). Monte Carlo: \(\mathrm{std}\cdot\sqrt{n}\) = 2.13, 2.01, 1.98, 1.91 for \(n\) = 25 to 1600, and \(\hat m\) is biased upward at small \(n\) (1.22 at \(n\) = 25). At 5 MHz, F#2 a ±0.1 map needs a 7 mm window: Nakagami images trade resolution for precision.

SimulationLine: 2/√n. Dots: Monte Carlo, 400 trials of n exponential samples, JS.
Counts independent cells, not pixels; \(m < 1\) (heavier tails) inflates the variance further; windowed estimates near boundaries mix two populations.
Speckle size is the FWHM of \(\lvert\mu\rvert^2\): about \(0.44\,c/B\) axially and \(0.84\,\lambda F_\#\) laterally, slightly smaller than the −6 dB PSF.
axial \(0.441\,c/B\)lateral \(0.835\,\lambda F_\#\) (rect), 1.40 (Hann)5 MHz, 60%, F#2: 0.23 × 0.51 mm

Siegert: for circular Gaussian fields the intensity autocovariance is \(\langle I\rangle^2\lvert\mu\rvert^2\), \(\mu = \mathcal{F}^{-1}\{\lvert H\rvert^2\}\) normalized (26.8). Gaussian pulse: \(\lvert\mu\rvert^2 = \exp(-2\pi^2 s_f^2\tau^2)\), FWHM \(4\ln 2/(\pi B) = 0.883/B\) in time (26.9). Rect aperture, two-way: \(\lvert H\rvert\) is a triangle and the FWHM of \(\lvert\mu\rvert^2\) is \(0.835\,\lambda F_\#\), against \(0.886\,\lambda F_\#\) for the two-way PSF at −6 dB (26.10). This sharpens Module 07's "speckle is the PSF size" into numbers.

Exact formula\(\lvert\mu\rvert^2\) from the transform of the squared coarray triangle; PSF = sinc² (dashed).
Focal plane, narrowband laterally, delta-correlated scatterers. Off focus the lateral size follows the local beam width; log images and scan conversion change the apparent size.
N independent looks cut intensity speckle contrast by √N; correlated looks count as \(N_{\mathrm{eff}} = N^2/\sum\rho_I\).
\(C = 1/\sqrt{N_{\mathrm{eff}}}\)\(N/(1 + (N-1)\rho_I)\)\(\rho_I = \lvert r\rvert^2\)

The variance of a sum is the sum of all covariances: \(\mathrm{var}\,S = \langle I\rangle^2\sum_{kl}\rho_{I,kl}\), mean \(N\langle I\rangle\) (26.11). Siegert makes the intensity correlation the squared field correlation, so a field correlation of 0.5 is only 0.25 in intensity. Nine looks at \(\rho_I\) = 0.25 are worth three. For envelope averaging use \(\rho_A \approx 0.915\lvert r\rvert^2\): \(\mathrm{SNR}_0 = 1.913\,N/\sqrt{\sum\rho_A}\) (26.12).

Exact formula
Equal mean intensities and registered looks; motion between looks blurs the compound; each look usually has worse resolution than the full aperture or band.
Coherent compounding never reduces speckle contrast: it makes speckle finer, not smoother.
\(\mathrm{SNR}_0\) = 1.91 for any \(N\)\(C_I = 1\)noise gain \(10\log_{10}N\)

A weighted sum of jointly circular Gaussian fields is circular Gaussian, so the envelope is Rayleigh again. What changes is the k-space support: wider support means smaller speckle cells, lower sidelobes and clutter, and incoherent noise down by \(10\log_{10}N\) relative to the signal. In the MATLAB lab the coherent \(\mathrm{SNR}_0\) stays at 1.90–1.94 for every \(N\) while the lateral speckle width shrinks from 25 to 1.4 pixels.

Exact formula
Holds for fully developed speckle. If the single-look image was noise- or clutter-dominated, coherent summation changes the measured contrast through those terms. Motion between angles makes the sum partly incoherent.
Two speckle patterns are as correlated as their k-space supports overlap.
\(r = \int H_1H_2^*/\lVert H_1\rVert\lVert H_2\rVert\)shift \(w/2\): \(\rho_I\) = 0.06plane waves: \(r = 1 - \lvert\Delta\sin\theta\rvert F_\#\)

Same delta-correlated scatterers through two linear systems: \(\mathrm{E}[z_1z_2^*] \propto \int H_1H_2^*\) (26.13). Translating a sub-aperture of width \(w\) by \(d\) shifts its triangular coarray by \(2d\): \(d = w/4\) gives \(\rho_I\) = 0.52, \(d = w/2\) gives 0.0625. Steering a plane wave shifts the receive rect by \(\sin\theta/\lambda\): at F# 1.5, 10° gives \(\rho_I\) = 0.55. The same rule is the van Cittert–Zernike theorem of Module 07 seen from the other side.

Exact formula
Delta-correlated scatterers, focal-zone narrowband geometry, no motion or aberration. Elevation and bandwidth add dimensions to the overlap.
Speckle reduction is resolution traded away: N disjoint sub-bands give √N lower contrast and N times the axial speckle size.
\(B \to B/N\) per bandaxial cell × \(N\)\(C \to 1/\sqrt{N}\)

Disjoint bands have zero k-space overlap, so the looks are independent, but each band's axial speckle and resolution scale as \(0.441c/(B/N)\). Spatial compounding with sub-apertures of width \(w < D\) pays laterally in the same way. Static tissue offers only one realization per geometry, so diversity in angle, aperture or frequency is the only way to average, and the lesion-detectability gain (\(C\sqrt{N_{\mathrm{cells}}}\)) partly cancels against the larger cells.

SchematicIdealized rectangular sub-bands of a 6 MHz band; real filters overlap and give fewer effective looks.
Partially overlapping bands give fewer effective looks; nonlinear despeckling filters change the statistics and invalidate CNR comparisons.
For an anechoic cyst in fully developed speckle, CNR tops out at 1.91 on the envelope and 1 on intensity; on the dB image your dynamic range sets it.
env 1.913intensity 1.000dB: ≈ (DR − 2.5)/5.57 = 10.3 at 60 dB

Inside: zero signal. Envelope: \(\mathrm{CNR} = \mu_{\mathrm{out}}/\sigma_{\mathrm{out}} = \mathrm{SNR}_0\) = 1.913; intensity: 1. On the dB image the inside sits at the display floor −DR with zero variance, while the background has mean −2.51 dB and std 5.57 dB: \(\mathrm{CNR} = (\mathrm{DR} - 2.51)/5.57\), i.e. 4.9 at 30 dB and 10.3 at 60 dB, set by the display, not the beamformer. A measured envelope CNR above 1.91 on a speckle background is sampling noise (1.955 in the MATLAB lab) or a non-Rayleigh background.

Approx. modeldB-image line uses the unclipped background moments, \((\mathrm{DR} - 2.51)/5.57\); the exact clipped value is in the visual-summary card.
Noise and clutter fill real cysts; partially developed or post-Rayleigh backgrounds can exceed 1.91; small ROIs scatter.
If a display map can change your metric, someone can tune it: gCNR is invariant to every strictly monotone transform, CNR is not.
gCNR = 1 − OVL\(P_{\mathrm{correct}} = (1 + \mathrm{gCNR})/2\)−20 dB cyst: CNR 2.60 → 3.03, gCNR 0.944

The overlap of two densities is unchanged when both are pushed through the same strictly monotone map (same Jacobian); means and variances are not. MATLAB lab, −20 dB cyst: log-display CNR 2.60 at 80 dB DR and 3.03 at 20 dB; power law 0.3: 2.61; linear envelope: 1.76. gCNR is 0.944 in every case (Rodriguez-Molares et al., IEEE TUFFC, 2020). Report gCNR next to any CNR, and the CNR's representation.

Exact formulaCNR of \(I^p\) for two exponential regions from Gamma-function moments; gCNR from (26.15).
Clipping is not injective and can only lower gCNR. gCNR ignores spatial structure (edge sharpness, resolution) and says nothing about lesions smaller than a few speckle cells.
Between two fully developed speckle regions gCNR depends only on their intensity ratio: 6 dB gives 0.47, 20 dB gives 0.94.
\(\mathrm{gCNR} = r^{-1/(r-1)}(1 - 1/r)\)3 dB: 0.2510 dB: 0.70

Two exponential densities with means 1 and \(r\) cross once at \(t = r\ln r/(r - 1)\); \(\mathrm{gCNR} = F_{\mathrm{in}}(t) - F_{\mathrm{out}}(t)\) gives the formula (26.15). Noise inside a hypoechoic cyst is circular Gaussian too, so the formula applies with the measured ratio. Consequence: if a new beamformer's gCNR equals this curve at its measured contrast, the gain is contrast alone; values above the curve mean the distributions changed shape (for example, non-Rayleigh output of a nonlinear beamformer).

Exact formula
Both regions Rayleigh; breaks for point targets, Rician or K statistics and for clipped displays.
gCNR is biased upward: two identical speckle regions of about 100 cells read about 0.1, not 0.
≈ 85 cells: 0.09–0.12grows with bin count

The histogram overlap of two finite samples falls below 1 even when the parent densities are identical, because the minimum of two noisy bin estimates is biased low. In the MATLAB lab two halves of the same background annulus (9400 pixels, about 85 independent cells each) give gCNR = 0.089 with 10 bins and 0.118 with 100 bins. Fix the bin rule, count cells rather than pixels, and report the null value from two regions of the same statistics.

Estimator dependent (bins, kernel density, quantile edges); spatially correlated pixels make the effective sample far smaller than the pixel count.
Lesion detectability grows as contrast times the square root of the number of speckle cells in the lesion.
\(d' \approx C\sqrt{N_{\mathrm{cells}}}\)3 mm, 1 dB: \(d' \approx 1.5\)5 mm, 1 dB: 2.5

Averaging intensity over \(N\) independent cells gives a relative std of \(1/\sqrt{N}\) (\(C_I = 1\) per cell), so a mean-intensity difference \(C\langle I\rangle\) stands out by \(C\sqrt{N}\) (26.18). With \(C = 1 - 10^{-0.1}\) = 0.206 (1 dB) and 7.5 cells/mm², a 3 mm lesion (53 cells) gives \(d' \approx 1.5\). Halving the speckle cell area gains \(\sqrt{2}\) at the same contrast; this is the speckle-limited contrast-detail curve (Smith, Wagner and colleagues, 1983).

Approx. modelLow-contrast ideal-observer approximation (26.18), 7.5 cells/mm².
Low contrast, lesion much larger than a speckle cell, ideal observer who knows the location; human observers do worse; edge blur matters for small lesions.
Threshold a Rayleigh background at \(-\ln P_{fa}\) times its mean power: 11.4 dB for a false-alarm rate of \(10^{-6}\).
\(P_{fa} = e^{-T/\langle I\rangle}\)\(10^{-3}\): 8.4 dB\(P_d\) = 0.87 at 13 dB SNR

A Rayleigh envelope means an exponential intensity: \(P(I > T) = e^{-T/\langle I\rangle}\), so \(10^{-6}\) needs \(T = 13.8\langle I\rangle\) = 11.4 dB (26.16). A bubble of power \(S\) on that background is Rician and \(P_d = Q_1(\sqrt{2S/\langle I\rangle},\ \sqrt{2T/\langle I\rangle})\): 0.25 at 10 dB, 0.87 at 13 dB, 0.997 at 15 dB. Residual clutter is heavier-tailed: a K background with \(\alpha = 1\) gives 2.1×10⁻³ at the same threshold.

Exact formulaRician tail probability by numerical quadrature.
Independent pixels and a known, stationary background power; local noise normalization is needed when \(\langle I\rangle\) varies with depth; tracking changes the effective false-alarm rate.

26.10Cheat sheet

ConceptFormulaUltrasound use
Rayleigh envelope\(\mathrm{SNR}_0 = \sqrt{\pi/(4-\pi)}\) = 1.913; \(\mathrm{E}R = 0.886\sqrt{\langle I\rangle}\)Fully developed speckle check
Exponential intensity\(C_I = 1\); \(P(I > T) = e^{-T/\langle I\rangle}\)Speckle contrast, detection thresholds
Log specklemean −2.51 dB, std 5.57 dB re \(10\log\langle I\rangle\)Contrast bias, display calibration
Nakagami \(m\)\(m = (\mathrm{E}R^2)^2/\mathrm{var}\,R^2 = 1/C_I^2\); \(\mathrm{std}\,\hat m \approx 2/\sqrt{n}\)Envelope-statistics QUS maps
Poisson + coherent\(m = (1+r)^2/(1+2r+\kappa/\bar N)\)Density vs coherent structure
Rician\(m = (1+r)^2/(1+2r)\), \(r = A^2/P\)Specular/point targets, bubble detection
K\(C_I^2 = 1 + 2/\alpha\); \(\alpha \approx 2N_c\) (equal scatterers)Sparse scatterers, heavy-tailed clutter
Scatterers per cell\(N_c = n(\int\lvert h\rvert^2)^2/\int\lvert h\rvert^4 = N_{-6\,\mathrm{dB}}/\ln 2\) (Gaussian PSF)"≥ 10 per cell" rule, stated precisely
Siegert relation\(\mathrm{cov}_I(\Delta) = \langle I\rangle^2\lvert\mu(\Delta)\rvert^2\), \(\mu = \mathcal{F}^{-1}\{\lvert H\rvert^2\}\)Speckle size, decorrelation in intensity
Speckle sizeaxial \(0.441\,c/B\); lateral \(0.835\,\lambda F_\#\) (rect), \(1.40\,\lambda F_\#\) (Hann)Cells per ROI, kernel sizes
Incoherent compounding\(N_{\mathrm{eff}} = N^2/\sum\rho_I\); \(C = 1/\sqrt{N_{\mathrm{eff}}}\)Spatial, angular, frequency compounding
Coherent compounding\(C_I = 1\), \(\mathrm{SNR}_0\) = 1.91 for any \(N\)Plane-wave compounding
k-space overlap\(r = \int H_1H_2^*/\lVert H_1\rVert\lVert H_2\rVert\); plane waves \(1 - \lvert\Delta\sin\theta\rvert F_\#\)Look decorrelation, tracking decorrelation
CNR (this module)\(\lvert\mu_{\mathrm{in}} - \mu_{\mathrm{out}}\rvert/\sqrt{\sigma_{\mathrm{in}}^2 + \sigma_{\mathrm{out}}^2}\) on the envelope; anechoic ceiling 1.91Lesion contrast reporting
gCNR\(1 - \int\min(p_{\mathrm{in}}, p_{\mathrm{out}})\); two Rayleigh regions \(r^{-1/(r-1)}(1 - 1/r)\)Display-invariant lesion separability
Detectability\(d' \approx C\sqrt{N_{\mathrm{cells}}}\)Contrast-detail analysis

26.11Interview and qualifying-exam questions

I1Derive the speckle SNR of 1.91. What do you need to assume? · screen

Many scatterers with independent uniform phases per cell: by the central limit theorem \(z\) is circular complex Gaussian with \(\mathrm{E}\lvert z\rvert^2 = \langle I\rangle\). Then \(R = \lvert z\rvert\) is Rayleigh, \(\mathrm{E}R = (\sqrt{\pi}/2)\sqrt{\langle I\rangle}\), \(\mathrm{E}R^2 = \langle I\rangle\), \(\mathrm{var}\,R = (1 - \pi/4)\langle I\rangle\), so \(\mathrm{SNR}_0 = \sqrt{\pi/4}/\sqrt{1 - \pi/4} = \sqrt{\pi/(4-\pi)}\) = 1.913 (26.1). Intensity is exponential with contrast 1.

Follow-up: what is it on the intensity, and on the log image? Intensity mean/std = 1; log image std 5.57 dB, independent of gain.

I2Why is the standard deviation of log-compressed speckle 5.57 dB whatever the gain, and what is its mean? · screen

For exponential \(I\), \(u = \ln(I/\langle I\rangle)\) has a fixed (Gumbel) distribution with mean \(-\gamma\) and variance \(\pi^2/6\). Gain multiplies \(I\), which adds a constant to \(u\). Times \(10/\ln 10\): std = \(4.343\cdot\pi/\sqrt{6}\) = 5.57 dB, mean = −2.51 dB re \(10\log_{10}\langle I\rangle\) (26.3).

Follow-up: so how should you compute lesion contrast from a dB image? From mean intensities (linearize first), or accept the −2.51 dB offset cancels only if both regions are Rayleigh.

I3You measure \(\mathrm{SNR}_0\) = 2.4 in a liver ROI. List the explanations in the order you would check them. · onsite

(1) Processing: was it computed on the linear envelope, before log compression, spatial filtering, persistence or scan conversion? Smoothing raises \(\mathrm{SNR}_0\). (2) Compounding: incoherent compounding raises it by \(\sqrt{N_{\mathrm{eff}}}\). (3) Physics: a coherent component (specular structures, regular spacing), Rician-like; \(\mathrm{SNR}_0\) = 2.4 matches a Rician \(r \approx 1.75\) (2.4 dB), i.e. \(m \approx 1.7\) (invert a measured \(\hat m\) with \(r = m - 1 + \sqrt{m^2 - m}\)). (4) Check the same pipeline on a phantom: it must read \(1.91 \pm 1.36/\sqrt{n}\).

Follow-up: how would you tell coherent structure from smoothing? Smoothing widens the autocovariance beyond \(\lvert\mu\rvert^2\) of the PSF; a coherent component does not.

I4What is the speckle size at 5 MHz, 60% bandwidth, F#2, and how many independent cells are in a 10 mm × 10 mm ROI? · onsite

Axial FWHM of \(\lvert\mu\rvert^2\) = \(0.441\,c/B\) = 0.441 × 1540/3×10⁶ = 0.226 mm. Lateral (rect, focal zone) \(0.835\,\lambda F_\#\) = 0.835 × 0.308 × 2 = 0.514 mm. Correlation area \(\iint\lvert\mu\rvert^2\) = 0.241 × 0.554 mm², i.e. 7.5 cells/mm², so about 750 cells in 1 cm². \(\mathrm{SNR}_0\) from that ROI scatters by about \(1.36/\sqrt{750}\) = 0.05.

Follow-up: what changes off focus? The lateral size follows the local beam width, so the count falls with distance from focus.

I5Why does coherent plane-wave compounding not reduce speckle, while spatial compounding does? · screen

Coherent compounding sums IQ: the result is a linear combination of circular Gaussian fields, hence circular Gaussian, so \(C_I = 1\) and \(\mathrm{SNR}_0\) = 1.91. It widens the k-space support, so speckle gets finer and resolution and clutter improve. Spatial compounding averages detected images from looks with low k-space overlap, so contrast falls as \(1/\sqrt{N_{\mathrm{eff}}}\) (26.11).

Follow-up: does coherent compounding change gCNR of a cyst? Yes, through contrast: lower sidelobes and clutter mean less energy in the cyst, so a higher intensity ratio; at a given ratio, gCNR is the same function of \(r\).

I6Three incoherently compounded plane-wave angles at −10°, 0°, +10°, receive F# 1.5. How many effective looks? · onsite

Field correlation \(r = 1 - \lvert\Delta\sin\theta\rvert\cdot F_\#\) (26.13): adjacent pairs \(\Delta\sin\theta\) = 0.174, \(r\) = 0.74, \(\rho_I\) = 0.55; outer pair \(\Delta\sin\theta\) = 0.347, \(r\) = 0.48, \(\rho_I\) = 0.23. \(\sum\rho_I\) = 3 + 2(0.55 + 0.55 + 0.23) = 5.65, so \(N_{\mathrm{eff}}\) = 9/5.65 = 1.6 and \(C_I\) = 0.79, not 0.58.

Follow-up: what angle span gives independent looks? \(\Delta\sin\theta \ge 1/F_\#\) = 0.67 between neighbours, about 40°, which is why speckle reduction needs large angles or sub-aperture translation of half the aperture or more.

I7A reviewer notes that your CNR rose 40% between two figures and the only change was the display dynamic range. Respond. · onsite

The reviewer is right that a dB-image CNR depends on the clipping floor: for an anechoic cyst it is ≈ \((\mathrm{DR} - 2.51)/5.57\). State the CNR definition (we use the linear envelope, pooled variance), recompute on the envelope, which has a ceiling of 1.91 for an anechoic cyst in fully developed speckle, and add gCNR, which is invariant to strictly monotone display maps.

Follow-up: can gCNR also be gamed? Not by the display map; it can be changed by ROI choice, bin rule and sample size (positive bias around 0.1 with about 100 cells), and by beamformers that change the output distribution's shape.

I8Define gCNR, give its operational meaning and its value for a −6 dB lesion in fully developed speckle. · research

\(\mathrm{gCNR} = 1 - \int\min(p_{\mathrm{in}}, p_{\mathrm{out}})\), the total-variation distance between the pixel-value distributions inside and outside. An ideal pixel classifier with equal priors is correct with probability \((1 + \mathrm{gCNR})/2\). For two exponential intensities with ratio \(r = 10^{0.6}\): \(\mathrm{gCNR} = r^{-1/(r-1)}(1 - 1/r)\) = 0.471 (classifier 74% correct).

Follow-up: how do you estimate it? Histograms with a fixed rule (quantile bins of the pooled data make it exactly display invariant), enough independent cells, and a null measurement on two same-statistics regions.

I9You want a Nakagami-m image of liver with ±0.05 precision. What window size, and what limits it? · research

\(\mathrm{std}\,\hat m \approx 2/\sqrt{n}\), so \(n\) ≈ 1600 independent cells. At 7.5 cells/mm² (5 MHz, 60%, F#2) that is about 210 mm², roughly a 15 mm window: far coarser than the B-mode. Limits: depth normalization (TGC, attenuation, diffraction) inside the window, mixtures at boundaries, small-\(n\) bias (\(\hat m\) is biased upward), and the lateral cell growing off focus.

Follow-up: how to do better? Compound independent looks or frequency bands, use 3D data, or use estimators with better efficiency (maximum likelihood) and report precision maps.

I10How should a ULM pipeline set the bubble detection threshold after SVD filtering? · research

Model the background as complex Gaussian: \(P_{fa} = e^{-T/\langle I\rangle}\), so \(T = 13.8\langle I\rangle\) (11.4 dB) for \(10^{-6}\) (26.16). Because \(\langle I\rangle\) varies with depth and region, normalize by a local noise estimate (CFAR) rather than a global threshold. Check the background histogram: residual clutter is heavier-tailed (K with \(\alpha = 1\) gives 2×10⁻³ at that threshold). Detection probability for a bubble at SNR \(S\) is Marcum Q: 0.87 at 13 dB.

Follow-up: why does a \(10^{-6}\) rate still matter? 10⁹ pixel tests per acquisition give about 1000 false detections, which tracking (track length, velocity limits) must remove.

I11How many scatterers per resolution cell do you need for fully developed speckle, and why is "10" an incomplete answer? · onsite

\(C_I^2 = 1 + \kappa/N_c\) with \(N_c\) counted in the correlation cell \((\int\lvert h\rvert^2)^2/\int\lvert h\rvert^4\), which for a Gaussian PSF is 1.44 times the −6 dB cell (26.6). Ten per −6 dB cell gives \(m\) = 0.935 (equal amplitudes) or 0.878 (Rayleigh amplitudes), \(\mathrm{SNR}_0\) ≈ 1.88 or 1.84. Incomplete because the answer depends on the cell definition, the amplitude distribution, the dimensionality (elevation) and on how precisely you can measure \(m\) (std \(2/\sqrt{n}\)).

Follow-up: can density alone be estimated from \(m\)? Only if there is no coherent component; otherwise use homodyned K.

26.12MATLAB 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.

Download the scripts: m26_rayleigh_rician_k.m · m26_nakagami_density.m · m26_speckle_size.m · m26_compounding_speckle.m · m26_cyst_cnr_gcnr.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m26_rayleigh_rician_k.mRayleigh, Rician and K envelope statistics vs closed forms (CI, m, SNR0, log mean and std) with pdf overlays · change r_dB, alphas · runtime ≈ 2 s 
%% m26_rayleigh_rician_k.m  --  envelope statistics: Rayleigh, Rician, K, and the log image
% Purpose : draw speckle IQ samples from three models and check their first-order
%           statistics against closed forms:
%             Rayleigh : circular complex Gaussian (fully developed speckle)
%             Rician   : coherent phasor A plus diffuse speckle of power s, r = A^2/s
%             K        : exponential speckle whose local mean fluctuates as Gamma(alpha)
%           For each: intensity contrast C_I = std(I)/mean(I), Nakagami m = 1/C_I^2,
%           envelope SNR0 = mean(R)/std(R), and the mean and std of the log image
%           10*log10(I/<I>) (Rayleigh: -2.51 dB and 5.57 dB).
% Change  : r_dB (coherent-to-diffuse ratios), alphas (K shape), Ns (samples).
% Expect  : Rayleigh C_I = 1, m = 1, SNR0 = 1.913, log mean -2.51 dB, log std 5.57 dB.
%           Rician  m = (1+r)^2/(1+2r): 1.33 (0 dB), 2.77 (6 dB), 5.76 (10 dB).
%           K       m = alpha/(alpha+2): 0.33 (alpha 1), 0.67 (alpha 4), 0.83 (alpha 10).
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
rng(1);
Ns = 400000;                                  % samples (independent resolution cells)
cg = @(n) (randn(1,n) + 1i*randn(1,n)) / sqrt(2);   % circular Gaussian, E|z|^2 = 1
stat = @(I) [std(I)/mean(I), mean(I)^2/var(I), mean(sqrt(I))/std(sqrt(I)), ...
             mean(10*log10(I/mean(I))), std(10*log10(I))];

%% 1. Rayleigh (fully developed speckle)
I_ray = abs(cg(Ns)).^2;
s = stat(I_ray);
fprintf('model            C_I     m      SNR0   logmean(dB) logstd(dB)   theory\n');
fprintf('Rayleigh       %6.3f %6.3f %6.3f   %7.2f   %7.2f      C_I=1 m=1 SNR0=%.3f mean=%.2f std=%.2f\n', ...
    s, sqrt(pi/(4-pi)), -10/log(10)*0.5772157, 10/log(10)*pi/sqrt(6));

%% 2. Rician: coherent component A on top of diffuse power s = 1
r_dB = [0 6 10];
I_ric = zeros(numel(r_dB), Ns);
for k = 1:numel(r_dB)
    r = 10^(r_dB(k)/10);
    z = sqrt(r) + cg(Ns);                     % A = sqrt(r), diffuse power 1
    I_ric(k,:) = abs(z).^2;
    s = stat(I_ric(k,:));
    fprintf('Rician r=%2d dB %6.3f %6.3f %6.3f   %7.2f   %7.2f      C_I=%.3f m=%.3f\n', ...
        r_dB(k), s, sqrt(1+2*r)/(1+r), (1+r)^2/(1+2*r));
end

%% 3. K distribution as a compound: I = x * E, x ~ Gamma(alpha, 1/alpha), E ~ Exp(1)
% Gamma variates without toolboxes: Marsaglia-Tsang for shape >= 1, boosted for shape < 1.
alphas = [1 4 10];
I_K = zeros(numel(alphas), Ns);
for k = 1:numel(alphas)
    a = alphas(k);  ab = a + (a < 1);         % boost shape below 1
    d = ab - 1/3;  c = 1/sqrt(9*d);
    x = zeros(1, Ns);  todo = 1:Ns;
    while ~isempty(todo)
        n = numel(todo);  g = randn(1, n);  v = (1 + c*g).^3;  u = rand(1, n);
        ok = v > 0 & log(u) < 0.5*g.^2 + d - d*v + d*log(max(v, eps));
        x(todo(ok)) = d * v(ok);  todo = todo(~ok);
    end
    if a < 1, x = x .* rand(1, Ns).^(1/a); end
    x = x / a;                                % mean 1, variance 1/alpha
    I_K(k,:) = x .* (-log(rand(1, Ns)));      % exponential speckle with fluctuating mean
    s = stat(I_K(k,:));
    fprintf('K alpha=%-4g   %6.3f %6.3f %6.3f   %7.2f   %7.2f      C_I=%.3f m=%.3f\n', ...
        a, s, sqrt(1+2/a), a/(a+2));
end

%% 4. Figures: envelope histograms with closed-form pdfs; log-image histograms
Rgrid = linspace(0.005, 4, 400);
figure('Position', [80 80 1150 380]);
subplot(1,3,1);                               % Rayleigh and Rician envelopes
cols = [0.2 0.4 0.8; 0.8 0.3 0.5; 0.7 0.5 0.1; 0.1 0.55 0.4];
edges = linspace(0, 6, 61);  ctr = (edges(1:end-1) + edges(2:end))/2;  bw = edges(2) - edges(1);
R = sqrt(I_ray);  cnt = histc(R, edges);  plot(ctr, cnt(1:end-1)/(Ns*bw), 'o', 'Color', cols(1,:)); hold on;
plot(Rgrid, 2*Rgrid.*exp(-Rgrid.^2), '-', 'Color', cols(1,:), 'LineWidth', 1.5);
Rg2 = linspace(0.005, 6, 500);
for k = 1:numel(r_dB)
    r = 10^(r_dB(k)/10);  A = sqrt(r);
    R = sqrt(I_ric(k,:));  cnt = histc(R, edges);
    plot(ctr, cnt(1:end-1)/(Ns*bw), 'o', 'Color', cols(k+1,:));
    % Rice pdf with diffuse power 1: 2R exp(-(R^2+A^2)) I0(2AR), via scaled besseli
    plot(Rg2, 2*Rg2.*exp(-(Rg2 - A).^2).*besseli(0, 2*A*Rg2, 1), '-', 'Color', cols(k+1,:), 'LineWidth', 1.5);
end
xlabel('envelope R (diffuse rms = 1)'); ylabel('pdf');
title('Rayleigh and Rician (r = 0, 6, 10 dB)'); xlim([0 6]);

subplot(1,3,2);                               % K envelopes (mean intensity 1) with Nakagami fit
edges = linspace(0, 4, 81);  ctr = (edges(1:end-1) + edges(2:end))/2;  bw = edges(2) - edges(1);
for k = 1:numel(alphas)
    a = alphas(k);  b = 2*sqrt(a);            % <I> = 1
    R = sqrt(I_K(k,:));  cnt = histc(R, edges);
    semilogy(ctr, cnt(1:end-1)/(Ns*bw), 'o', 'Color', cols(k+1,:)); hold on;
    pK = 2*b/gamma(a) * (b*Rgrid/2).^a .* besselk(a-1, b*Rgrid);
    semilogy(Rgrid, pK, '-', 'Color', cols(k+1,:), 'LineWidth', 1.5);
    m = a/(a+2);                              % Nakagami with the same m and Omega = 1
    semilogy(Rgrid, 2*m^m/gamma(m) * Rgrid.^(2*m-1) .* exp(-m*Rgrid.^2), '--', 'Color', cols(k+1,:));
end
semilogy(Rgrid, 2*Rgrid.*exp(-Rgrid.^2), 'k-', 'LineWidth', 1);
ylim([1e-4 3]); xlabel('envelope R (<R^2> = 1)'); ylabel('pdf (log scale)');
title('K: alpha = 1, 4, 10 (solid), Nakagami (dashed)');

subplot(1,3,3);                               % log-compressed speckle
edges = -45:1:15;  ctr = (edges(1:end-1) + edges(2:end))/2;
D = 10*log10(I_ray);  cnt = histc(D, edges);  plot(ctr, cnt(1:end-1)/Ns, 'o', 'Color', cols(1,:)); hold on;
Dg = linspace(-45, 15, 400);  Ig = 10.^(Dg/10);
plot(Dg, log(10)/10 * Ig .* exp(-Ig), '-', 'Color', cols(1,:), 'LineWidth', 1.5);
D = 10*log10(I_ric(3,:)/mean(I_ric(3,:)));  cnt = histc(D, edges);  plot(ctr, cnt(1:end-1)/Ns, 's', 'Color', cols(4,:));
D = 10*log10(I_K(1,:));  cnt = histc(D, edges);  plot(ctr, cnt(1:end-1)/Ns, '^', 'Color', cols(2,:));
xlabel('10 log_{10}(I/<I>)  (dB)'); ylabel('fraction per dB');
title('log image: Rayleigh (line), Rician 10 dB, K alpha=1'); xlim([-45 15]);

What to look for. Rayleigh: \(C_I\) 1.001, \(m\) 0.998, \(\mathrm{SNR}_0\) 1.913, log mean −2.50 dB, log std 5.56 dB (theory 1, 1, 1.913, −2.51, 5.57). Rician \(m\) = 1.335, 2.762, 5.744 at \(r\) = 0, 6, 10 dB (theory 1.333, 2.768, 5.762) with log std 5.09, 3.29, 2.00 dB. K: \(m\) = 0.332, 0.667, 0.825 for \(\alpha\) = 1, 4, 10 (theory \(\alpha/(\alpha+2)\)). In the middle panel the Nakagami pdf with the same \(m\) misses the K tail beyond \(R \approx 2.5\): matching two moments does not match the tail.

m26_nakagami_density.mNakagami m and SNR0 vs scatterers per cell (Poisson vs fixed count) and the 2/√n spread of m̂ · change Nlist, nlist · runtime ≈ 5 s 
%% m26_nakagami_density.m  --  how many scatterers make speckle "fully developed"?
% Purpose : sum unit phasors with uniform random phases in each resolution cell and
%           measure the envelope statistics as the number of scatterers per cell grows.
%             Poisson count, mean N  : C_I^2 = 1 + 1/N  ->  Nakagami m = N/(N+1)  (K-like, m < 1)
%             exactly N per cell     : C_I^2 = 1 - 1/N  ->  m = N/(N-1)          (sub-Rayleigh, m > 1)
%           The Poisson result also holds for PSF-weighted scatterers if N = n*Vc, where
%           n is the number density and Vc = (int |h|^2)^2 / int |h|^4 the correlation area.
%           Part 2 checks the spread of the moment estimator m^ = (E[R^2])^2 / Var(R^2)
%           from n independent cells (std about 2/sqrt(n) at m = 1).
% Change  : Nlist (mean scatterers per cell), Ncell, nlist (cells per estimate).
% Expect  : Poisson N = 10 -> m = 0.909, SNR0 = 1.86 (97 % of 1.913);
%           fixed N = 10 -> m = 1.111; std(m^) * sqrt(n) close to 2.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
rng(1);
Ncell = 40000;                                 % independent resolution cells per point
Nlist = [1 2 3 5 10 20 50 100];                % mean number of scatterers per cell
poissrnd0 = @(lam, n) sum(cumsum(-log(rand(ceil(lam + 7*sqrt(lam) + 12), n)) / lam, 1) <= 1, 1); % Poisson via exponential gaps
nak  = @(R) mean(R.^2)^2 / var(R.^2);          % Nakagami m, moment estimator
snr0 = @(R) mean(R) / std(R);                  % envelope SNR (Rayleigh: 1.913)

%% 1. Statistics vs mean scatterer count
mP = zeros(size(Nlist));  sP = mP;  mF = mP;  sF = mP;
fprintf('  N    Poisson: m  (theory)  SNR0  | fixed N: m  (theory)  SNR0\n');
for k = 1:numel(Nlist)
    N = Nlist(k);
    cnt = poissrnd0(N, Ncell);                % Poisson number of scatterers in each cell
    maxc = max(cnt);
    ph = 2*pi*rand(maxc, Ncell);
    use = bsxfun(@le, (1:maxc)', cnt);        % keep the first cnt(j) phasors of column j
    RP = abs(sum(exp(1i*ph) .* use, 1));
    RF = abs(sum(exp(1i*2*pi*rand(N, Ncell)), 1));   % exactly N phasors per cell
    mP(k) = nak(RP);  sP(k) = snr0(RP);  mF(k) = nak(RF);  sF(k) = snr0(RF);
    if N == 1                                 % one unit phasor: R = 1 exactly, Var(R^2) = 0
        fprintf('%4d   %8.3f  (%5.3f)  %5.3f  |   R = 1 exactly (no speckle)\n', N, mP(k), N/(N+1), sP(k));
    else
        fprintf('%4d   %8.3f  (%5.3f)  %5.3f  | %8.3f  (%5.3f)  %5.3f\n', N, mP(k), N/(N+1), sP(k), mF(k), N/(N-1), sF(k));
    end
end
fprintf('Rayleigh limit: m = 1, SNR0 = %.3f\n', sqrt(pi/(4-pi)));

%% 2. Spread of the Nakagami estimate from n independent Rayleigh cells
nlist = [25 100 400 1600];  Ntrial = 4000;
fprintf('\n   n   mean(m^)   std(m^)   std*sqrt(n)\n');
for n = nlist
    R = abs(randn(n, Ntrial) + 1i*randn(n, Ntrial));     % Rayleigh envelopes
    I = R.^2;
    mh = mean(I, 1).^2 ./ var(I, 1, 1);                   % per-trial moment estimate
    fprintf('%5d   %7.3f   %7.3f   %7.2f\n', n, mean(mh), std(mh), std(mh)*sqrt(n));
end

%% 3. Figure
figure('Position', [100 100 1000 380]);
subplot(1,2,1);
Ng = logspace(0, 2, 100);
semilogx(Ng, Ng./(Ng+1), 'b-', 'LineWidth', 1.5); hold on;
semilogx(Nlist, mP, 'bo', 'MarkerFaceColor', 'b');
semilogx(Ng(Ng > 1.5), Ng(Ng > 1.5)./(Ng(Ng > 1.5)-1), 'r-', 'LineWidth', 1.5);
semilogx(Nlist(2:end), mF(2:end), 'rs', 'MarkerFaceColor', 'r');
plot([1 100], [1 1], 'k--');
ylim([0.4 2.2]); xlabel('mean scatterers per resolution cell N'); ylabel('Nakagami m');
title('Poisson count (blue) vs fixed count (red)'); grid on;
subplot(1,2,2);
semilogx(Nlist, sP, 'bo-', 'MarkerFaceColor', 'b'); hold on;
semilogx(Nlist(2:end), sF(2:end), 'rs-', 'MarkerFaceColor', 'r');
plot([1 100], sqrt(pi/(4-pi))*[1 1], 'k--');
xlabel('mean scatterers per resolution cell N'); ylabel('envelope SNR_0 = mean/std');
title('SNR_0 approaches 1.913 from both sides'); grid on; ylim([1.2 2.3]);

What to look for. Poisson counts approach Rayleigh from below: \(m\) = 0.504, 0.843, 0.899, 0.962 at \(N\) = 1, 5, 10, 20 (theory \(N/(N+1)\)), \(\mathrm{SNR}_0\) 1.853 at \(N\) = 10. A fixed count approaches from above (\(m = N/(N-1)\): 1.108 at \(N\) = 10), a reminder that "few scatterers" is not automatically K-like. \(\mathrm{std}(\hat m)\cdot\sqrt{n}\) = 2.13, 2.01, 1.98, 1.91 for \(n\) = 25 … 1600, and mean \(\hat m\) = 1.22 at \(n\) = 25 (small-sample bias).

m26_speckle_size.mSpeckle intensity-autocovariance widths vs f0, bandwidth, F#, apodization against |μ|2 theory · change cases · runtime ≈ 2 s 
%% m26_speckle_size.m  --  speckle size from the PSF: autocovariance widths vs f0, bandwidth, F#
% Purpose : build fully developed speckle IQ images in k-space (white circular Gaussian
%           scatterer field times the pulse-echo transfer function H), then measure the
%           FWHM of the intensity autocovariance axially and laterally and compare with
%             axial   (Gaussian pulse, B = -6 dB two-way bandwidth): 0.8825/B in time,
%                     i.e. 0.441 c/B in depth, equal to the pulse-envelope FWHM
%             lateral (focal plane, CW, rect tx and rx apertures, two-way PSF sinc^2):
%                     0.835 lambda F#   (two-way -6 dB beam width is 0.886 lambda F#)
%           Theory = FWHM of |mu|^2, where mu = inverse FT of |H|^2 (Siegert relation).
%           Also prints the envelope and log-image ACF widths, which are narrower.
% Change  : cases (f0, fractional bandwidth, F#, apodization), grid sizes.
% Expect  : 5 MHz, 60 %, F#2, rect: axial 0.227 mm, lateral 0.514 mm (theory); measured
%           within a few percent; log-image axial width about 0.77 x the intensity width.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
rng(1);
c  = 1540;                                     % sound speed (m/s)
Nz = 256;  Nx = 256;  dz = 20e-6;  dx = 40e-6; % grid: 5.1 mm x 10.2 mm
uz = [0:Nz/2-1, -Nz/2:-1]' / (Nz*dz);          % axial spatial frequency (cycles/m), FFT order
ux = [0:Nx/2-1, -Nx/2:-1]  / (Nx*dx);          % lateral spatial frequency (cycles/m)
nreal = 4;                                     % independent images averaged per case
cases = [5 0.6 2 0;  3 0.6 2 0;  7.5 0.6 2 0;  5 0.3 2 0;  5 0.9 2 0;  5 0.6 1 0;  5 0.6 3 0;  5 0.6 2 1];
% columns: f0 (MHz), fractional -6 dB bandwidth, F#, apodization (0 rect, 1 Hann)

% FWHM of a sampled, peak-normalized symmetric profile p(1) = 1 (index 1 = zero lag)
fwhm = @(p) 2*((find(p < 0.5, 1) - 2) + (p(find(p < 0.5, 1) - 1) - 0.5) / (p(find(p < 0.5, 1) - 1) - p(find(p < 0.5, 1))));
acf2 = @(I) real(ifft2(abs(fft2(I - mean(I(:)))).^2));     % circular autocovariance

% Lateral transfer function: autocorrelation of the apodization (coarray), sampled on ux.
% The aperture coordinate is scaled so that the coarray half-width is 1/(lambda F#).
coarray = @(u, U, apod) interp1(linspace(-1, 1, 2001), ...
    conv(apod(linspace(-0.5, 0.5, 1001)), apod(linspace(-0.5, 0.5, 1001))) / ...
    sum(apod(linspace(-0.5, 0.5, 1001)).^2), u / U, 'linear', 0);
rectw = @(t) ones(size(t));
hannw = @(t) 0.5 + 0.5*cos(2*pi*t);

fprintf(' f0   BW   F#  apod |  axial FWHM (mm) th / meas   | lateral FWHM (mm) th / meas  (th/lamF#) | env/I  log/I (axial) | speckles/mm^2\n');
for k = 1:size(cases, 1)
    f0 = cases(k,1)*1e6;  B = cases(k,2)*f0;  Fn = cases(k,3);  lam = c/f0;
    if cases(k,4) == 0, apod = rectw; an = 'rect'; else, apod = hannw; an = 'hann'; end
    sf = B / (2*sqrt(2*log(2)));               % |H(f)| = exp(-f^2/(2 sf^2)), -6 dB full width = B
    Hz = exp(-(c*uz/2).^2 / (2*sf^2));         % axial: temporal f = c*uz/2 (round trip)
    Hx = coarray(ux, 1/(lam*Fn), apod);        % lateral: coarray, half-width 1/(lambda F#)
    H  = Hz * Hx;                              % separable pulse-echo transfer function
    % theory: |mu|^2 with mu = IFFT(|H|^2), on a 16x finer 1D grid
    Lf = 4096;  uzf = [0:Lf/2-1, -Lf/2:-1]'/(Lf*dz/4);  muz = abs(ifft(exp(-(c*uzf/2).^2/sf^2)));
    muz = (muz/muz(1)).^2;  thz = fwhm(muz) * dz/4;
    uxf = [0:Lf/2-1, -Lf/2:-1]/(Lf*dx/4);  mux = abs(ifft(coarray(uxf, 1/(lam*Fn), apod).^2));
    mux = (mux/mux(1)).^2;  thx = fwhm(mux) * dx/4;
    Vc = 1 / (sum(muz)*dz/4 * sum(mux)*dx/4);  % 1 / speckle correlation area (int int |mu|^2), per m^2
    % simulation
    AI = 0;  AA = 0;  AL = 0;
    for r = 1:nreal
        W = (randn(Nz, Nx) + 1i*randn(Nz, Nx)) / sqrt(2);
        z = ifft2(fft2(W) .* H);  I = abs(z).^2;
        AI = AI + acf2(I);  AA = AA + acf2(sqrt(I));  AL = AL + acf2(log(I));
    end
    pz = AI(:,1)/AI(1,1);  px = AI(1,:)'/AI(1,1);
    mz = fwhm(pz)*dz;  mx = fwhm(px)*dx;
    ea = fwhm(AA(:,1)/AA(1,1))*dz / mz;  la = fwhm(AL(:,1)/AL(1,1))*dz / mz;
    fprintf('%4.1f %4.1f %3.0f  %s |      %6.3f / %6.3f         |      %6.3f / %6.3f   (%5.3f)    |  %4.2f   %4.2f        | %6.2f\n', ...
        f0/1e6, B/f0, Fn, an, thz*1e3, mz*1e3, thx*1e3, mx*1e3, thx/(lam*Fn), ea, la, Vc*1e-6);
    if k == 1, Ishow = I; pz1 = pz; px1 = px; mzT = muz; mxT = mux; end
end
fprintf('axial theory 0.8825/B * c/2 = 0.441 c/B; lateral rect theory 0.835 lambda F#\n');

%% Figure: one speckle image and its autocovariance profiles vs theory
figure('Position', [80 80 1100 380]);
subplot(1,3,1);
imagesc((0:Nx-1)*dx*1e3, (0:Nz-1)*dz*1e3, 10*log10(Ishow/mean(Ishow(:))), [-30 10]);
colormap(gray); axis image; colorbar; xlabel('lateral (mm)'); ylabel('depth (mm)');
title('speckle, 5 MHz, 60 %, F#2 (dB)');
subplot(1,3,2);
lagz = (0:Nz/2)*dz*1e3;  plot(lagz, pz1(1:Nz/2+1), 'bo'); hold on;
lf = (0:Lf/2)*dz/4*1e3;  plot(lf, mzT(1:Lf/2+1), 'b-', 'LineWidth', 1.5);
plot([0 1], [0.5 0.5], 'k--'); xlim([0 0.6]); grid on;
xlabel('axial lag (mm)'); ylabel('intensity autocovariance'); title('axial: sim (o) vs |mu|^2');
subplot(1,3,3);
lagx = (0:Nx/2)*dx*1e3;  plot(lagx, px1(1:Nx/2+1), 'ro'); hold on;
lf = (0:Lf/2)*dx/4*1e3;  plot(lf, mxT(1:Lf/2+1), 'r-', 'LineWidth', 1.5);
plot([0 2], [0.5 0.5], 'k--'); xlim([0 1.5]); grid on;
xlabel('lateral lag (mm)'); ylabel('intensity autocovariance'); title('lateral: sim (o) vs |mu|^2');

What to look for. 5 MHz, 60%, F#2, rect: axial 0.227 mm theory and measured, lateral 0.514 mm theory vs 0.496 measured (one-realization scatter), \(0.834\,\lambda F_\#\); 7.49 cells/mm². Halving the bandwidth (30%) doubles the axial width (0.453 mm); F# 1 and 3 give 0.257 and 0.771 mm laterally; Hann gives \(1.401\,\lambda F_\#\). The envelope ACF is 0.95–0.98 and the log-image ACF 0.76–0.80 of the intensity width.

m26_compounding_speckle.mIncoherent vs coherent compounding of translated sub-aperture looks: CI vs √(ΣρI)/N, envelope SNR0, coherent speckle size · change Nlist, dwlist · runtime ≈ 2 s 
%% m26_compounding_speckle.m  --  incoherent vs coherent compounding of speckle looks
% Purpose : N looks of the same scatterers from sub-apertures of width w whose centres are
%           translated by d (all focused at the same point). In k-space each look keeps the
%           same axial band and a lateral coarray triangle shifted by 2d (in aperture units).
%           Field correlation between looks = overlap of the two triangles:
%             rho_F(s) = 1 - 1.5 s^2 + 0.75 s^3 (s <= 1),  (2 - s)^3 / 4 (1 < s < 2),  s = 2 d / w
%           Intensity correlation rho_I = rho_F^2. Summing intensities of N looks gives
%             C_I = sqrt(sum_ij rho_I(i,j)) / N,   N_eff = N^2 / sum_ij rho_I(i,j)
%           Coherent summation (IQ) widens the k-space support: finer speckle, but it is
%           still circular Gaussian, so envelope SNR0 stays 1.91 for any N.
% Change  : Nlist, dwlist (d/w), grid.
% Expect  : d/w >= 0.5 makes looks almost independent (rho_I <= 0.06): C_I ~ 1/sqrt(N);
%           coherent SNR0 = 1.91 +- 0.03 for every case; coherent lateral speckle shrinks.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
rng(1);
Nz = 128;  Nx = 512;                           % grid (pixels); lateral unit = 1 pixel
ux = [0:Nx/2-1, -Nx/2:-1] / Nx;                % lateral frequency (cycles/pixel)
uz = [0:Nz/2-1, -Nz/2:-1]' / Nz;
U  = 1/32;                                     % coarray half-width of one look (cycles/pixel)
Hz = exp(-uz.^2 / (2*0.08^2));                 % axial band (Gaussian)
tri = @(u, c0) max(0, 1 - abs(u - c0)/U);
rhoF = @(s) (abs(s) <= 1).*(1 - 1.5*s.^2 + 0.75*abs(s).^3) + (abs(s) > 1 & abs(s) < 2).*(2 - abs(s)).^3/4;
acfx = @(I) mean(real(ifft(abs(fft(bsxfun(@minus, I, mean(I, 2)), [], 2)).^2, [], 2)), 1);
fw = @(p) 2*((find(p/p(1) < 0.5, 1) - 2) + (p(find(p/p(1) < 0.5, 1) - 1)/p(1) - 0.5) / ...
          ((p(find(p/p(1) < 0.5, 1) - 1) - p(find(p/p(1) < 0.5, 1)))/p(1)));
W = (randn(Nz, Nx) + 1i*randn(Nz, Nx)) / sqrt(2);
S = fft2(W);
single_look = ifft2(S .* (Hz * tri(ux, 0)));
w1 = fw(acfx(abs(single_look).^2));

Nlist = [1 3 5 9];  dwlist = [0.1 0.25 0.5 1.0];
fprintf('single look: SNR0 %.3f, lateral speckle FWHM %.1f px\n', mean(abs(single_look(:)))/std(abs(single_look(:))), w1);
fprintf('  N  d/w  | incoh. intensity C_I  meas / theory  N_eff | incoh. env SNR0 | coherent SNR0  lat. FWHM (px)\n');
res = zeros(numel(Nlist), numel(dwlist), 3);
for a = 1:numel(Nlist)
    N = Nlist(a);
    for b = 1:numel(dwlist)
        s = 2*dwlist(b);  ctr = ((0:N-1) - (N-1)/2) * s * U;
        Isum = 0;  Asum = 0;  Zsum = 0;
        for k = 1:N
            zk = ifft2(S .* (Hz * tri(ux, ctr(k))));
            Isum = Isum + abs(zk).^2;  Asum = Asum + abs(zk);  Zsum = Zsum + zk;
        end
        [ii, jj] = meshgrid(1:N);
        rhoI = rhoF((ii - jj)*s).^2;
        Cth = sqrt(sum(rhoI(:)))/N;  Neff = N^2/sum(rhoI(:));
        Cm = std(Isum(:))/mean(Isum(:));
        snrA = mean(Asum(:))/std(Asum(:));
        snrC = mean(abs(Zsum(:)))/std(abs(Zsum(:)));
        wc = fw(acfx(abs(Zsum).^2));
        res(a, b, :) = [Cm Cth snrC];
        fprintf('%3d  %4.2f |      %5.3f / %5.3f    %5.2f |     %5.3f       |    %5.3f       %5.1f\n', ...
            N, dwlist(b), Cm, Cth, Neff, snrA, snrC, wc);
        if N == 5 && dwlist(b) == 0.5, Iinc = Asum; Icoh = abs(Zsum); end
    end
end
fprintf('limit for independent looks: C_I = 1/sqrt(N); envelope SNR0 -> 1.91*sqrt(N)\n');

%% Figure
figure('Position', [60 60 1150 420]);
show = @(A) 20*log10(A/mean(A(:)));
xr = 1:256;
subplot(2,3,1); imagesc(show(abs(single_look(:, xr))), [-25 8]); colormap(gray); axis image; title('single look'); xlabel('lateral (px)'); ylabel('depth (px)');
subplot(2,3,2); imagesc(show(Iinc(:, xr)), [-25 8]); axis image; title('incoherent, N = 5, d/w = 0.5'); xlabel('lateral (px)');
subplot(2,3,3); imagesc(show(Icoh(:, xr)), [-25 8]); axis image; title('coherent, N = 5, d/w = 0.5'); xlabel('lateral (px)');
subplot(2,3,[4 5 6]);
cols = 'brkm';
for a = 1:numel(Nlist)
    plot(dwlist, squeeze(res(a,:,1)), [cols(a) 'o'], 'MarkerFaceColor', cols(a)); hold on;
    dg = linspace(0, 1, 101);  th = zeros(size(dg));
    for q = 1:numel(dg)
        [ii, jj] = meshgrid(1:Nlist(a));  r = rhoF((ii - jj)*2*dg(q)).^2;  th(q) = sqrt(sum(r(:)))/Nlist(a);
    end
    plot(dg, th, [cols(a) '-']);
end
xlabel('look translation d / sub-aperture width w'); ylabel('intensity speckle contrast C_I');
title('incoherent compounding, N = 1, 3, 5, 9: simulation (o) vs sqrt(sum rho_I)/N'); grid on; ylim([0 1.1]);

What to look for. Incoherent intensity contrast follows the overlap theory: \(N\) = 5 at \(d/w\) = 0.1, 0.25, 0.5, 1.0 gives 0.818, 0.623, 0.463, 0.440 vs 0.848, 0.617, 0.469, 0.447 (\(N_{\mathrm{eff}}\) 1.39 to 5.0). The single look reads 0.960 instead of 1 because one 128 × 512 image holds only a few hundred independent cells. Coherent \(\mathrm{SNR}_0\) stays between 1.90 and 1.94 for every case while the lateral speckle width falls from 25 to 1.4 pixels: finer, not smoother.

m26_cyst_cnr_gcnr.mCyst phantom: contrast, CNR on envelope / intensity / dB image (two denominators) and gCNR vs cyst contrast and display dynamic range; gCNR null test · change cset_dB, DRs · runtime ≈ 2 s 
%% m26_cyst_cnr_gcnr.m  --  cyst phantom: contrast, CNR (four definitions) and gCNR
% Purpose : simulate a speckle phantom with a circular cyst (5 MHz, 60 % bandwidth, F#2,
%           Hann apodization), then measure image quality as reviewers see it reported.
%   Contrast   C   = 10 log10(<I_in>/<I_out>)  (dB, mean intensities)
%   CNR (ours)     = |mu_in - mu_out| / sqrt(sigma_in^2 + sigma_out^2), on the ENVELOPE
%   variants       : same formula on intensity, on the dB image, and on the dB image
%                    with sigma_out only in the denominator
%   gCNR           = 1 - OVL,  OVL = sum over bins of min(p_in, p_out)
%                    (Rodriguez-Molares et al., IEEE TUFFC 2020). Bins = quantiles of the
%                    pooled data, so the estimate is exactly invariant to strictly monotone
%                    transforms (log, power law, gain).
%   Rayleigh theory for intensity ratio r = <I_out>/<I_in>:  gCNR = r^(-1/(r-1)) (1 - 1/r)
%   CNR_env theory = (sqrt(pi)/2)(1 - q) / sqrt((1 - pi/4)(1 + q^2)), q = 1/sqrt(r)
% Change  : cyst contrasts cset_dB, radius Rc, dynamic ranges DRs, apodization.
% Expect  : CNR_env saturates near 1.91 for an anechoic cyst, CNR on dB images grows with
%           contrast and with clipping; gCNR is identical for log, power and linear maps
%           and can only drop when the display clips (loss of information).
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
rng(1);
c = 1540;  f0 = 5e6;  lam = c/f0;  B = 0.6*f0;  Fn = 2;
Nz = 256;  Nx = 256;  dz = 40e-6;  dx = 50e-6;      % 10.2 mm x 12.8 mm
uz = [0:Nz/2-1, -Nz/2:-1]'/(Nz*dz);  ux = [0:Nx/2-1, -Nx/2:-1]/(Nx*dx);
sf = B/(2*sqrt(2*log(2)));  Hz = exp(-(c*uz/2).^2/(2*sf^2));
t = linspace(-0.5, 0.5, 1001);                       % aperture coordinate (units of width)
coarray = @(apod) interp1(linspace(-1, 1, 2001), conv(apod, apod)/max(conv(apod, apod)), ux*lam*Fn, 'linear', 0);
H_hann = Hz * coarray(0.5 + 0.5*cos(2*pi*t));
H_rect = Hz * coarray(ones(size(t)));
[X, Z] = meshgrid(((0:Nx-1) - Nx/2)*dx, ((0:Nz-1) - Nz/2)*dz);
Rc = 2.5e-3;  rr = sqrt(X.^2 + Z.^2);
inR  = rr < 0.8*Rc;                                  % ROI inside, inset from the edge
outR = rr > 1.3*Rc & rr < 1.9*Rc;                    % background annulus

% ---- metric helpers (anonymous only) ----
cnr   = @(a, b) abs(mean(a) - mean(b)) / sqrt(var(a) + var(b));
qedge = @(p, nb) unique(p(max(1, round(linspace(1, numel(p), nb + 1)))));
fold  = @(h) [h(1:end-2); h(end-1) + h(end)];        % histc puts x == last edge in its own bin
ovl   = @(a, b, e) sum(min(fold(histc(a(:), e))/numel(a), fold(histc(b(:), e))/numel(b)));
gcnrn = @(a, b, nb) 1 - ovl(a, b, qedge(sort([a(:); b(:)]), nb));
gcnr  = @(a, b) gcnrn(a, b, 100);
gth   = @(r) r.^(-1./(r - 1)) .* (1 - 1./r);
cth   = @(r) sqrt(pi)/2*(1 - 1./sqrt(r)) ./ sqrt((1 - pi/4)*(1 + 1./r));
W = (randn(Nz, Nx) + 1i*randn(Nz, Nx))/sqrt(2);     % one scatterer realization for all cases
noise_dB = -50;                                      % thermal noise power re. background speckle
Nn = (randn(Nz, Nx) + 1i*randn(Nz, Nx))/sqrt(2);
Ps = mean(abs(H_hann(:)).^2);                        % speckle power per pixel for unit reflectivity
Nn = Nn * sqrt(Ps * 10^(noise_dB/10));

%% 1. Metrics vs cyst contrast (log display, no clipping)
cset_dB = [-Inf -20 -10 -6 -3];
fprintf('set C | meas C (dB) | CNR env (theory) | CNR int | CNR dB | CNR dB, sigma_out | gCNR (theory at meas C)\n');
for k = 1:numel(cset_dB)
    refl = ones(Nz, Nx);  refl(rr < Rc) = 10^(cset_dB(k)/20);
    A = abs(ifft2(fft2(W .* refl) .* H_hann) + Nn);  I = A.^2;
    L = 10*log10(I / mean(I(outR)));                 % dB image, background at 0 dB
    Cm = 10*log10(mean(I(inR))/mean(I(outR)));  r = 10^(-Cm/10);
    fprintf('%5.0f |   %7.2f   |  %5.3f (%5.3f)  |  %5.3f  | %5.3f  |      %5.3f        |  %5.3f (%5.3f)\n', ...
        cset_dB(k), Cm, cnr(A(inR), A(outR)), cth(r), cnr(I(inR), I(outR)), cnr(L(inR), L(outR)), ...
        abs(mean(L(inR)) - mean(L(outR)))/std(L(outR)), gcnr(A(inR), A(outR)), gth(r));
    if k == 2, A20 = A; end
    if k == 1, A0 = A; end
end
Ar = abs(ifft2(fft2(W .* (rr >= Rc)) .* H_rect) + Nn);
fprintf('anechoic cyst, rect apodization: C = %.1f dB (Hann above: wider mainlobe tails leak more near the edge,\n   rect sidelobes leak more far from it; noise floor %d dB)\n', ...
    10*log10(mean(Ar(inR).^2)/mean(Ar(outR).^2)), noise_dB);

%% 2. Same image, different display maps and dynamic ranges (cyst at -20 dB)
A = A20;  An = A / sqrt(mean(A(outR).^2));           % background rms = 1
DRs = [80 50 40 30 20];
fprintf('\n-20 dB cyst. Display v = T(max(An, 10^(-DR/20))): CNR and gCNR of the DISPLAYED values\n');
fprintf(' DR (dB) | log: CNR  gCNR | power 0.3: CNR  gCNR | linear: CNR  gCNR\n');
for DR = DRs
    Ac = max(An, 10^(-DR/20));
    Dl = 20*log10(Ac);  Dp = Ac.^0.3;  Dn = Ac;
    fprintf('  %4d   |   %5.3f  %5.3f |      %5.3f  %5.3f    |    %5.3f  %5.3f\n', DR, ...
        cnr(Dl(inR), Dl(outR)), gcnr(Dl(inR), Dl(outR)), cnr(Dp(inR), Dp(outR)), gcnr(Dp(inR), Dp(outR)), ...
        cnr(Dn(inR), Dn(outR)), gcnr(Dn(inR), Dn(outR)));
end

%% 3. Estimator sanity check: two halves of the same background must give gCNR ~ 0
half1 = outR & X < 0;  half2 = outR & X >= 0;
mu = abs(ifft2(abs(H_hann).^2));  mu = mu/mu(1,1);  Ac = sum(mu(:).^2)*dx*dz;   % speckle correlation area (m^2)
fprintf('\nsame-statistics regions (%d and %d pixels, but only about %.0f independent speckle cells each):\n', ...
    sum(half1(:)), sum(half2(:)), sum(half1(:))*dx*dz/Ac);
fprintf('   gCNR should be 0; the excess is sampling noise set by the cell count, not by the pixel count\n');
for nb = [10 20 50 100]
    fprintf('   %3d bins: gCNR(same) = %.3f   gCNR(-20 dB cyst) = %.3f\n', nb, gcnrn(A(half1), A(half2), nb), gcnrn(A(inR), A(outR), nb));
end

%% Figures
figure('Position', [60 60 1150 400]);
ax = ((0:Nx-1) - Nx/2)*dx*1e3;  az = ((0:Nz-1) - Nz/2)*dz*1e3;
subplot(1,3,1); imagesc(ax, az, 20*log10(A0/sqrt(mean(A0(outR).^2))), [-50 10]); colormap(gray); axis image; colorbar;
xlabel('lateral (mm)'); ylabel('depth (mm)'); title('anechoic cyst, 60 dB display');
subplot(1,3,2); imagesc(ax, az, 20*log10(max(An, 10^(-30/20))), [-30 10]); axis image; colorbar;
xlabel('lateral (mm)'); title('-20 dB cyst, 40 dB display');
subplot(1,3,3);
Dl = 20*log10(An);  edges = -60:1:15;  ctr = edges(1:end-1) + 0.5;
hi = histc(Dl(inR), edges);  ho = histc(Dl(outR), edges);
plot(ctr, hi(1:end-1)/sum(inR(:)), 'r-', ctr, ho(1:end-1)/sum(outR(:)), 'b-', 'LineWidth', 1.5); hold on;
area(ctr, min(hi(1:end-1)/sum(inR(:)), ho(1:end-1)/sum(outR(:))), 'FaceColor', [0.7 0.7 0.7]);
xlabel('displayed value (dB)'); ylabel('fraction per dB'); grid on;
title(sprintf('-20 dB cyst: inside (red), outside (blue), overlap (gray); gCNR = %.2f', gcnr(A(inR), A(outR))));

What to look for. Measured \(\mathrm{CNR}_{\mathrm{env}}\) tracks the Rayleigh theory (−20 dB cyst: 1.764 vs 1.704; anechoic 1.955, above the 1.913 ceiling by sampling noise) and gCNR matches \(r^{-1/(r-1)}(1 - 1/r)\) (0.944 vs 0.941 at −19.6 dB). For the −20 dB cyst the log-display CNR rises from 2.60 to 3.03 as DR falls from 80 to 20 dB, the power-law and linear maps give 2.61 and 1.76, while gCNR = 0.944 in every column. Two halves of the same background read gCNR 0.09–0.12: the estimator's null value with about 85 cells per region.

26.13Self-check

Q1A phantom ROI gives envelope mean 41.2 and std 21.6 (linear units). Is the speckle fully developed?

\(\mathrm{SNR}_0\) = 41.2/21.6 = 1.91: consistent with Rayleigh. Whether 1.91 vs, say, 1.85 can be distinguished depends on the number of independent cells (std ≈ \(1.36/\sqrt{n}\)).

Q2Intensity-ratio contrast of a lesion is −6 dB. What do you get if you instead subtract mean dB values, (a) both regions Rayleigh, (b) lesion Rician with r = 10 dB?

(a) −6.0 dB: the −2.51 dB offsets cancel. (b) The lesion's log mean is −0.41 dB re its own mean intensity, so the mean-dB difference is −6 + 2.51 − 0.41 = −3.9 dB: a 2.1 dB bias.

Q3Nakagami m̂ = 0.80 from 100 independent cells. What can you conclude?

\(\mathrm{std}\,\hat m \approx 2/\sqrt{100}\) = 0.2, so 0.80 is within about one standard deviation of 1: not significantly pre-Rayleigh. With 1600 cells (std 0.05) the same value would be 4σ below 1.

Q4At 7.5 MHz, 80% bandwidth, F#1.5, rect apodization, what are the axial and lateral speckle sizes?

\(B\) = 6 MHz: axial 0.441 × 1540/6×10⁶ = 0.113 mm. \(\lambda\) = 0.205 mm, \(\lambda F_\#\) = 0.308 mm, lateral 0.835 × 0.308 = 0.257 mm.

Q5Four looks with pairwise field correlation 0.6 are averaged as intensities. What is the speckle contrast?

\(\rho_I\) = 0.36. \(N_{\mathrm{eff}} = 4/(1 + 3\times 0.36)\) = 1.92, contrast \(1/\sqrt{1.92}\) = 0.72 (independent looks would give 0.50).

Q6An anechoic cyst displayed with 40 dB dynamic range: what is the dB-image CNR, and what is the envelope CNR?

dB image: (40 − 2.51)/5.57 = 6.7 (exact clipped computation: 6.74). Envelope: 1.91, independent of the display.

Q7What threshold (re mean background intensity) gives a per-pixel false-alarm rate of \(10^{-4}\) on complex Gaussian noise?

\(T = -\ln 10^{-4}\) = 9.21\(\langle I\rangle\) = 9.6 dB.

26.14Figures from the literature

Published figures with real or carefully simulated data behind the statistics and metrics of this module. Citations and figure numbers were checked against the papers; images are not reproduced, follow the links.

A display curve that inflates CNR without changing what can be seen

A. Rodriguez-Molares, O. M. H. Rindal, J. D'hooge, S.-E. Måsøy, A. Austeng, M. A. Lediju Bell and H. Torp, "The generalized contrast-to-noise ratio: a formal definition for lesion detectability," IEEE TUFFC 67(4):745–759, 2020, Fig. 1 (open access, CC BY 4.0). doi:10.1109/TUFFC.2019.2956855

Fig. 1: the same DAS image before and after an S-shaped compression curve ("S-DAS"): the estimated CNR changes although the information is identical. That is the argument of Section 26.6 and rule 14 in one picture: CNR depends on the pixel-value map, the overlap of the two pdfs does not.

Inside and outside pdfs, and C, CNR and gCNR across beamformers in silico, in vitro and in vivo

Same paper, Figs. 4, 5, 7, 8 and 9.

Fig. 4: the inside (circle) and outside (square) ROIs used for every case, the counterpart of the drawing above Demo 26.5. Fig. 5: estimated envelope pdfs inside and outside a cyst for DAS at a fixed SNR; their overlap is \(1 - \mathrm{gCNR}\) (26.14). Figs. 7–9: contrast, CNR and gCNR for a set of beamformers on simulated, phantom and in vivo data. Compare how the three metrics rank the beamformers: where a contrast or CNR gain is not matched by gCNR, the output distributions changed shape rather than separating, the situation rule 15 warns about.

Coherent compounding on a standard phantom: finer, not smoother

H. Liebgott, A. Rodriguez-Molares, F. Cervenansky, J. A. Jensen and O. Bernard, "Plane-wave imaging challenge in medical ultrasound," Proc. IEEE International Ultrasonics Symposium, 2016, Figs. 2 and 3. doi:10.1109/ULTSYM.2016.7728908

Fig. 2: simulated and experimental (CIRS 040GSE phantom) contrast and resolution images from a single 0° plane wave; Fig. 3: the same with 75 coherently compounded plane waves. Compare the cyst edges and the background texture between the two figures: by rule 10 compounding should sharpen the cysts and lower clutter while the speckle grain gets finer, not smoother. The challenge's CNR is computed from image gray levels and expressed in dB, so its value depends on the display map (rule 14).

Envelope-statistics parameters as tissue features

M. L. Oelze and J. Mamou, "Review of quantitative ultrasound: envelope statistics and backscatter coefficient imaging and contributions to diagnostic ultrasound," IEEE TUFFC 63(2):336–351, 2016, Fig. 3. doi:10.1109/TUFFC.2015.2513958

Fig. 3: three tumour types (rat fibroadenoma, mouse carcinoma, mouse sarcoma) plotted in a feature space of ESD and the homodyned-K parameters \(\mu\) (scatterers per resolution cell) and \(k\) (coherent-to-diffuse ratio). It shows the two departures of Section 26.2 (density and coherence) estimated separately on real tissue, the thing a single Nakagami \(m\) cannot do.

26.15Go deeper

  • J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications, Roberts & Company, 2007. First- and second-order speckle statistics, Siegert relation, sums of speckle patterns.
  • C. B. Burckhardt, "Speckle in ultrasound B-mode scans," IEEE Trans. Sonics Ultrason., 1978.
  • R. F. Wagner, S. W. Smith, J. M. Sandrik, H. Lopez, "Statistics of speckle in ultrasound B-scans," IEEE Trans. Sonics Ultrason. 30(3):156–163, 1983.
  • S. W. Smith, R. F. Wagner, J. M. Sandrik, H. Lopez, "Low contrast detectability and contrast/detail analysis in medical ultrasound," IEEE Trans. Sonics Ultrason., 1983.
  • B. J. Oosterveld, J. M. Thijssen, W. A. Verhoef, "Texture of B-mode echograms: 3-D simulations and experiments of the effects of diffraction and scatterer density," Ultrasonic Imaging, 1985.
  • G. E. Trahey, S. W. Smith, O. T. von Ramm, "Speckle pattern correlation with lateral aperture translation: experimental results and implications for spatial compounding," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 1986.
  • E. Jakeman, P. N. Pusey, "A model for non-Rayleigh sea echo," IEEE Trans. Antennas Propag., 1976.
  • V. Dutt, J. F. Greenleaf, "Ultrasound echo envelope analysis using a homodyned K distribution signal model," Ultrasonic Imaging, 1994.
  • P. M. Shankar, "A general statistical model for ultrasonic backscattering from tissues," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 2000.
  • F. Destrempes, G. Cloutier, "A critical review and uniformized representation of statistical distributions modeling the ultrasound echo envelope," Ultrasound Med. Biol., 2010.
  • A. Rodriguez-Molares et al., "The generalized contrast-to-noise ratio: a formal definition for lesion detectability," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 67(4):745–759, 2020.
  • G. Montaldo, M. Tanter, J. Bercoff, N. Benech, M. Fink, "Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 2009.