Part IV · Module 25

Adaptive and coherence-based beamforming: MV, CF, GCF, PCF, DMAS, SLSC, FXPF

After the delays, every pixel is a vector of 64 channel samples. DAS adds them with fixed weights. Every method in this module looks at that vector first and then decides what to output. Here we derive each one, put all of them on one simulated dataset with honest metrics, and show which improvements are real and which come from dynamic-range stretching.

You will be able to

Derive the MV (Capon) weights and explain why ultrasound needs subarray smoothing, temporal averaging and diagonal loading, and what ESBMV changes
Relate CF, the Wiener postfilter, GCF, PCF and SCF to one aperture-domain picture, and predict their values for a point, speckle, noise and off-axis clutter
Explain DMAS as a \(2f_0\) sum-coarray beamformer, SLSC as coherence imaging and FXPF as linear prediction across the aperture
Run a fair comparison: matched data, stated parameters, FWHM, PSL, CR, CNR, gCNR, speckle SNR and cost, and spot the dynamic-range trap
Predict which methods lose most under aberration, and why coherence is also a focus metric

Assumes

Module 07 · Random vectors, covariance & correlation: covariance, VCZ coherence, CF basics, RMB and loading
Module 23 · Beam patterns, k-space and compounding: DAS, apodization, plane-wave compounding and the coarray
Optional further reading: Module 26 (speckle statistics; CR, CNR and gCNR are defined here where used)

Exact formulaclosed form plotted as is (identities, narrowband array factors, iid phase-error loss)SimulationJS in the page: plane-wave channel data of a speckle phantom synthesized from the delay law, all ten beamformers and metrics run on it; Monte Carlo trials; FFTsApprox. modelnarrowband or expected-covariance physics (VCZ triangle, single source in white noise, operation counts); assumptions under the badgeSchematicdrawing for intuition, not computedNo Field II or k-Wave output is shown; the MATLAB scripts in 25.12 re-simulate everything in Octave/MATLAB.

25.1One pixel, one aperture vector: where each method acts

Fix a pixel \(\mathbf{r} = (x, z)\). For every transmit angle θ and element i, read the IQ channel sample at the round-trip delay, rotate it by the carrier phase of that delay, and sum over angles. The result is the aperture vector \(\mathbf{s}(\mathbf{r}) \in \mathbb{C}^N\), with N = 64 here. Its real part is the delayed RF sample on each element, which is what DMAS and SLSC use.

\[ \begin{aligned} s_i(\mathbf{r}) &= \sum_\theta \mathrm{IQ}_{i,\theta}\big(\tau_{\theta,i}(\mathbf{r})\big)\, e^{\,j2\pi f_0 \tau_{\theta,i}(\mathbf{r})} \\ c\,\tau_{\theta,i}(\mathbf{r}) &= x\sin\theta + z\cos\theta \\ &\quad + \sqrt{(x - x_i)^{2} + z^{2}} \\ y_{\mathrm{DAS}}(\mathbf{r}) &= \mathbf{w}^{H}\mathbf{s} = \frac{1}{N}\sum_i s_i, \quad \mathbf{w} = \frac{\mathbf{1}}{N} \end{aligned} \]

(25.1)

\(s_i(\mathbf{r})\): delayed, phase-rotated, angle-compounded sample of element i at pixel \(\mathbf{r}\)
\(\mathrm{IQ}_{i,\theta}(t)\): IQ channel data of element i for transmit angle θ, interpolated at time t
\(\tau_{\theta,i}(\mathbf{r})\): round-trip delay: plane-wave transmit plus receive path to element \(x_i\), µs
\(f_0\): centre frequency, 5 MHz
\(N\): number of receive elements, 64
\(\mathbf{w}\): receive weights; DAS uses \(\mathbf{1}/N\)
\(x,\ z\): lateral and axial pixel coordinates, mm
\(c\): sound speed, 1540 m/s
\(\theta\): plane-wave transmit angle (summed over all angles), rad
\(j\): imaginary unit
\(\mathbf{s},\ y_{\mathrm{DAS}}\): vector of the \(s_i\) over the aperture; DAS output
\((\cdot)^{H}\): conjugate transpose
\(\mathbf{1}\): vector of ones

After the delays a focused echo is the same on every element, so the steering vector is \(\mathbf{a} = \mathbf{1}\). Everything off focus (off-axis echoes, clutter, noise) is a departure from a constant vector: a phase ramp, a random phase, or an amplitude pattern.

The ten methods differ only in what they do with \(\mathbf{s}\) before or instead of that sum (Figure 25.A):

Family	Methods	Output	What it estimates from s	Linear in the data?
Channel-data filter	FXPF, ADMIRE	DAS of filtered channels	The part of the aperture data a model can predict (plane waves, modelled scatterers)	No: coefficients are fitted to the data
Adaptive weights	MV (Capon), ESBMV	\(\mathbf{w}(\mathbf{s})^{H}\mathbf{s}\)	Covariance \(\mathbf{R}\) from neighbouring samples and subarrays; weights that null everything except the focus	No, but the output is still an amplitude with signal phase
Coherence weighting	CF, Wiener postfilter, GCF, PCF, SCF	\(W(\mathbf{s})\, y_{\mathrm{DAS}}\), \(W \in [0, 1]\)	How constant the aperture vector is, from energy, spatial spectrum or phase	No: a data-dependent gain
Nonlinear summation	DMAS, F-DMAS	\(\sum_{i<j}\operatorname{sign}(s_i s_j)\sqrt{\lvert s_i s_j\rvert}\), band-passed at \(2f_0\)	Pairwise agreement between channels	No: second order in the data, square root restores amplitude units
Coherence imaging	SLSC	\(\sum_{m=1}^{M} R(m)\), a correlation	Normalized spatial coherence at short element lags	No, and the output is not an amplitude

Schematic Figure 25.A · The aperture-domain data cube and where each beamformer acts: channel data (θ, i, t) are delayed and compounded into an element-by-depth slab per image column; each method reads a different part of it (dimensions of this module's dataset, not to scale)

Demo 25.1

The aperture at one pixel: what CF, GCF, PCF, SLSC and DMAS each see

SimulationAperture data at one pixel synthesized in JS for the chosen scene (point, speckle, noise, off-axis target), nine compounded plane waves, optional Gaussian-correlated phase screen; coherence curve over a ±0.15 mm axial kernel.

\(\mathrm{Re}\, s_i\)phase, ±σ bandGCF band \(\lvert k\rvert \le M_0\)\(R(m)\), SLSC lags 1…MVCZ envelope for this transmit sequence

Pixel content

Channel SNR re speckle 20 dB

Phase screen rms 0 ns

GCF cut-off M₀ 1

SLSC maximum lag M 8

Model: pixel at (0, 30) mm, 64 elements, 0.3 mm pitch, 5 MHz, 60 % bandwidth, nine plane waves from −4° to 4° coherently compounded; the coherence curve uses a ±0.15 mm axial kernel. Try: point → every element agrees, CF ≈ 1, the spectrum is one bar at k = 0 and \(R(m) \approx 1\) at all lags. Speckle → the spectrum spreads over \(\lvert k\rvert \lesssim 4\) and \(R(m)\) follows the VCZ envelope, which falls to zero near lag 16 because nine angles at 1° synthesize a 4.7 mm transmit aperture. Noise → phases fill ±π (σ → 1.81 rad), CF → 1/N, \(R(m \ge 1) \approx 0\). Off-axis → a clean phase ramp: the spectrum is one bar at k ≈ 4, so GCF rejects it until \(M_0 \ge 4\), while \(R(m)\) oscillates.

Intuition

Read the aperture spectrum as a map of arrival directions: bin k is a wave arriving from lateral offset \(\Delta x \approx k\lambda z/(Np)\), 0.48 mm per bin at 30 mm. DAS keeps bin 0 and nothing else. CF is the fraction of energy in bin 0. GCF widens the window to \(\lvert k\rvert \le M_0\). MV and ESBMV build weights that null whichever bins carry interference. FXPF keeps whatever few discrete bins (plane waves) explain the data best.

Methods card · DAS (delay and sum)64 el. · uniform · 9 angles

What it does	Sums the delayed, compounded aperture vector with equal weights, \(\lvert\sum_i s_i\rvert/N\).
Parameters in Demo 25.6	64 elements, rectangular apodization, no aperture growth (F/1.2–F/1.8); 9 plane waves, −4° to 4°.
Source of each value	Classical baseline; uniform full aperture: ours.
Why this value	Same receive aperture for every method, so differences come from the method; uniform weights give the narrowest mainlobe.
Known weaknesses	Sidelobes (−13.3 dB one-way) and off-axis clutter fill the cyst; FWHM 0.48 mm, PSL −19.4 dB, pair unresolved.
Compute cost / pixel	N = 64 adds (10²), after the shared delay and compounding (about 60 µs/pixel). Ours: 0.6 µs/pixel; your run: run Demo 25.6.
Tuned value	default already optimal under the rule; Hann apodization raises held-out gCNR to 0.942 but widens the point to 0.74 mm, beyond the 0.53 mm limit (protocol)

25.2Minimum variance (Capon) and the three fixes ultrasound needs

Module 07 introduced MVDR and the snapshot problem. Here is the beamformer as it is actually implemented in ultrasound (Synnevåg, Austeng & Holm, IEEE TUFFC 2007 and 2009). Minimize the output power subject to unit gain toward the focus:

\[ \begin{gathered} \min_{\mathbf{w}}\ \mathbf{w}^{H}\mathbf{R}\mathbf{w} \quad \text{s.t.} \quad \mathbf{w}^{H}\mathbf{a} = 1 \\ \Rightarrow\ \ \mathbf{w} = \frac{\mathbf{R}^{-1}\mathbf{a}}{\mathbf{a}^{H}\mathbf{R}^{-1}\mathbf{a}}, \qquad \mathbf{a} = \mathbf{1} \\ P_{\mathrm{MV}} = \frac{1}{\mathbf{a}^{H}\mathbf{R}^{-1}\mathbf{a}} \end{gathered} \]

(25.2)

\(\mathbf{w}\): adaptive receive weights (N, or L with subarrays)
\(\mathbf{R}\): covariance of the aperture vector, \(\mathrm{E}[\mathbf{s}\mathbf{s}^H]\)
\(\mathbf{a}\): steering vector toward the focus; all ones after the delays
\(P_{\mathrm{MV}}\): minimum output power
\((\cdot)^{H},\ (\cdot)^{-1}\): conjugate transpose, matrix inverse
\(\mathbf{1}\): vector of ones
\(\mathrm{E}[\cdot]\): expectation
\(\min,\ \text{s.t.}\): minimize over \(\mathbf{w}\), subject to the distortionless constraint \(\mathbf{w}^H\mathbf{a} = 1\)

Lagrangian \(\mathbf{w}^H\mathbf{R}\mathbf{w} + \lambda(1 - \mathbf{w}^H\mathbf{a}) + \lambda^*(1 - \mathbf{a}^H\mathbf{w})\); setting the gradient with respect to \(\mathbf{w}^*\) to zero gives \(\mathbf{R}\mathbf{w} = \lambda\mathbf{a}\) (Capon, Proc. IEEE 1969). With \(\mathbf{R} = \sigma^2\mathbf{I}\) the weights are \(\mathbf{1}/N\): DAS is MV in white noise.

Three things make the textbook formula fail on echo data, and each has a standard fix:

\[ \begin{aligned} \hat{\mathbf{R}} &= \frac{1}{(2K+1)P}\sum_{k=-K}^{K}\sum_{l=0}^{N-L} \\ &\qquad \mathbf{s}_l(t+k)\,\mathbf{s}_l^{H}(t+k) + \varepsilon\mathbf{I} \\ \varepsilon &= \Delta\cdot\mathrm{tr}(\hat{\mathbf{R}})/L \\ y_{\mathrm{MV}}(t) &= \frac{1}{P}\sum_l \mathbf{w}^{H}\mathbf{s}_l(t) \\ P &= N - L + 1 \\ \mathbf{s}_l &= [s_l, \dots, s_{l+L-1}]^{T} \end{aligned} \]

(25.3)

\(\mathbf{s}_l\): subarray of L consecutive elements starting at element l
\(L\): subarray length (adaptive aperture), here 32
\(P\): number of subarrays, \(N - L + 1 = 33\)
\(K\): half-length of the temporal average, samples (here 2 samples of 75 µm)
\(\Delta, \varepsilon\): loading as a fraction of the mean element power, and its absolute value
\(\hat{\mathbf{R}}\): estimated covariance matrix, \(L\times L\)
\(t,\ k\): depth sample index and temporal-average offset, samples
\(l\): subarray start index
\(\mathbf{I},\ \mathrm{tr}\): identity matrix and trace
\(y_{\mathrm{MV}},\ \mathbf{w}\): MV output and weights from (25.2), computed with \(\hat{\mathbf{R}}\)
\(N\): number of receive elements, as in (25.1)
\((\cdot)^{T}\): transpose

(1) Coherent signals: the echo of interest and the echoes that should be suppressed are copies of one transmit pulse, so they are fully correlated. Forward subarray smoothing (L-element subarrays, \(L \le N/2\) usually) decorrelates them (Shan, Wax & Kailath, IEEE Trans. ASSP 1985). Forward-backward averaging, \(\hat{\mathbf{R}}_{\mathrm{FB}} = \tfrac12(\hat{\mathbf{R}} + \mathbf{J}\hat{\mathbf{R}}^{T}\mathbf{J})\) with \(\mathbf{J}\) the exchange (flip) matrix, doubles the effective snapshot count for the symmetric steering vector \(\mathbf{a} = \mathbf{1}\) and is a common addition (Asl & Mahloojifar, IEEE TUFFC 2011); this page uses forward smoothing only. (2) Few snapshots: temporal averaging over ±K samples, of the order of the pulse length, raises the snapshot count to \((2K+1)(N-L+1)\). (3) Model errors and noise: diagonal loading ε, here \(\Delta = 10^{-2}\) of the mean element power (Cox, Zeskind & Owen, IEEE Trans. ASSP 1987). Our engine uses L = 32, K = 2 samples of 75 µm (±0.19 µs), Δ = 0.01.

Why smoothing is not optional. Let the desired echo be s and a coherent interferer be βs with steering \(\mathbf{v}\). With one snapshot model \(\mathbf{x} = s(\mathbf{a} + \beta\mathbf{v})\), \(\mathbf{R} = \sigma^2\mathbf{u}\mathbf{u}^H\) with \(\mathbf{u} = \mathbf{a} + \beta\mathbf{v}\). Minimizing \(\mathbf{w}^H\mathbf{R}\mathbf{w}\) drives \(\mathbf{w}^H\mathbf{u} \to 0\), and with \(\mathbf{w}^H\mathbf{a} = 1\) that forces \(\beta\,\mathbf{w}^H\mathbf{v} = -1\): the interferer is used to cancel the signal (Widrow et al., IEEE Trans. Antennas Propag. 1982). Subarray l sees the interferer with an extra phase \(e^{\,j\psi l}\), so averaging over \(P = N - L + 1\) subarrays multiplies the cross term by the Dirichlet factor:

\[ \begin{gathered} D_P(\psi) = \frac{1}{P}\sum_{l=0}^{P-1} e^{\,j\psi l} \\ \lvert D_P\rvert \ll 1\ \ \text{once}\ \ P\psi \gtrsim 2\pi \\ \psi = \frac{2\pi p\,\Delta x}{\lambda z} \end{gathered} \]

(25.4)

\(\psi\): phase step per element of an interferer at lateral offset Δx in the focal plane at depth z, rad
\(p\): element pitch, 0.3 mm
\(D_P(\psi)\): residual coherence between desired echo and interferer after smoothing
\(P\): number of subarrays, as in (25.3)
\(l\): subarray index
\(\Delta x\): lateral offset of the interferer, mm
\(z,\ \lambda\): depth and wavelength, mm
\(j\): imaginary unit

Methods card · MV (Capon)L = 32 · K = ±2 · Δ = 0.01

What it does	Capon weights from a smoothed, loaded covariance; output \(\mathbf{w}^H\bar{\mathbf{s}}\) (25.3).
Parameters in Demo 25.6	L = 32 (33 subarrays); K = ±2 rows of 75 µm; Δ = 0.01, ε = Δ·tr(R̂)/L; forward smoothing only.
Source of each value	Method, L ≤ N/2, K ≈ pulse length: Synnevåg et al. 2007 doi:10.1109/TUFFC.2007.431; Synnevåg et al. 2009 doi:10.1109/TUFFC.2009.1263. Loading: Cox et al. 1987 doi:10.1109/TASSP.1987.1165054. Values: ours.
Why this value	L = N/2: best resolution in the usual range; ±0.15 mm ≈ one pulse; 1 % loading stabilizes the solve without pulling MV toward DAS.
Known weaknesses	Lone-point width measures SNR, not resolution; points vanish between pixels; speckle SNR 1.53 (DAS 1.86); K = 0 collapses speckle.
Compute cost / pixel	≈ 1.6·10⁴ complex MACs. Ours: 158 µs/pixel; your run: run Demo 25.6.
Tuned value	tuned: L = 40, K = ±4, Δ = 10⁻⁴ (gCNR 0.978 → 0.982) (protocol)

\[ \begin{aligned} \hat{\mathbf{R}} &= \sum_i \lambda_i\,\mathbf{e}_i\mathbf{e}_i^{H} \\ \mathbf{E}_s &= \left[\,\mathbf{e}_i : \lambda_i \ge \delta\lambda_1\,\right] \\ \mathbf{w}_{\mathrm{ESB}} &= \mathbf{E}_s\mathbf{E}_s^{H}\,\mathbf{w}_{\mathrm{MV}} \end{aligned} \]

(25.5)

\(\lambda_i, \mathbf{e}_i\): eigenvalues (descending) and eigenvectors of the unloaded \(\hat{\mathbf{R}}\)
\(\mathbf{E}_s\): "signal subspace": eigenvectors within a factor δ of the largest
\(\delta\): relative threshold, here 0.5
\(\hat{\mathbf{R}}\): as in (25.3)
\(\lambda_1\): largest eigenvalue
\(\mathbf{w}_{\mathrm{MV}},\ \mathbf{w}_{\mathrm{ESB}}\): MV weights from (25.2) and their projection onto the signal subspace
\((\cdot)^{H}\): conjugate transpose

ESBMV: Asl & Mahloojifar (IEEE TUFFC 2010). Projecting the weights onto the dominant eigenvectors removes the part of \(\mathbf{w}\) that only fits noise. δ is a free parameter: δ → 0 returns MV. In our speckle with δ = 0.5, on average 1.3 of 32 eigenvectors are kept; the engine finds the four dominant ones by warm-started subspace iteration (checked against a full eigendecomposition: 0.01 % mean error in speckle). The projection does not preserve \(\mathbf{w}^H\mathbf{a} = 1\); Section 25.6 shows what that does next to a bright scatterer.

Methods card · ESBMV (eigenspace MV)δ = 0.5 · MV settings

What it does	Projects the MV weights onto the dominant eigenvectors, \(\mathbf{E}_s\mathbf{E}_s^H\mathbf{w}_{\mathrm{MV}}\) (25.5).
Parameters in Demo 25.6	δ = 0.5 (keep λ_i ≥ 0.5 λ₁ of the unloaded R̂), at most the 4 largest eigenpairs; L, K, Δ as MV.
Source of each value	Method, δ free: Asl & Mahloojifar 2010 doi:10.1109/TUFFC.2010.1706. δ and the 4-vector cap: ours.
Why this value	Mid-range of the Demo 25.7b sweep; speckle keeps 1.3 of 32 eigenvectors on average, so the cap rarely binds.
Known weaknesses	Drops w^Ha = 1: black band beside the bright scatterer (−52 dB re DAS); CR 70.5 dB is mostly stretch (gCNR 0.988 vs 0.958); speckle SNR 1.22.
Compute cost / pixel	≈ 4·10⁴ (MV + subspace iteration). Ours: 362 (incl. MV) µs/pixel; your run: run Demo 25.6.
Tuned value	tuned: L = 48, δ = 0.3 (K, Δ as MV default) (gCNR 0.988 → 0.991) (protocol)

Demo 25.2

MV weights against an interferer: beam pattern, smoothing and loading

Approx. modelNarrowband, focal plane, expected smoothed covariance of one desired source plus one interferer in white noise (channel SNR 20 dB); weights by a matrix solve in JS.

DAS (uniform)MV, effective full-aperture responseinterferer response

Interferer

Lateral offset 1.5 mm

Interferer re signal 0 dB

Subarray length L 32

Loading Δ 10^−2

Model: focal plane at z = 30 mm, 64 elements, 0.3 mm pitch, 5 MHz narrowband; an interferer at offset Δx gives a phase step \(\psi = 2\pi p\Delta x/(\lambda z)\) per element (25.4); expected smoothed covariance with channel SNR 20 dB. The effective response is \(\mathbf{w}^H\mathbf{v}_L(x)\cdot D_P(\psi(x))\). Try: coherent, L = 64 (no smoothing): the signal gain falls to about −66 dB, signal cancellation. Drop to L = 48 and it recovers to −0.3 dB. Move the interferer to 0.3 mm: now L must be smaller before the subarrays decorrelate it. Switch to incoherent: cancellation disappears at any L; the null deepens with INR. Push Δ to 10: the pattern becomes DAS with a trapezoidal taper (box L convolved with box N − L + 1), white-noise gain 16.9 dB instead of 18.1.

Pitfall: the MV point is an SNR meter

For one point in white noise the MV output has a closed form, \(\lvert y(x)\rvert = \sigma^2\lvert\mathbf{v}^H\mathbf{a}\rvert/Q\) with \(Q = N(\sigma^2 + N) - \lvert\mathbf{v}^H\mathbf{a}\rvert^2\) with signal power 1 per channel, whose −6 dB width is \(0.914/\sqrt{N\cdot\mathrm{SNR}}\) times the DAS width (it matches a direct solve to 1 % at 10, 20 and 30 dB: 0.036, 0.011 and 0.0036 of the DAS width). A lone point therefore reports SNR and loading, not resolution: loading caps the effective SNR near 1/Δ. Judge resolution with two points (the pair dip in Demo 25.6) and sample laterally finer than the MV peak, or bright points vanish between pixels (−23 dB at half a 0.15 mm pixel in our dataset).

25.3Coherence weighting: CF, the Wiener postfilter, GCF, PCF and SCF

All five multiply the DAS output by a number in [0, 1] computed from the same aperture vector. Module 07 gave CF's values for a focused transmit (point 1, speckle 2/3, noise 1/N). Here we connect the five and compute them for the plane-wave sequence we image with.

\[ \begin{aligned} \mathrm{CF} &= \frac{\lvert\sum_i s_i\rvert^{2}}{N\sum_i \lvert s_i\rvert^{2}} \\ &= \frac{\lvert\bar{y}\rvert^{2}}{\lvert\bar{y}\rvert^{2} + \hat{\sigma}^{2}} \\ \bar{y} &= \frac{1}{N}\sum_i s_i \\ \hat{\sigma}^{2} &= \frac{1}{N}\sum_i \lvert s_i - \bar{y}\rvert^{2} \\ H_{\mathrm{Wiener}} &= \frac{\lvert\bar{y}\rvert^{2}}{\lvert\bar{y}\rvert^{2} + \hat{\sigma}^{2}/N} \\ &= \frac{N\cdot\mathrm{CF}}{(N - 1)\,\mathrm{CF} + 1} \end{aligned} \]

(25.6)

\(\mathrm{CF}\): coherence factor, in [0, 1]; its expected value for white noise is 1/N
\(\bar{y}\): aperture mean, equal to the DAS output
\(\hat{\sigma}^2\): spread of the channels about the mean (channel "noise" estimate)
\(H_{\mathrm{Wiener}}\): Wiener postfilter gain applied to \(y_{\mathrm{DAS}}\)
\(s_i,\ N\): as in (25.1)
\(\lvert\cdot\rvert\): magnitude

The identity \((1/N)\sum\lvert s_i\rvert^2 = \lvert\bar{y}\rvert^2 + \hat{\sigma}^2\) turns CF into a signal-to-(signal + noise) ratio. The Wiener postfilter of Nilsen & Holm (IEEE TUFFC 2010) has the same estimates but divides the channel noise by N, because DAS averages N independent noise samples. CF therefore overestimates the output noise by a factor N and suppresses far harder: for pure noise CF = 1/N (−36 dB on an amplitude image), \(H = N/(2N - 1) \approx 0.5\) (−6 dB).

Methods card · CF (coherence factor)no free parameter

What it does	Weights DAS by \(\lvert\sum_i s_i\rvert^2/(N\sum_i\lvert s_i\rvert^2)\) (25.6).
Parameters in Demo 25.6	None; 64 channels, one sample per pixel.
Source of each value	Definition: Mallart & Fink 1994 doi:10.1121/1.410562.
Why this value	Parameter-free reference for the coherence weights.
Known weaknesses	Noise → 1/N, speckle → 0.25: CR 52.3 dB is stretch; gCNR 0.935 (DAS 0.958); speckle SNR 0.91 (dark holes); aberration lowers it further.
Compute cost / pixel	≈ 2N (10²). Ours: 1.0 µs/pixel; your run: run Demo 25.6.
Tuned value	tuned: no CF weight (exponent 0, i.e. DAS) (gCNR 0.935 → 0.958) (protocol)

Methods card · Wiener postfilterno free parameter

What it does	Weights DAS by \(H = N\,\mathrm{CF}/((N-1)\mathrm{CF}+1)\): CF with the output noise divided by N (25.6).
Parameters in Demo 25.6	None; ȳ, σ̂² from one aperture sample per pixel, uniform weights.
Source of each value	H: Nilsen & Holm 2010 doi:10.1109/TUFFC.2010.1553. Single-sample DAS estimator: ours (simplest form).
Why this value	Shows the noise-correct version of CF on the same data.
Known weaknesses	Noise only −6 dB (H ≈ 0.5); coherent clutter passes (H ≈ 1); noisy one-sample estimate. CR 31.1 dB, gCNR 0.943.
Compute cost / pixel	as CF (10²). Ours: 0.3 µs/pixel; your run: run Demo 25.6.
Tuned value	default already optimal under the rule (no free parameter in our implementation) (protocol)

\[ \begin{aligned} \mathrm{GCF} &= \frac{\sum_{\lvert k\rvert \le M_0} \lvert S(k)\rvert^{2}}{\sum_k \lvert S(k)\rvert^{2}} \\ S(k) &= \sum_i s_i\, e^{-j2\pi k i/N} \end{aligned} \]

(25.7)

\(S(k)\): aperture spectrum (DFT across elements); bin k is a wave from lateral offset \(k\lambda z/(Np)\)
\(M_0\): cut-off bin; \(M_0 = 0\) gives CF by Parseval
\(s_i,\ N\): as in (25.1)
\(i\): element index
\(j\): imaginary unit
\(z,\ p,\ \lambda\): depth, element pitch and wavelength

\[ \begin{aligned} \mathrm{PCF} &= \max\!\Big(0,\ 1 - \frac{\gamma}{\sigma_0}\,\sigma_{\min}\Big) \\ \sigma_{\min} &= \min\big(\sigma(\varphi), \sigma(\varphi_A)\big) \\ \sigma_0 &= \frac{\pi}{\sqrt{3}}, \qquad \varphi_A = \varphi - \pi\operatorname{sign}(\varphi) \\ \mathrm{SCF} &= \big\lvert 1 - \sqrt{1 - \bar{b}^{2}}\,\big\rvert^{p} \\ \bar{b} &= \frac{1}{N}\sum_i \operatorname{sign}(\mathrm{Re}\, s_i) \end{aligned} \]

(25.8)

\(\sigma(\varphi)\): standard deviation of the channel phases across the aperture, rad
\(\sigma_{\min}\): the smaller of the two phase spreads
\(\varphi_A\): auxiliary phase that handles wrapping at ±π
\(\gamma\): PCF sensitivity, here 1
\(\bar{b}\): mean of the channel signs (one bit per channel)
\(p\): SCF exponent
\(\varphi\): phase of the channel samples, rad
\(\mathrm{PCF},\ \mathrm{SCF}\): phase and sign coherence factors, in [0, 1]
\(\sigma_0\): std of phases uniform on \([-\pi, \pi]\), \(\pi/\sqrt{3}\), rad
\(\operatorname{sign},\ \mathrm{Re}\): sign function, real part
\(\min,\ \max\): minimum and maximum

GCF: Li & Li (IEEE TUFFC 2003); it accepts slow phase variation (speckle, mild aberration) that CF penalizes. PCF and SCF: Camacho, Parrilla & Fritsch (IEEE TUFFC 2009). \(\sigma_0 = 2\pi/\sqrt{12}\) is the standard deviation of a phase uniform on (−π, π], so PCF reaches 0 for noise at γ = 1. SCF needs one bit per channel: the standard deviation of ±1 signs is \(\sqrt{1 - \bar{b}^2}\).

For our sequence (nine angles over ±4°), the VCZ coherence of speckle falls to zero near lag \(m_{\mathrm{tx}} = N_\theta z\,\Delta(\sin\theta)/p \approx 15.7\), so the expected speckle CF is \(\sum_m (N - \lvert m\rvert)C(m)/N^2 = 0.225\), far below the 2/3 of a fully focused transmit; the simulated per-pixel mean is 0.25. The noise floor stays at 1/N. CF therefore pushes speckle down by about 13 dB and noise by 36 dB: a 23 dB stretch of the dynamic range before anything is resolved better. GCF with \(M_0 = 1\) keeps 60 % of the expected speckle aperture energy and rejects an off-axis wavefront from 2 mm, which lands at \(k = Np\,\Delta x/(\lambda z) = 4.2\).

Methods card · GCF (generalized CF)M₀ = 1

What it does	Weights DAS by the share of aperture-spectrum energy in bins \(\lvert k\rvert \le M_0\) (25.7).
Parameters in Demo 25.6	M₀ = 1: bins −1, 0, 1 of a 64-point DFT (spacing ≈ 0.48 mm at 30 mm); one sample per pixel.
Source of each value	Method: Li & Li 2003 doi:10.1109/TUFFC.2003.1182117. M₀: ours.
Why this value	Smallest M₀ ≠ CF: keeps 60 % of speckle energy (CF 22 %), still rejects a wavefront from 2 mm (k ≈ 4.2).
Known weaknesses	Arbitrary cut-off; no mainlobe narrowing (FWHM 0.43 mm); aberration spreads energy past M₀. gCNR 0.975, speckle SNR 1.50.
Compute cost / pixel	≈ 3N (10²). Ours: 2.4 µs/pixel; your run: run Demo 25.6.
Tuned value	default kept under the rule: M₀ = 2 scores 0.0006 higher on the held-out data, within the 0.005 tie band (protocol)

Methods card · PCF (phase coherence factor)γ = 1

What it does	Weights DAS by \(\max(0, 1 - \gamma\sigma_{\min}/\sigma_0)\), the spread of channel phases (25.8).
Parameters in Demo 25.6	γ = 1; σ₀ = π/√3; auxiliary phase φ − π·sign φ; 64 complex samples, one per pixel.
Source of each value	Method, σ₀, auxiliary phase: Camacho et al. 2009 doi:10.1109/TUFFC.2009.1128. γ: ours.
Why this value	At γ = 1, uniformly random (noise) phases get weight 0.
Known weaknesses	Clipping to 0 gives dark holes (speckle SNR 0.99); γ-critical (1.5 zeros most tissue); CR 41.9 dB but gCNR 0.914, below DAS.
Compute cost / pixel	N phases + sums (10²–10³ flops). Ours: 3.3 µs/pixel; your run: run Demo 25.6.
Tuned value	tuned: γ = 0.25 (gCNR 0.914 → 0.957) (protocol)

Ultrasound link

CF-type weights were introduced as focusing-quality measures (Mallart & Fink, JASA 1994; Hollman, Rigby & O'Donnell, IEEE Ultrasonics Symposium 1999). The same number that darkens a cyst also scores focus: averaged over a region, CF or GCF peaks at the right sound speed and drops with aberration (Module 07, Module 14). Used as a per-pixel gain, the weight inherits the speckle variance of its own estimate: the CF of speckle is broadly distributed (card above), so CF images show dark holes in tissue and a speckle SNR well below 1.91.

25.4Nonlinear summation and coherence imaging: DMAS and SLSC

DMAS (Matrone, Savoia, Caliano & Magenes, IEEE TMI 2015) multiplies every pair of delayed RF channels, restores amplitude units with a signed square root, sums, and band-passes at \(2f_0\) (F-DMAS):

\[ \begin{aligned} y_{\mathrm{DMAS}}(t) &= \sum_{i<j} \operatorname{sign}(s_i s_j)\sqrt{\lvert s_i s_j\rvert} \\ &= \frac{1}{2}\Big[\Big(\sum_i \hat{s}_i\Big)^{2} - \sum_i \hat{s}_i^{2}\Big] \\ \hat{s}_i &= \operatorname{sign}(s_i)\sqrt{\lvert s_i\rvert} \end{aligned} \]

(25.9)

\(s_i\): delayed real RF sample of element i (Re of the aperture vector)
\(\hat{s}_i\): signed square root of \(s_i\)
\(y_{\mathrm{DMAS}}\): DMAS output before the \(2f_0\) band-pass
\(\sum_{i<j}\): sum over all distinct element pairs
\(\operatorname{sign}\): sign function
\(t\): depth (fast time), µs

The second form costs O(N) instead of \(N(N - 1)/2 = 2016\) pair products for N = 64. \(s_i\) is real RF here, sampled every 18.75 µm ≈ λ/16 in depth (41 MHz equivalent) so that the \(2f_0\) band is not aliased. For an aligned echo every pair term equals \(\lvert \mathrm{rf}(t)\rvert\), whose spectrum has a DC term, a \(2f_0\) band and weaker \(4f_0\): the band-pass (here Gaussian, centred at \(2f_0\) with \(\sigma = \sqrt{2}\,\sigma_f\)) keeps only the \(2f_0\) part.

Why it narrows the beam: in a narrowband picture the \(2f_0\) term of the pair (i, j) carries the phase \(\varphi_i + \varphi_j\), as if a virtual element at \((x_i + x_j)/2\) received at twice the frequency. The pairs fill the same aperture with a triangular weight (the sum coarray, compare Module 23). A triangle at λ/2 has a −6 dB width 0.74 times that of a uniform aperture at λ and −26.5 dB sidelobes instead of −13.3 dB (computed for N = 64). Look at the algebra and the reason is plain:

\[ \sum_{i<j} e^{\,j(\varphi_i + \varphi_j)} = \tfrac{1}{2}\Big[\Big(\sum_i e^{\,j\varphi_i}\Big)^{2} - \sum_i e^{\,j2\varphi_i}\Big] \]

(25.10)

\(\varphi_i\): narrowband phase of element i after the delays
\(j\): imaginary unit
\(\sum_{i<j}\): sum over all distinct element pairs

So the narrowband DMAS beam is the DAS beam squared (minus a small self term).

Its −6 dB width is the DAS −3 dB width (0.886/1.207 ≈ 0.74), and its sidelobes are DAS's doubled in dB. For an isolated point that is exactly what a power law p = 2 on the DAS image would show; DMAS differs because the signed square root keeps the output linear in scatterer amplitude, so speckle and targets of different strength are not stretched like DAS². On our broadband plane-wave data the gain is smaller: FWHM 0.43 mm vs 0.48 mm, peak sidelobe −27 vs −19 dB, because the transmit part of the PSF is shared and the square-root compression mixes harmonics.

Methods card · F-DMAS (filtered delay multiply and sum)2f₀ Gaussian band-pass

What it does	Sums signed-root pair products (O(N) form, (25.9)), band-passes at 2f₀, envelope.
Parameters in Demo 25.6	RF = Re s on 18.75 µm rows (41 MHz); Gaussian band-pass at 10 MHz, σ = √2 σ_f = 1.80 MHz; output / [N(N−1)/2].
Source of each value	Method, 2f₀ band-pass: Matrone et al. 2015 doi:10.1109/TMI.2014.2371235. Filter shape and width: ours.
Why this value	√2 σ_f is the width of the 2f₀ band of rf² for our Gaussian pulse; the fine grid keeps 2f₀ unaliased.
Known weaknesses	Small gain on broadband compounded data (FWHM 0.43 vs 0.48 mm); the signed root mixes harmonics; speckle SNR 1.24.
Compute cost / pixel	≈ 2N × 4 RF rows + FFT filter (10³). Ours: 5.0 µs/pixel; your run: run Demo 25.6.
Tuned value	default already optimal under the rule (no free parameter in our implementation) (protocol)

SLSC (Lediju, Trahey, Byram & Dahl, IEEE TUFFC 2011) displays coherence instead of amplitude:

\[ \begin{aligned} R(m) &= \frac{1}{N - m}\sum_{i=1}^{N-m} \rho_i(m) \\ \rho_i(m) &= \frac{\sum_n s_i(n)\, s_{i+m}(n)}{\sqrt{\sum_n s_i^{2}(n)\,\sum_n s_{i+m}^{2}(n)}} \\ R_{\mathrm{sl}} &= \sum_{m=1}^{M} R(m) \end{aligned} \]

(25.11)

\(R(m)\): normalized spatial coherence at element lag m
\(\rho_i(m)\): normalized correlation of elements i and i + m over the kernel
\(n\): sample index over an axial kernel (here ±0.15 mm, 17 RF samples)
\(M\): maximum lag, here 8 (12.5 % of N)
\(R_{\mathrm{sl}}\): short-lag spatial coherence, the displayed value
\(i\): element index, \(1\dots N\)
\(N\): number of receive elements
\(s_i(n)\): RF sample \(n\) of element \(i\)

M is commonly 5–30 % of N. A point gives \(R(m) \approx 1\), so \(R_{\mathrm{sl}} \approx M\); fully developed speckle follows the VCZ curve (\(R_{\mathrm{sl}} \approx \sum(1 - m/m_{\mathrm{tx}}) = 5.7\) for M = 8; 6.0 in the noise-free simulation); noise and spatially incoherent clutter give \(R(m \ge 1) \approx 0\). Bell, Dahl & Trahey (IEEE TUFFC 2015) analysed SLSC resolution and brightness versus M.

Two consequences follow from the definition. SLSC is normalized per channel pair, so a bright off-axis scatterer's sidelobe in a dark region is still incoherent there and stays dark, and weak but coherent tissue stays bright: that is the clutter robustness. And a point target is wide in SLSC (about 1.9 mm here), because coherence stays near 1 wherever the point dominates the noise, even far down its DAS sidelobes. Robust variants (for example lag-weighted sums that emphasize the most reliable lags) trade some of the brightness for resolution.

Methods card · SLSC (short-lag spatial coherence)M = 8 · kernel ±0.15 mm

What it does	Displays summed RF coherence over lags 1…M (25.11), not amplitude.
Parameters in Demo 25.6	M = 8 (12.5 % of N); kernel ±8 RF rows of 18.75 µm (±0.15 mm); A = max(R_sl, 0), linear.
Source of each value	Method, M = 5–30 % of N: Lediju et al. 2011 doi:10.1109/TUFFC.2011.1957. M, kernel: ours.
Why this value	gCNR peaks near M = 4–8 here; kernel ≈ one wavelength (0.31 mm).
Known weaknesses	Not amplitude: point FWHM 1.87 mm, pair unresolved, CR 15.6 dB; CR and CNR not comparable with the others.
Compute cost / pixel	N·M·(2K_s+1) ≈ 8.7·10³ (10⁴). Ours: 25.8 µs/pixel; your run: run Demo 25.6.
Tuned value	tuned: M = 32, kernel ±0.30 mm (gCNR 0.991 → 0.963); the default point FWHM (1.87 mm) fails the 0.53 mm limit (protocol)

25.5Channel-data filters: f-x prediction (FXPF) and the ADMIRE idea

FXPF borrows f-x deconvolution from seismic noise attenuation (Canales, SEG Expanded Abstracts, 1984) and applies it to the delayed channel data before DAS (Shin & Huang, Proc. SPIE 9790, 2016; Shin, Huang & Yen, IEEE TUFFC 65(10), 2018). In a short depth window, take the FFT of each channel along time. At each frequency the aperture is a complex sequence \(X_1(f), \dots, X_N(f)\). A sum of q plane waves across the aperture is exactly predictable by a linear recursion of order q (the Prony / linear-prediction property). Fit that recursion by least squares, replace the data by its prediction, invert the FFT and beamform:

\[ \begin{aligned} X_i(f) &\approx \sum_{k=1}^{p} a_k(f)\,X_{i-k}(f) \\ X_i(f) &\approx \sum_{k=1}^{p} b_k(f)\,X_{i+k}(f) \\ \hat{\mathbf{a}} &= \arg\min_{\mathbf{a}} \sum_i \Big\lvert X_i - \sum_k a_k X_{i-k}\Big\rvert^{2} \\ &= (\mathbf{A}^{H}\mathbf{A} + \mu\mathbf{I})^{-1}\mathbf{A}^{H}\mathbf{x} \end{aligned} \]

(25.12)

\(X_i(f)\): FFT along depth of channel i in one window, at frequency f
\(a_k, b_k\): forward (first line) and backward (second line) prediction coefficients across the aperture
\(p\): prediction (AR) order, here 4
\(\mathbf{A}, \mathbf{x}\): regressor matrix of shifted channels and target vector
\(\mu\): regularization, \(10^{-3}\) of the mean regressor power
\(i,\ f\): channel index and frequency
\(\hat{\mathbf{a}}\): estimated forward prediction coefficients
\(\arg\min\): value of the argument minimizing the cost
\(\mathbf{I},\ (\cdot)^{H}\): identity matrix, conjugate transpose

The cleaned channel \(\hat{X}_i\) is the mean of the available predictions. Our implementation: Hann windows of W = 16 output samples (1.2 mm) with 50 % overlap, p = 4, \(\mu = 10^{-3}\) of the mean regressor power, forward and backward predictions averaged (only one is available within p elements of either edge).

What it removes and what it cannot. The focused echo is constant across the aperture and survives exactly. An off-axis echo is a linear phase ramp: also predictable, so FXPF keeps coherent off-axis clutter and leaves its suppression to the DAS sum. What it removes is spatially incoherent energy: electronic noise, and reverberation or multipath clutter that does not look like a few plane waves. A second, less advertised effect: when the channel SNR is low the fitted coefficients shrink toward zero (regression dilution), so weak echoes are attenuated too. FXPF is therefore an SNR-dependent gain, like the Wiener postfilter, but built from prediction rather than energy ratios.

Methods card · FXPF (f-x prediction filter)p = 4 · W = 16 · μ = 10⁻³

What it does	Replaces each channel spectrum by linear prediction across the aperture in short depth windows, then DAS (25.12).
Parameters in Demo 25.6	p = 4; W = 16 rows (1.2 mm), Hann, 50 % overlap; μ = 10⁻³ × mean regressor power; forward and backward averaged.
Source of each value	Method: Shin & Huang 2016; Shin, Huang & Yen 2018 doi:10.1109/TUFFC.2018.2860962. p, W, μ: ours.
Why this value	p ≪ N predicts a few plane waves exactly: noise −14 dB, echo −2 dB at 0 dB SNR (Demo 25.5); W spans several pulses.
Known weaknesses	Keeps coherent off-axis clutter (CR +1.4 dB at −20 dB noise); regression dilution weakens low-SNR echoes.
Compute cost / pixel	≈ 6·10³ (10³–10⁴). Ours: 76.9 µs/pixel; your run: run Demo 25.6.
Tuned value	tuned: p = 3, W = 32 (2.4 mm) (gCNR 0.962 → 0.966) (protocol)

Demo 25.5

f-x prediction at one frequency bin

Simulation64-channel focused echo plus white complex noise, forward-backward AR(p) least-squares prediction in JS; lower plot is a Monte Carlo (60 trials per SNR for the echo, 1000 for noise).

noisy dataprediction \(\hat{X}\)noise-free signallower plot: solid = echo gain, dashed = noise-only gain

AR order p 4

Channel SNR 0 dB

add an off-axis wavefront (2.5 mm)

Model: 64 channels, unit-amplitude focused echo plus white complex noise; the lower plot is a Monte Carlo of the DAS-output gain on the echo (solid, 60 trials per SNR) and on noise alone (dashed, 1000 trials). Try: at 30 dB the prediction is the signal. At 0 dB, p = 4 removes most of the noise and costs the echo about 2 dB. Add the off-axis wavefront: it is predicted and kept, which is why FXPF does not remove off-axis clutter. Set p = 1: about 21 dB of noise removed, but the echo now loses about 6 dB at 0 dB SNR. Converged Monte Carlo (4000 trials) gives noise-only gains of −20.9, −17.5, −14.1 and −10.4 dB (±0.3 dB) for p = 1, 2, 4, 8, 2–3 dB below the forward-only projection estimate \(10\log_{10}(p/(N - p))\) = −18.0, −14.9, −11.8, −8.5 dB; the plotted 1000-trial values agree within about 0.5 dB.

Model-based suppression: ADMIRE

FXPF asks whether the aperture data are predictable; ADMIRE (aperture domain model image reconstruction; Byram, Dei, Tierney & Dumont, IEEE TUFFC 62(11), 2015) asks where each part of them came from. In every short depth window and at every frequency bin of the pulse band it writes the delayed aperture vector as a sum of modelled wavefronts: echoes of scatterers inside a small region of interest (ROI) around the pixel, and echoes of scatterers outside it, at other lateral positions and with extra path length (off-axis and multipath sources). An elastic-net regression picks a sparse set of those wavefronts; the channel data are rebuilt from the ROI wavefronts only, transformed back to time and summed (DAS). Because the model knows the direction and timing of each wavefront, it can remove coherent off-axis clutter that FXPF keeps. Everything below follows the authors' released code (VU-BEAM-Lab/ADMIRE on GitHub, Apache-2.0: model space, window and bin rules, α, \(c_\lambda\), solver), applied to this module's plane-wave data; where we had to adapt it, the table says so.

\[ \begin{aligned} \mathbf{y}(f) &= \mathbf{X}_{\mathrm{ROI}}\,\mathbf{\beta}_{\mathrm{ROI}} + \mathbf{X}_{\mathrm{out}}\,\mathbf{\beta}_{\mathrm{out}} + \mathbf{\varepsilon} \\ X_{ej} &= \frac{A_{ej}\, e^{-jk\zeta_{ej}}}{\lVert \mathbf{A}_{j} \rVert}, \quad k = \frac{2\pi f\kappa}{c} \\ \sqrt{u_e^{2} + \zeta_{ej}^{2}/4} + \frac{\zeta_{ej}}{2} &= \sqrt{(u_e - x_j)^{2} + z_j^{2}} + d_j + 2z_c - z_j \\ A_{ej} &= \sqrt{G_{ej}/\max_{e} G_{ej}}, \quad G_{ej} = \int_{p^-_{ej}}^{p^+_{ej}} e^{-\left(\frac{\zeta - \zeta_{ej}}{c\,\sigma_t}\right)^{2}} d\zeta \end{aligned} \]

(25.13)

\(\mathbf{y}(f)\): aperture vector: one STFT bin (frequency f) of the delayed channels in one depth window, \(N = 64\) complex values
\(\mathbf{X}_{\mathrm{ROI}}, \mathbf{X}_{\mathrm{out}}\): model matrices: columns are modelled aperture signals of scatterers inside / outside the ROI
\(\mathbf{\beta}_{\mathrm{ROI}}, \mathbf{\beta}_{\mathrm{out}}\): complex coefficients (one per modelled scatterer and path)
\(\mathbf{\varepsilon}\): misfit (noise and whatever the model cannot represent)
\(X_{ej}\): element e of column j, normalized so that each column has unit norm
\(e, j\): element index (1…64) and predictor index
\(A_{ej}, \mathbf{A}_j\): amplitude of predictor j on element e (column vector \(\mathbf{A}_j\)): square root of the fraction of the pulse energy that falls inside the window
\(k, \kappa\): wavenumber at the bin, rad/mm; per-bin calibration factor (1.10, 0.95, 0.88 here)
\(f, c\): bin frequency (3.85, 5.13, 6.42 MHz) and sound speed, 1.54 mm/µs
\(\zeta_{ej}\): two-way position (mm) at which the echo of predictor j appears on element e after dynamic receive focusing on the image line
\(u_e\): lateral position of element e relative to the image line, mm
\((x_j, z_j)\): scatterer position: lateral offset from the line and depth, mm
\(d_j\): extra path length (distance offset), mm: sets the arrival time; models multipath and the unknown transmit path
\(z_c\): depth of the window centre, mm
\(G_{ej}\): energy of the Gaussian pulse envelope inside the window
\(p^{\pm}_{ej}\): integration limits: \(\zeta_{ej} \pm h\) clipped to the window \([2z_{\mathrm{first}}, 2z_{\mathrm{last}}]\), with \(h\) = window length in depth units (0.30 mm); \(G_{ej} = 0\) if the pulse misses the window
\(\sigma_t\): rms duration of the pulse envelope, 0.125 µs
\(\zeta\): integration variable, mm

The first line is the regression; the next three generate one column: a scatterer at \((x_j, z_j)\), its wavefront delayed like the data (the left side of the third line is this page's receive delay law with a 0° transmit), its phase at the bin, and its amplitude from how much of its pulse lands in the window. The left side of the third line is the round-trip path the beamformer assigns to the sample at two-way position ζ on element e; setting it equal to the actual path of the echo gives the position ζ at which that echo appears in the delayed data.

\[ \begin{aligned} \hat{\mathbf{\beta}} &= \arg\min_{\mathbf{\beta}}\ \frac{1}{2N_o}\big\lVert \tilde{\mathbf{y}} - \tilde{\mathbf{X}}\mathbf{\beta} \big\rVert_2^{2} + \lambda\Big(\alpha\lVert\mathbf{\beta}\rVert_1 + \frac{1 - \alpha}{2}\lVert\mathbf{\beta}\rVert_2^{2}\Big) \\ \tilde{\mathbf{y}} &= \begin{bmatrix} \mathrm{Re}\,\mathbf{y} \\ \mathrm{Im}\,\mathbf{y} \end{bmatrix}, \quad \tilde{\mathbf{X}} = \begin{bmatrix} \mathrm{Re}\,\mathbf{X} & -\mathrm{Im}\,\mathbf{X} \\ \mathrm{Im}\,\mathbf{X} & \mathrm{Re}\,\mathbf{X} \end{bmatrix}, \quad \lambda = c_\lambda\, \mathrm{rms}(\tilde{\mathbf{y}}) \\ \hat{\beta}_j &\leftarrow \frac{S\big(\tilde{\mathbf{x}}_j^{T}\mathbf{r}_{(-j)}/N_o,\ \lambda\alpha\big)}{1/N_o + \lambda(1 - \alpha)}, \qquad \hat{\mathbf{y}}_{\mathrm{dec}} = \mathbf{X}_{\mathrm{ROI}}\,\hat{\mathbf{\beta}}_{\mathrm{ROI}} \end{aligned} \]

(25.14)

\(\hat{\mathbf{\beta}}\): fitted real coefficient vector (real and imaginary parts of the complex coefficients stacked)
\(\tilde{\mathbf{y}}, \tilde{\mathbf{X}}\): real-valued forms of \(\mathbf{y}\) and \(\mathbf{X}\) (each complex coefficient becomes two real unknowns)
\(N_o\): number of real observations, \(2N = 128\)
\(\lambda\): regularization strength, proportional to the rms of the data in this window and bin
\(\alpha\): L1/L2 mix, 0.9
\(c_\lambda\): regularization scale, 0.0189 (released default)
\(\lVert\cdot\rVert_1, \lVert\cdot\rVert_2\): sum of absolute values; Euclidean norm
\(\mathrm{rms}(\tilde{\mathbf{y}})\): root mean square of the 128 real observations
\(S(a, t)\): soft threshold: \(\operatorname{sign}(a)\max(\lvert a\rvert - t, 0)\)
\(\tilde{\mathbf{x}}_j\): column j of \(\tilde{\mathbf{X}}\) (unit norm)
\(\mathbf{r}_{(-j)}\): residual \(\tilde{\mathbf{y}} - \tilde{\mathbf{X}}\hat{\mathbf{\beta}}\) without the contribution of column j
\(\hat{\mathbf{y}}_{\mathrm{dec}}\): decluttered aperture vector at this bin: the ROI part of the fit
\(\arg\min\): value of the argument minimizing the cost
\(\mathrm{Re}, \mathrm{Im}\): real and imaginary parts

The third line is one step of cyclic coordinate descent, the solver of the released C code (y is first divided by its standard deviation). After the fit, \(\hat{\mathbf{y}}_{\mathrm{dec}}\) is placed in its bin, all other bins are set to zero, the window is inverse-transformed and summed over elements.

Why the regression must be regularized. With the released model space each fit has about 31 000 complex predictors (62 000 real unknowns) for 128 real numbers: hundreds of times underdetermined. The L1 term makes the answer sparse (here 30–200 nonzero coefficients per fit); the small L2 term spreads weight across near-identical neighbouring columns instead of picking one at random. Because λ scales with the rms of each window's data, the method is invariant to gain: what decides the outcome is how the energy divides between ROI and outer wavefronts, not the absolute level.

What it keeps and what it drops. A scatterer inside the resolution cell produces a nearly flat wavefront after focusing; ROI columns match it (model–data correlation ≥ 0.999 for isolated points in our tests) and it survives. An off-axis scatterer produces a tilted, partly truncated wavefront; outer columns take it and it is discarded. In an anechoic cyst the window contains only off-axis clutter and noise, so the fit often uses no ROI column at all and the pixel is exactly zero. Electronic noise matches neither family well: part of it ends up in the ROI part, part in the residual.

Item	Value used here	Source
Input	delayed, 9-angle compounded aperture data of Demo 25.6, 64 elements, complex samples on the image's 0.075 mm axial grid	adaptation: compound first, then ADMIRE on the synthetically focused data, as the authors did for plane waves (Schlunk, Dei & Byram, IEEE TUFFC 2021, in vivo; Schlunk et al., Ultrasonics 2023). The released code targets focused, walked apertures.
STFT window	rectangular, \(\lceil \mathrm{FWHM}_{\mathrm{pulse}}\cdot f_s\rceil\) = ⌈0.294 µs × 10.27 MHz⌉ = 4 samples (0.30 mm), zero-padded to 8, no overlap	released rule (rectangular, overlap 0)
Bins fitted	\(f_0 \pm 0.6\,\mathrm{BW}\,f_0\) = 3.2–6.8 MHz → 3.85, 5.13, 6.42 MHz	released rule; our data are analytic, so no mirrored negative band
ROI model	x ∈ ±0.5 res_lat in steps of 0.0179 res_lat; z ∈ z_c ± 0.5 res_ax in steps of 0.143 res_ax; d from −0.8 to 0.4 mm in steps of 0.0485 λ	released defaults (res_lat from the aperture transform, 0.38–0.78 mm here; res_ax = 2 res_lat)
Outer model	x within ±10.6 mm (half aperture + 1 mm) in steps of 1.42 res_lat; z = z_c − 0.5 res_ax and z_c + 0.21 res_ax; d from −8 to 3.2 mm in steps of 0.121 λ; outside the ROI ellipse enlarged by 0.5 mm	released defaults, except the depth range: the README's documented off-axis-only option instead of 0…1.05 z_c (this simulation has no reverberation, and the default range gives 85 000–495 000 predictors per fit, measured)
Wavenumber calibration κ	1.10, 0.95, 0.88 for the three bins	ours, fitted the way the release's per-probe factors are (best model–data correlation for isolated point scatterers); the release's own 5 MHz factors (1.29, 0.995, 0.85, …) follow the same trend
Element factor, calibration offsets	none (A_es = 1, offsets 0)	ours: the simulation has point elements and exact delays
Regression	α = 0.9, \(c_\lambda\) = 0.0189, coordinate descent from zero, tolerance 10⁻⁵, no ICA, all 64 elements (no aperture growth)	released defaults (tolerance: README's value without ICA); full aperture so that every panel of Demo 25.6 uses the same receive aperture
Predictors per fit	full ≈ 31 000 complex; reduced (live button) ≈ 900: ROI x, z, d sampled 8×, 4×, 8× more coarsely, outer d 8× more coarsely	measured; reduced model is ours

Methods card · ADMIRE (aperture domain model image reconstruction)α = 0.9 · λ = 0.0189·rms

What it does	Fits each window and frequency bin with modelled ROI and outer wavefronts (elastic net), keeps the ROI part, then DAS (25.14).
Parameters in Demo 25.6	α = 0.9; λ = 0.0189·rms(ỹ); tol 10⁻⁵; 4-sample (0.30 mm) windows; bins 3.85, 5.13, 6.42 MHz; κ = 1.10, 0.95, 0.88; ≈ 31 000 predictors.
Source of each value	Method: Byram et al. 2015 doi:10.1109/TUFFC.2015.007004. α, c_λ, tol, window, model space: authors' code (its default c_λ = 0.0189; Schlunk et al. 2021 doi:10.1109/TUFFC.2020.3012165 used 0.0189/2 and 0.0189/10 with α = 0.9). κ, outer depths: ours.
Why this value	Released defaults wherever they apply; κ fitted to this simulator as the release fits each probe.
Known weaknesses	Cost; dark band beside a bright scatterer; exact zeros in the cyst make CR arbitrary (67 % of cyst pixels here; ∞ when all are 0); mainlobe not narrowed; model mismatch loses signal.
Compute cost / pixel	≈ 10⁷–10⁸ (31 000 predictors × 64 channels per sweep, 3 bins). Ours: ≈ 0.05 s per fit, ≈ 0.04 s per pixel offline (Node.js; about 10 min per image); live reduced model 76 s per image (headless Chromium, our test).
Tuned value	tuned on the reduced model only: c_λ = 0.0189/4 (reduced-model gCNR 0.982 → 0.976 on the displayed data, 0.958 → 0.962 on the held-out data; the default point FWHM 0.5271 mm fails the 0.5267 mm limit) (protocol)

Demo 25.5b

ADMIRE at one window: data, fitted ROI and outer wavefronts, reconstruction

SimulationThe Demo 25.6 phantom at −20 dB channel noise, no phase screen (same simulation code and seeds), delayed and compounded on one image line; model (25.13) and elastic net (25.14) solved live in JS with the module's ADMIRE engine. "Full" is the released model space (≈ 31 000 predictors per fit; building the model and fitting takes 0.1–0.3 s), "reduced" the live-button model (≈ 900), "no outer model" keeps only the ROI columns. The depth profile runs the same fit on the 9 windows around the selected one (3 bins each; with the full model this takes several seconds).

data / DASROI part (kept) / ADMIREouter part (removed)residual

Window (image line, depth)

Frequency bin

Model size

Regularization c_λ 0.0189

computing…

Plots: top, the real part of the aperture vector at the chosen bin (rotated so that the data's mean phase is zero) with the fitted ROI and outer parts; middle, where the selected predictors sit (lateral offset from the image line, depth; marker area ∝ fitted energy summed over distance offsets; dashed box = ROI grid, ±0.5 res_lat × ±0.5 res_ax, drawn 2 px wider so that it stays visible); bottom, DAS and ADMIRE amplitude along the line over ±1.2 mm around the window (shaded), both relative to the DAS maximum in that range. Try: in speckle the fitted part is almost all ROI (outer part near zero, the rest residual); in the cyst centre the ROI part is zero or negligible and the pixel is black; at the point target the ROI part carries the echo. Remove the outer model and the cyst fills in (depth profile: most rows become nonzero, about 12 dB below DAS at the centre; a single window and bin can still come out empty): with nothing else to explain it, clutter is fitted by ROI columns and kept. Raise \(c_\lambda\) and fewer, smaller coefficients survive, so the ROI part shrinks toward zero; lower it and more of the clutter is fitted as well. Beside the bright scatterer, its off-axis energy dominates λ and the weak local speckle is shrunk away: the dark band visible in the ADMIRE panel of Demo 25.6.

Demo 25.5c

ADMIRE across conditions: 8 noise levels × 3 phase screens

SimulationFull model, precomputed offline in Node.js by running this page's own simulation, delay, ADMIRE and metric code with the same seeds as Demo 25.6 (shared-recording noise; channel noise −40…+30 dB in 10 dB steps × receive phase screen 0, 20, 40 ns rms); stored as 8-bit dB images (display 60 dB here) plus metrics computed at full precision. DAS values come from the same offline run and match the page's own DAS. Metrics use the ROIs and definitions of Demo 25.6; the pair dip is n/a at +30 dB / 0 ns, +20 dB / 40 ns and +30 dB / 40 ns (no resolvable pair).

ADMIRE (full model)DASfraction of ADMIRE cyst pixels exactly 0solid 0 ns, dashed 20 ns, dotted 40 nsunder each tile: gCNR ADMIRE / DAS

click a tile to open that setting in Demo 25.6

What it shows (full model, c_λ = 0.0189): with no phase screen ADMIRE holds gCNR 0.98–0.99 against 0.96 for DAS from −40 to −10 dB channel noise, keeps 0.94 against 0.82 at 0 dB, and both fall together once noise dominates (0.64 / 0.43 at +10 dB, 0.23 / 0.21 at +30 dB). The fraction of cyst pixels set exactly to zero falls from 67 % at −20 dB to 30 % at 0 dB and 21 % at +10 dB and beyond: noise matches neither family of wavefronts, so part of it is fitted by ROI columns. Speckle SNR stays near 1.82 down to −10 dB (DAS 1.86) and drops to 1.37 at +10 dB and 0.70 at +30 dB, because the same shrinkage that zeros the cyst also zeros weak speckle. A phase screen costs ADMIRE as much as DAS: at −20 dB, gCNR 0.984 → 0.967 → 0.885 for 0, 20, 40 ns (DAS 0.958 → 0.922 → 0.834); the model assumes exact delays, so an aberrated wavefront is split between ROI and outer columns, and the point widens (0.50 → 0.66 mm at 40 ns, DAS 0.48 → 0.61 mm). Click a tile to open that noise and screen setting in Demo 25.6 and compare the other nine methods on the same data.

25.6One dataset, ten beamformers and ADMIRE: images and metrics

Everything below runs in your browser on one simulated acquisition, so every method sees exactly the same channel data. Array and sequence: 64 elements, 0.3 mm pitch, 5 MHz, 60 % bandwidth (Gaussian pulse), IQ sampled at 20 MHz, nine plane waves from −4° to 4° in 1° steps (transmit grating lobe at \(\lambda/\Delta(\sin\theta) = 17.6\) mm, outside the phantom), coherent compounding per channel, rectangular receive apodization over all 64 elements. Phantom: fully developed speckle (80 scatterers/mm², z = 26–34.6 mm), an anechoic cyst of radius 2.5 mm at (0, 30) mm, a bright scatterer at (4.2, 30.5) mm, a point at (−3, 24) mm and a pair at x = 2.25 and 2.75 mm (0.5 mm apart) in an echo-free band. Noise: band-limited complex channel noise, independent per angle and element, at −20 dB re the mean speckle channel power by default (slider −40 to +30 dB). Optional aberration: a receive-path phase screen (Gaussian-correlated element delays, 3 mm correlation length). Grid: 81 × 161 pixels (0.15 × 0.075 mm); DMAS and SLSC use RF on a 4× finer axial grid.

Parameters: MV L = 32, K = ±2 samples, \(\Delta = 10^{-2}\); ESBMV δ = 0.5; GCF \(M_0 = 1\); PCF γ = 1; Wiener postfilter per sample; F-DMAS Gaussian \(2f_0\) band-pass; SLSC M = 8, kernel ±0.15 mm; FXPF p = 4, W = 16 samples, \(\mu = 10^{-3}\); ADMIRE: released model space, α = 0.9, \(c_\lambda = 0.0189\), 4-sample windows, 3 bins (25.5, table there). These are the defaults; the "Settings" selector switches to the values tuned by one common rule (25.7; tuned column of the table below).

Metrics (all computed on each method's own output, A ≥ 0): FWHM = lateral −6 dB width of the point at x = −3 mm; PSL = highest sidelobe beyond the first minima within ±2.5 mm; pair dip = depth of the valley between the 0.5 mm pair relative to the weaker peak (0 dB = not resolved; "n/a" = no resolvable pair: the weaker peak or the valley is exactly 0, or the two equal-amplitude peaks differ by more than 12 dB, which happens when noise dominates); \(\mathrm{CR} = 10\log_{10}(\langle A^2\rangle_{\mathrm{bg}}/\langle A^2\rangle_{\mathrm{cyst}})\) on the linear output; \(\mathrm{CNR} = \lvert\mu_{\mathrm{bg}} - \mu_{\mathrm{cyst}}\rvert/\sqrt{\sigma_{\mathrm{bg}}^2 + \sigma_{\mathrm{cyst}}^2}\) on the displayed dB image clipped at −DR; gCNR = 1 − OVL with 100 quantile bins of the pooled samples (exactly invariant to strictly monotone maps), on the linear output and on the clipped display; speckle SNR = μ/σ of A in the background ROI (1.91 for Rayleigh). Cyst ROI: disc of radius 2 mm; background ROI: x = −5.8…−3.4 mm, z = 28…32 mm. SLSC values are coherences, displayed linearly, and enter the metrics as \(A = \max(R_{\mathrm{sl}}, 0)\).

Methods and parameters at a glance

Method	Key parameters (values in Demo 25.6)	Source	Tuned value	Main weakness
DAS	uniform, 64 el., 9 angles	classical baseline	default already optimal	sidelobe clutter fills the cyst
MV	L = 32, K = ±2, Δ = 0.01	Synnevåg 2007/2009; values ours	L = 40, K = ±4, Δ = 10⁻⁴ (gCNR 0.978 → 0.982)	lone-point width measures SNR; low speckle SNR
ESBMV	δ = 0.5 (≤ 4 eigvecs) + MV	Asl & Mahloojifar 2010; δ ours	L = 48, δ = 0.3 (K, Δ as MV default) (gCNR 0.988 → 0.991)	breaks w^Ha = 1: band beside bright scatterer
CF	none	Mallart & Fink 1994	no CF weight (exponent 0, i.e. DAS) (gCNR 0.935 → 0.958)	dynamic-range stretch, dark holes
GCF	M₀ = 1	Li & Li 2003; M₀ ours	default kept (M₀ = 2 is +0.0006, within the 0.005 tie band)	arbitrary cut-off; no resolution gain
PCF	γ = 1, σ₀ = π/√3	Camacho 2009; γ ours	γ = 0.25 (gCNR 0.914 → 0.957)	γ-critical clipping to 0
SLSC	M = 8, kernel ±0.15 mm	Lediju 2011; M ours	M = 32, kernel ±0.30 mm (gCNR 0.991 → 0.963)	coherence, not amplitude; wide point
F-DMAS	Gaussian 2f₀ BP, σ = 1.80 MHz	Matrone 2015; filter ours	no free parameter	small gain on broadband data
FXPF	p = 4, W = 16 (1.2 mm), μ = 10⁻³	Shin 2016/2018; values ours	p = 3, W = 32 (2.4 mm) (gCNR 0.962 → 0.966)	keeps coherent off-axis clutter
Wiener	none (per sample)	Nilsen & Holm 2010	no free parameter	weak on noise (−6 dB), none on clutter
ADMIRE	α = 0.9, λ = 0.0189·rms, 4-sample windows, 3 bins	Byram 2015, authors' code; κ ours	c_λ = 0.0189/4 (gCNR 0.982 → 0.976, reduced model)	cost; zeros make CR meaningless

Each method name opens its methods card in Sections 25.1–25.5 (full citations with DOI, why each value, compute cost). "Ours" = our choice, with the reason in the card. Card timings: "ours" = headless Chromium in our test; "your run" fills in from the Demo 25.6 table.

Scope of this comparison

This is an educational comparison on one simulated phantom (a public dataset will follow). It is not a ranking and not clinical evidence: every number below depends on the parameters and on what one chooses to optimize. Each method runs with published or stated defaults (cards above), and is also shown tuned by one common rule: maximize gCNR subject to point FWHM ≤ 1.1 × DAS (the "Settings" selector in Demo 25.6; protocol and trade-off curves in 25.7). If you know a better setting for a method, please send it through the "Report an error" link at the top of the page.

Demo 25.6

Side-by-side comparison on one simulated dataset

SimulationPlane-wave IQ channel data of the phantom synthesized in JS from the delay law (Gaussian pulse, 9 angles, 64 elements, band-limited noise); all ten beamformers and the metrics run on it in the page (own code, not Field II). The image-SNR readout beamforms (DAS) the stored noise-only part and the noise-free part of the same data over the background ROI. The noise-texture selector and "New noise frame" change only the DAS panel. Its three modes: Plane-wave: the simulated data (frame 1) or, for later frames, a fresh channel-noise realization (same generator, filter and power, new seed) beamformed with DAS over the image and added to the stored noise-free DAS image. Line-by-line: emulated in the image domain, not by simulating focused transmits: each image column (one scan line) gets an independent complex Gaussian trace filtered along depth by the same Gaussian pulse (t = 2z/c), carried at f₀ like the echoes and scaled to the measured image noise power, added to the plane-wave noise-free DAS image; only the noise statistics of a line-by-line scan are modelled, not its transmit focusing. Per-pixel (non-physical): independent complex Gaussian noise per pixel at the measured image noise power. Line-by-line for all beamformers (checkbox, new run): each image line (column, and each lateral-profile column) gets its own independent channel-noise recording (same generator, power and pulse band, independent per element and angle, new seed per line) added to the noise-free channel data, and all ten beamformers, the metrics and the image-SNR readout are computed line by line from it; the transmit focusing of a real line-by-line scanner is not simulated (the signal is still plane-wave compounded), and the 25.7b/c sweeps keep the shared plane-wave noise.

ADMIRE panel: loading the precomputed full-model results…

x (mm) horizontal, z (mm) verticalcyst ROIbackground ROI

Display dynamic range 60 dB

Channel noise re speckle −20 dB

image SNR: run to measure

Noise texture in the DAS panel (DAS display only, instant)

noise frame: run first

Line-by-line noise for all beamformers (slower, ~13 s): reruns all ten beamformers

off: all panels use the shared plane-wave noise

Phase screen rms 0 ns

show ROIs

Settings

default settings

starts automatically when the module opens

Run time: about 9 s for the images and the fine lateral profiles in a desktop browser (the image-SNR readout adds well under 0.1 s); the sweeps in 25.7 then take about 15 s more in the background. The dynamic-range slider only redraws. Channel noise runs from −40 to +30 dB re the mean speckle channel power in 5 dB steps; noise and phase screen need a new run. The noise-texture selector and "New noise frame" recompute only the DAS panel (the table and the other panels keep the run's noise): line-by-line and per-pixel frames take a few milliseconds after a one-off 1–2 s for the noise-free DAS image; a new plane-wave frame regenerates all 576 channel-noise records and re-beamforms them over the image, about 2 s. With "Line-by-line noise for all beamformers" ticked, a run takes about 13 s instead of 9 s (in our test browser: about 3.7 s of per-line noise generation over the touched sample windows for 81 image lines and 416 profile lines, the rest beamforming as before); the DAS-only selector is then disabled because every panel already shows the line-by-line result, and the 25.7b/c sweeps still use the shared plane-wave noise. Settings: "Tuned" reruns every method with the setting chosen by the common rule of 25.7 (MV, ESBMV, CF, PCF, SLSC and FXPF change; DAS, GCF, Wiener and F-DMAS keep their defaults); the run takes longer (16.8 s against 14.9 s for the default run in our headless test, because the tuned MV, ESBMV and SLSC cost more per pixel), the 25.7b/c sweeps and the precomputed full-model ADMIRE panel stay at the defaults, and "Run reduced ADMIRE live" then uses the tuned c_λ. Measured in headless Chromium (our test machine), the live reduced ADMIRE run takes 76 s with the default c_λ and 116 s with the tuned c_λ (a smaller λ leaves more nonzero coefficients, so the solver needs more sweeps). In line-by-line mode the FWHM and pair-dip columns of the table measure noise, not resolution: every fine profile line has its own noise, so at high noise the profiles turn jagged (CF FWHM 0.28 → 0.07 mm at +10 dB). Read resolution from the shared-recording run.

Method	FWHM (mm)	PSL (dB)	pair dip (dB)	CR (dB)	CNR @ 60 dB	gCNR	gCNR @ 60 dB	speckle SNR	time (ms; per pixel)

What to look for (default settings; your browser reproduces these numbers exactly). Resolution: MV and ESBMV shrink the point to 0.03 mm and resolve the pair with a 12 dB and 9 dB dip, PCF and F-DMAS also open the pair (about −11 dB), DAS, GCF, Wiener and FXPF do not (−1.3 dB). In the MV and ESBMV images the right point of the pair looks dimmer: it sits at x = 2.75 mm, 0.05 mm from the nearest pixel column, and the MV peak is narrower than that. Contrast: CF, PCF and ESBMV report CR of 42–71 dB but their gCNR (0.91–0.99) is close to DAS's 0.96; their speckle SNR falls to 0.9–1.2 (dark holes in tissue). SLSC has the lowest CR (16 dB) and the highest gCNR (0.99) with a speckle SNR near 10, because it images coherence, not brightness. FXPF changes little here (CR +1.4 dB): at −20 dB noise the cyst is filled mostly by coherent off-axis clutter from the surrounding speckle, which prediction keeps. Raise the noise to −5 dB and FXPF lifts CR from 17.9 to 25.6 dB and gCNR from 0.91 to 0.96 while leaving the resolution and the speckle SNR (1.83) untouched. Cost: in our test browser the MV solve takes about 85 µs per pixel and ESBMV about 200 µs, against 0.3 µs for the DAS sum and about 60 µs for the delayed, compounded aperture data that every method shares (the table shows your machine's numbers).

Channel SNR, image SNR, and the texture of noise

The slider sets the noise per channel, relative to the mean speckle power in one element and one angle. The image sees something else. DAS sums 64 elements × 9 angles = 576 independent noise records, so the noise power grows by about 576 (27.6 dB; 568 after the linear interpolation of the 20 MHz IQ samples, which loses 1.3 %; 536 or 27.3 dB in our run, one realization of a ROI average whose spread is about ±12 % between noise draws). Speckle adds partly coherently across elements and angles: in our run its power grows by about 13 100 (41.2 dB; other scatterer realizations give 7 000–13 000). The image speckle-to-noise ratio is therefore \(\mathrm{SNR}_{\mathrm{img}} \approx \mathrm{SNR}_{\mathrm{ch}} + 10\log_{10}(G_{\mathrm{s}}/G_{\mathrm{n}})\) with \(G_{\mathrm{s}}/G_{\mathrm{n}} \approx 24\): the readout gives +13.9 dB at 0 dB channel SNR, +33.9 dB at the default −20 dB, and −16.1 dB at +30 dB (no phase screen; an aberrating screen lowers \(G_{\mathrm{s}}\)). That is why the 0 dB image still looks like tissue.

Why noise can look like speckle although it never touched the sound field. Speckle-like texture is not an acoustic phenomenon: it appears whenever a random complex field passes through a band-limited linear system, and the beamformer is one. Each beamformed pixel is a weighted sum of channel samples, so the noise covariance between two pixels is

\[ R_n(\mathbf{r}_1, \mathbf{r}_2) = \mathrm{E}\big[n(\mathbf{r}_1)\, n^{*}(\mathbf{r}_2)\big] = \sigma_n^2 \sum_{e} h_e(\mathbf{r}_1)\, h_e^{*}(\mathbf{r}_2) \]

\(n(\mathbf{r})\): beamformed (DAS) noise at pixel \(\mathbf{r}\), complex
\(\mathbf{r}_1, \mathbf{r}_2\): two pixel positions, mm
\(\mathrm{E}\): expectation over noise realizations
\(\sigma_n^2\): power of the white noise driving each channel
\(e\): index over all recorded channel samples (element, transmit, time)
\(h_e(\mathbf{r})\): weight with which channel sample e enters pixel \(\mathbf{r}\): receive filter, delay interpolation, apodization and phase rotation combined
\({}^{*}\): complex conjugate

Neighbouring pixels draw on nearly the same channel samples with nearly the same weights, so the noise stays correlated over about one resolution cell; its amplitude scales \(R_n\) but does not change its shape.

The result is a circular complex Gaussian field with grains of about the PSF size and a Rayleigh envelope. In our default run the noise-only DAS envelope in the background ROI has a one-pixel correlation of 0.71 laterally and 0.71 axially against 0.77 and 0.76 for the noise-free speckle, and μ/σ = 1.74 (1.86 for the speckle; 1.91 for an ideal Rayleigh variable).

This depends on the acquisition. It holds here because every pixel is beamformed from the same recordings (plane-wave or synthetic-aperture imaging). In conventional line-by-line focused imaging each image line comes from its own transmit and recording, so pixels on different lines share no channel samples and \(R_n\) between lines is zero; along a line the noise stays correlated over about a pulse length. The texture is then fine and nearly salt-and-pepper laterally, streaky along depth (the scan converter's interpolation between lines softens it to the line spacing). The selector shows both, plus the non-physical case. At +30 dB channel noise (image SNR −16.1 dB, background ROI essentially noise) the one-pixel envelope correlations lateral/axial are 0.73/0.70 for plane-wave noise, 0.04/0.71 for the line-by-line emulation and 0.00/−0.01 for per-pixel noise, all with μ/σ of 1.7–1.9. Per-pixel salt and pepper appears in real systems only when the noise enters after the band-limited imaging chain (detector or display electronics, quantization, image-domain processing), or when the pixel grid is much coarser than the PSF.

Line-by-line for all beamformers. Nonlinear beamformers change the texture further: CF and MV re-weight noise pixel by pixel. With the line-by-line option for all beamformers at +10 dB every panel turns laterally fine and streaky along depth, but the region metrics barely move for most methods (CR / gCNR shared → line-by-line: CF 14.8 → 13.5 dB / 0.40 → 0.43, PCF 10.0 → 9.3 dB / 0.37 → 0.39, SLSC 15.2 → 15.3 dB / 0.80 → 0.78, MV 4.8 → 4.2 dB / 0.35 → 0.34); ESBMV loses most (14.3 → 10.7 dB, 0.72 → 0.63) because its covariance and subspace now change from line to line, and the profile metrics become noise-limited (each fine profile line has its own noise, so FWHM and the pair dip read jagged noise: CF FWHM 0.28 → 0.07 mm, SLSC 1.31 → 0.03 mm; at −20 dB only the far sidelobes change, MV PSL −42 → −20 dB).

The test that always works: time. In every mode the noise is a new realization in each frame, while the speckle of stationary tissue is the same from frame to frame. Press "New noise frame" at +10 dB: the background grain changes while the point targets, the cyst and much of the speckle pattern stay put. On a scanner, hold the probe still on a phantom and compare consecutive frames, or average N frames: noise power falls by N, stationary speckle does not.

Pitfall: one dataset, one parameter set

The ranking above is for one noise level, one sequence and one parameter set per method. Change the noise to −5 dB and MV's gCNR drops below DAS (0.83 vs 0.91) while FXPF's rises above it; add a 30 ns screen and ESBMV's advantage shrinks (25.7). A beamformer paper that does not show the parameter sweep and the noise dependence has shown one point of a curve.

The dark streak beside the bright scatterer in ESBMV

The ESBMV image shows a black horizontal band at the depth of the bright scatterer (4.2, 30.5) mm, reaching from the cyst rim to the scatterer, with a thin bright tail around the scatterer itself. It is not a numerical artefact: recomputing those pixels with a full 32 × 32 eigendecomposition instead of the warm-started subspace iteration gives the same output to 0.1 dB. It is what (25.5) does when one bright off-axis scatterer dominates the local covariance. Within about 2.4 mm of a scatterer that stands some 30 dB above the surrounding speckle in the DAS image, its wavefront (a phase ramp across the subarray) is the largest eigenvector by a wide margin (\(\lambda_2/\lambda_1 \approx 0.1\) at z = 30.5 mm), so with δ = 0.5 the signal subspace \(\mathbf{E}_s\) holds that one vector only. That vector is nearly orthogonal to \(\mathbf{a} = \mathbf{1}\), and the projection drops the distortionless constraint: \(\lvert\mathbf{w}_{\mathrm{ESB}}^H\mathbf{a}\rvert \le 0.002\) (−54 dB) across the band. The local speckle, which arrives along \(\mathbf{a}\), is erased; what passes is the scatterer itself, scaled by MV's imperfect null on it, which forms the tail. The band is as tall as the pulse plus the ±K temporal average.

Numbers (default run, band x = 1.8–3.2 mm, z = 30.35–30.8 mm): ESBMV sits 51.8 dB below DAS there against 9.2 dB below DAS in the background speckle ROI; MV in the same band sits only 6.3 dB below DAS. Lowering δ shrinks the hole: −47.6, −35.7, −19.9 and −11.5 dB re DAS at δ = 0.2, 0.1, 0.05 and 0.01. The band touches the rim of the 2 mm cyst ROI, so it also adds a little to ESBMV's CR. Remedies that keep the eigen-idea: a lower or noise-referenced threshold (count eigenvalues above the estimated noise floor rather than above δλ₁), or a subspace that must retain the steering direction; the second is our suggestion and was not run here.

25.7Fair comparison: dynamic range, parameters and aberration

The dynamic-range trap (Rindal, Austeng, Fatemi & Rodriguez-Molares, IEEE TUFFC 2019). A displayed B-mode image is \(D = \max(20\log_{10}(A/A_{\max}), -\mathrm{DR})\). Any weight \(W \in [0, 1]\) adds \(20\log_{10}W\) pixel by pixel, and W is small exactly where the signal is incoherent, so the cyst sinks further than the tissue. A pure power law \(A \to A^p\) already does this: it multiplies every dB value by p, so a contrast measured as a difference of mean dB values scales exactly by p (the power-ratio CR used in our table roughly so: 24.7 → 46.4 dB for DAS² here) and the displayed CNR changes, while the pixel ordering, and therefore gCNR, is untouched.

\[ \begin{aligned} D_W(\mathbf{r}) &= D_{\mathrm{DAS}}(\mathbf{r}) + 20\log_{10}W(\mathbf{r}) \\ \mathrm{CR}_{\mathrm{dB}}(A^{p}) &= p\cdot\mathrm{CR}_{\mathrm{dB}}(A) \\ \mathrm{gCNR}\big(T(A)\big) &= \mathrm{gCNR}(A) \end{aligned} \]

(25.15)

\(D_W, D_{\mathrm{DAS}}\): displayed dB images (before clipping) with and without the weight W
\(A\): linear output amplitude
\(\mathrm{CR}_{\mathrm{dB}}\): difference of the mean dB values of tissue and cyst, before clipping
\(T\): any strictly increasing map (power law, log compression, S-curve)
\(\mathbf{r}\): pixel position
\(W(\mathbf{r})\): adaptive weight map (e.g. CF or PCF), dimensionless
\(p\): exponent applied to the weight
\(\mathrm{gCNR}\): generalized contrast-to-noise ratio, in [0, 1]
\(\log_{10}\): base-10 logarithm

Clipping at −DR is monotone but not strictly, so it can only lower gCNR (ties merge the two distributions at the floor). CF behaves roughly like DAS with a pixel-dependent power law: its CR in Demo 25.6 (52 dB) is close to that of DAS² (46 dB), but its gCNR is lower than DAS's because its noisy weights reorder tissue pixels.

Demo 25.7a

CR and CNR move with the display; gCNR does not

SimulationMetrics recomputed in JS from the Demo 25.6 images and ROIs for each display dynamic range.

Methods

Power-law stretch of DAS, p 2

Uses the images of Demo 25.6 and its ROIs. Try: compare CF with \(\mathrm{DAS}^p\) at p = 2: a similar linear CR (52 and 46 dB), and both need about 80 dB of display before it shows, yet DAS² keeps DAS's gCNR exactly (0.958) while CF's is lower (0.935). Displayed CR and CNR keep rising with DR until the cyst's own values are no longer clipped, then saturate; gCNR is flat until the floor starts to clip the tissue itself: below about 35 dB for DAS, but already below about 70 dB for CF, whose weights push many speckle pixels toward the floor.

Parameter sensitivity. Each method has one or two knobs, and the published images are usually at the best setting for the claim being made. The sweep below recomputes the metrics on the ROI and profile pixels for each value with everything else at default.

Demo 25.7b

Parameter sensitivity on the same data

SimulationEach beamformer re-run in JS on the Demo 25.6 channel data for every parameter value (ROI and profile pixels only).

Parameter

waiting for the comparison run

What the sweeps show: MV: FWHM falls from 0.48 mm (L = 1, DAS) to 0.03 mm (L = 32), CR peaks at L = 16 (31.7 dB) and speckle SNR drops above L = 16 (1.53 at 32). Without temporal averaging (K = 0) MV speckle collapses (speckle SNR 1.0, gCNR 0.86); K = 2 restores 1.53 and 0.98. Loading moves MV toward DAS: Δ = 10 gives FWHM 0.40 mm. ESBMV's δ changes CR from 55 to 72 dB at constant gCNR (0.99): pure stretch. GCF \(M_0 = 0\) is CF (gCNR 0.94); \(M_0\) = 1–4 gives gCNR 0.975 with speckle SNR 1.5–1.8. PCF with γ ≥ 1.5 zeroes most speckle (gCNR 0.38). SLSC: larger M narrows the point (5.3 → 0.63 mm for M = 1 → 24) and lowers speckle SNR (22 → 2.3). FXPF is insensitive to p and W at this noise level.

Equal tuning: one rule for every method

Defaults are open to the objection that a method was left untuned. So every method is also tuned by one rule, the same for all: maximize gCNR subject to point FWHM ≤ 1.1 × DAS FWHM, at −20 dB channel noise and no phase screen. gCNR is the objective because the display cannot move it (Demo 25.7a); the FWHM ceiling stops a method from buying contrast with resolution. The "Settings" selector of Demo 25.6 runs either set, and the methods cards and the summary table give the tuned values.

Tuning protocol

Objective: gCNR of the linear output, cyst ROI against background ROI (definitions of 25.6). Constraint: lateral −6 dB FWHM of the point at (−3, 24) mm ≤ 1.1 × the FWHM of default DAS (rectangular, 64 elements) on the same data, 1.1 × 0.4788 = 0.5267 mm. Settings whose half-maximum crossings fall outside the ±2.5 mm profile (SLSC with M = 1, three settings) have no valid FWHM and are excluded.

Data: a held-out realization of the same phantom, with new scatterer positions and amplitudes (seed 1007) and new channel noise (seed 1108); the displayed run uses seeds 7 and 108. ROI and profile pixels only, as in Demo 25.7b. The chosen setting is then evaluated once on the displayed data (full image, the metrics of the Demo 25.6 table). Everything ran offline in Node with this page's own engine code.

Grids: coarse, 293 settings in total, listed in the table below (MV L, K and Δ jointly; ESBMV δ with L; CF exponent including "none"; GCF M₀; PCF γ; SLSC M with the axial kernel; FXPF p with W; DAS rectangular or Hann apodization). Wiener and F-DMAS have no free parameter in our implementation; FXPF μ, the angles and the aperture stay fixed. Selection: the feasible setting with the highest gCNR (exact ties: the narrower point); the default is kept if it is feasible and within 0.005 of that maximum. For scale, the same DAS gives gCNR 0.928 on the held-out data and 0.958 on the displayed data: realization-to-realization differences are much larger than 0.005, so gains of a few thousandths below are within that spread.

ADMIRE: only c_λ ∈ {0.0189/4, 0.0189/2, 0.0189, 0.0189·2} is searched, with α fixed at the published 0.9, and only with the reduced model (≈ 900 predictors per fit). Tuning the full model (≈ 31 000 predictors per fit; its 24-setting grid in Demo 25.5c alone took 3.8 CPU-hours) over a parameter grid was too costly here, so the precomputed full-model panel keeps c_λ = 0.0189.

Other objectives would change the winner. Maximizing CR rewards dynamic-range stretch (ESBMV's default reaches 70.5 dB on the displayed data), minimizing FWHM rewards MV's SNR-dependent point width (25.2), and a different ceiling, noise level or phase screen moves every choice. The trade-off curves below show each method's whole range so that another rule can be applied by eye.

Demo 25.7d

Trade-off curves: point FWHM against gCNR for every setting tried

SimulationComputed offline in Node with this page's own engine (the M16E and ADM code of this page): 293 settings of 11 methods on the held-out realization (scatterer seed 1007, noise seed 1108), −20 dB channel noise, no phase screen, ROI and profile pixels as in Demo 25.7b; ADMIRE with the reduced model. The values are embedded in the page, not recomputed in your browser.

Methods

gCNR axis

What it shows: every setting tried, as point FWHM (log axis) against gCNR on the held-out data. Up and to the left is better on both. The line through each method's dots is its Pareto front (no other setting of that method is both narrower and higher), the dashed vertical line is the ceiling 1.1 × DAS FWHM = 0.527 mm, the open marker is the default and the filled marker the tuned setting (both coincide where the default was kept). Read off: DAS with Hann apodization gains gCNR (0.942) only by widening the point to 0.74 mm; CF's gCNR falls steadily as the exponent grows (0.928 without the weight, 0.892 at exponent 2, while the point narrows from 0.479 to 0.233 mm), so under this rule its best setting is no weight; PCF collapses above γ = 1 (0.599, 0.314, 0.129 at γ = 1.25, 1.5, 2); SLSC's best gCNR settings all have points wider than 0.6 mm, and only M = 32 with the ±0.30 mm kernel meets the ceiling. The ranking depends on the objective: without the ceiling SLSC and ESBMV tie at the top of the held-out data (0.991); with it ESBMV leads; a rule on CR, on the narrowest point or on a fixed display range would reorder the methods again.

Method	Grid searched (settings)	Default: FWHM (mm) / gCNR, held-out	Tuned setting	Tuned: FWHM (mm) / gCNR, held-out	gCNR on the displayed data, default → tuned	Ceiling binding?
DAS	apodization {rectangular, Hann} (2)	0.479 / 0.928	= default	—	0.958	yes
MV	L {8, 16, 24, 32, 40, 48} × K {0…4} × Δ {10⁻⁴ … 10} (decades) (180)	0.029 / 0.978	L = 40, K = ±4, Δ = 10⁻⁴	0.018 / 0.988	0.978 → 0.982	no
ESBMV	L {16, 24, 32, 40, 48} × δ {0.01, 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 0.9} (40)	0.029 / 0.982	L = 48, δ = 0.3 (K, Δ as MV default)	0.018 / 0.991	0.988 → 0.991	no
CF	exponent {0 (none), 0.25, 0.5, 1, 1.5, 2} (6)	0.296 / 0.897	no CF weight (exponent 0, i.e. DAS)	0.479 / 0.928	0.935 → 0.958	no
GCF	M₀ {0, 1, 2, 3, 4, 6, 8} (7)	0.433 / 0.960	= default	—	0.975	no
PCF	γ {0.25, 0.5, 0.75, 1, 1.25, 1.5, 2} (7)	0.302 / 0.884	γ = 0.25	0.434 / 0.929	0.914 → 0.957	no
Wiener	none (no parameter) (1)	0.468 / 0.906	= default	—	0.943	no
F-DMAS	none (band-pass fixed at √2 σ_f) (1)	0.427 / 0.950	= default	—	0.972	no
SLSC	M {1, 2, 4, 6, 8, 12, 16, 24, 32} × kernel ±{0.075, 0.15, 0.30} mm (27)	1.871 / 0.981 ✗	M = 32, kernel ±0.30 mm	0.474 / 0.964	0.991 → 0.963	yes
FXPF	p {1, 2, 3, 4, 6, 8} × W {8, 16, 32} samples (18)	0.478 / 0.938	p = 3, W = 32 (2.4 mm)	0.478 / 0.943	0.962 → 0.966	no
ADMIRE	c_λ {0.0189/4, 0.0189/2, 0.0189, 0.0189·2}, reduced model, α = 0.9 (4)	0.5271 / 0.958 ✗	c_λ = 0.0189/4	0.495 / 0.962	0.982 → 0.976	yes

"Binding" = the highest-gCNR setting of the grid fails the FWHM ceiling, so the constraint decided the choice. ✗ = the default fails the ceiling. ADMIRE rows use the reduced model on both data sets; the gCNR of the precomputed full model on the displayed data is in the Demo 25.6 table. ADMIRE's default misses the ceiling by 0.0004 mm (0.5271 against 0.5267 mm), a margin far below any resolution difference that matters, so its "tuned" value is a formality of the rule; the tuned c_λ is better on the held-out data and worse on the displayed data (0.982 → 0.976), and SLSC must give up gCNR (0.991 → 0.963) to meet the ceiling: the rule is a constraint, not a promise of improvement.

Aberration robustness. A near-field phase screen adds a random delay \(\delta_i\) to each element. For a point with independent Gaussian phase errors of rms \(\sigma_\varphi = 2\pi f_0\sigma_t\), the DAS power and the CF-weighted amplitude are

\[ \begin{aligned} \mathrm{E}\lvert\bar{y}\rvert^{2} &= \frac{1 + (N - 1)\,e^{-\sigma_\varphi^{2}}}{N} \\ \sigma_\varphi &= 2\pi f_0\sigma_t \\ \mathrm{CF}\cdot\lvert\bar{y}\rvert &= \lvert\bar{y}\rvert^{3} \end{aligned} \]

(25.16)

\(\sigma_t\): rms element delay error, ns
\(\sigma_\varphi\): rms phase error at \(f_0\), rad
\(\mathrm{E}\): expectation
\(\bar{y},\ \mathrm{CF},\ N\): as in (25.6) and (25.1)
\(f_0\): centre frequency, MHz

For equal-amplitude channels CF equals \(\lvert\bar{y}\rvert^2\) per realization, so the CF-weighted amplitude falls three times as fast in dB as DAS. Adaptive weights fitted to an aberrated covariance can also lose their resolution abruptly, and coherence images degrade gracefully because they normalize each channel pair.

Demo 25.7c

Metrics versus phase-screen strength

SimulationDemo 25.6 channel data re-beamformed in JS with Gaussian-correlated receive phase screens of increasing rms.

Metric

Screen correlation length

waiting for the comparison run

What it shows (3 mm correlation, 0 → 60 ns rms): DAS gCNR 0.96 → 0.72 and point brightness −11.6 dB. CF loses −31.8 dB of point brightness, close to three times DAS, and its gCNR falls to 0.50. GCF holds up better than CF (0.69). MV's 0.03 mm point widens to 1.1 mm by 50 ns. ESBMV shares the top with SLSC at 0–20 ns (gCNR ≥ 0.98) and is among the worst beyond 40 ns (0.60 at 60 ns). SLSC degrades most gracefully (gCNR 0.99 → 0.87, brightness −7.9 dB): normalized short-lag coherence is the least sensitive to a smooth screen. Shorter correlation lengths hurt more at the same rms.

Ultrasound link: targets B and C

The same curves are the basis of coherence-based sound-speed and aberration estimation: a region's mean CF, GCF or SLSC falls monotonically with the residual delay error, so maximizing it over an assumed speed or over a phase law corrects the focus (Module 07 for the VCZ upper bound, Module 14 for the optimization). Two practical consequences: estimate the focus quality with a method that degrades smoothly (SLSC or region-averaged CF), and apply the adaptive beamformer after the aberration correction, not before, because MV and ESBMV amplify a residual screen into resolution loss.

25.8Visual summary

Fourteen pictures that carry the module. Each card is a question to answer in five seconds, the answer with its number, and a figure computed from the formula (N = 64, 0.3 mm pitch, 5 MHz, z = 30 mm unless stated); the badge says whether it is a closed form, a simplified model or a simulation.

25.9How experts think about it

An expert reduces every method to what it does with one aperture vector: re-weight it (MV, ESBMV), score its constancy and scale the DAS output (CF family), combine it nonlinearly (DMAS), report its normalized correlation (SLSC), or clean it before summing (FXPF, ADMIRE). Each choice has a predictable signature in four numbers: the width of a lone point (often an SNR measurement), the valley between two points (actual resolution), gCNR (actual separability) and speckle SNR (what the weights did to tissue). Then they ask two questions of any claimed improvement: would a power law on DAS give the same CR, and does it survive a phase screen.

MV re-solves the weights at every pixel: unit gain toward the focus, minimum power from everything else.

\(\mathbf{w} = \mathbf{R}^{-1}\mathbf{a}/(\mathbf{a}^H\mathbf{R}^{-1}\mathbf{a})\)\(\mathbf{R} = \sigma^2\mathbf{I}\) → DASideal R, 2 mm, INR 10 dB: −121 vs −29 dB

Minimize \(\mathbf{w}^H\mathbf{R}\mathbf{w}\) subject to \(\mathbf{w}^H\mathbf{a} = 1\); the Lagrangian gives \(\mathbf{R}\mathbf{w} \propto \mathbf{a}\) (25.2). After the delays \(\mathbf{a} = \mathbf{1}\), so MV is DAS whose weights are recomputed from the local covariance. With the exact covariance of an incoherent interferer 2 mm off-axis at 10 dB, MV puts a −121 dB null on it while DAS leaves the −29 dB array-factor sidelobe there. Real covariance estimates from a few samples give far shallower nulls; the principle is the same.

Approx. modelNarrowband focal-plane array factor; MV weights from the expected covariance (source, interferer, white noise), matrix solve in JS.

Assumes the steering vector is exact (correct delays and sound speed) and enough independent snapshots; both fail in tissue without the fixes below.

Echoes are coherent copies of one pulse, so MV without subarray smoothing cancels the signal it is constrained to pass.

L = N: −66 dBL = 48: −0.3 dBneeds \((N - L + 1)\psi \gtrsim 2\pi\)

With a coherent interferer, \(\mathbf{R} \propto \mathbf{u}\mathbf{u}^H\) with \(\mathbf{u} = \mathbf{a} + \beta\mathbf{v}\). MV drives \(\mathbf{w}^H\mathbf{u}\) to zero, and \(\mathbf{w}^H\mathbf{a} = 1\) then forces \(\beta\,\mathbf{w}^H\mathbf{v} = -1\): the interferer is subtracted from the signal. Averaging the covariance over \(P = N - L + 1\) subarrays multiplies the cross term by the Dirichlet factor \(D_P(\psi)\) (25.4), which decorrelates them. Interferer at 1.5 mm, equal amplitude: −66 dB at L = 64, −0.3 dB at L = 48. A closer interferer has a smaller ψ and needs more subarrays.

Approx. modelNarrowband, expected forward-smoothed covariance of a coherent equal-amplitude pair, channel SNR 20 dB, loading 10⁻²; matrix solve in JS.

Forward smoothing assumes a shift-invariant (near-Toeplitz) covariance across the aperture; strong apodization, element dropout or near-field geometry break it.

The subarray length is the resolution knob: L = 1 is DAS, L = N/2 is the usual ceiling.

2 sources 0.5 mm apart: L = 16 unresolvedL = 32: −10.9 dB dipimage FWHM 0.48 → 0.03 mm

The adaptive aperture is L elements, so the Capon spectrum resolves like an L-element array: two incoherent sources 0.5 mm apart at 30 mm are unresolved for L = 16, dip −10.9 dB at L = 32 and −26.8 dB at L = 64 (exact covariance, 20 dB SNR). With L = 1 the covariance is a scalar and w = 1: DAS exactly. Above N/2 there are fewer subarrays to average, and the speckle suffers (simulated speckle SNR 1.86 at L = 16, 1.53 at 32, 1.43 at 48).

Approx. modelNarrowband Capon spectrum \(1/(\mathbf{v}^H\mathbf{R}^{-1}\mathbf{v})\) with the exact covariance of two incoherent sources, 20 dB SNR, loading 10⁻²; matrix solve in JS.

Two equal incoherent sources in white noise; coherent pairs, unequal amplitudes and speckle backgrounds change the numbers.

Diagonal loading slides MV toward DAS; quote Δ as a fraction of the mean element power.

\(\varepsilon = \Delta\cdot\mathrm{tr}(\mathbf{R})/L\)\(\Delta = 10^{-2}\): 1 %Δ = 10: FWHM 0.40 mm

Loading adds \(\varepsilon\mathbf{I}\), a pretend white-noise floor, so MV behaves as if the SNR were at most about 1/Δ. For a lone point the −6 dB width scales as \(\sqrt{\sigma^2 + \varepsilon}/\sqrt{N}\), so loading sets the width once ε exceeds the true noise. As Δ → ∞ the weights become 1/L and the smoothed output is DAS with a trapezoidal taper (box L ∗ box N − L + 1): white-noise gain 16.9 dB instead of 18.1 dB for L = 32. In the image, FWHM goes 0.021 → 0.029 → 0.40 mm for \(\Delta = 10^{-4}, 10^{-2}, 10\).

Approx. modelLone-point width law \(0.914\sqrt{(\sigma^2 + \varepsilon)/N}\) relative to DAS, narrowband, L = N; capped at the DAS width.

Some papers define \(\varepsilon = \Delta\cdot\mathrm{tr}(\mathbf{R})\) without the 1/L; compare the loading-to-power ratio, not the symbol.

A lone point under MV is as narrow as the SNR allows: its width measures SNR, not resolution.

width ratio \(0.914/\sqrt{N\cdot\mathrm{SNR}}\)30 dB: 0.0036 × DASsimulated 0.03 vs 0.48 mm

For one source in white noise, Woodbury gives \(\lvert y(x)\rvert = \sigma^2\lvert\mathbf{v}^H\mathbf{a}\rvert/Q\), \(Q = N(\sigma^2 + N) - \lvert\mathbf{v}^H\mathbf{a}\rvert^2\). Near the peak \(\lvert\mathbf{v}^H\mathbf{a}\rvert/N \approx 1 - \psi^2(N^2 - 1)/24\), so the output halves when \(2N(1 - \lvert\mathbf{v}^H\mathbf{a}\rvert/N) = \sigma^2\), giving a −6 dB width of \(0.914/\sqrt{N\cdot\mathrm{SNR}}\) times the DAS width (checked against a direct solve). Hence the striking MV point images. Two points separated by less than the DAS width are the honest test: the pair dip in Demo 25.6 is −12 dB.

Approx. modelClosed-form MV output for one source in white noise, narrowband, L = N, no loading; plotted as is at three SNRs.

Narrowband, exact covariance, L = N; subarrays, temporal averaging and bandwidth widen it, loading caps it.

MV peaks are narrower than your pixels: sample the lateral grid finely or bright points flicker.

half a 0.15 mm pixel: MV −23 dBDAS −0.5 dBCF −1.4 dB

Delay-and-sum images are band-limited, so λ/2 lateral sampling loses little. MV and ESBMV outputs are not: in our dataset a point 0.075 mm from the nearest pixel column reads −23 dB (MV) and −25 dB (ESBMV) below its peak (the narrowband L = 32 model in the figure gives −30 dB), and the pair point at 2.75 mm, 0.05 mm off the grid, looks dimmer than its twin in Demo 25.6. Interpolating the MV image afterwards cannot recover it; evaluate the beamformer on a finer grid or at the scatterer.

Approx. modelNarrowband lone point, L = 32 forward-smoothed covariance, loading 10⁻², matrix solve in JS; dots at the 0.15 mm pixel centres.

Matters for points and ULM-like sparse targets; speckle regions are dominated by the transmit PSF and loading, and change less.

Temporal averaging over about a pulse length is what keeps MV speckle looking like speckle.

K = 0: speckle SNR 1.02, gCNR 0.86K = ±2: 1.53, 0.98DAS: 1.86, 0.96

With K = 0 each pixel's covariance comes from one time sample: N − L + 1 = 33 subarrays of the same snapshot, a rank-deficient estimate that changes completely from sample to sample, so the weights and the output brightness fluctuate on top of the speckle. Averaging ±K samples (here ±0.15 mm, about one pulse length) stabilizes the estimate. Synnevåg, Austeng & Holm (IEEE TUFFC 2009) made this point on speckle statistics; the sweep in Demo 25.7b reproduces it.

Longer averaging blurs axially and mixes depths with different interference; K beyond a pulse length buys little (speckle SNR 1.71 at K = 4).

ESBMV keeps only the eigenvectors that carry signal: δ near 1 keeps one, δ → 0 returns MV.

\(\mathbf{w}_{\mathrm{ESB}} = \mathbf{E}_s\mathbf{E}_s^H\mathbf{w}\)speckle keeps 1.3 of 32 at δ = 0.5δ 0.1 → 0.9: CR 55 → 72 dB, gCNR 0.99

A point's subarray covariance is rank one plus noise; speckle's eigenvalues decay smoothly (the VCZ coherence is a Toeplitz triangle); noise is flat. Projecting the MV weights on the eigenvectors with \(\lambda \ge \delta\lambda_1\) removes weight components that only fit noise (25.5). In the sweep, δ moves CR by 17 dB while gCNR stays at 0.99: the threshold mostly stretches the dynamic range. Next to a dominant bright scatterer the kept vector is that scatterer's wavefront and the local tissue is erased (the dark streak in Demo 25.6, 52 dB below DAS).

Approx. modelEigenvalues of L = 32 Toeplitz covariances: rank-one point, VCZ triangle (\(m_{\mathrm{tx}}\) = 15.7) for speckle, identity for noise, each plus 1 % white noise.

Costs an eigendecomposition per pixel (or an iterative approximation); the threshold has no universal value, breaks the distortionless constraint next to bright off-axis scatterers, and degrades fastest under aberration (gCNR 0.60 at 60 ns).

CF is the Wiener postfilter with the noise overestimated by a factor N.

\(\mathrm{CF} = \lvert\bar{y}\rvert^2/(\lvert\bar{y}\rvert^2 + \hat{\sigma}^2)\)\(H = \lvert\bar{y}\rvert^2/(\lvert\bar{y}\rvert^2 + \hat{\sigma}^2/N)\)noise: CF −36 dB, H −6 dB

Since \((1/N)\sum\lvert s_i\rvert^2 = \lvert\bar{y}\rvert^2 + \hat{\sigma}^2\), CF is a signal over signal-plus-channel-noise ratio. The DAS output averages N channels, so its noise power is \(\hat{\sigma}^2/N\); the Wiener gain with that noise is \(H = N\cdot\mathrm{CF}/((N - 1)\mathrm{CF} + 1)\) (25.6) (Nilsen & Holm, IEEE TUFFC 2010). For a point at 0 dB channel SNR, CF = 0.51 and H = 0.98; in pure noise CF = 1/64 and H = 0.50.

Approx. modelExpected CF of a point in white noise taken as a ratio of means, \((N\cdot\mathrm{SNR} + 1)/(N(\mathrm{SNR} + 1))\); H from the exact identity.

Both use one-sample estimates unless you average over an axial kernel; the per-pixel values are noisy and darken speckle (simulated speckle SNR 0.91 for CF, 1.67 for Wiener).

GCF counts the low aperture spatial frequencies: \(M_0 = 0\) is CF, each extra bin admits waves from 0.5 mm further off-axis.

bin k ↔ \(\Delta x = k\lambda z/(Np)\)2 mm → k = 4.2\(M_0 = 1\) keeps 60 % of speckle

By Parseval, \(\lvert S(0)\rvert^2/\sum\lvert S(k)\rvert^2\) is CF (25.7). A wave from lateral offset Δx is a phase ramp of \(Np\,\Delta x/(\lambda z)\) cycles across the aperture, 0.48 mm per bin at 30 mm. Speckle spreads its aperture energy over \(\lvert k\rvert \lesssim N/m_{\mathrm{tx}}\), so CF (\(M_0 = 0\)) keeps only 22 % of it and GCF with \(M_0 = 1\) keeps 60 %, while off-axis clutter from beyond about \(M_0 + 1\) bins is still rejected. In the image: CF speckle SNR 0.91, GCF(1) 1.50, gCNR 0.94 vs 0.97.

Approx. modelExpected speckle aperture spectrum from the VCZ triangle (\(m_{\mathrm{tx}}\) = 15.7); off-axis wave as a Dirichlet kernel at k = 4.2.

Assumes a linear-phase picture of off-axis energy (focal zone); aberration spreads every bin and erodes the separation.

PCF and SCF ignore amplitude: a phase spread of \(\pi/\sqrt{3}\) means noise, zero means a focused echo.

\(\sigma_0 = \pi/\sqrt{3} = 1.81\) radSCF: 1 bit per channelγ ≥ 1.5: gCNR 0.38

A phase uniform on (−π, π] has standard deviation \(2\pi/\sqrt{12} = \pi/\sqrt{3}\), so \(\mathrm{PCF} = 1 - \gamma\sigma(\varphi)/\sigma_0\) reaches zero for noise at γ = 1 (25.8) (Camacho, Parrilla & Fritsch, IEEE TUFFC 2009). Speckle phases spread too, which is why PCF images have dark speckle (speckle SNR 0.99 here) and why γ is critical: at γ = 1.5 most tissue pixels are clipped to zero and gCNR collapses to 0.38 while CR, with the zeroed cyst pixels floored at \(10^{-12}\) of the maximum, reads over 200 dB.

Exact formula

Phase is only defined where there is signal; wrapping needs the auxiliary phase; strongly amplitude-modulated channels (apodization, element directivity) are treated as if equal.

For a point with phase errors, the CF-weighted brightness falls exactly three times as fast in dB as DAS.

\(\mathrm{CF}\cdot\lvert\bar{y}\rvert = \lvert\bar{y}\rvert^3\)iid: DAS \(-4.34\,\sigma_\varphi^2\) dB60 ns screen: −11.6 vs −31.8 dB

For unit-amplitude channels with phase errors, \(\mathrm{CF} = \lvert\sum e^{\,j\varphi}\rvert^2/N^2 = \lvert\bar{y}\rvert^2\), so the weighted amplitude is \(\lvert\bar{y}\rvert^3\): three times the DAS loss in dB, per realization. With independent Gaussian errors \(\mathrm{E}\lvert\bar{y}\rvert^2 = (1 + (N - 1)e^{-\sigma^2})/N \approx e^{-\sigma^2}\) (25.16), i.e. \(-4.34\,\sigma^2\) dB for DAS and \(-13\,\sigma^2\) dB for CF. Our correlated 60 ns screen gives −11.6 dB (DAS) and −31.8 dB (CF).

Exact formula

Exact for a single point and phase-only errors; amplitude errors and speckle backgrounds add terms.

DMAS is a \(2f_0\) beamformer on the sum coarray: lower sidelobes and a slightly narrower main lobe, at O(N) cost.

\(\sum_{i<j}\hat{s}_i\hat{s}_j = \tfrac12[(\sum\hat{s})^2 - \sum\hat{s}^2]\)narrowband: −26.5 dB, 0.74 × widthsimulated: −27 vs −19 dB

The \(2f_0\) term of pair (i, j) has phase \(\varphi_i + \varphi_j\): a virtual element at \((x_i + x_j)/2\) receiving at λ/2 with triangular multiplicity. A triangle at λ/2 has −26.5 dB sidelobes and 0.74 times the −6 dB width of a uniform aperture at λ; equivalently, the narrowband DMAS beam is the DAS beam squared (25.10), so for a lone point the gain equals a p = 2 power law. On broadband plane-wave data the measured gains are smaller (FWHM 0.43 vs 0.48 mm) because the transmit PSF is common and the signed square root mixes harmonics. The identity turns 2016 pair products into 64 square roots.

Exact formula

Narrowband picture only to first order; DMAS output is nonlinear in the scatterer amplitudes, so speckle statistics change (speckle SNR 1.24).

The DMAS sum has a DC term; without the \(2f_0\) band-pass you image the DC.

aligned echo: each pair = \(\lvert\mathrm{rf}(t)\rvert\)DC + \(2f_0\) + \(4f_0\) …RF every ≈ λ/16: 41 MHz

For an aligned echo every pair term is \(\operatorname{sign}(\mathrm{rf}^2)\sqrt{\mathrm{rf}^2} = \lvert\mathrm{rf}(t)\rvert\), a full-wave rectified pulse: a large DC (envelope) component, a band at \(2f_0\) (energy 0.22 of the DC band) and weaker even harmonics. The low-frequency part carries the difference phases \(\varphi_i - \varphi_j\), whose pair sum \((\lvert\sum e^{\,j\varphi}\rvert^2 - N)/2\) is the DAS intensity pattern: DAS on a squared scale, nothing the \(2f_0\) part does not already give. The "F" in F-DMAS is the band-pass that keeps \(2f_0\). The products need RF sampled well above \(4f_0\) so that nothing folds into the pass band: we form them on an 18.75 µm (≈ λ/16) depth grid, 41 MHz equivalent.

SimulationFFT in JS of the full-wave rectified Gaussian pulse (5 MHz, 60 % band); dashed: the Gaussian band-pass at \(2f_0\).

The filter shape changes the axial response and the noise level of F-DMAS; papers rarely state it.

SLSC brightness in speckle is set by the transmit aperture: coherence hits zero at the lag equal to the (synthesized) transmit width.

\(m_{\mathrm{tx}} = N_\theta z\,\Delta(\sin\theta)/p = 15.7\)\(R_{\mathrm{sl}}(8)\) = 5.7 VCZ, 6.0 simulatedpoint ≈ M, noise ≈ 0

VCZ: the receive coherence of speckle is the autocorrelation of the transmit aperture function. Nine plane waves at 1° synthesize a 4.7 mm aperture at 30 mm, so \(R(m)\) falls to zero near lag 15.7 and \(R_{\mathrm{sl}} = \sum_{m \le M}(1 - m/m_{\mathrm{tx}}) = 5.7\) for M = 8 (25.11). A full focused transmit (\(m_{\mathrm{tx}} = 64\)) would give 7.4. Choose M inside the coherent range: larger M narrows points (5.3 → 0.63 mm for M = 1 → 24) and lowers speckle SNR (22 → 2.3).

Approx. modelVCZ triangle \(R(m) = 1 - m/m_{\mathrm{tx}}\) for three transmit widths; fully developed speckle in the focal zone, no noise.

Fully developed speckle in the focal zone; noise scales every lag ≥ 1 by SNR/(1 + SNR) (Module 07); coarse angle steps add transmit-grating-lobe oscillations to \(R(m)\).

FXPF keeps whatever is predictable across the aperture: noise loses about \(p/(N - p)\) of its power, weak echoes are shrunk too.

p = 4: noise −14 dBecho at 0 dB SNR: −2 dBoff-axis waves: kept

q plane waves across the aperture obey an order-q linear recursion, so the least-squares prediction (25.12) reproduces them; white noise is reproduced only through its projection on p regressors, about \(p/(N - p)\) of its power (−11.8 dB for p = 4), and a further 2 dB after forward-backward averaging (−14.1 dB in Monte Carlo). Noisy regressors bias the fitted coefficients toward zero, so a 0 dB echo loses about 2 dB (p = 4) or 6 dB (p = 1). In the image at −5 dB noise: CR 17.9 → 25.6 dB, gCNR 0.91 → 0.96.

SimulationMonte Carlo in JS: 64-channel echo plus white noise, forward-backward AR(p) LS prediction, DAS-output gain on echo (solid) and noise (dashed).

Coherent off-axis clutter is predictable and survives; reverberation that looks like a few plane waves survives too; windows shorter than a few pulse lengths make the per-frequency fits noisy.

ADMIRE explains the aperture data with modelled wavefronts and keeps only those from inside the resolution cell: off-axis clutter is not attenuated, it is assigned elsewhere and deleted.

α = 0.9, \(c_\lambda\) = 0.0189≈ 31 000 predictors for 128 real data per fitDemo 25.6, −20 dB: gCNR 0.984 vs DAS 0.958; 67 % of cyst pixels exactly 0

Per window and bin, (25.13) stacks the aperture signatures of in-ROI and out-of-ROI scatterers and (25.14) picks a sparse combination; the image is DAS of \(\mathbf{X}_{\mathrm{ROI}}\hat{\mathbf{\beta}}_{\mathrm{ROI}}\). A point with five +9.5 dB off-axis points around it keeps its amplitude within 0.7 dB, while the clutter DAS leaves at the point (−33.6 dB) becomes exactly zero (validation run, full model). In Demo 25.6 (full model, −20 dB, no screen) it gives gCNR 0.984 against 0.958 for DAS with speckle SNR 1.82 (DAS 1.86) and a point FWHM of 0.496 mm against 0.479 mm: contrast without a resolution gain.

The model must match the physics: our per-bin wavenumber calibration (as in the released code), the 0° receive delay law and the absence of an element factor are assumptions; energy 0.3–1 mm off the line has no outer predictor and stays; near a bright scatterer λ (∝ window rms) shrinks weak speckle to zero; cost here ≈ 0.04 s per pixel (0.05 s per fit) in Node.js, about 10 min per image.

Every [0, 1] weight stretches the dynamic range: raise the image to the power p and its dB contrast scales by p; gCNR does not move.

\(\mathrm{CR}_{\mathrm{dB}}(A^p) = p\cdot\mathrm{CR}_{\mathrm{dB}}\)\(\mathrm{gCNR}(T(A)) = \mathrm{gCNR}(A)\)CF: CR 52 dB, gCNR 0.94 < DAS 0.96

A power law multiplies every dB value by p, so the cyst-to-tissue difference of mean dB values scales exactly by p (a power-ratio CR roughly: 24.7 → 46.4 dB for DAS²) while the ordering of pixels, and hence the overlap of the two histograms, is unchanged (25.15). CF and ESBMV act like data-dependent power laws (Rindal et al., IEEE TUFFC 2019): CF doubles DAS's CR (24.7 → 52.3 dB) yet has a lower gCNR, because its noisy weights also scramble the tissue ordering. Report gCNR, and compare CR only at matched dynamic range.

Exact formula

gCNR estimates need a fixed bin rule and enough independent speckle cells; clipping at the display floor can only lower it.

Compare at matched data, noise and dynamic range, show the parameter sweep, and quote the cost per pixel.

MV ≈ 16 000 MACs ≈ 250 × the DAS sumESBMV ≈ 2–3 × MVCF family, DMAS ≈ 1–10 ×

Per pixel, DAS costs N complex multiply-adds; the CF family 2–4N; DMAS O(N) with the square-sum identity; SLSC N·M·kernel; FXPF a few thousand; MV \((2K + 1)(L(N - L + 1) + L^2)\) \(+\, L^3/6 \approx 16\,000\) for L = 32, K = 2 when the subarray-averaged covariance is built by running sums along its diagonals, as this page's engine does (the direct sum of P outer products costs \((2K + 1)(N - L + 1)L^2 \approx 170\,000\)); ESBMV adds the eigen-analysis. In this page's JavaScript the MV solve takes about 85 µs per pixel and ESBMV about 200 µs, against 0.3 µs for the DAS sum and about 60 µs for the delayed aperture data all methods share. The PICMUS challenge (Liebgott et al., IEEE IUS 2016) standardized data and metrics for exactly this reason.

Approx. modelOperation counts per pixel (complex MACs) from the formulas in the text; memory traffic ignored.

Operation counts ignore memory traffic, which dominates on GPUs; FXPF and SLSC costs scale with window and kernel sizes.

25.10Cheat sheet

Concept	Formula	Ultrasound use
Aperture vector	\(\mathbf{s}(\mathbf{r}) \in \mathbb{C}^N\) after delays; focus ⇒ \(\mathbf{a} = \mathbf{1}\); \(y_{\mathrm{DAS}} = \frac{1}{N}\sum_i s_i\)	common input of every method here
Off-axis wave	\(\psi = 2\pi p\,\Delta x/(\lambda z)\) per element; bin \(k = Np\,\Delta x/(\lambda z)\) (0.48 mm per bin at 30 mm)	read clutter direction off \(S(k)\)
MV (Capon)	\(\mathbf{w} = \frac{\mathbf{R}^{-1}\mathbf{a}}{\mathbf{a}^H\mathbf{R}^{-1}\mathbf{a}},\ \ P = \frac{1}{\mathbf{a}^H\mathbf{R}^{-1}\mathbf{a}}\)	adaptive apodization per pixel
MV estimator	L-element subarrays (\(L \le N/2\)), ±K samples (≈ 1 pulse), \(\varepsilon = \Delta\cdot\mathrm{tr}(\hat{\mathbf{R}})/L\)	Synnevåg-style implementation; L = 32, K = 2, \(\Delta = 10^{-2}\)
Smoothing factor	\(D_P(\psi) = \frac{1}{P}\sum_l e^{\,j\psi l},\ \ P = N - L + 1\)	decorrelates coherent echoes; needs \(P\psi \gtrsim 2\pi\)
MV lone-point width	\(0.914/\sqrt{N\cdot\mathrm{SNR}}\) × DAS width	reads SNR or loading, not resolution
ESBMV	\(\mathbf{w}_{\mathrm{ESB}} = \mathbf{E}_s\mathbf{E}_s^H\mathbf{w},\ \ \lambda_i \ge \delta\lambda_1\)	stronger contrast, mostly dynamic-range stretch; dark streaks next to bright scatterers
CF	\(\frac{\lvert\sum s\rvert^2}{N\sum\lvert s\rvert^2} = \frac{\lvert\bar{y}\rvert^2}{\lvert\bar{y}\rvert^2 + \hat{\sigma}^2}\): point 1, noise 1/N, speckle \(\sum(N - \lvert m\rvert)C(m)/N^2\)	focus metric; weight on DAS
Wiener postfilter	\(H = \frac{\lvert\bar{y}\rvert^2}{\lvert\bar{y}\rvert^2 + \hat{\sigma}^2/N} = \frac{N\cdot\mathrm{CF}}{(N - 1)\mathrm{CF} + 1}\)	noise suppression without CF's over-darkening
GCF	\(\sum_{\lvert k\rvert \le M_0}\lvert S(k)\rvert^2/\sum\lvert S(k)\rvert^2\), \(M_0 = 0\) ⇒ CF	tolerant of speckle and mild aberration
PCF / SCF	\(\max(0, 1 - \gamma\sigma(\varphi)/\sigma_0),\ \sigma_0 = \pi/\sqrt{3}\); \(\lvert 1 - \sqrt{1 - \bar{b}^2}\rvert^p\)	amplitude-free coherence, 1-bit hardware
DMAS	\(\sum_{i<j}\hat{s}_i\hat{s}_j = \frac12[(\sum\hat{s})^2 - \sum\hat{s}^2],\ \hat{s} = \operatorname{sign}(s)\sqrt{\lvert s\rvert}\), band-pass at \(2f_0\)	sum-coarray beam (≈ DAS beam squared) at O(N)
SLSC	\(R(m)\) = normalized lag-m correlation; \(R_{\mathrm{sl}} = \sum_{m=1}^{M}R(m)\), M ≈ 5–30 % of N	clutter-robust coherence image
VCZ speckle coherence	\(R(m) \approx 1 - m/m_{\mathrm{tx}}\); plane waves \(m_{\mathrm{tx}} = N_\theta z\,\Delta(\sin\theta)/p\)	expected SLSC and CF of speckle
FXPF	\(X_i(f) \approx \sum_{k=1}^{p}a_k X_{i-k}(f)\), LS fit, forward + backward	removes spatially incoherent noise before DAS
ADMIRE	\(\hat{\mathbf{\beta}} = \arg\min \frac{1}{2N_o}\lVert\tilde{\mathbf{y}} - \tilde{\mathbf{X}}\mathbf{\beta}\rVert^2 + \lambda\big(\alpha\lVert\mathbf{\beta}\rVert_1 + \frac{1-\alpha}{2}\lVert\mathbf{\beta}\rVert^2\big)\), \(\lambda = c_\lambda\,\mathrm{rms}(\tilde{\mathbf{y}})\); keep \(\mathbf{X}_{\mathrm{ROI}}\hat{\mathbf{\beta}}_{\mathrm{ROI}}\)	removes coherent off-axis (and multipath) clutter before DAS; α = 0.9, \(c_\lambda\) = 0.0189
Aberration loss	\(\mathrm{E}\lvert\bar{y}\rvert^2 = \frac{1 + (N - 1)e^{-\sigma_\varphi^2}}{N}\); \(\mathrm{CF}\cdot\lvert\bar{y}\rvert = \lvert\bar{y}\rvert^3\) for a point	robustness ranking, autofocus metric
Dynamic-range trap	\(D_W = D_{\mathrm{DAS}} + 20\log_{10}W\); \(\mathrm{CR}_{\mathrm{dB}}(A^p) = p\cdot\mathrm{CR}_{\mathrm{dB}}\); \(\mathrm{gCNR}(T(A)) = \mathrm{gCNR}(A)\)	report gCNR; compare CR at matched DR

25.11Interview and qualifying-exam questions

I1Derive the MV weights. When is MV identical to DAS? · screen

Minimize \(\mathbf{w}^H\mathbf{R}\mathbf{w}\) subject to \(\mathbf{w}^H\mathbf{a} = 1\). The Lagrangian gradient gives \(\mathbf{R}\mathbf{w} \propto \mathbf{a}\), and the constraint fixes the scale: \(\mathbf{w} = \mathbf{R}^{-1}\mathbf{a}/(\mathbf{a}^H\mathbf{R}^{-1}\mathbf{a})\), output power \(1/(\mathbf{a}^H\mathbf{R}^{-1}\mathbf{a})\) (25.2). After the delays \(\mathbf{a} = \mathbf{1}\), so with \(\mathbf{R} = \sigma^2\mathbf{I}\) (white noise) the weights are 1/N: DAS. The same holds for L = 1 (scalar covariance) and for Δ → ∞ (then the smoothed MV is DAS with a trapezoidal taper). Follow-up: "What does the steering vector assume?" Correct delays: a sound-speed error makes \(\mathbf{a}\) wrong and MV suppresses the signal itself.

I2Why does textbook MVDR fail on ultrasound echoes, and what are the three standard fixes? · onsite

(1) The interfering echoes are coherent copies of the same transmit pulse, so MV uses them to cancel the desired signal: with an equal-amplitude echo 1.5 mm away the gain is −66 dB at L = N. Forward subarray smoothing (\(L \le N/2\)) decorrelates them through the factor \(D_P(\psi)\); at L = 48 the gain is back to −0.3 dB. (2) One time sample is one snapshot: averaging ±K samples over about a pulse length stabilizes \(\hat{\mathbf{R}}\) (speckle SNR 1.02 → 1.53 in the lab). (3) Model errors and noise: diagonal loading \(\varepsilon = \Delta\cdot\mathrm{tr}(\hat{\mathbf{R}})/L\) with \(\Delta \approx 10^{-2}\) to \(10^{-1}\) (25.3). Follow-up: "What is the cost of smoothing?" The adaptive aperture shrinks to L elements, so resolution and null depth drop.

I3A paper shows an MV point 0.03 mm wide against 0.48 mm for DAS. Is the resolution 16 times better? · onsite

No. For a lone point in white noise the MV width is \(0.914/\sqrt{N\cdot\mathrm{SNR}}\) times the DAS width: it measures SNR (capped by loading near 1/Δ), not the ability to separate two scatterers. The honest test is two points closer than the DAS width: in the lab, a 0.5 mm pair shows a −12 dB dip with MV against −1.3 dB for DAS, which is a real gain, but nowhere near 16×. Also ask about lateral sampling: an MV peak narrower than the pixel pitch reads −23 dB half a pixel off the grid. Follow-up: "What would you ask the authors to show?" Pair dip versus separation, loading and subarray sweeps, and gCNR.

I4What CF do you expect for a point, for speckle and for noise? Why is our speckle CF 0.23 and not 2/3? · screen

Point at focus: 1. Noise: 1/N (0.0156 for 64 elements). Speckle: the ratio of means is \(\sum(N - \lvert m\rvert)C(m)/N^2\) with C the VCZ coherence. A fully focused transmit with the full receive aperture gives \(C(m) = 1 - \lvert m\rvert/N\) and CF ≈ 2/3. Nine plane waves at 1° synthesize only a 4.7 mm transmit aperture at 30 mm, so C falls to zero at \(m_{\mathrm{tx}} = 15.7\) lags and the expected CF is 0.225 (simulated 0.23). Follow-up: "So what does CF do to the speckle brightness?" It lowers speckle by about 13 dB and noise by 36 dB here, a 23 dB stretch.

I5How are CF and the Wiener postfilter related? · onsite

Write \((1/N)\sum\lvert s_i\rvert^2 = \lvert\bar{y}\rvert^2 + \hat{\sigma}^2\). Then \(\mathrm{CF} = \lvert\bar{y}\rvert^2/(\lvert\bar{y}\rvert^2 + \hat{\sigma}^2)\): signal over signal plus channel noise. The MMSE gain for the DAS output needs the output noise \(\hat{\sigma}^2/N\), so \(H = N\cdot\mathrm{CF}/((N - 1)\mathrm{CF} + 1)\) (25.6) (Nilsen and Holm, 2010). For noise CF = 1/N (−36 dB) but H ≈ 0.5 (−6 dB). Follow-up: "Is H optimal per pixel?" Only with good estimates; single-sample estimates are noisy, so average over an axial kernel or subarrays.

I6Explain GCF's \(M_0\) physically. How would you choose it? · onsite

\(S(k)\) is the DFT of the aperture vector; bin k is a wave arriving from lateral offset \(\Delta x = k\lambda z/(Np)\), 0.48 mm per bin at 30 mm. GCF keeps \(\lvert k\rvert \le M_0\), so \(M_0\) sets how far from the focus incoming energy still counts as signal. Speckle occupies \(\lvert k\rvert \lesssim N/m_{\mathrm{tx}} \approx 4\) here: \(M_0 = 0\) (CF) keeps 22 % of it, \(M_0 = 1\) keeps 60 %. Choose \(M_0\) small enough to exclude the clutter directions you care about and large enough to pass speckle and residual aberration; the lab finds gCNR flat at 0.975 for \(M_0\) = 1–4. Follow-up: "What does aberration do to \(S(k)\)?" It spreads every bin, so GCF degrades more gently than CF but still degrades.

I7DMAS costs \(N(N - 1)/2\) products per sample. Is that true, and where does its gain come from? · onsite

No: \(\sum_{i<j}\hat{s}_i\hat{s}_j = \frac12[(\sum\hat{s})^2 - \sum\hat{s}^2]\) (25.9), so it is O(N) (64 square roots instead of 2016 products). The \(2f_0\) part of each pair carries \(\varphi_i + \varphi_j\), a virtual element at the mid-point receiving at λ/2: the sum coarray, a triangle. In narrowband terms the beam is the DAS beam squared: −26.5 dB sidelobes and 0.74× width, the same as a p = 2 power law for an isolated point. The signed square root keeps amplitude linearity, which is what separates it from DAS². Follow-up: "Why the band-pass?" The pair sum also has a DC (difference-coarray) part that only reproduces the DAS intensity; it must be removed, and the products need RF sampled well above \(4f_0\).

I8Why is a point target wide in SLSC, and what sets the speckle brightness? · research

SLSC normalizes each channel pair, so it reports coherence, not amplitude. Wherever a point dominates the local noise, even far down its DAS sidelobes, the aperture signal is a coherent wavefront and \(R(m) \approx 1\): the point stays bright over about 1.9 mm here (M = 8). Larger M includes lags where the off-focus wavefront decorrelates, so points narrow (5.3 mm at M = 1 to 0.63 mm at M = 24). In speckle \(R(m)\) follows the VCZ curve, the autocorrelation of the transmit aperture, so \(R_{\mathrm{sl}} \approx \sum(1 - m/m_{\mathrm{tx}}) = 5.7\) for M = 8, and noise multiplies every lag ≥ 1 by SNR/(1 + SNR). Follow-up: "Why does SLSC beat DAS in clutter?" Reverberation and off-axis clutter in a dark region are incoherent there and contribute R ≈ 0 regardless of brightness.

I9A reviewer sees CF with CR 52 dB against DAS 25 dB. What do you say? · onsite

CR on the linear output is not a fair comparison for any [0, 1] weight. The weight adds \(20\log_{10}W\) per pixel and is smallest where the signal is incoherent, which stretches the dynamic range: DAS raised to the power p multiplies a dB-domain contrast by p (the power-ratio CR goes from 24.7 to 46.4 dB at p = 2) and leaves gCNR unchanged (25.15). On the same data CF's gCNR is 0.94 against DAS's 0.96, and its speckle SNR falls to 0.91. Report gCNR (Rindal et al., 2019), and CR or CNR only at a matched displayed dynamic range. Follow-up: "Is there any real CF benefit?" Suppression of incoherent noise and clutter where it is the dominant contrast limit, which gCNR will then show.

I10FXPF and ADMIRE both filter channel data before DAS. What can each remove? · research

FXPF fits a linear predictor across the aperture at each temporal frequency. A focused echo and any off-axis plane wave are predictable, so both survive; what goes is spatially incoherent energy (noise, diffuse reverberation): about −14 dB on white noise at p = 4. Its coefficients shrink at low SNR, so weak echoes lose too (−2 dB at 0 dB channel SNR). ADMIRE fits a dictionary of modelled scatterer signatures, inside and outside the region of interest, with elastic-net regularization and reconstructs only the in-region part, so it can remove coherent off-axis and multipath clutter, at a much higher cost. Follow-up: "Why put them before aberration correction?" They clean the data that the coherence or delay estimators use, without changing the delays.

I11You want to estimate sound speed or a phase screen from coherence. Which coherence measure, and where does MV go in the pipeline? · research

Use a region-averaged measure that falls smoothly and monotonically with the residual delay error: mean CF or GCF over a speckle region, or SLSC (here SLSC's gCNR fell only from 0.99 to 0.87 at 60 ns rms). For a point with phase errors the CF-weighted brightness falls three times as fast in dB as DAS (\(\lvert\bar{y}\rvert^3\)), which makes CF sensitive but also noisy as a per-pixel weight. Maximize over the assumed speed or the phase law, as in Module 14. Apply MV or ESBMV after the correction: they amplify residual delay errors into resolution loss (MV's 0.03 mm point widens to 1.1 mm at 50 ns). Follow-up: "Why not use MV's output power as the focus metric?" It depends on loading and SNR as much as on focus.

I12An ESBMV image shows a black horizontal band next to a bright point target that DAS and MV do not show. Explain and fix. · research

Near the bright scatterer its wavefront dominates the subarray covariance, so with a relative threshold \(\lambda_i \ge \delta\lambda_1\) the signal subspace holds only that eigenvector. It is nearly orthogonal to \(\mathbf{a}\), and the projection \(\mathbf{E}_s\mathbf{E}_s^H\mathbf{w}_{\mathrm{MV}}\) drops the unit-gain constraint (\(\lvert\mathbf{w}_{\mathrm{ESB}}^H\mathbf{a}\rvert \le 0.002\) in our band): the local tissue vanishes over the pulse length plus the temporal average, 52 dB below DAS in Demo 25.6. Fixes: a lower δ (−36 dB at δ = 0.1, −11.5 dB at 0.01, at the price of less contrast stretch), or a threshold referenced to the noise floor instead of \(\lambda_1\). Follow-up: "Does it affect the metrics?" Only where the band meets a ROI; here it touches the cyst rim and inflates CR slightly.

I13A colleague adds noise at 0 dB SNR to every channel of a 64-element, 9-angle plane-wave dataset and says the noise "did not work" because the DAS image still looks clean. What do you tell them? · onsite

Per-channel SNR is not image SNR. Independent noise adds incoherently, a power gain of \(64 \times 9 = 576\) (27.6 dB), while speckle adds partly coherently; in Demo 25.6 speckle gains about 41 dB, so 0 dB per channel becomes +12 to +14 dB speckle-to-noise in the image, depending on the speckle realization (+13.9 dB in the displayed run). In this plane-wave dataset what remains also looks like speckle: every pixel is a weighted sum of the same band-limited channel samples, so neighbouring pixels share most of their noise (\(R_n(\mathbf{r}_1,\mathbf{r}_2) \propto \sum_e h_e(\mathbf{r}_1)h_e^{*}(\mathbf{r}_2)\)) and the noise becomes a complex Gaussian field with grains of about the PSF size and a Rayleigh envelope. To check, beamform the noise-only data and compare ROI powers, look in an echo-free region, or compare two frames: noise changes from frame to frame, the speckle of stationary tissue does not. Follow-up: "Thermal noise never interacts with the sound field, so how can it look like speckle?" Speckle-like texture comes from any random complex field passing through a band-limited linear system, not from acoustics; the beamformer is that system. "Would a focused line-by-line scan look the same?" No: each line has its own transmit and recording, so noise is independent from line to line and correlated only over about a pulse length along depth, giving a fine, streaky texture that is nearly salt and pepper laterally. True per-pixel salt and pepper needs noise entering after the band-limited chain (post-detection electronics, quantization, image-domain processing) or a pixel grid that undersamples the PSF.

I14A paper reports that ADMIRE improves cyst contrast by 20 dB over DAS. What do you check before you believe it, and what does ADMIRE actually decide per pixel? · research

Per depth window and frequency bin it fits the delayed aperture vector with modelled wavefronts of scatterers inside a small ROI around the line and outside it, by an elastic net (α = 0.9, \(\lambda = 0.0189\,\mathrm{rms}(\mathbf{y})\) in the released code), and keeps only the ROI part (25.14). It is underdetermined by a factor of hundreds (about 31 000 complex predictors for 64 complex data with the released model space), so sparsity decides. Checks: (1) CR on the linear output is meaningless if the cyst is set to exactly zero; in our full-model run 67 % of the cyst pixels are exactly zero at −20 dB channel noise (CR 35.7 dB against 24.7 dB for DAS, set by the remaining third), so ask for gCNR and speckle SNR. (2) Is the comparison at matched dynamic range? (3) Was the model calibrated to the system (the released code carries per-probe wavenumber factors)? (4) What happens to speckle near bright scatterers (λ scales with the window's rms, so weak tissue next to a bright target is shrunk away; this motivated iADMIRE, Schlunk et al. 2021)? (5) Cost: about 0.05 s per fit (10 min per 81 × 161 image) in our Node.js run, real time only with GPU and ICA-compressed models (Khan et al. 2021). Follow-up: "Why can it remove off-axis clutter when FXPF cannot?" FXPF keeps anything predictable, and an off-axis plane wave is; ADMIRE has an explicit outer model that claims that wavefront, so it is dropped.

25.12MATLAB lab

Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8 and simulates its own data. The scripts use MATLAB's random generator, so their phantoms differ from the browser's and the numbers differ by a few dB in CR; the rankings agree. Paste into the editor and press Run.

Download the scripts: m25_compare_all.m · m25_mv_esbmv.m · m25_coherence_factors.m · m25_dmas.m · m25_slsc.m · m25_fxpf.m · m25_admire.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m25_compare_all.mAll ten beamformers on one plane-wave dataset: image grid and metric table with gCNR at 60 and 40 dB · change noise_dB, screen_ns and every method parameter · runtime ≈ 23 s

%% m25_compare_all.m  -  ten beamformers on one simulated plane-wave dataset
% Purpose : simulate 64-channel IQ data for 9 compounded plane waves (point targets,
%           a 0.5 mm pair, an anechoic cyst in speckle, a bright scatterer, channel
%           noise, optional receive phase screen), beamform with DAS, MV, ESBMV, CF,
%           GCF, PCF, Wiener, F-DMAS, SLSC and FXPF, and print a metric table:
%           FWHM, PSL, pair dip, CR, CNR and gCNR at two display ranges, speckle SNR, time.
% Change  : noise_dB (channel noise re speckle), screen_ns (rms phase screen),
%           L, Delta, K (MV), delta (ESBMV), M0 (GCF), M (SLSC), p, W (FXPF).
% Expect  : MV/ESBMV FWHM ~0.03 mm and a deep pair dip; CF/PCF/ESBMV CR far above
%           DAS but gCNR close to DAS; SLSC lowest CR, highest gCNR; runtime ~20 s.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);

%% Acquisition and phantom
c = 1540; f0 = 5e6; fs = 20e6; lam = c/f0;            % m/s, Hz, Hz, m
N = 64; pitch = 0.3e-3; xe = ((0:N-1) - (N-1)/2)*pitch; % elements (m)
bw = 0.6; sigf = bw*f0/(2*sqrt(2*log(2))); sigt = 1/(2*pi*sigf);   % Gaussian pulse
ang = (-4:4)*pi/180; Na = numel(ang);                 % plane waves, 1 deg steps
t0 = 25e-6; Nt = 580; tt = t0 + (0:Nt-1)'/fs;         % IQ time axis (s)
noise_dB = -20; screen_ns = 0;                        % channel noise, phase screen
% speckle (80 /mm^2) with an anechoic cyst, three points, one bright scatterer
area = [-6.8 6.8 26 34.6]*1e-3; dens = 80e6;          % m, scatterers per m^2
nsc = round((area(2)-area(1))*(area(4)-area(3))*dens);
xs = area(1) + (area(2)-area(1))*rand(nsc,1); zs = area(3) + (area(4)-area(3))*rand(nsc,1);
keep = (xs.^2 + (zs-30e-3).^2) > (2.5e-3)^2; xs = xs(keep); zs = zs(keep); as = randn(numel(xs),1);
xs = [xs; -3e-3; 2.25e-3; 2.75e-3; 4.2e-3]; zs = [zs; 24e-3; 24e-3; 24e-3; 30.5e-3]; as = [as; 15; 15; 15; 60];

%% Channel data: iq(t) = a g(t - tau) exp(-j 2 pi f0 tau), accumulated sample by sample
tic; IQ = zeros(Nt, N, Na); K4 = ceil(4*sigt*fs); E = repmat(1:N, numel(xs), 1);
for a = 1:Na
  tau = (xs*sin(ang(a)) + zs*cos(ang(a)))/c + sqrt((xs - xe).^2 + zs.^2)/c;
  ph = as.*exp(-1j*2*pi*f0*tau); n0 = (tau - t0)*fs; accr = zeros(Nt*N,1); acci = accr;
  for k = -K4:K4
    n = floor(n0) + k; g = exp(-((n - n0)/fs).^2/(2*sigt^2)); ok = n >= 0 & n < Nt;
    idx = n(ok) + 1 + (E(ok)-1)*Nt; v = g(ok).*ph(ok);
    accr = accr + accumarray(idx, real(v), [Nt*N 1]); acci = acci + accumarray(idx, imag(v), [Nt*N 1]);
  end
  IQ(:,:,a) = reshape(accr + 1j*acci, Nt, N);
end
band = tt > 2*27.5e-3/c & tt < 2*33e-3/c; Pref = mean(mean(abs(IQ(band,:,5)).^2));
taps = exp(-((-K4:K4)'/fs).^2/(2*sigt^2)); sn = sqrt(Pref*10^(noise_dB/10)/sum(taps.^2));
for a = 1:Na
  w = (randn(Nt+2*K4, N) + 1j*randn(Nt+2*K4, N))/sqrt(2);
  IQ(:,:,a) = IQ(:,:,a) + sn*conv2(w, taps, 'valid');
end
scr = zeros(1,N);                                      % receive phase screen (s)
if screen_ns > 0
  wv = randn(1, N+40); ker = exp(-(-20:20).^2/(2*(1.5e-3/pitch)^2)); scr = conv(wv, ker, 'valid');
  scr = scr - mean(scr); scr = scr/sqrt(mean(scr.^2))*screen_ns*1e-9;
end
fprintf('simulated %d scatterers, 9 angles x 64 channels in %.1f s\n', numel(xs), toc);

%% Parameters
L = 32; Delta = 1e-2; Kt = 2; delta = 0.5; M0 = 1; gam = 1; M = 8; Ks = 8; p = 4; W = 16; mu = 1e-3;
sub = 4; dz = 0.075e-3; dzf = dz/sub; fsz = c/(2*dzf);  % output and RF axial steps, RF rate (Hz)
names = {'DAS','MV','ESBMV','CF','GCF','PCF','Wiener','F-DMAS','SLSC','FXPF'}; nm = numel(names);
T = zeros(1, nm);                                     % accumulated runtime per method (s)

%% Two passes: pass 1 = image grid, pass 2 = lateral profile through the point at (-3, 24) mm
for pass = 1:2
  if pass == 1
    xg = (-6:0.15:6)*1e-3; z0 = 22e-3; nzo = 161; rows = 1:nzo;
  else
    xg = [(-5.48:0.04:-3.64) (-3.6:0.02:-3.16) (-3.15:0.004:-2.85) (-2.84:0.02:-2.4) (-2.36:0.04:-0.52) (1.5:0.02:3.5)]*1e-3;
    z0 = 22.8e-3; nzo = 33; rows = 17;
  end
  zf = z0 + (0:(nzo-1)*sub)'*dzf; nzf = numel(zf); nx = numel(xg);
  % aperture data: S(row, col, element), delayed, phase-corrected, compounded (linear interpolation of IQ)
  [Zp, Xp] = ndgrid(zf, xg); Xp = Xp(:); Zp = Zp(:); np_ = numel(Xp); S = zeros(np_, N);
  Ttx = (Xp*sin(ang) + Zp*cos(ang))/c; Ptx = exp(1j*2*pi*f0*Ttx);        % np x Na
  for e = 1:N
    trx = sqrt((Xp - xe(e)).^2 + Zp.^2)/c; acc = zeros(np_, 1);
    for a = 1:Na
      fn = (Ttx(:,a) + trx - scr(e) - t0)*fs; n = floor(fn); fr = fn - n; ok = n >= 0 & n < Nt-1;
      col = IQ(:,e,a); v = zeros(np_, 1); v(ok) = col(n(ok)+1).*(1 - fr(ok)) + col(n(ok)+2).*fr(ok);
      acc = acc + v.*Ptx(:,a);
    end
    S(:,e) = acc.*exp(1j*2*pi*f0*(trx - scr(e)));
  end
  S = reshape(S, nzf, nx, N); So = S(1:sub:end,:,:);  % output-row samples
  out = zeros(nzo, nx, nm);
  % --- DAS and coherence weights (per sample)
  tic; ys = sum(So,3); out(:,:,1) = abs(ys)/N; T(1) = T(1) + toc;
  tic; P = sum(abs(So).^2,3); CF = abs(ys).^2./(N*P + eps); out(:,:,4) = CF.*abs(ys)/N; T(4) = T(4) + toc;
  tic; Sk = fft(So, [], 3); low = abs(Sk(:,:,1)).^2;
  for k = 1:M0, low = low + abs(Sk(:,:,1+k)).^2 + abs(Sk(:,:,N+1-k)).^2; end
  out(:,:,5) = low./(N*P + eps).*abs(ys)/N; T(5) = T(5) + toc;
  tic; phs = angle(So); pha = phs - pi*sign(phs); sg = min(std(phs,1,3), std(pha,1,3));
  out(:,:,6) = max(0, 1 - gam*sg/(pi/sqrt(3))).*abs(ys)/N; T(6) = T(6) + toc;
  tic; out(:,:,7) = (N*CF./((N-1)*CF + 1)).*abs(ys)/N; T(7) = T(7) + toc;
  % --- F-DMAS on RF = real part, fine axial grid, Gaussian band-pass at 2 f0
  tic; R = real(S); sh = sign(R).*sqrt(abs(R)); yd = (sum(sh,3).^2 - sum(abs(R),3))/2;
  nf = 2^nextpow2(nzf + 64); f = [0:nf/2, -nf/2+1:-1]'*fsz/nf; H = 2*exp(-(f - 2*f0).^2/(2*(sqrt(2)*sigf)^2)).*(f > 0);
  env = abs(ifft(fft(yd, nf, 1).*H, [], 1)); out(:,:,8) = env(1:sub:nzf, :)/(N*(N-1)/2); T(8) = T(8) + toc;
  % --- SLSC on RF, axial kernel +-Ks fine samples (kernel sums by cumulative sums along depth)
  tic; ro = (rows-1)*sub + 1; hi = min(nzf, ro + Ks); lo = max(1, ro - Ks);
  Cs = cumsum(cat(1, zeros(1, nx, N), R.^2), 1); En = Cs(hi+1,:,:) - Cs(lo,:,:); sl = zeros(numel(rows), nx);
  for m = 1:M
    Cs = cumsum(cat(1, zeros(1, nx, N-m), R(:,:,1:N-m).*R(:,:,1+m:N)), 1); num = Cs(hi+1,:,:) - Cs(lo,:,:);
    sl = sl + mean(num./sqrt(En(:,:,1:N-m).*En(:,:,1+m:N) + eps), 3);
  end
  out(rows, :, 9) = max(sl, 0); T(9) = T(9) + toc;
  % --- MV and ESBMV: subarrays L, temporal averaging +-Kt output samples, loading Delta
  Pn = N - L + 1; idx = (0:L-1)' + (1:Pn); tmv = 0; tes = 0;
  for ix = 1:nx
    for r = rows
      t1 = tic; X = [];
      for k = -Kt:Kt
        if r+k >= 1 && r+k <= nzo, x = squeeze(So(r+k, ix, :)); X = [X, x(idx)]; end %#ok<AGROW>
      end
      Rm = (X*X')/size(X,2); Rl = Rm + Delta*real(trace(Rm))/L*eye(L);
      w = Rl\ones(L,1); w = w/(ones(1,L)*w); x = squeeze(So(r, ix, :)); b = mean(x(idx), 2);
      out(r, ix, 2) = abs(w'*b); tmv = tmv + toc(t1);
      t2 = tic; [V, D] = eig((Rm + Rm')/2); lam_ = real(diag(D)); [lam_, o] = sort(lam_, 'descend'); V = V(:,o);
      Es = V(:, lam_ >= delta*lam_(1)); we = Es*(Es'*w); out(r, ix, 3) = abs(we'*b); tes = tes + toc(t2);
    end
  end
  T(2) = T(2) + tmv; T(3) = T(3) + tmv + tes;
  % --- FXPF: Hann windows (W samples, hop W/2); per bin, forward and backward AR(p) least squares across the
  %     aperture, solved for all columns at once (vectorized elimination); predictions averaged; then DAS
  tic; win = 0.5 - 0.5*cos(2*pi*(0:W-1)'/W); hop = W/2; Y = zeros(nzo, nx, N);
  for s0 = -hop:hop:nzo-1
    rr = s0 + (1:W); ok = rr >= 1 & rr <= nzo; seg = zeros(W, nx, N); seg(ok,:,:) = So(rr(ok),:,:);
    F = fft(seg.*win, [], 1); G = zeros(W, nx, N);
    for kb = 1:W
      Xb = reshape(F(kb,:,:), nx, N); pr = zeros(nx, N); cnt = zeros(1, N);
      for dirn = [1 -1]
        if dirn > 0, tgt = p+1:N; else, tgt = 1:N-p; end
        A = zeros(nx, numel(tgt), p); for q = 1:p, A(:,:,q) = Xb(:, tgt - dirn*q); end
        Gm = zeros(nx, p, p); h = zeros(nx, p);
        for q1 = 1:p, h(:,q1) = sum(conj(A(:,:,q1)).*Xb(:,tgt), 2); for q2 = 1:p, Gm(:,q1,q2) = sum(conj(A(:,:,q1)).*A(:,:,q2), 2); end, end
        tr = zeros(nx,1); for q = 1:p, tr = tr + real(Gm(:,q,q)); end
        for q = 1:p, Gm(:,q,q) = Gm(:,q,q) + mu*tr/p + 1e-30; end
        for k = 1:p                                   % Gaussian elimination, vectorized over columns
          for i = k+1:p
            fct = Gm(:,i,k)./Gm(:,k,k); for j = k:p, Gm(:,i,j) = Gm(:,i,j) - fct.*Gm(:,k,j); end, h(:,i) = h(:,i) - fct.*h(:,k);
          end
        end
        cf_ = zeros(nx, p);
        for k = p:-1:1, sres = h(:,k); for j = k+1:p, sres = sres - Gm(:,k,j).*cf_(:,j); end, cf_(:,k) = sres./Gm(:,k,k); end
        P_ = zeros(nx, numel(tgt)); for q = 1:p, P_ = P_ + cf_(:,q).*A(:,:,q); end
        pr(:, tgt) = pr(:, tgt) + P_; cnt(tgt) = cnt(tgt) + 1;
      end
      G(kb,:,:) = reshape(pr./cnt, [1 nx N]);
    end
    g = ifft(G, [], 1); Y(rr(ok),:,:) = Y(rr(ok),:,:) + g(ok,:,:);
  end
  out(:,:,10) = abs(sum(Y,3))/N; T(10) = T(10) + toc;
  if pass == 1, img = out; xi = xg; zi = z0 + (0:nzo-1)*dz; else, prof = squeeze(out(17,:,:)); xp = xg; end
  fprintf('pass %d: %d x %d pixels, cumulative times (s): %s\n', pass, nzo, nx, sprintf('%.1f ', T));
end

%% Metrics
[ZI, XI] = ndgrid(zi, xi); cyst = (XI.^2 + (ZI-30e-3).^2) < (2e-3)^2;
bgr = XI >= -5.8e-3 & XI <= -3.4e-3 & ZI >= 28e-3 & ZI <= 32e-3;
fprintf('\n%-7s %6s %6s %6s %6s %7s %7s %6s %7s %7s %6s %7s\n', 'method', 'FWHM', 'PSL', 'dip', 'CR', 'CNR60', 'CNR40', 'gCNR', 'gCNR60', 'gCNR40', 'sSNR', 'time');
for m = 1:nm
  A = max(img(:,:,m), 0); mx = max(A(:)); ci = max(A(cyst), 1e-12*mx); bi = max(A(bgr), 1e-12*mx);
  res = zeros(1, 5); DRs = [Inf 60 40];
  for d = 1:3
    vc = max(20*log10(ci/mx), -DRs(d)); vb = max(20*log10(bi/mx), -DRs(d));
    v = [vc; vb]; lab = [zeros(numel(vc),1); ones(numel(vb),1)]; [vs, o] = sort(v); lab = lab(o); n = numel(v);
    bins = floor((0:n-1)'*100/n); [~, first, ic] = unique(vs, 'first'); bins = bins(first(ic));
    hc = accumarray(bins(lab == 0)+1, 1, [100 1])/numel(vc); hb = accumarray(bins(lab == 1)+1, 1, [100 1])/numel(vb);
    res(d) = 1 - sum(min(hc, hb));
    if d > 1, res(3+d-1) = abs(mean(vb) - mean(vc))/sqrt(var(vb,1) + var(vc,1)); end
  end
  CR = 10*log10(mean(bi.^2)/mean(ci.^2)); sSNR = mean(bi)/std(bi,1);
  sel = xp < 0; pr = prof(sel, m); xq = xp(sel); [pk, ip] = max(pr.*(abs(xq(:) + 3e-3) < 0.6e-3)); half = pk/2;
  il = ip; while il > 1 && pr(il) > half, il = il - 1; end
  ir = ip; while ir < numel(pr) && pr(ir) > half, ir = ir + 1; end
  xl = xq(il) + (half - pr(il))/(pr(il+1) - pr(il))*(xq(il+1) - xq(il));
  xr = xq(ir-1) + (half - pr(ir-1))/(pr(ir) - pr(ir-1))*(xq(ir) - xq(ir-1));
  ml = il; while ml > 1 && pr(ml-1) <= pr(ml), ml = ml - 1; end
  mr = ir; while mr < numel(pr) && pr(mr+1) <= pr(mr), mr = mr + 1; end
  side = pr; side(ml:mr) = 0; side(abs(xq(:) - xq(ip)) >= 2.5e-3) = 0; PSL = 20*log10(max(max(side), 1e-9*pk)/pk);
  q = prof(~sel, m); xx = xp(~sel); p1 = max(q(abs(xx - 2.25e-3) < 0.1e-3)); p2 = max(q(abs(xx - 2.75e-3) < 0.1e-3));
  dip = 20*log10(min(q(xx > 2.3e-3 & xx < 2.7e-3))/min(p1, p2));
  fprintf('%-7s %6.3f %6.1f %6.1f %6.1f %7.2f %7.2f %6.3f %7.3f %7.3f %6.2f %6.1fs\n', names{m}, (xr-xl)*1e3, PSL, dip, CR, res(4), res(5), res(1), res(2), res(3), sSNR, T(m));
end

%% Figure: 2 x 5 image grid, 60 dB (SLSC linear)
figure('Position', [50 50 1300 560]);
for m = 1:nm
  subplot(2, 5, m); A = max(img(:,:,m), 0); A = A/max(A(:));
  if m == 9, imagesc(xi*1e3, zi*1e3, A, [0 1]); else, imagesc(xi*1e3, zi*1e3, 20*log10(A + 1e-12), [-60 0]); end
  axis image; colormap(gray); title(names{m}); xlabel('x (mm)'); ylabel('z (mm)');
end

What to look for: DAS FWHM 0.479 mm, MV and ESBMV 0.031 mm with pair dips of −12.2 and −8.7 dB; CF CR 49.9 dB with gCNR 0.914 against DAS 23.4 dB and 0.944; ESBMV CR 63.9 dB, gCNR 0.992 on the linear output but 0.807 at a 40 dB display; SLSC the lowest CR (15.2 dB) and the highest gCNR (0.994); MV costs about 9 s and ESBMV about 16 s in Octave (machine-dependent) against well under 3 s for each of the rest. This is the reproducible counterpart of Demo 25.6.

m25_mv_esbmv.mMV beam pattern with an interferer, signal cancellation versus L, the lone-point width law, then MV (L, K) and ESBMV images · change L, Kt, Delta, delta · runtime ≈ 21 s

%% m25_mv_esbmv.m  -  minimum variance (Capon) and eigenspace-based MV for ultrasound
% Purpose : (A) narrowband focal-plane analysis: MV vs DAS beam pattern with an
%           interferer, signal cancellation of a coherent echo versus subarray length L,
%           and the lone-point width law 0.914/sqrt(N*SNR); (B) broadband plane-wave
%           images of a point pair and a cyst in speckle with DAS, MV (L = 16, 32;
%           temporal averaging K = 0, 2) and ESBMV, with FWHM, pair dip, CR, gCNR,
%           speckle SNR and run time.
% Change  : dx_int, INR_dB (part A); L, Kt, Delta, delta, noise_dB (part B).
% Expect  : coherent echo cancelled (about -66 dB) at L = N, recovered (-0.3 dB) at
%           L = 48; MV point FWHM ~0.03 mm vs 0.48 mm DAS; K = 0 ruins speckle SNR.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8. Runtime ~10 s.
clear; close all; rng(2);
c = 1540; f0 = 5e6; lam = c/f0; N = 64; pitch = 0.3e-3; zF = 30e-3;
psi = @(dx) 2*pi*pitch*dx/(lam*zF);                  % phase step per element for a lateral offset dx at zF

%% A1. Beam pattern of MV weights (exact covariance, incoherent interferer)
dx_int = 2e-3; INR_dB = 10; SNR_dB = 20; L = 32; Delta = 1e-2; P = N - L + 1;
v = @(dx, n) exp(1j*psi(dx)*(0:n-1)');
R = ones(L) + 10^(INR_dB/10)*(v(dx_int,L)*v(dx_int,L)') + 10^(-SNR_dB/10)*eye(L);
R = R + Delta*real(trace(R))/L*eye(L); w = R\ones(L,1); w = w/sum(w);
xs = linspace(-8e-3, 8e-3, 1601); Bmv = zeros(size(xs)); Bdas = Bmv;
for k = 1:numel(xs)
  Bmv(k) = abs(w'*v(xs(k),L))*abs(mean(exp(1j*psi(xs(k))*(0:P-1))));   % subarray weights x Dirichlet of the smoothing
  Bdas(k) = abs(sum(v(xs(k),N)))/N;
end
[~, ki] = min(abs(xs - dx_int));
fprintf('A1 interferer at %.1f mm, %d dB: MV response %.1f dB, DAS %.1f dB (exact covariance)\n', dx_int*1e3, INR_dB, 20*log10(Bmv(ki)), 20*log10(Bdas(ki)));

%% A2. Signal cancellation: coherent interferer (same pulse), desired-signal gain vs L
dxc = 1.5e-3; beta = 1; u = 1 + beta*exp(1j*psi(dxc)*(0:N-1)'); Ls = 1:N; gain = zeros(size(Ls));
for iL = Ls
  P = N - iL + 1; R = zeros(iL);
  for l = 1:P, uu = u(l:l+iL-1); R = R + uu*uu'; end
  R = R/P + 10^(-SNR_dB/10)*eye(iL); R = R + Delta*real(trace(R))/iL*eye(iL);
  w = R\ones(iL,1); w = w/sum(w);
  gain(iL) = abs(1 + beta*(w'*exp(1j*psi(dxc)*(0:iL-1)'))*mean(exp(1j*psi(dxc)*(0:P-1))));
end
fprintf('A2 coherent echo at %.1f mm: desired-signal gain %.1f dB (L = 64), %.1f dB (L = 48), %.1f dB (L = 32), %.1f dB (L = 1)\n', dxc*1e3, 20*log10(gain([64 48 32 1])));

%% A3. Lone point: MV -6 dB width vs SNR (closed form vs direct solve)
xx = linspace(-0.3e-3, 0.3e-3, 6001); das = abs(sum(exp(1j*psi(linspace(-1e-3,1e-3,20001))'*(0:N-1)), 2))/N;
xd = linspace(-1e-3, 1e-3, 20001); wd = max(xd(das >= 0.5)) - min(xd(das >= 0.5));
for snr = [10 20 30]
  s2 = 10^(-snr/10); va = abs(sum(exp(1j*psi(xx)'*(0:N-1)), 2)); y = s2*va./(N*(s2 + N) - va.^2);
  Ri = inv(ones(N) + s2*eye(N)); vt = exp(1j*psi(0.05e-3)*(0:N-1)'); wt = Ri*vt/(vt'*Ri*vt); yd = abs(wt'*ones(N,1));
  [~, k5] = min(abs(xx - 0.05e-3)); ii = find(y >= 0.5); i1 = ii(1); i2 = ii(end);  % interpolate the -6 dB crossings (grid step 0.1 um)
  w6 = interp1(y([i2 i2+1]), xx([i2 i2+1]), 0.5) - interp1(y([i1-1 i1]), xx([i1-1 i1]), 0.5);
  fprintf('A3 SNR %2d dB: MV -6 dB width %.4f mm = %.4f x DAS (law 0.914/sqrt(N SNR) = %.4f); closed form vs solve at 0.05 mm: %.2e\n', snr, w6*1e3, w6/wd, 0.914/sqrt(N*10^(snr/10)), abs(y(k5) - yd));
end

%% B. Broadband plane-wave data: pair at z = 24 mm, cyst in speckle, channel noise
fs = 20e6; bw = 0.6; sigf = bw*f0/(2*sqrt(2*log(2))); sigt = 1/(2*pi*sigf); ang = (-4:4)*pi/180; Na = numel(ang);
xe = ((0:N-1) - (N-1)/2)*pitch; t0 = 25e-6; Nt = 560; tt = t0 + (0:Nt-1)'/fs; noise_dB = -20;
area = [-5.5 5.5 26 34]*1e-3; nsc = round(11*8*80); xsc = area(1) + 11e-3*rand(nsc,1); zsc = area(3) + 8e-3*rand(nsc,1);
keep = (xsc.^2 + (zsc-30e-3).^2) > (2.5e-3)^2; xsc = [xsc(keep); -0.25e-3; 0.25e-3; -2.5e-3]; zsc = [zsc(keep); 24e-3; 24e-3; 24e-3];
asc = [randn(sum(keep),1); 15; 15; 15];
IQ = zeros(Nt, N, Na); K4 = ceil(4*sigt*fs); E = repmat(1:N, numel(xsc), 1); tic;
for a = 1:Na
  tau = (xsc*sin(ang(a)) + zsc*cos(ang(a)))/c + sqrt((xsc - xe).^2 + zsc.^2)/c; ph = asc.*exp(-1j*2*pi*f0*tau); n0 = (tau - t0)*fs;
  accr = zeros(Nt*N,1); acci = accr;
  for k = -K4:K4
    n = floor(n0) + k; g = exp(-((n - n0)/fs).^2/(2*sigt^2)); ok = n >= 0 & n < Nt; idx = n(ok) + 1 + (E(ok)-1)*Nt; vv = g(ok).*ph(ok);
    accr = accr + accumarray(idx, real(vv), [Nt*N 1]); acci = acci + accumarray(idx, imag(vv), [Nt*N 1]);
  end
  IQ(:,:,a) = reshape(accr + 1j*acci, Nt, N);
end
Pref = mean(mean(abs(IQ(tt > 2*27.5e-3/c & tt < 2*33e-3/c, :, 5)).^2)); taps = exp(-((-K4:K4)'/fs).^2/(2*sigt^2));
for a = 1:Na, IQ(:,:,a) = IQ(:,:,a) + sqrt(Pref*10^(noise_dB/10)/sum(taps.^2))*conv2((randn(Nt+2*K4,N) + 1j*randn(Nt+2*K4,N))/sqrt(2), taps, 'valid'); end
fprintf('B  simulated %d scatterers in %.1f s\n', numel(xsc), toc);
dz = 0.075e-3; Kt = 2; delta = 0.5; Delta = 1e-2;
grids = {(-4:0.1:4)*1e-3, 23e-3 + (0:133)*dz; [(-0.8:0.005:0.8) (-2.8:0.005:-2.2)]*1e-3, 24e-3 + (-Kt:Kt)*dz};
for gi = 1:2
  xg = grids{gi,1}; zg = grids{gi,2}; [Zp, Xp] = ndgrid(zg, xg); Xp = Xp(:); Zp = Zp(:); S = zeros(numel(Xp), N);
  Ptx = exp(1j*2*pi*f0*(Xp*sin(ang) + Zp*cos(ang))/c); Ttx = (Xp*sin(ang) + Zp*cos(ang))/c;
  for e = 1:N
    trx = sqrt((Xp - xe(e)).^2 + Zp.^2)/c; acc = 0;
    for a = 1:Na
      fn = (Ttx(:,a) + trx - t0)*fs; n = floor(fn); fr = fn - n; ok = n >= 0 & n < Nt-1; col = IQ(:,e,a); vv = zeros(size(fn));
      vv(ok) = col(n(ok)+1).*(1-fr(ok)) + col(n(ok)+2).*fr(ok); acc = acc + vv.*Ptx(:,a);
    end
    S(:,e) = acc.*exp(1j*2*pi*f0*trx);
  end
  S = reshape(S, numel(zg), numel(xg), N); nz = numel(zg); nx = numel(xg);
  if gi == 1, rows = 1:nz; else, rows = Kt+1; end
  out = zeros(nz, nx, 5); tm = zeros(1,5);                % DAS, MV L32 K2, MV L32 K0, MV L16 K2, ESBMV L32 K2
  out(:,:,1) = abs(sum(S,3))/N;
  cfg = [32 Kt; 32 0; 16 Kt];
  for ix = 1:nx
    for r = rows
      for ci = 1:3
        t1 = tic; L = cfg(ci,1); K = cfg(ci,2); P = N - L + 1; id = (0:L-1)' + (1:P); X = [];
        for k = -K:K, if r+k >= 1 && r+k <= nz, x = squeeze(S(r+k, ix, :)); X = [X, x(id)]; end, end %#ok<AGROW>
        Rm = (X*X')/size(X,2); w = (Rm + Delta*real(trace(Rm))/L*eye(L))\ones(L,1); w = w/sum(w);
        x = squeeze(S(r, ix, :)); b = mean(x(id), 2); out(r, ix, 1+ci) = abs(w'*b); tm(1+ci) = tm(1+ci) + toc(t1);
        if ci == 1
          t2 = tic; [V, D] = eig((Rm + Rm')/2); [lv, o] = sort(real(diag(D)), 'descend'); Es = V(:, o(lv >= delta*lv(1)));
          out(r, ix, 5) = abs((Es*(Es'*w))'*b); tm(5) = tm(5) + toc(t2);
        end
      end
    end
  end
  if gi == 1, img = out; xi = xg; zi = zg; T = tm; T(5) = T(5) + T(2); else, prof = squeeze(out(Kt+1, :, :)); xp = xg; end
end

%% Metrics
names = {'DAS', 'MV L32 K2', 'MV L32 K0', 'MV L16 K2', 'ESBMV'};
[ZI, XI] = ndgrid(zi, xi); cyst = (XI.^2 + (ZI-30e-3).^2) < (2e-3)^2; bgr = abs(XI) >= 3e-3 & abs(XI) <= 3.9e-3 & ZI >= 28e-3 & ZI <= 32e-3;
fprintf('\n%-10s %7s %7s %6s %6s %6s %7s\n', 'method', 'FWHM', 'dip', 'CR', 'gCNR', 'sSNR', 'time');
for m = 1:5
  A = img(:,:,m); ci = A(cyst); bi = A(bgr); v = [ci; bi]; lab = [zeros(numel(ci),1); ones(numel(bi),1)]; [~, o] = sort(v); lab = lab(o);
  bins = floor((0:numel(v)-1)'*100/numel(v)); g = 1 - sum(min(accumarray(bins(lab==0)+1,1,[100 1])/numel(ci), accumarray(bins(lab==1)+1,1,[100 1])/numel(bi)));
  pr = prof(:,m); sel = xp < -2e-3; q = pr(sel); xq = xp(sel); [pk, ip] = max(q); above = xq(q >= pk/2);
  pp = pr(~sel); xx = xp(~sel); p1 = max(pp(abs(xx + 0.25e-3) < 0.05e-3)); p2 = max(pp(abs(xx - 0.25e-3) < 0.05e-3));
  dip = 20*log10(min(pp(abs(xx) < 0.15e-3))/min(p1, p2));
  fprintf('%-10s %6.3fmm %6.1fdB %6.1f %6.3f %6.2f %6.1fs\n', names{m}, (max(above) - min(above))*1e3, dip, 10*log10(mean(bi.^2)/mean(ci.^2)), g, mean(bi)/std(bi,1), T(m));
end

%% Figures
figure('Position', [50 50 1200 420]);
subplot(1,2,1); plot(xs*1e3, 20*log10(Bdas), 'k', xs*1e3, 20*log10(Bmv + 1e-12), 'b'); hold on; plot([1 1]*dx_int*1e3, [-80 0], 'r--');
ylim([-80 2]); xlabel('lateral position at z = 30 mm (mm)'); ylabel('beam pattern (dB)'); legend('DAS', 'MV (L = 32)', 'interferer'); title('A1: MV places a null on the interferer');
subplot(1,2,2); plot(Ls, 20*log10(gain), 'b', 'LineWidth', 1.5); xlabel('subarray length L'); ylabel('desired-signal gain (dB)'); title('A2: coherent echo cancels without smoothing'); grid on;
figure('Position', [50 50 1300 330]);
for m = 1:5
  subplot(1,5,m); A = img(:,:,m)/max(max(img(:,:,m))); imagesc(xi*1e3, zi*1e3, 20*log10(A + 1e-12), [-60 0]); axis image; colormap(gray);
  title(names{m}); xlabel('x (mm)'); ylabel('z (mm)');
end

What to look for: interferer response −120.8 dB (MV) against −28.8 dB (DAS); desired-signal gain −66.0 dB at L = 64 and −0.3 dB at L = 48; lone-point widths 0.0361, 0.0114 and 0.0036 × DAS at 10, 20, 30 dB, matching \(0.914/\sqrt{N\cdot\mathrm{SNR}}\); in the image, K = 0 drops the speckle SNR to 0.98 and gCNR to 0.86, while K = 2 keeps 1.67 and 0.98; ESBMV CR 67 dB at gCNR 0.995.

m25_coherence_factors.mCF, Wiener, GCF, PCF and SCF: identities, expected values in speckle and noise, images and metrics · change M0, gam, noise level · runtime ≈ 8 s

%% m25_coherence_factors.m  -  CF, Wiener postfilter, GCF, PCF and SCF on one dataset
% Purpose : compute the five coherence weights from the same plane-wave aperture data
%           (cyst in speckle, bright off-axis scatterer, channel noise), verify the
%           identity H = N CF/((N-1) CF + 1), compare the speckle CF with the VCZ
%           prediction for this transmit sequence, and report CR, gCNR, speckle SNR.
% Change  : noise_dB, M0 list, gam (PCF), pSCF.
% Expect  : speckle CF ratio of means ~0.23 (VCZ triangle 0.225), CF in noise 1/N;
%           CR rises CF > PCF > GCF > Wiener > DAS while gCNR does not follow.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8. Runtime ~6 s.
clear; close all; rng(3);
c = 1540; f0 = 5e6; fs = 20e6; lam = c/f0; N = 64; pitch = 0.3e-3; xe = ((0:N-1) - (N-1)/2)*pitch;
bw = 0.6; sigf = bw*f0/(2*sqrt(2*log(2))); sigt = 1/(2*pi*sigf); ang = (-4:4)*pi/180; Na = numel(ang);
t0 = 25e-6; Nt = 580; tt = t0 + (0:Nt-1)'/fs; noise_dB = -15;

%% Simulation: speckle with an anechoic cyst (r = 2.5 mm) and a bright scatterer 4.2 mm off-axis
nsc = round(13.6*8.6*80); xs = (rand(nsc,1)*13.6 - 6.8)*1e-3; zs = (26 + 8.6*rand(nsc,1))*1e-3;
keep = (xs.^2 + (zs-30e-3).^2) > (2.5e-3)^2; xs = [xs(keep); 4.2e-3]; zs = [zs(keep); 30.5e-3]; as = [randn(sum(keep),1); 60];
IQ = zeros(Nt, N, Na); K4 = ceil(4*sigt*fs); E = repmat(1:N, numel(xs), 1);
for a = 1:Na
  tau = (xs*sin(ang(a)) + zs*cos(ang(a)))/c + sqrt((xs - xe).^2 + zs.^2)/c; ph = as.*exp(-1j*2*pi*f0*tau); n0 = (tau - t0)*fs;
  accr = zeros(Nt*N,1); acci = accr;
  for k = -K4:K4
    n = floor(n0) + k; g = exp(-((n - n0)/fs).^2/(2*sigt^2)); ok = n >= 0 & n < Nt; idx = n(ok) + 1 + (E(ok)-1)*Nt; v = g(ok).*ph(ok);
    accr = accr + accumarray(idx, real(v), [Nt*N 1]); acci = acci + accumarray(idx, imag(v), [Nt*N 1]);
  end
  IQ(:,:,a) = reshape(accr + 1j*acci, Nt, N);
end
Pref = mean(mean(abs(IQ(tt > 2*27.5e-3/c & tt < 2*33e-3/c, :, 5)).^2)); taps = exp(-((-K4:K4)'/fs).^2/(2*sigt^2));
for a = 1:Na, IQ(:,:,a) = IQ(:,:,a) + sqrt(Pref*10^(noise_dB/10)/sum(taps.^2))*conv2((randn(Nt+2*K4,N) + 1j*randn(Nt+2*K4,N))/sqrt(2), taps, 'valid'); end

%% Aperture data on the image grid
xg = (-6:0.15:6)*1e-3; zg = (26:0.075:34)*1e-3; [Zp, Xp] = ndgrid(zg, xg); Xp = Xp(:); Zp = Zp(:); S = zeros(numel(Xp), N);
Ttx = (Xp*sin(ang) + Zp*cos(ang))/c; Ptx = exp(1j*2*pi*f0*Ttx);
for e = 1:N
  trx = sqrt((Xp - xe(e)).^2 + Zp.^2)/c; acc = 0;
  for a = 1:Na
    fn = (Ttx(:,a) + trx - t0)*fs; n = floor(fn); fr = fn - n; ok = n >= 0 & n < Nt-1; col = IQ(:,e,a); v = zeros(size(fn));
    v(ok) = col(n(ok)+1).*(1-fr(ok)) + col(n(ok)+2).*fr(ok); acc = acc + v.*Ptx(:,a);
  end
  S(:,e) = acc.*exp(1j*2*pi*f0*trx);
end
nz = numel(zg); nx = numel(xg);

%% Weights
ybar = mean(S, 2); P = mean(abs(S).^2, 2); sig2 = mean(abs(S - ybar).^2, 2);
CF = abs(sum(S,2)).^2./(N*sum(abs(S).^2,2) + eps);                  % |sum s|^2/(N sum |s|^2)
Hw = abs(ybar).^2./(abs(ybar).^2 + sig2/N + eps);                    % Wiener postfilter from its definition
fprintf('identity check: max |H - N CF/((N-1)CF+1)| = %.2e,  max |P - (|ybar|^2 + sig2)| = %.2e\n', max(abs(Hw - N*CF./((N-1)*CF + 1))), max(abs(P - abs(ybar).^2 - sig2))/max(P));
Sk = fft(S, [], 2); M0s = [0 1 3]; GCF = zeros(numel(Xp), numel(M0s));
for i = 1:numel(M0s), low = abs(Sk(:,1)).^2; for k = 1:M0s(i), low = low + abs(Sk(:,1+k)).^2 + abs(Sk(:,N+1-k)).^2; end, GCF(:,i) = low./(sum(abs(Sk).^2,2) + eps); end
fprintf('GCF(M0 = 0) equals CF: max difference %.2e\n', max(abs(GCF(:,1) - CF)));
gam = 1; phs = angle(S); pha = phs - pi*sign(phs); PCF = max(0, 1 - gam*min(std(phs,1,2), std(pha,1,2))/(pi/sqrt(3)));
pSCF = 1; bbar = mean(sign(real(S)), 2); SCF = abs(1 - sqrt(max(0, 1 - bbar.^2))).^pSCF;

%% Statistics: speckle CF vs VCZ, noise CF vs 1/N
[ZI, XI] = ndgrid(zg, xg); cyst = (XI.^2 + (ZI-30e-3).^2) < (2e-3)^2; bgr = XI >= -5.8e-3 & XI <= -3.4e-3 & ZI >= 28e-3 & ZI <= 32e-3;
mtx = Na*30e-3*(sin(ang(2)) - sin(ang(1)))/pitch; m = -(N-1):(N-1); Ctri = max(0, 1 - abs(m)/mtx);
Sb = S(bgr(:), :);
fprintf('speckle: CF ratio of means %.3f (VCZ triangle, m_tx = %.1f: %.3f), mean of per-pixel CF %.3f\n', mean(abs(sum(Sb,2)).^2)/(N*mean(sum(abs(Sb).^2,2))), mtx, sum((N - abs(m)).*Ctri)/N^2, mean(CF(bgr(:))));
Wn = (randn(20000, N) + 1j*randn(20000, N))/sqrt(2); cfn = abs(sum(Wn,2)).^2./(N*sum(abs(Wn).^2,2));
fprintf('white noise: mean CF %.4f (1/N = %.4f), mean Wiener H %.3f (N/(2N-1) = %.3f at CF = 1/N), mean PCF %.3f\n', mean(cfn), 1/N, mean(N*cfn./((N-1)*cfn + 1)), N/(2*N-1), mean(max(0, 1 - std(angle(Wn),1,2)/(pi/sqrt(3)))));

%% Image metrics (amplitude = weight x |DAS|)
das = abs(ybar); imgs = [das, CF.*das, Hw.*das, GCF(:,2).*das, GCF(:,3).*das, PCF.*das, SCF.*das];
names = {'DAS', 'CF', 'Wiener', 'GCF M0=1', 'GCF M0=3', 'PCF', 'SCF'};
fprintf('\n%-9s %6s %6s %6s %8s\n', 'method', 'CR', 'gCNR', 'sSNR', 'w(spk)');
for i = 1:numel(names)
  A = imgs(:,i); ci = max(A(cyst(:)), 1e-12*max(A)); bi = A(bgr(:)); v = [ci; bi]; lab = [zeros(numel(ci),1); ones(numel(bi),1)];
  [vs, o] = sort(v); lab = lab(o); bins = floor((0:numel(v)-1)'*100/numel(v)); [~, first, ic] = unique(vs, 'first'); bins = bins(first(ic));
  g = 1 - sum(min(accumarray(bins(lab==0)+1,1,[100 1])/numel(ci), accumarray(bins(lab==1)+1,1,[100 1])/numel(bi)));
  fprintf('%-9s %6.1f %6.3f %6.2f %8.3f\n', names{i}, 10*log10(mean(bi.^2)/mean(ci.^2)), g, mean(bi)/std(bi,1), mean(A(bgr(:))./das(bgr(:))));
end

%% Figures
figure('Position', [50 50 1300 330]);
for i = 1:numel(names)
  subplot(1, numel(names), i); A = reshape(imgs(:,i), nz, nx); A = A/max(A(:));
  imagesc(xg*1e3, zg*1e3, 20*log10(A + 1e-12), [-60 0]); axis image; colormap(gray); title(names{i}); xlabel('x (mm)'); ylabel('z (mm)');
end
figure; edges = linspace(0, 1, 41); hs = histc(CF(bgr(:)), edges); hn = histc(cfn, edges);
bar(edges + 0.0125, [hs/sum(hs) hn/sum(hn)]*40, 1); xlim([0 1]); legend('speckle (background ROI)', 'white noise'); xlabel('CF'); ylabel('pdf'); title('CF is broadly distributed in speckle');

What to look for: the Wiener identity and \(\mathrm{GCF}(M_0 = 0) = \mathrm{CF}\) hold to \(10^{-15}\); speckle CF 0.233 against the VCZ prediction 0.225; noise CF 0.0158 against 1/N = 0.0156; CF raises CR from 21 to 46 dB but lowers gCNR (0.914 against 0.939) and speckle SNR (0.82); GCF with \(M_0 = 3\) keeps speckle SNR 1.64 and the best gCNR; SCF stretches most and has the lowest gCNR.

m25_dmas.mO(N) identity, the DC + 2f₀ spectrum, the narrowband sum-coarray beam, DAS / DMAS / F-DMAS images · change noise_dB, the band-pass width · runtime ≈ 12 s

%% m25_dmas.m  -  delay-multiply-and-sum (DMAS / F-DMAS)
% Purpose : verify the O(N) identity sum_{i<j} s_i s_j = ((sum s)^2 - sum s^2)/2 for the
%           signed-square-root products, show the DC + 2f0 spectrum of the DMAS sum,
%           compute the narrowband sum-coarray beam (DMAS at 2f0) against DAS, and
%           compare DAS, unfiltered DMAS and F-DMAS images (points, cyst in speckle).
% Change  : noise_dB, the band-pass width (sigh), N.
% Expect  : identity error ~1e-13; narrowband PSL -13.3 dB (DAS) vs -26.5 dB (DMAS),
%           width ratio 0.74; images: F-DMAS PSL ~-27 dB vs ~-19 dB DAS.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8. Runtime ~6 s.
clear; close all; rng(4);
c = 1540; f0 = 5e6; fs = 20e6; lam = c/f0; N = 64; pitch = 0.3e-3; xe = ((0:N-1) - (N-1)/2)*pitch;
bw = 0.6; sigf = bw*f0/(2*sqrt(2*log(2))); sigt = 1/(2*pi*sigf); ang = (-4:4)*pi/180; Na = numel(ang);

%% 1. The O(N) identity (random RF aperture vectors)
s = randn(1000, N); sh = sign(s).*sqrt(abs(s)); slow = zeros(1000,1);
for i = 1:N-1, for j = i+1:N, slow = slow + sign(s(:,i).*s(:,j)).*sqrt(abs(s(:,i).*s(:,j))); end, end
fast = (sum(sh,2).^2 - sum(abs(s),2))/2;
fprintf('identity: max |pairwise - O(N)| = %.2e (relative %.1e); %d pair products replaced by %d square roots\n', max(abs(slow - fast)), max(abs(slow - fast))/max(abs(slow)), N*(N-1)/2, N);

%% 2. Spectrum of the DMAS sum for an aligned echo: every pair term equals |rf(t)|
fsr = 200e6; t = (-2048:2047)'/fsr; rf = exp(-t.^2/(2*sigt^2)).*cos(2*pi*f0*t); Y = abs(fft(abs(rf))); fsp = (0:4095)'*fsr/4096;
e0 = sum(Y(fsp < 2.5e6).^2); e2 = sum(Y(fsp > 7.5e6 & fsp < 12.5e6).^2); e4 = sum(Y(fsp > 17.5e6 & fsp < 22.5e6).^2);
fprintf('aligned echo, |rf| spectrum: energy DC band : 2f0 band : 4f0 band = 1 : %.3f : %.4f (positive frequencies)\n', e2/e0, e4/e0);

%% 3. Narrowband sum coarray at z = 30 mm: 2f0 term of pair (i,j) has phase phi_i + phi_j
x = linspace(-8e-3, 8e-3, 32001); q = 2*pi*pitch*x/(lam*30e-3); AF = exp(1j*q'*(0:N-1)); a1 = sum(AF,2); a2 = sum(AF.^2,2);
dasNB = abs(a1)/N; dmNB = abs((a1.^2 - a2)/2)/(N*(N-1)/2);
wid = @(y) max(x(y >= 0.5)) - min(x(y >= 0.5));
i0 = 16001; il = i0; while dasNB(il-1) <= dasNB(il), il = il - 1; end
jl = i0; while dmNB(jl-1) <= dmNB(jl), jl = jl - 1; end
fprintf('narrowband: DAS FWHM %.3f mm PSL %.1f dB; DMAS(2f0) FWHM %.3f mm PSL %.1f dB; width ratio %.3f\n', wid(dasNB)*1e3, 20*log10(max(dasNB(1:il))), wid(dmNB)*1e3, 20*log10(max(dmNB(1:jl))), wid(dmNB)/wid(dasNB));

%% 4. Plane-wave simulation: points at z = 24 mm, cyst in speckle, noise
t0 = 25e-6; Nt = 580; tt = t0 + (0:Nt-1)'/fs; noise_dB = -20;
nsc = round(13.6*8.6*80); xs = (rand(nsc,1)*13.6 - 6.8)*1e-3; zs = (26 + 8.6*rand(nsc,1))*1e-3;
keep = (xs.^2 + (zs-30e-3).^2) > (2.5e-3)^2; xs = [xs(keep); -3e-3; 2.25e-3; 2.75e-3]; zs = [zs(keep); 24e-3; 24e-3; 24e-3]; as = [randn(sum(keep),1); 15; 15; 15];
IQ = zeros(Nt, N, Na); K4 = ceil(4*sigt*fs); E = repmat(1:N, numel(xs), 1);
for a = 1:Na
  tau = (xs*sin(ang(a)) + zs*cos(ang(a)))/c + sqrt((xs - xe).^2 + zs.^2)/c; ph = as.*exp(-1j*2*pi*f0*tau); n0 = (tau - t0)*fs;
  accr = zeros(Nt*N,1); acci = accr;
  for k = -K4:K4
    n = floor(n0) + k; g = exp(-((n - n0)/fs).^2/(2*sigt^2)); ok = n >= 0 & n < Nt; idx = n(ok) + 1 + (E(ok)-1)*Nt; v = g(ok).*ph(ok);
    accr = accr + accumarray(idx, real(v), [Nt*N 1]); acci = acci + accumarray(idx, imag(v), [Nt*N 1]);
  end
  IQ(:,:,a) = reshape(accr + 1j*acci, Nt, N);
end
Pref = mean(mean(abs(IQ(tt > 2*27.5e-3/c & tt < 2*33e-3/c, :, 5)).^2)); taps = exp(-((-K4:K4)'/fs).^2/(2*sigt^2));
for a = 1:Na, IQ(:,:,a) = IQ(:,:,a) + sqrt(Pref*10^(noise_dB/10)/sum(taps.^2))*conv2((randn(Nt+2*K4,N) + 1j*randn(Nt+2*K4,N))/sqrt(2), taps, 'valid'); end

%% 5. RF aperture data on a fine axial grid (lambda/16), then DAS, DMAS, F-DMAS
sub = 4; dz = 0.075e-3; dzf = dz/sub; fsz = c/(2*dzf); xg = (-6:0.05:6)*1e-3; zf = (22e-3:dzf:34e-3)';
[Zp, Xp] = ndgrid(zf, xg); Xp = Xp(:); Zp = Zp(:); S = zeros(numel(Xp), N); Ttx = (Xp*sin(ang) + Zp*cos(ang))/c; Ptx = exp(1j*2*pi*f0*Ttx);
for e = 1:N
  trx = sqrt((Xp - xe(e)).^2 + Zp.^2)/c; acc = 0;
  for a = 1:Na
    fn = (Ttx(:,a) + trx - t0)*fs; n = floor(fn); fr = fn - n; ok = n >= 0 & n < Nt-1; col = IQ(:,e,a); v = zeros(size(fn));
    v(ok) = col(n(ok)+1).*(1-fr(ok)) + col(n(ok)+2).*fr(ok); acc = acc + v.*Ptx(:,a);
  end
  S(:,e) = acc.*exp(1j*2*pi*f0*trx);
end
nzf = numel(zf); nx = numel(xg); R = reshape(real(S), nzf, nx, N); Sc = reshape(S, nzf, nx, N);
tic; das = abs(sum(Sc,3))/N; tdas = toc;
tic; shh = sign(R).*sqrt(abs(R)); yd = (sum(shh,3).^2 - sum(abs(R),3))/2;
nf = 2^nextpow2(nzf + 64); f = [0:nf/2, -nf/2+1:-1]'*fsz/nf; sigh = sqrt(2)*sigf;
Yd = fft(yd, nf, 1); envF = abs(ifft(Yd.*(2*exp(-(f - 2*f0).^2/(2*sigh^2)).*(f > 0)), [], 1)); envF = envF(1:nzf,:);
envU = abs(ifft(Yd.*(2*(f > 0) + (f == 0)), [], 1)); envU = envU(1:nzf,:); tdm = toc;   % unfiltered: analytic-signal envelope
imgs = {das(1:sub:end,:), envU(1:sub:end,:), envF(1:sub:end,:)}; names = {'DAS', 'DMAS, no filter', 'F-DMAS'}; zo = zf(1:sub:end);
fprintf('time per image: DAS %.2f s, DMAS + band-pass %.2f s (%d x %d RF samples)\n', tdas, tdm, nzf, nx);
[ZI, XI] = ndgrid(zo, xg); cyst = (XI.^2 + (ZI-30e-3).^2) < (2e-3)^2; bgr = XI >= -5.8e-3 & XI <= -3.4e-3 & ZI >= 28e-3 & ZI <= 32e-3;
[~, rp] = min(abs(zo - 24e-3)); xsel = xg >= -5.5e-3 & xg <= -0.5e-3; xq = xg(xsel);
fprintf('\n%-16s %6s %6s %6s %6s %6s\n', 'image', 'FWHM', 'PSL', 'CR', 'gCNR', 'sSNR');
for i = 1:3
  A = imgs{i}; pr = A(rp, xsel); [pk, ip] = max(pr); half = pk/2;
  il = ip; while pr(il) > half, il = il - 1; end; ir = ip; while pr(ir) > half, ir = ir + 1; end
  xl = xq(il) + (half - pr(il))/(pr(il+1) - pr(il))*(xq(il+1) - xq(il)); xr = xq(ir-1) + (half - pr(ir-1))/(pr(ir) - pr(ir-1))*(xq(ir) - xq(ir-1));
  ml = il; while ml > 1 && pr(ml-1) <= pr(ml), ml = ml - 1; end; mr = ir; while mr < numel(pr) && pr(mr+1) <= pr(mr), mr = mr + 1; end
  side = pr; side(ml:mr) = 0; ci = A(cyst); bi = A(bgr); v = [ci; bi]; lab = [zeros(numel(ci),1); ones(numel(bi),1)]; [~, o] = sort(v); lab = lab(o);
  bins = floor((0:numel(v)-1)'*100/numel(v)); g = 1 - sum(min(accumarray(bins(lab==0)+1,1,[100 1])/numel(ci), accumarray(bins(lab==1)+1,1,[100 1])/numel(bi)));
  fprintf('%-16s %5.3fmm %6.1f %6.1f %6.3f %6.2f\n', names{i}, (xr - xl)*1e3, 20*log10(max(side)/pk), 10*log10(mean(bi.^2)/mean(ci.^2)), g, mean(bi)/std(bi,1));
end
fprintf('(FWHM from the 0.05 mm image columns through the point at x = -3 mm, z = 24 mm)\n');

%% Figures
figure('Position', [50 50 1200 330]);
for i = 1:3
  subplot(1,4,i); A = imgs{i}/max(imgs{i}(:)); imagesc(xg*1e3, zo*1e3, 20*log10(A + 1e-12), [-60 0]); axis image; colormap(gray); title(names{i}); xlabel('x (mm)'); ylabel('z (mm)');
end
subplot(1,4,4); plot(fsp(1:500)/1e6, Y(1:500)/Y(1), 'b'); hold on; plot(fsp(1:500)/1e6, 0.5*exp(-(fsp(1:500) - 2*f0).^2/(2*sigh^2)), 'r--');
xlabel('frequency (MHz)'); ylabel('|spectrum| (norm.)'); title('DMAS sum of an aligned echo'); legend('|rf| spectrum', 'band-pass');
figure; plot(x*1e3, 20*log10(dasNB + 1e-9), 'k', x*1e3, 20*log10(dmNB + 1e-9), 'b'); ylim([-60 0]); xlim([-4 4]);
xlabel('lateral position (mm)'); ylabel('beam (dB)'); legend('DAS', 'DMAS at 2f_0 (sum coarray)'); title('Narrowband: DMAS = triangular aperture at \lambda/2');

What to look for: identity error \(5 \times 10^{-13}\); \(2f_0\)-to-DC energy ratio 0.221; narrowband FWHM 0.580 against 0.429 mm (ratio 0.740) and PSL −13.3 against −26.5 dB; F-DMAS PSL −27.2 dB against −18.9 dB for DAS, CR 35.6 against 23.1 dB, gCNR 0.961 against 0.943; without the band-pass gCNR collapses to 0.778.

m25_slsc.mSpatial coherence of speckle, point and noise against the VCZ triangle, noise scaling, SLSC images and the M sweep · change noise_dB, M_list, Ks · runtime ≈ 13 s

%% m25_slsc.m  -  short-lag spatial coherence (SLSC) imaging
% Purpose : simulate 64-channel plane-wave data (9 angles, -4..4 deg) of speckle with an
%           anechoic cyst and a point target, measure the spatial coherence R(m) of the
%           delayed RF in speckle, at the point and in the cyst, compare speckle R(m)
%           with the van Cittert-Zernike (VCZ) triangle of the synthesized transmit
%           aperture, show the noise scaling SNR/(1+SNR), and sweep the maximum lag M.
% Change  : noise_dB (channel noise re speckle), M_list, Ks (axial kernel, fine samples).
% Expect  : speckle R(m) close to 1 - m/m_tx with m_tx = 15.7 lags; R_sl(M=8) ~ 5.7 (VCZ);
%           at 0 dB channel SNR the lags m >= 1 halve; point R(m) ~ 1; cyst R(m >= 1) ~ 0;
%           larger M -> narrower point, lower speckle SNR. Runtime ~15 s in Octave.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(2);

%% Acquisition, phantom, channel data (signal and unit noise kept separate: beamforming is linear)
c = 1540; f0 = 5e6; fs = 20e6; lam = c/f0; N = 64; pitch = 0.3e-3;        % m/s, Hz, Hz, m, -, m
xe = ((0:N-1) - (N-1)/2)*pitch; bw = 0.6; sigf = bw*f0/(2*sqrt(2*log(2))); sigt = 1/(2*pi*sigf);
ang = (-4:4)*pi/180; Na = numel(ang); t0 = 25e-6; Nt = 580; tt = t0 + (0:Nt-1)'/fs;
noise_dB = -20; Ks = 8; M_list = [1 2 4 8 12 16 24];
nsc = round(13.6*8.6*80); xs = -6.8e-3 + 13.6e-3*rand(nsc,1); zs = 26e-3 + 8.6e-3*rand(nsc,1);
keep = (xs.^2 + (zs-30e-3).^2) > (2.5e-3)^2; xs = [xs(keep); -3e-3]; zs = [zs(keep); 24e-3]; as = [randn(sum(keep),1); 15];
tic; IQ = zeros(Nt, N, Na); K4 = ceil(4*sigt*fs); E = repmat(1:N, numel(xs), 1);
for a = 1:Na
  tau = (xs*sin(ang(a)) + zs*cos(ang(a)))/c + sqrt((xs - xe).^2 + zs.^2)/c; ph = as.*exp(-1j*2*pi*f0*tau); n0 = (tau - t0)*fs;
  accr = zeros(Nt*N,1); acci = accr;
  for k = -K4:K4
    n = floor(n0) + k; g = exp(-((n - n0)/fs).^2/(2*sigt^2)); ok = n >= 0 & n < Nt; idx = n(ok) + 1 + (E(ok)-1)*Nt; v = g(ok).*ph(ok);
    accr = accr + accumarray(idx, real(v), [Nt*N 1]); acci = acci + accumarray(idx, imag(v), [Nt*N 1]);
  end
  IQ(:,:,a) = reshape(accr + 1j*acci, Nt, N);
end
Pref = mean(mean(abs(IQ(tt > 2*27.5e-3/c & tt < 2*33e-3/c, :, 5)).^2));   % mean speckle channel power
taps = exp(-((-K4:K4)'/fs).^2/(2*sigt^2)); NQ = zeros(Nt, N, Na);
for a = 1:Na, NQ(:,:,a) = sqrt(Pref/sum(taps.^2))*conv2((randn(Nt+2*K4,N) + 1j*randn(Nt+2*K4,N))/sqrt(2), taps, 'valid'); end
fprintf('simulated %d scatterers in %.1f s\n', numel(xs), toc);

%% Delayed, compounded aperture data on a fine axial grid (dz/4, RF rate 41 MHz)
sub = 4; dz = 0.075e-3; dzf = dz/sub; xg = (-6:0.15:6)*1e-3; zf = (22e-3:dzf:34e-3)'; nzf = numel(zf); nx = numel(xg);
[Zp, Xp] = ndgrid(zf, xg); Xp = Xp(:); Zp = Zp(:); Ttx = (Xp*sin(ang) + Zp*cos(ang))/c; Ptx = exp(1j*2*pi*f0*Ttx);
Ssig = zeros(numel(Xp), N); Snoi = Ssig;
for e = 1:N
  trx = sqrt((Xp - xe(e)).^2 + Zp.^2)/c; as_ = 0; an_ = 0;
  for a = 1:Na
    fn = (Ttx(:,a) + trx - t0)*fs; n = floor(fn); fr = fn - n; ok = n >= 0 & n < Nt-1;
    v = zeros(size(fn)); col = IQ(:,e,a); v(ok) = col(n(ok)+1).*(1-fr(ok)) + col(n(ok)+2).*fr(ok); as_ = as_ + v.*Ptx(:,a);
    v = zeros(size(fn)); col = NQ(:,e,a); v(ok) = col(n(ok)+1).*(1-fr(ok)) + col(n(ok)+2).*fr(ok); an_ = an_ + v.*Ptx(:,a);
  end
  Ssig(:,e) = as_.*exp(1j*2*pi*f0*trx); Snoi(:,e) = an_.*exp(1j*2*pi*f0*trx);
end
Ssig = reshape(Ssig, nzf, nx, N); Snoi = reshape(Snoi, nzf, nx, N);
ro = (1:sub:nzf)'; nzo = numel(ro); zo = zf(ro); lo = max(1, ro - Ks); hi = min(nzf, ro + Ks);
[ZI, XI] = ndgrid(zo, xg); cyst = (XI.^2 + (ZI-30e-3).^2) < (2e-3)^2; bgr = XI >= -5.8e-3 & XI <= -3.4e-3 & ZI >= 28e-3 & ZI <= 32e-3;

%% 1. Coherence curves R(m), m = 0..40, in speckle (noise-free and 0 dB), at the point, in the cyst
m_tx = Na*30e-3*(sin(ang(2)) - sin(ang(1)))/pitch;                        % VCZ: synthesized transmit width in lags
fprintf('VCZ: 9 angles x 1 deg at z = 30 mm synthesize %.2f mm = %.1f element lags\n', m_tx*pitch*1e3, m_tx);
cases = {Inf, 0}; Rc = zeros(41, 4); ipt = find(abs(xg + 3e-3) < 1e-9); [~, rpt] = min(abs(zo - 24e-3));
for ic = 1:2
  snr = cases{ic}; if snr == Inf, R = real(Ssig); else, R = real(Ssig + 10^(-snr/20)*Snoi); end
  Cs = cumsum(cat(1, zeros(1, nx, N), R.^2), 1); En = Cs(hi+1,:,:) - Cs(lo,:,:);
  for m = 0:40
    Cs = cumsum(cat(1, zeros(1, nx, N-m), R(:,:,1:N-m).*R(:,:,1+m:N)), 1); num = Cs(hi+1,:,:) - Cs(lo,:,:);
    Rm = mean(num./sqrt(En(:,:,1:N-m).*En(:,:,1+m:N) + eps), 3);
    Rc(m+1, ic) = mean(Rm(bgr));
    if ic == 1, Rc(m+1, 3) = Rm(rpt, ipt); end
    if ic == 2, Rc(m+1, 4) = mean(Rm(cyst)); end
  end
end
vcz = max(0, 1 - (0:40)'/m_tx);
fprintf('speckle R(1), R(4), R(8), R(16): sim %.2f %.2f %.2f %.2f | VCZ %.2f %.2f %.2f %.2f\n', Rc([2 5 9 17],1), vcz([2 5 9 17]));
fprintf('R_sl(M = 8) in speckle: sim %.2f, VCZ %.2f; at 0 dB channel SNR: sim %.2f, predicted x SNR/(1+SNR) = %.2f\n', ...
  sum(Rc(2:9,1)), sum(vcz(2:9)), sum(Rc(2:9,2)), 0.5*sum(Rc(2:9,1)));
fprintf('point (-3, 24) mm: R(1) = %.3f, R_sl(8) = %.2f (ideal 8);  cyst at 0 dB SNR: mean R(1..8) = %.3f\n', Rc(2,3), sum(Rc(2:9,3)), mean(Rc(2:9,4)));

%% 2. Images: DAS and SLSC at the default noise level, sweep of the maximum lag M
R = real(Ssig + 10^(noise_dB/20)*Snoi); Sc = Ssig(ro,:,:) + 10^(noise_dB/20)*Snoi(ro,:,:);
das = abs(sum(Sc,3))/N; Cs = cumsum(cat(1, zeros(1, nx, N), R.^2), 1); En = Cs(hi+1,:,:) - Cs(lo,:,:);
Rlag = zeros(nzo, nx, max(M_list));
for m = 1:max(M_list)
  Cs = cumsum(cat(1, zeros(1, nx, N-m), R(:,:,1:N-m).*R(:,:,1+m:N)), 1); num = Cs(hi+1,:,:) - Cs(lo,:,:);
  Rlag(:,:,m) = mean(num./sqrt(En(:,:,1:N-m).*En(:,:,1+m:N) + eps), 3);
end
gcnr = @(a, b) 1 - sum(min(histc(a, linspace(min([a;b]), max([a;b]), 101))/numel(a), histc(b, linspace(min([a;b]), max([a;b]), 101))/numel(b)));
fw = @(pr) sum(pr >= max(pr)/2)*0.15;                                     % -6 dB width on the 0.15 mm grid (mm)
met = zeros(numel(M_list)+1, 4);
A = das; met(1,:) = [fw(A(rpt, abs(xg + 3e-3) < 2.6e-3)), 10*log10(mean(A(bgr).^2)/mean(A(cyst).^2)), gcnr(A(cyst), A(bgr)), mean(A(bgr))/std(A(bgr),1)];
for k = 1:numel(M_list)
  A = max(sum(Rlag(:,:,1:M_list(k)), 3), 0); ci = max(A(cyst), 1e-6); bi = A(bgr);
  met(k+1,:) = [fw(A(rpt, abs(xg + 3e-3) < 2.6e-3)), 10*log10(mean(bi.^2)/mean(ci.^2)), gcnr(ci, bi), mean(bi)/std(bi,1)];
end
fprintf('\n%-10s %10s %8s %7s %8s\n', 'image', 'point -6dB', 'CR(dB)', 'gCNR', 'sSNR');
fprintf('%-10s %8.2fmm %8.1f %7.3f %8.2f\n', 'DAS', met(1,:));
for k = 1:numel(M_list), fprintf('SLSC M=%-3d %8.2fmm %8.1f %7.3f %8.2f\n', M_list(k), met(k+1,:)); end
fprintf('(point width is quantized to the 0.15 mm pixel grid; SLSC enters the metrics as max(R_sl, 0))\n');

%% Figures
figure('Position', [50 50 1200 380]);
subplot(1,3,1); imagesc(xg*1e3, zo*1e3, 20*log10(das/max(das(:)) + 1e-12), [-60 0]); axis image; colormap(gray); title('DAS (60 dB)'); xlabel('x (mm)'); ylabel('z (mm)');
A = max(sum(Rlag(:,:,1:8), 3), 0); subplot(1,3,2); imagesc(xg*1e3, zo*1e3, A/max(A(:)), [0 1]); axis image; title('SLSC, M = 8 (linear)'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(1,3,3); plot(0:40, Rc(:,1), 'b', 0:40, Rc(:,2), 'b--', 0:40, vcz, 'k:', 0:40, Rc(:,3), 'r', 0:40, Rc(:,4), 'm', 'LineWidth', 1.3);
xlabel('element lag m'); ylabel('R(m)'); title('Spatial coherence'); legend('speckle', 'speckle, 0 dB SNR', 'VCZ triangle', 'point', 'cyst (noise)'); ylim([-0.2 1.05]); grid on;
figure('Position', [50 50 1100 300]);
lab = {'point -6 dB width (mm)', 'CR (dB)', 'gCNR', 'speckle SNR'};
for q = 1:4, subplot(1,4,q); plot(M_list, met(2:end,q), 'ro-', M_list, met(1,q)*ones(size(M_list)), 'k--'); xlabel('maximum lag M'); title(lab{q}); grid on; end

What to look for: \(m_{\mathrm{tx}} = 15.7\) lags; speckle R(1), R(4), R(8) = 0.94, 0.81, 0.52 against 0.94, 0.74, 0.49 (VCZ); \(R_{\mathrm{sl}}(8)\) 6.0 against 5.7; at 0 dB per-angle channel SNR \(R_{\mathrm{sl}}\) falls to 3.4 (compounding nine angles raises the effective SNR slightly above 0 dB); point \(R_{\mathrm{sl}}\) = 8.00; the M sweep narrows the point from 4.95 to 0.75 mm and lowers speckle SNR from 29 to 2.1, with gCNR peaking near M = 4–8.

m25_fxpf.mPlane waves are linearly predictable across the aperture; Monte Carlo of the FXPF gain on noise and on echoes; DAS against FXPF at −5 dB noise · change p, W, mu, noise_dB · runtime ≈ 11 s

%% m25_fxpf.m  -  f-x prediction filtering (FXPF) of delayed channel data
% Purpose : (1) show that q plane waves across the aperture obey an order-q linear
%           recursion (exact prediction for p >= q); (2) Monte Carlo of the FXPF gain on
%           white noise and on a focused echo versus channel SNR (regression dilution);
%           (3) image a cyst/point phantom with DAS and with FXPF (Hann windows of W
%           samples, 50 % overlap, per-bin forward + backward AR(p) least squares, DAS).
% Change  : p, W, mu, noise_dB; add more plane waves in section 1.
% Expect  : section 1 residual ~1e-15 once p >= q; noise-only gain about -14 dB for p = 4
%           (forward-only projection estimate 10 log10(p/(N-p)) = -11.8 dB); a 0 dB echo
%           loses ~2 dB (p = 4) and ~6 dB (p = 1); at -5 dB channel noise FXPF raises CR
%           and gCNR while the point width and speckle SNR stay. Runtime ~15 s in Octave.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(3);
N = 64; c = 1540; f0 = 5e6; lam = c/f0; pitch = 0.3e-3; z = 30e-3; mu = 1e-3;
psi = @(dx) 2*pi*pitch*dx/(lam*z);                      % phase step per element for a wave from lateral offset dx

%% 1. Linear prediction across the aperture: residual of a forward AR(p) fit to q plane waves
fprintf('relative LS residual of a forward AR(p) fit, noise-free, q plane waves (offsets 0, 1.3, -2.1 mm)\n');
offs = [0 1.3e-3 -2.1e-3]; amps = [1 0.7 0.5];
for q = 1:3
  x = zeros(N,1); for k = 1:q, x = x + amps(k)*exp(1j*(psi(offs(k))*(0:N-1)' + k)); end
  r = zeros(1,4);
  for p = 1:4
    A = zeros(N-p, p); for k = 1:p, A(:,k) = x(p+1-k:N-k); end; b = x(p+1:N); a = A\b; r(p) = norm(b - A*a)/norm(b);
  end
  fprintf('  q = %d: p = 1..4 -> %s\n', q, sprintf('%9.1e', r));
end

%% 2. Monte Carlo: FXPF prediction (forward + backward, loading mu) on a constant echo plus white noise
p_list = [1 2 4 8]; snr_list = -20:5:30; T = 200; Tn = 2000;   % trials per echo SNR, trials for noise only (heavy-tailed: needs more)
gs = zeros(numel(p_list), numel(snr_list)); gn = zeros(numel(p_list), 1);
for ip = 1:numel(p_list)
  p = p_list(ip);
  for is = 0:numel(snr_list)                             % is = 0: noise only
    if is == 0, Aamp = 0; else, Aamp = 10^(snr_list(is)/20); end
    acc = 0; pn = 0; pn0 = 0;
    for t = 1:(T + (is == 0)*(Tn - T))
      nz = (randn(N,1) + 1j*randn(N,1))/sqrt(2); x = Aamp + nz; y = zeros(N,1); cnt = zeros(N,1);
      for dirn = [1 -1]
        if dirn > 0, tgt = (p+1:N)'; else, tgt = (1:N-p)'; end
        A = zeros(numel(tgt), p); for k = 1:p, A(:,k) = x(tgt - dirn*k); end
        G = A'*A; a = (G + mu*real(trace(G))/p*eye(p))\(A'*x(tgt)); y(tgt) = y(tgt) + A*a; cnt(tgt) = cnt(tgt) + 1;
      end
      y = y./cnt;
      if is == 0, pn = pn + abs(sum(y))^2; pn0 = pn0 + abs(sum(nz))^2; else, acc = acc + sum(y); end
    end
    if is == 0, gn(ip) = 10*log10(pn/pn0); else, gs(ip, is) = 20*log10(abs(acc/T)/(N*Aamp)); end
  end
end
fprintf('\nnoise-only DAS-output power gain (dB), %d trials, and forward-only estimate 10log10(p/(N-p)):\n', Tn);
for ip = 1:numel(p_list), fprintf('  p = %d: %6.1f dB   (estimate %6.1f dB)\n', p_list(ip), gn(ip), 10*log10(p_list(ip)/(N-p_list(ip)))); end
fprintf('echo amplitude gain (dB) at channel SNR %s dB:\n', sprintf('%4d', snr_list));
for ip = 1:numel(p_list), fprintf('  p = %d: %s\n', p_list(ip), sprintf('%6.1f', gs(ip,:))); end

%% 3. Image: plane-wave phantom (9 angles), DAS vs FXPF at -5 dB channel noise
fs = 20e6; xe = ((0:N-1) - (N-1)/2)*pitch; bw = 0.6; sigf = bw*f0/(2*sqrt(2*log(2))); sigt = 1/(2*pi*sigf);
ang = (-4:4)*pi/180; Na = numel(ang); t0 = 25e-6; Nt = 580; tt = t0 + (0:Nt-1)'/fs; noise_dB = -5; p = 4; W = 16;
nsc = round(13.6*8.6*80); xs = -6.8e-3 + 13.6e-3*rand(nsc,1); zs = 26e-3 + 8.6e-3*rand(nsc,1);
keep = (xs.^2 + (zs-30e-3).^2) > (2.5e-3)^2; xs = [xs(keep); -3e-3]; zs = [zs(keep); 24e-3]; as = [randn(sum(keep),1); 15];
IQ = zeros(Nt, N, Na); K4 = ceil(4*sigt*fs); E = repmat(1:N, numel(xs), 1);
for a = 1:Na
  tau = (xs*sin(ang(a)) + zs*cos(ang(a)))/c + sqrt((xs - xe).^2 + zs.^2)/c; ph = as.*exp(-1j*2*pi*f0*tau); n0 = (tau - t0)*fs;
  accr = zeros(Nt*N,1); acci = accr;
  for k = -K4:K4
    n = floor(n0) + k; g = exp(-((n - n0)/fs).^2/(2*sigt^2)); ok = n >= 0 & n < Nt; idx = n(ok) + 1 + (E(ok)-1)*Nt; v = g(ok).*ph(ok);
    accr = accr + accumarray(idx, real(v), [Nt*N 1]); acci = acci + accumarray(idx, imag(v), [Nt*N 1]);
  end
  IQ(:,:,a) = reshape(accr + 1j*acci, Nt, N);
end
Pref = mean(mean(abs(IQ(tt > 2*27.5e-3/c & tt < 2*33e-3/c, :, 5)).^2)); taps = exp(-((-K4:K4)'/fs).^2/(2*sigt^2));
for a = 1:Na, IQ(:,:,a) = IQ(:,:,a) + sqrt(Pref*10^(noise_dB/10)/sum(taps.^2))*conv2((randn(Nt+2*K4,N) + 1j*randn(Nt+2*K4,N))/sqrt(2), taps, 'valid'); end
dz = 0.075e-3; xg = (-6:0.15:6)*1e-3; zo = (22e-3:dz:34e-3)'; nzo = numel(zo); nx = numel(xg);
[Zp, Xp] = ndgrid(zo, xg); Xp = Xp(:); Zp = Zp(:); Ttx = (Xp*sin(ang) + Zp*cos(ang))/c; Ptx = exp(1j*2*pi*f0*Ttx); S = zeros(numel(Xp), N);
for e = 1:N
  trx = sqrt((Xp - xe(e)).^2 + Zp.^2)/c; acc = 0;
  for a = 1:Na
    fn = (Ttx(:,a) + trx - t0)*fs; n = floor(fn); fr = fn - n; ok = n >= 0 & n < Nt-1; col = IQ(:,e,a); v = zeros(size(fn));
    v(ok) = col(n(ok)+1).*(1-fr(ok)) + col(n(ok)+2).*fr(ok); acc = acc + v.*Ptx(:,a);
  end
  S(:,e) = acc.*exp(1j*2*pi*f0*trx);
end
S = reshape(S, nzo, nx, N);
% FXPF: per window, FFT along depth, per bin vectorized forward/backward LS over all columns, IFFT, overlap-add
tic; win = 0.5 - 0.5*cos(2*pi*(0:W-1)'/W); hop = W/2; Y = zeros(nzo, nx, N);   % periodic Hann: 50 % overlap sums to 1
for s0 = -hop:hop:nzo-1
  rr = s0 + (1:W); ok = rr >= 1 & rr <= nzo; seg = zeros(W, nx, N); seg(ok,:,:) = S(rr(ok),:,:);
  F = fft(seg.*win, [], 1); G = zeros(W, nx, N);
  for kb = 1:W
    Xb = reshape(F(kb,:,:), nx, N); pr = zeros(nx, N); cnt = zeros(1, N);
    for dirn = [1 -1]
      if dirn > 0, tgt = p+1:N; else, tgt = 1:N-p; end
      A = zeros(nx, numel(tgt), p); for q = 1:p, A(:,:,q) = Xb(:, tgt - dirn*q); end
      Gm = zeros(nx, p, p); h = zeros(nx, p);
      for q1 = 1:p, h(:,q1) = sum(conj(A(:,:,q1)).*Xb(:,tgt), 2); for q2 = 1:p, Gm(:,q1,q2) = sum(conj(A(:,:,q1)).*A(:,:,q2), 2); end, end
      tr = zeros(nx,1); for q = 1:p, tr = tr + real(Gm(:,q,q)); end
      for q = 1:p, Gm(:,q,q) = Gm(:,q,q) + mu*tr/p + 1e-30; end
      for k = 1:p, for i = k+1:p, fct = Gm(:,i,k)./Gm(:,k,k); for j = k:p, Gm(:,i,j) = Gm(:,i,j) - fct.*Gm(:,k,j); end, h(:,i) = h(:,i) - fct.*h(:,k); end, end
      cf_ = zeros(nx, p); for k = p:-1:1, sres = h(:,k); for j = k+1:p, sres = sres - Gm(:,k,j).*cf_(:,j); end, cf_(:,k) = sres./Gm(:,k,k); end
      P_ = zeros(nx, numel(tgt)); for q = 1:p, P_ = P_ + cf_(:,q).*A(:,:,q); end
      pr(:, tgt) = pr(:, tgt) + P_; cnt(tgt) = cnt(tgt) + 1;
    end
    G(kb,:,:) = reshape(pr./cnt, [1 nx N]);
  end
  g = ifft(G, [], 1); Y(rr(ok),:,:) = Y(rr(ok),:,:) + g(ok,:,:);
  if s0 == 72, Fchk = F; Gchk = G; end
end
tfx = toc;
imgs = {abs(sum(S,3))/N, abs(sum(Y,3))/N}; names = {'DAS', sprintf('FXPF p=%d W=%d', p, W)};
[ZI, XI] = ndgrid(zo, xg); cyst = (XI.^2 + (ZI-30e-3).^2) < (2e-3)^2; bgr = XI >= -5.8e-3 & XI <= -3.4e-3 & ZI >= 28e-3 & ZI <= 32e-3;
[~, rpt] = min(abs(zo - 24e-3)); sel = abs(xg + 3e-3) < 1.5e-3;
gcnr = @(a, b) 1 - sum(min(histc(a, linspace(min([a;b]), max([a;b]), 101))/numel(a), histc(b, linspace(min([a;b]), max([a;b]), 101))/numel(b)));
fprintf('\nimage at %d dB channel noise (FXPF %.1f s for %d x %d pixels)\n%-14s %8s %7s %7s %7s\n', noise_dB, tfx, nzo, nx, 'image', 'CR(dB)', 'gCNR', 'sSNR', 'pt -6dB');
for i = 1:2
  A = imgs{i}; pr = A(rpt, sel);
  fprintf('%-14s %8.1f %7.3f %7.2f %5.2fmm\n', names{i}, 10*log10(mean(A(bgr).^2)/mean(A(cyst).^2)), gcnr(A(cyst), A(bgr)), mean(A(bgr))/std(A(bgr),1), sum(pr >= max(pr)/2)*0.15);
end
% cross-check: the vectorized solver equals the per-vector solver of section 2 on one aperture vector
x = squeeze(Fchk(3, 40, :)); y = zeros(N,1); cnt = zeros(N,1);
for dirn = [1 -1]
  if dirn > 0, tgt = (p+1:N)'; else, tgt = (1:N-p)'; end
  A = zeros(numel(tgt), p); for k = 1:p, A(:,k) = x(tgt - dirn*k); end
  Gq = A'*A; a = (Gq + mu*real(trace(Gq))/p*eye(p))\(A'*x(tgt)); y(tgt) = y(tgt) + A*a; cnt(tgt) = cnt(tgt) + 1;
end
fprintf('vectorized vs per-vector prediction (window at row 73, bin 3, column 40): max rel. difference %.1e\n', max(abs(y./cnt - squeeze(Gchk(3,40,:))))/max(abs(y)));

%% Figures
figure('Position', [50 50 1150 330]);
subplot(1,3,1); hold on; cols = 'brmk';
for ip = 1:numel(p_list), plot(snr_list, gs(ip,:), [cols(ip) '-o']); plot(snr_list([1 end]), gn(ip)*[1 1], [cols(ip) ':']); end
xlabel('channel SNR (dB)'); ylabel('DAS-output gain (dB)'); title('FXPF: echo (solid) and noise (dotted)'); legend('p=1','','p=2','','p=4','','p=8',''); grid on;
for i = 1:2
  subplot(1,3,i+1); A = imgs{i}/max(imgs{i}(:)); imagesc(xg*1e3, zo*1e3, 20*log10(A + 1e-12), [-60 0]); axis image; colormap(gray);
  title([names{i} ', ' num2str(noise_dB) ' dB noise']); xlabel('x (mm)'); ylabel('z (mm)');
end

What to look for: the forward AR residual drops to \(10^{-15}\) as soon as p ≥ q; noise-only gains −20.8, −17.6, −14.4, −10.1 dB for p = 1, 2, 4, 8 (2000 trials; estimate \(10\log_{10}(p/(N - p))\) = −18.0 … −8.5 dB); a 0 dB echo loses 5.8 dB at p = 1 and 1.7 dB at p = 4; in the image FXPF lifts CR from 19.1 to 26.7 dB and gCNR from 0.931 to 0.965 with the point width (0.45 mm) and speckle SNR (≈ 1.9–2.0) unchanged; the vectorized and per-vector solvers agree to \(10^{-15}\).

m25_admire.mADMIRE on one line: STFT windows, ROI and off-axis wavefront models, elastic net (FISTA), keep the ROI part; target kept, off-axis clutter removed · change clutAmp, snrDB, cLam, useOuter · runtime ≈ 20–45 s

%% m25_admire.m  --  ADMIRE (aperture domain model image reconstruction) on one image line
% Purpose : a compact version of the ADMIRE pipeline of Module 25 (Byram et al., IEEE TUFFC 2015;
%           model space and regularization from the authors' released code, VU-BEAM-Lab/ADMIRE):
%           delayed aperture data -> short-time FFT along depth -> per window and frequency, fit the
%           aperture vector with modelled wavefronts of scatterers inside (ROI) and outside (off-axis)
%           the resolution cell (elastic net) -> keep the ROI part -> inverse FFT -> DAS.
% Scene   : one line at x = 0 through a point target at (0, 30) mm, 120 off-axis scatterers at
%           1.5 <= |x| <= 9 mm (z = 28-32 mm) and white channel noise; 64 elements, 0.3 mm pitch,
%           5 MHz, 60 % bandwidth, one 0 deg plane wave, dynamic receive focusing.
% Differences from the page: one plane wave instead of nine, the reduced model space of the page's
%           live button (about 900 predictors), and the elastic net solved by FISTA (vectorized)
%           instead of coordinate descent; the objective, alpha and lambda rule are the same.
% Change  : clutAmp (clutter strength), snrDB, cLam (lambda scale), useOuter (0 = no off-axis model).
% Expected: printed target loss (ADMIRE vs DAS) and clutter level relative to the target, both
%           methods; figure with the depth profiles and the aperture fit at the target window.
% Runtime : about 10-20 s in Octave 8.
rng(1);
c = 1.54; f0 = 5; BW = 0.6; N = 64; pitch = 0.3;            % mm/us, MHz, -, -, mm
xe = ((0:N-1) - (N-1)/2) * pitch;                           % element x (mm)
sigf = BW*f0/(2*sqrt(2*log(2))); sigt = 1/(2*pi*sigf);      % pulse: Gaussian envelope (us)
dz = 0.075; z = (27.5:dz:32.5)'; nz = numel(z);             % image rows on the line x = 0 (mm)
clutAmp = 0.25; nClut = 120; snrDB = 20; cLam = 0.0189; alpha = 0.9; useOuter = 1;

%% Delayed aperture data on the line x = 0: s(n,e) = sum_s a_s g(dt) exp(j 2 pi f0 dt), dt = tau_pix - tau_s
xc = (1.5 + 7.5*rand(nClut,1)) .* sign(rand(nClut,1) - 0.5); zcl = 28 + 4*rand(nClut,1);
ac = clutAmp * (randn(nClut,1) + 1i*randn(nClut,1)) / sqrt(2);
tauPix = (z + sqrt(xe.^2 + z.^2)) / c;                      % nz x N, 0 deg transmit + receive
mk = @(xs, zs, as) reshape(sum(reshape(as,1,1,[]) .* exp(-(tauPix - reshape((zs + sqrt((xs - xe).^2 + zs.^2))', 1, N, []) / c).^2 ...
     / (2*sigt^2)) .* exp(1i*2*pi*f0*(tauPix - reshape((zs + sqrt((xs - xe).^2 + zs.^2))', 1, N, []) / c)), 3), nz, N);
sT = mk(0, 30, 1);                                          % target only
sC = zeros(nz, N); for k = 1:30:nClut, idx = k:min(k+29, nClut); sC = sC + mk(xc(idx), zcl(idx), ac(idx)); end
sigP = mean(abs(sT(abs(z-30) < 0.15, :)).^2, 1); nsd = sqrt(mean(sigP) * 10^(-snrDB/10) / 2);
sN = nsd * (randn(nz, N) + 1i*randn(nz, N));               % white noise, SNR re target peak per channel
scenes = {sT + sC + sN, sC + sN, sT};                       % A: target + clutter + noise, B: clutter + noise, C: target

%% STFT and model parameters (released rule: window = ceil(FWHM * fs), zero-pad x2, bins in f0 +- 0.6 BW f0)
fsa = c/(2*dz); Lw = ceil(sqrt(8*log(2))*sigt*fsa); Lp = 2*Lw;
fk = (0:Lp-1)*fsa/Lp; bins = find(fk >= f0*(1-0.6*BW) & fk <= f0*(1+0.6*BW));
kcal = [1.10 0.95 0.88];                                    % per-bin wavenumber calibration (page, fitted to this pulse)
h = Lw/fsa*c/2; st = c*sigt; nw = floor(nz/Lw);
gridv = @(a, s, b) a + (0:floor((b - a)/s + 1e-9))*s;      % MATLAB a:s:b without round-off surprises
outD = cell(3,1); fitShow = [];
tic;
for w = 1:nw
  rows = (w-1)*Lw + (1:Lw); zf = z(rows(1)); zc = mean(z(rows)); z2a = 2*zf; z2b = 2*z(rows(end));
  S = cell(3,1); for q = 1:3, X = fft(scenes{q}(rows, :), Lp, 1); S{q} = X; end
  R = cell(3,1); for q = 1:3, R{q} = zeros(Lp, N); end
  for b = 1:numel(bins)
    f = fk(bins(b)); kw = 2*pi*f*kcal(b)/c; lam = 2*pi/kw;
    resLat = 1.2067/(N*pitch)*lam*zc; resAx = 2*resLat;    % FWHM of |sinc| aperture transform (release computes it by FFT)
    % predictor grid: ROI and outer (off-axis) scatterers, offsets d (reduced sampling of the page's live model)
    [xr, zr, dr] = ndgrid(gridv(-0.5, 0.1432, 0.5)*resLat, zc + gridv(-0.5, 0.572, 0.5)*resAx, gridv(-0.8, 0.388*lam, 0.4));
    L = pitch*N/2 + 1; [xq, zq, dq] = ndgrid(gridv(-L, 1.4228*resLat, L), zc + gridv(-0.5, 0.7114, 0.5)*resAx, gridv(-8, 0.9688*lam, 3.2));
    keepO = (xq(:)/(resLat + 0.5)).^2 + ((zq(:) - zc)/(resAx + 0.5)).^2 > 1;
    xp = [xr(:); xq(keepO)]; zp = [zr(:); zq(keepO)]; dp = [dr(:); dq(keepO)]; nRoi = numel(xr);
    if ~useOuter, xp = xp(1:nRoi); zp = zp(1:nRoi); dp = dp(1:nRoi); end
    d0 = sqrt((xe - xp).^2 + zp.^2) + dp + 2*zc - zp;     % P x N total path
    zd = (d0.^2 - xe.^2) ./ d0;                             % two-way position after dynamic receive focusing
    pmin = max(zd - h, z2a); pmax = min(zd + h, z2b);
    g = erf((zd - pmin)/st) - erf((zd - pmax)/st); g(zd - h > z2b | zd + h < z2a) = 0;
    ok = max(g, [], 2) > 0; g = g(ok, :); zd = zd(ok, :); isRoi = (1:numel(ok))' <= nRoi; isRoi = isRoi(ok);
    A = sqrt(g ./ max(g, [], 2)) .* exp(-1i*kw*zd);        % modelled aperture signals (rows)
    A = (A ./ sqrt(sum(abs(A).^2, 2))).';                   % N x P, unit-norm columns
    Xr = [real(A) -imag(A); imag(A) real(A)]; P = size(A, 2); No = 2*N;
    Lip = norm(Xr)^2/No;                                    % Lipschitz constant of the smooth part
    for q = 1:3
      y = S{q}(bins(b), :).'; yv = [real(y); imag(y)]; sd = std(yv, 1);
      if sd == 0, continue; end
      ys = yv/sd; lm = cLam*sqrt(mean(yv.^2))/sd; tstep = 1/(Lip + lm*(1-alpha));
      B = zeros(2*P, 1); V = B; tk = 1;                     % FISTA on (1/2No)||y - Xb||^2 + lm(alpha|b|_1 + (1-alpha)/2 |b|^2)
      for it = 1:400
        G = V - tstep*(Xr'*(Xr*V - ys)/No + lm*(1-alpha)*V);
        Bn = sign(G).*max(abs(G) - tstep*lm*alpha, 0);
        tn = (1 + sqrt(1 + 4*tk^2))/2; V = Bn + (tk - 1)/tn*(Bn - B);
        if norm(Bn - B) < 1e-7*max(1, norm(B)), B = Bn; break; end
        B = Bn; tk = tn;
      end
      B = B*sd; bc = B(1:P) + 1i*B(P+1:end);
      R{q}(bins(b), :) = (A(:, isRoi) * bc(isRoi)).';      % keep only the in-ROI part
      if q == 1 && abs(zc - 30) < Lw*dz/2 && b == 2, fitShow = struct('y', y, 'roi', A(:, isRoi)*bc(isRoi), 'out', A(:, ~isRoi)*bc(~isRoi), 'P', P, 'nRoi', sum(isRoi)); end
    end
  end
  for q = 1:3, r = ifft(R{q}, [], 1); outD{q}(rows, 1) = abs(sum(r(1:Lw, :), 2))/N; end
end
tRun = toc;
das = cellfun(@(s) abs(sum(s, 2))/N, scenes, 'UniformOutput', false);
nr = nw*Lw; zz = z(1:nr);

%% Numbers: target kept? clutter removed?
iT = abs(zz - 30) < 0.04;                                   % rows at the target depth
iC = abs(zz - 30) > 0.6;                                    % rows away from the target (clutter only)
pk = @(v) max(v(iT));
dasT = pk(das{3}(1:nr)); admT = pk(outD{3});               % target alone
fprintf('Window %d samples (%.2f mm), bins %s MHz, last fit %d predictors (%d ROI), %.1f s\n', Lw, Lw*dz, mat2str(round(fk(bins)*100)/100), fitShow.P, fitShow.nRoi, tRun);
fprintf('Target only: ADMIRE/DAS peak = %.2f dB (amplitude loss from the elastic-net shrinkage and the 3-bin band)\n', 20*log10(admT/dasT));
fprintf('Target + clutter + noise at the target: DAS %.2f dB, ADMIRE %.2f dB re its own target-only peak\n', 20*log10(pk(das{1}(1:nr))/dasT), 20*log10(pk(outD{1})/admT));
cD = 10*log10(mean(das{2}(iC).^2)/dasT^2); pA = mean(outD{2}(iC).^2);
if pA > 0, sA = sprintf('%.1f dB', 10*log10(pA/admT^2)); else, sA = 'exactly 0 (no ROI coefficient survives)'; end
fprintf('Clutter + noise only, rows > 0.6 mm from the target: DAS %.1f dB re target peak, ADMIRE %s\n', cD, sA);
fprintf('Clutter-only rows set exactly to zero by ADMIRE: %.0f %%; in the target window it keeps %.0f %% of the rows nonzero\n', 100*mean(outD{2}(iC) == 0), 100*mean(outD{1}(iT) > 0));

%% Figure
figure(1); clf;
subplot(1,2,1);
plot(zz, 20*log10(das{1}(1:nr)/dasT + 1e-6), 'k', zz, 20*log10(outD{1}/admT + 1e-6), 'r', ...
     zz, 20*log10(das{2}(1:nr)/dasT + 1e-6), 'k:', zz, 20*log10(outD{2}/admT + 1e-6), 'r:', 'LineWidth', 1.2);
ylim([-70 5]); xlabel('depth z (mm)'); ylabel('amplitude (dB re target-only peak)');
legend('DAS, target+clutter', 'ADMIRE, target+clutter', 'DAS, clutter only', 'ADMIRE, clutter only', 'Location', 'southwest');
title('Line x = 0, depth profiles');
subplot(1,2,2);
plot(xe, real(fitShow.y), 'k.', xe, real(fitShow.roi), 'r', xe, real(fitShow.out), 'b', 'LineWidth', 1.2);
xlabel('element x (mm)'); ylabel('Re Y(f) at 5.13 MHz (a.u.)');
legend('data', 'ROI part (kept)', 'outer part (removed)', 'Location', 'best');
title('Aperture fit, target window');

What to look for (Octave 8): 4-sample windows (0.30 mm), bins 3.85, 5.13, 6.42 MHz, about 810 predictors per fit (≈ 106 ROI). Alone, the target comes out 3.3 dB below DAS (shrinkage and the 3-bin band; the display normalizes it away). With 120 off-axis scatterers and noise added, the target pixel changes by −0.3 dB for ADMIRE and −0.01 dB for DAS, while away from the target DAS shows clutter at −32.9 dB re the target and ADMIRE returns exactly zero on every clutter-only row: no ROI coefficient survives. Set useOuter = 0 and the clutter returns, fitted by ROI columns. The script uses one plane wave, the reduced model space and FISTA instead of coordinate descent; the objective, α and the λ rule are the page's.

25.13Self-check

Q1A point echo sits in white noise at 0 dB channel SNR, N = 64. What are CF and the Wiener gain H?

With signal power 1 and noise power \(r^{-1} = 1\) per channel, \(\mathrm{E}\lvert\sum s\rvert^2 \approx N^2 + N\) and \(\mathrm{E}\sum\lvert s\rvert^2 = 2N\), so \(\mathrm{CF} = (N + 1)/(2N) = 0.508\). \(H = N\cdot\mathrm{CF}/((N - 1)\mathrm{CF} + 1)\) = 32.5/33.0 = 0.985. CF costs 5.9 dB of output amplitude, H 0.1 dB.

Q2At z = 30 mm, 5 MHz, 0.3 mm pitch and N = 64, which aperture bin does an echo from 3 mm off-axis occupy? Does GCF with \(M_0 = 2\) reject it?

\(k = Np\,\Delta x/(\lambda z)\) = 64 × 0.3 × 3/(0.308 × 30) = 6.2. Yes: it lies outside \(\lvert k\rvert \le 2\) (most of its energy is in bins 6 and 7).

Q3Fifteen plane waves at 0.5° spacing image a speckle region at 40 mm with a 0.3 mm pitch. Where does speckle coherence reach zero, and what is \(R_{\mathrm{sl}}\) for M = 8?

\(m_{\mathrm{tx}} = N_\theta z\,\Delta(\sin\theta)/p\) = 15 × 40 × 0.00873/0.3 = 17.5 lags. \(R_{\mathrm{sl}} = \sum_{m=1}^{8}(1 - m/17.5)\) = 8 − 36/17.5 = 5.9.

Q4How many complex multiply-adds per pixel does MV need for N = 64, L = 16, K = ±1 (covariance plus a Cholesky solve)?

\((2K + 1)(L(N - L + 1) + L^2) + L^3/6\) = 3 × (16 × 49 + 256) + 683 = 3120 + 683 ≈ 3800, about 60× DAS's 64. With L = 32, K = 2 it is about 15 900.

Q5A phase screen adds independent delays of 30 ns rms at 5 MHz. Estimate the point-brightness loss for DAS and for CF-weighted DAS (N = 64).

\(\sigma_\varphi = 2\pi \times 5\ \mathrm{MHz} \times 30\ \mathrm{ns} = 0.94\) rad. \(\mathrm{E}\lvert\bar{y}\rvert^2 = (1 + 63e^{-0.888})/64 = 0.42\), so DAS loses 3.8 dB (\(-4.34\,\sigma_\varphi^2 = -3.9\) dB in the large-N limit). \(\mathrm{CF}\cdot\lvert\bar{y}\rvert = \lvert\bar{y}\rvert^3\) loses three times as much: about 11.3 dB.

Q6In a DAS image the tissue ROI averages −40 dB and the cyst −65 dB (unclipped), and gCNR = 0.958. You display DAS² instead. What are the dB contrast and gCNR, and what does that say about weighted beamformers?

\(A^2\) doubles every dB value: −80 and −130 dB, a 50 dB contrast instead of 25. Pixel ordering is unchanged, so gCNR = 0.958. (A power-ratio CR scales only approximately: 24.7 → 46.4 dB on the lab image.) Any method whose CR gain over DAS can be reproduced by a power law has shown a dynamic-range change, not better separability.

Q7The phases of a pixel's aperture vector have a standard deviation of 0.9 rad. What is PCF for γ = 1 and γ = 1.5?

\(\sigma_0 = \pi/\sqrt{3} = 1.814\) rad. γ = 1: 1 − 0.9/1.814 = 0.50. γ = 1.5: 1 − 1.35/1.814 = 0.26. In dB on the amplitude image, −6.0 and −11.8 dB.

Q8In Demo 25.6 the channel noise is set to +30 dB re speckle. Estimate the image speckle-to-noise ratio for DAS from the readout at 0 dB (+13.9 dB), and predict what the background ROI looks like.

DAS is linear, so with no phase screen the image SNR follows the channel SNR one for one: \(13.9 - 30 = -16.1\) dB, which the readout confirms. Noise then carries about 40 times the speckle power in the ROI, so the background is essentially beamformed noise: the cyst almost filled (DAS CR 0.2 dB, gCNR 0.21) and, because this plane-wave dataset forms every pixel from the same recordings, speckle-like grain at the PSF scale with a near-Rayleigh envelope (DAS speckle SNR 1.69 in the table), not salt and pepper. Two differences remain visible. In a line-by-line focused scan the same noise would be independent from line to line, fine laterally and streaky along depth (selector, line-by-line). And unlike speckle it changes with every frame: press "New noise frame" and the whole background re-draws, which a stationary tissue pattern would not do.

Q9Demo 25.6's image has 81 lines of 161 rows. With 4-row STFT windows, 3 bins and the released ADMIRE model space (≈ 31 000 complex predictors per fit), how many elastic-net fits does one image need, how many real unknowns and equations does each have, and why can the answer still be unique?

Each line is tiled by ⌈161/4⌉ = 41 windows (we compute 164 rows), so 81 × 41 × 3 = 9963 fits per image (plus 416 profile lines × 3 bins for the resolution metrics). Each fit has 2 × 31 000 ≈ 62 000 real unknowns and 2 × 64 = 128 real equations, about 480 times underdetermined. Least squares alone would have infinitely many exact solutions; the elastic net's L2 term \(\lambda(1-\alpha)\lVert\mathbf{\beta}\rVert^2/2\) makes the cost strictly convex, so the minimizer is unique, and the L1 term makes it sparse (tens to a few hundred nonzero coefficients). In our offline run a single image with the full model costs about 10 minutes of CPU time in Node.js (8 min for the 81 image lines, 1.3 min for the 416 profile lines; 3.8 CPU-hours for the 24 settings of Demo 25.5c), which is why Demo 25.6 shows it precomputed.

25.14Figures from the literature

Published figures with measured or carefully simulated data for the methods of this module. Links only (no reproduced images); figure numbers checked against the cited version.

SLSC against B-mode, in a phantom and in a human thyroid

M. A. Lediju, G. E. Trahey, B. C. Byram and J. J. Dahl, "Short-lag spatial coherence of backscattered echoes: imaging characteristics," IEEE TUFFC 58(7), 1377–1388, 2011. Fig. 1 (coherence curves: point, speckle, in vivo thyroid), Fig. 7 (4 mm anechoic lesions, SLSC at three M), Fig. 10 (thyroid cyst in vivo). doi:10.1109/TUFFC.2011.1957

Fig. 1 is the measured counterpart of Demo 25.1's coherence panel: R(m) near 1 for a point, a VCZ-like decay for speckle, and the in vivo curve falling faster (clutter and aberration). Fig. 10 shows a thyroid cyst at 1.5 cm in B-mode, in SLSC at increasing Q (their M as a percentage of N) and in spatially compounded B-mode: the cyst clears in SLSC, as in our Demo 25.6, and points broaden.

SLSC point width and brightness versus M, noise and clutter

M. A. Lediju Bell, J. J. Dahl and G. E. Trahey, "Resolution and brightness characteristics of short-lag spatial coherence (SLSC) images," IEEE TUFFC 62(7), 1265–1276, 2015. Fig. 7 (lateral point width versus M, and texture size), Fig. 9 (phantom B-mode and SLSC with and without added noise). doi:10.1109/TUFFC.2014.006909

Fig. 7(a) plots lateral point-target width against the maximum lag M for SLSC and against the number of receive elements for B-mode: the trend our M sweep reports (5.3 → 0.63 mm for M = 1 → 24). Fig. 9 shows experimental phantom images in B-mode and SLSC with and without added channel noise (point targets boxed), and Fig. 10 tracks tissue and point brightness against the noise-to-signal ratio: the per-pair normalization at work.

DMAS beam patterns and an in vivo carotid

G. Matrone, A. S. Savoia, G. Caliano and G. Magenes, "The delay multiply and sum beamforming algorithm in ultrasound B-mode medical imaging," IEEE Trans. Med. Imaging 34(4), 940–949, 2015. Fig. 3 (spectra before filtering, axial and lateral profiles), Fig. 6 (simulated cyst phantom), Fig. 7 (in vivo carotid). doi:10.1109/TMI.2014.2371235

Fig. 3(a) shows the DMAS spectrum before filtering with its baseband and \(2f_0\) parts (our rule card on the DC term); Fig. 3(c) compares lateral two-way beam patterns of DAS, F-DMAS and a triangle-apodized DAS, the sum-coarray picture of (25.10). Fig. 7 is a carotid artery imaged by DAS and F-DMAS over 70 dB, and Fig. 8 gives the two pixel histograms; compare them, then ask how much of the difference a power law on DAS would reproduce.

CNR moves with the display, gCNR does not: eight beamformers compared

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 (open access, CC BY 4.0). Fig. 1 (S-curve compression changes CNR, not gCNR), Figs. 7–9 (DAS, S-DAS, CF, PCF, GCF, DMAS and SLSC in silico, in vitro and in a carotid). doi:10.1109/TUFFC.2019.2956855

Fig. 1 is (25.15) on one image: an S-shaped grey map raises CNR from 0.89 to 2.29 while gCNR stays at 96.9 %. Figs. 7–9 put most of this module's coherence beamformers side by side: CNR scatters widely across methods while gCNR moves little, and S-DAS (pure dynamic-range manipulation) matches DAS in gCNR exactly. The dynamic-range companion paper is Rindal et al., IEEE TUFFC 66(7), 1198–1208, 2019 (doi:10.1109/TUFFC.2019.2911267).

25.15Go deeper

J. Capon, "High-resolution frequency-wavenumber spectrum analysis," Proc. IEEE 57(8), 1969. The MV (Capon) beamformer.
J.-F. Synnevåg, A. Austeng and S. Holm, "Adaptive beamforming applied to medical ultrasound imaging," IEEE TUFFC 54(8), 1606–1613, 2007 (doi); and "Benefits of minimum-variance beamforming in medical ultrasound imaging," IEEE TUFFC 56(9), 1868–1879, 2009 (doi). Subarray smoothing, temporal averaging, loading and speckle statistics.
B. M. Asl and A. Mahloojifar, "Eigenspace-based minimum variance beamforming applied to medical ultrasound imaging," IEEE TUFFC 57(11), 2010.
B. M. Asl and A. Mahloojifar, "Contrast enhancement and robustness improvement of adaptive ultrasound imaging using forward-backward minimum-variance beamforming," IEEE TUFFC 58(4), 858–867, 2011. doi:10.1109/TUFFC.2011.1880.
T.-J. Shan, M. Wax and T. Kailath, "On spatial smoothing for direction-of-arrival estimation of coherent signals," IEEE Trans. ASSP 33(4), 1985.
B. Widrow, K. M. Duvall, R. P. Gooch and W. C. Newman, "Signal cancellation phenomena in adaptive antennas: causes and cures," IEEE Trans. Antennas Propag. 30(3), 1982.
H. Cox, R. M. Zeskind and M. M. Owen, "Robust adaptive beamforming," IEEE Trans. ASSP 35(10), 1987. Diagonal loading.
C.-I. C. Nilsen and S. Holm, "Wiener beamforming and the coherence factor in ultrasound imaging," IEEE TUFFC 57(6), 2010.
R. Mallart and M. Fink, "Adaptive focusing in scattering media through sound-speed inhomogeneities: the van Cittert Zernike approach and focusing criterion," J. Acoust. Soc. Am. 96(6), 1994.
K. W. Hollman, K. W. Rigby and M. O'Donnell, "Coherence factor of speckle from a multi-row probe," Proc. IEEE Ultrasonics Symposium, 1999.
P.-C. Li and M.-L. Li, "Adaptive imaging using the generalized coherence factor," IEEE TUFFC 50(2), 128–141, 2003. doi:10.1109/TUFFC.2003.1182117.
J. Camacho, M. Parrilla and C. Fritsch, "Phase coherence imaging," IEEE TUFFC 56(5), 2009. PCF and SCF.
G. Matrone, A. S. Savoia, G. Caliano and G. Magenes, "The delay multiply and sum beamforming algorithm in ultrasound B-mode medical imaging," IEEE Trans. Med. Imaging 34(4), 940–949, 2015. doi:10.1109/TMI.2014.2371235.
M. A. Lediju, G. E. Trahey, B. C. Byram and J. J. Dahl, "Short-lag spatial coherence of backscattered echoes: imaging characteristics," IEEE TUFFC 58(7), 1377–1388, 2011 (doi); M. A. Lediju Bell, J. J. Dahl and G. E. Trahey, "Resolution and brightness characteristics of short-lag spatial coherence (SLSC) images," IEEE TUFFC 62(7), 1265–1276, 2015 (doi).
J. Shin and L. Huang, "Frequency-space prediction filtering for acoustic clutter and random noise attenuation in ultrasound imaging," Proc. SPIE 9790, 97901O, 2016 (doi); J. Shin, L. Huang and J. T. Yen, "Spatial prediction filtering for medical ultrasound in aberration and random noise," IEEE TUFFC 65(10), 1845–1856, 2018 (doi). R. L. Canales, "Random noise reduction," SEG Expanded Abstracts, 1984: the seismic f-x origin.
B. Byram, K. Dei, J. Tierney and D. Dumont, "A model and regularization scheme for ultrasonic beamforming clutter reduction," IEEE TUFFC 62(11), 1913–1927, 2015. doi:10.1109/TUFFC.2015.007004. ADMIRE.
K. Dei and B. Byram, "The impact of model-based clutter suppression on cluttered, aberrated wavefronts," IEEE TUFFC 64(10), 1450–1464, 2017. doi:10.1109/TUFFC.2017.2729944.
C. Khan, K. Dei, S. Schlunk, K. Ozgun and B. Byram, "A real-time, GPU-based implementation of aperture domain model image reconstruction," IEEE TUFFC 68(6), 2101–2116, 2021. doi:10.1109/TUFFC.2021.3056334. Code: github.com/VU-BEAM-Lab/ADMIRE (Apache-2.0), the source of this module's ADMIRE parameters.
S. Schlunk, K. Dei and B. Byram, "Iterative model-based beamforming for high dynamic range applications," IEEE TUFFC 68(3), 482–493, 2021. doi:10.1109/TUFFC.2020.3012165. iADMIRE; ADMIRE on compounded plane-wave data.
O. M. H. Rindal, A. Austeng, A. Fatemi and A. Rodriguez-Molares, "The effect of dynamic range alterations in the estimation of contrast," IEEE TUFFC 66(7), 1198–1208, 2019. doi:10.1109/TUFFC.2019.2911267.
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. doi:10.1109/TUFFC.2019.2956855.
H. Liebgott, A. Rodriguez-Molares, F. Cervenansky, J. A. Jensen and O. Bernard, "Plane-wave imaging challenge in medical ultrasound," Proc. IEEE IUS, 2016 (PICMUS).

← 24 Fourier beamforming 26 Speckle statistics →