Part V · Module 32

Sound-speed estimation (capstone)

Every beamformer assumes a sound speed, and every quantitative ultrasound number inherits the error of that assumption. Estimating the speed from the echoes is a delay-estimation problem (Module 11) sitting inside an inverse problem (Module 13), solved with the beamformer itself as the forward model (Modules 07, 14, 25). This capstone builds the estimators from the global focus-optimal speed, through the average-to-local inversion of layered media and echo-mode tomography (CUTE), to differentiable beamforming, and closes with the error budget that decides whether a liver measurement can separate healthy from steatotic tissue.

You will be able to

Compute the focusing tolerance on sound speed for a given F-number, depth and frequency, and say what a 1 % error does to phase, depth and coherence
Explain why a focus-quality maximizer returns an average speed (the RMS speed of seismology, not the harmonic mean), and invert averages to local speeds with Dix's formula, straight-ray differencing or regularized least squares, with their noise amplification
Derive the CUTE measurement (echo phase shifts between steered acquisitions), its wrapping limit, the phase jitter of correlated speckle, and the difference between transmit-only and common-mid-angle data
Set up limited-angle sound-speed tomography as \(\mathbf{d} = \mathbf{A}\,\Delta\mathbf{s} + \mathbf{n}\), read its resolution matrix, and separate noise, regularization and forward-model errors
Recognize differentiable-beamforming autofocus as nonlinear tomography, predict where its losses have local minima, and place deep-learning and full-wave methods on the same map
Write an error budget (CRLB, calibration, layer bias, model error) for a clinical liver sound-speed measurement and compare it with the 10–40 m/s effects that matter

Assumes

Module 07 · Random vectors, covariance & correlation: coherence and the van Cittert–Zernike theorem; Module 11 · Correlation & time-delay estimation: delay estimation and its jitter
Module 12 · Estimation theory and the Cramér–Rao bound; Module 13 · Inverse problems and regularization: Tikhonov, resolution matrices and Section 13.8 on limited-angle SoS tomography; Module 14 · Optimization & automatic differentiation: autofocus (Section 14.5) and cycle skipping (Section 14.4)
Module 25 · Adaptive and coherence-based beamforming: coherence factor and SLSC; Module 23 · Beam patterns, k-space and compounding: plane-wave compounding

Exact formulaclosed form plotted as is, or published numbers re-plottedSimulationsolvers that run here in JS: ray-traced channel delays, linear inversions, Monte CarloApprox. modelsimplified physics; assumptions under the badgeSchematicdrawing for intuition, not to scaleDefaults: \(c_0\) = 1540 m/s assumed by the beamformer, \(f_0\) = 5 MHz, pitch 0.3 mm, 128-element linear array (38.4 mm).

32.1Why sound speed matters: focus, geometry and a biomarker

Soft tissues span roughly 1450 m/s (fat) to 1600 m/s (muscle, healthy liver), a −6 % to +4 % band around the 1540 m/s that every scanner assumes. A wrong speed does three separate things. It defocuses: the receive delays have the wrong curvature, so the aperture sums out of phase. It misplaces: depth is computed as \(\hat c\,t/2\). And it hides a biomarker: fat lowers the sound speed of liver, so the speed itself is a quantitative measure of steatosis, and every downstream QUS estimate (attenuation, backscatter, Module 27) is computed from images formed with the assumed value.

Start with the defocus. Place the pixel by its two-way time, as every scanner does, so the on-axis delay is always right. A receive element at lateral offset \(x\) then sees a residual delay set by the error in the wavefront curvature. With \(\hat c = c + \Delta c\), depth \(\hat z = \hat c t_a\) (\(t_a\) the one-way apex time) and \(x \ll z\):

\[ \delta\tau(x) = \frac{\sqrt{z^{2}+x^{2}}-z}{c} - \frac{\sqrt{\hat z^{2}+x^{2}}-\hat z}{\hat c} \;\approx\; \frac{\Delta c}{c}\,\frac{x^{2}}{z\,c}, \qquad \left\lvert\frac{\Delta c}{c}\right\rvert < \frac{\lambda F^{2}}{2z} \;\Longleftrightarrow\; \lvert\Delta c\rvert < \frac{c^{2}F^{2}}{2 f_0 z} \]

(32.1)

\(\delta\tau(x)\): receive delay error of the element at offset \(x\) after the pixel is placed at the correct two-way time, s
\(c,\ \hat c,\ \Delta c\): true, assumed and error of the sound speed, m/s
\(z,\ \hat z\): true and displayed depth of the pixel, m (\(\hat z/z = \hat c/c\))
\(F\): receive F-number \(z/D\), with \(D = 2X\) the active aperture
\(\lambda,\ f_0\): wavelength \(c/f_0\) and centre frequency

The tolerance is the λ/8 (π/4 phase) criterion at the aperture edge \(x = X\). At 5 MHz and 30 mm it is 7.9 m/s at F/1, 32 m/s at F/2 and 71 m/s at F/3: a large aperture is both what needs the right speed and what can measure it. The curvature error is twice that of a pixel held at fixed depth, because the time-tied pixel moves by \(\Delta c/c\) as well.

The geometry and the parameterization follow from the same travel times. Displayed depth scales as \(\hat z = (\hat c/c)\,z\): 1 % of speed is 0.3 mm at 30 mm, larger than the axial resolution at 5 MHz. And every travel time is linear in slowness \(s = 1/c\), not in speed, which is why all the inversions below are written for the slowness deviation:

\[ \Delta s = \frac{1}{c} - \frac{1}{c_0}, \qquad \Delta c \approx -c_0^{2}\,\Delta s, \qquad 1\ \mathrm{ns/mm} \;\leftrightarrow\; 2.37\ \mathrm{m/s}\ \text{at}\ c_0 = 1540\ \mathrm{m/s} \]

(32.2)

\(\Delta s\): slowness deviation from the beamforming slowness, ns/mm (10⁻⁶ s/m)
\(c_0\): sound speed assumed by the beamformer
\(c\): local sound speed of the medium, m/s

A +40 m/s liver region is −16.4 ns/mm; 10 mm of it advances a one-way path by 164 ns, 0.8 of a period at 5 MHz (329 ns for the round trip of an echo from below it). The linearization error is \(\Delta c^{2}/c\): 1.0 m/s at 40 m/s, 2.5 % of the correction.

The biomarker. Imbault et al. (Phys. Med. Biol., 2017) estimated local liver speed in 17 patients and found 1492–1604 m/s across MRI proton-density fat fractions of 2–17 %, with an area under the ROC curve of 0.952 against biopsy and cut-offs of 1541 m/s (PDFF 5 %) and 1555 m/s (biopsy 10 %). Those numbers set the scale for everything in this module: a useful estimator must resolve about 10 m/s in a liver that sits under 10–30 mm of fat and muscle at 1450–1580 m/s. The overlying layers bias any average by tens of m/s (Section 32.3), so the hard part of liver sound speed is not precision, it is locality.

Family	Measures	Output	Key assumption	Where in this module
Focus-quality maximization (CF, coherence, brightness, sharpness)	focus metric versus assumed \(\hat c\)	one average speed per region	homogeneous medium above the region	32.2, Demo 32.1
Wavefront curvature (Anderson and Trahey, 1998)	arrival-time profile across the aperture	average speed to a depth	hyperbolic moveout	32.2, 32.7
Average-to-local inversion (Jakovljevic et al., 2018)	averages at many depths	local speed per layer or depth	layered medium, straight vertical rays	32.3, Demo 32.2
Echo-mode tomography: CUTE and spatial-domain methods	echo phase shifts between steered acquisitions	2D map of \(\Delta s\)	straight rays, linearized, limited angle	32.4–32.5, Demos 32.3–32.4
Differentiable beamforming autofocus	gradient of a focus or phase loss	2D map	differentiable forward model (usually straight rays)	32.6, Demo 32.5
Learned and full-wave inversion	channel data to map, or wave-equation migration velocity analysis	2D map	training distribution, or the wave equation	32.6

32.2Global sound speed by focus-quality maximization, and which average it returns

The oldest estimators sweep \(\hat c\), beamform, and keep the speed that maximizes a focus metric. Module 07 showed why that works: the van Cittert–Zernike theorem makes the aperture-domain coherence of speckle peak when the delays are right, and Module 14 turned the sweep into gradient ascent. The usual metrics, all evaluated on delayed channel signals \(s_i\) of a pixel or a region tied to two-way time:

\[ \mathrm{CF}(\hat c) = \frac{\left\lvert \sum_{i=1}^{N} s_i(\hat c) \right\rvert^{2}}{N \sum_{i=1}^{N} \lvert s_i(\hat c)\rvert^{2}}, \qquad R_1(\hat c) = \frac{\mathrm{Re}\sum_i s_i\, s_{i+1}^{*}}{\sqrt{\sum_i \lvert s_i\rvert^{2} \sum_i \lvert s_{i+1}\rvert^{2}}}, \qquad B(\hat c) = \sum_{\mathrm{ROI}} \Bigl\lvert \sum_i s_i \Bigr\rvert^{2} \]

(32.3)

\(s_i(\hat c)\): complex (IQ) signal of element \(i\) after delaying with speed \(\hat c\)
\(\mathrm{CF}\): coherence factor, 1 for perfectly aligned equal signals, ≈ 1/N for noise (Module 25)
\(R_1\): lag-one spatial coherence
\(B\): speckle brightness (energy) of a region
\(\hat c\): trial (assumed) sound speed used for the delays, m/s
\(i,\ N\): element index and number of elements
\(\mathrm{Re},\ (\cdot)^{*}\): real part, complex conjugate
\(\mathrm{ROI}\): region of interest (set of pixels summed)

For a quadratic residual phase with edge value \(\Phi\) across a uniform aperture, \(\mathrm{CF} = \lvert\langle e^{j\Phi u^{2}}\rangle\rvert^{2}\): 0.95 at \(\Phi = \pi/4\), 0.5 at 2.73 rad, 0.39 at π. The CF peak is therefore broad (half-width about 67 m/s for a 19.2 mm aperture at 30 mm, F/1.6), and its centre, not its width, carries the estimate.

Anderson and Trahey (J. Acoust. Soc. Am., 1998) estimated the speed more directly from the arrival-time profile of the echoes across the aperture: in a homogeneous medium the profile is a hyperbola whose curvature fixes \(\hat c\) once the apex time is known. Both families answer the same question, "which homogeneous speed best explains the wavefront?", and both are exact only for a homogeneous medium. Under a fat layer they return an average. Which one?

Write the vertical one-way time through layers of thickness \(d_i\) and speed \(c_i\) as \(T = \sum_i \Delta T_i\), \(\Delta T_i = d_i/c_i\), and the time weights \(w_i = \Delta T_i/T\). Two averages are natural:

\[ \begin{gathered} V_h = \frac{z}{T} = \sum_i w_i\, c_i, \qquad V_{\mathrm{rms}}^{2} = \sum_i w_i\, c_i^{2} \\ V_{\mathrm{rms}} - V_h \approx \frac{\operatorname{Var}_w(c)}{2V_h} \;\overset{\text{two layers}}{=}\; \frac{w_1 w_2\,(c_2 - c_1)^{2}}{2V_h} \end{gathered} \]

(32.4)

\(V_h\): harmonic (time-average) speed: depth over one-way time, m/s
\(V_{\mathrm{rms}}\): root-mean-square speed, time-weighted, m/s
\(w_i\): fraction of the vertical travel time spent in layer \(i\)
\(\operatorname{Var}_w(c)\): time-weighted variance of the layer speeds
\(z,\ T\): depth (m) and one-way travel time (s)
\(c_i\): sound speed of layer \(i\), m/s
\(w_1,\ w_2,\ c_1,\ c_2\): time fractions and speeds of the two layers

For small offsets the moveout of a reflection under flat layers is \(t^{2}(x) \approx t_0^{2} + x^{2}/V_{\mathrm{rms}}^{2}\) (Dix, Geophysics, 1955), so a curvature fit or a focus maximizer returns \(V_{\mathrm{rms}}\), not \(z/T\). For 10 mm of fat (1450) over liver (1580) at 40 mm: \(V_h\) = 1545.4, \(V_{\mathrm{rms}}\) = 1546.4 m/s. The difference (1.1 m/s) is small; the difference from the liver itself (33 m/s) is not.

m32_global_sos_focus.m checks this on simulated channel data with Snell refraction: the CF-optimal speed lands 0.6–1.0 m/s above \(V_{\mathrm{rms}}\) at every depth from 15 to 50 mm (1547.1 m/s at 40 mm), while the liver below is at 1580. Demo 32.1 runs the same physics on a point target and shows the residual phase across the aperture: at the optimum the quadratic part vanishes and only the non-hyperbolic (quartic) moveout of the layered medium is left.

Demo 32.1

Focus quality versus assumed speed under a fat layer: average, not local

SimulationPoint target on axis under a flat fat layer. True receive times by Snell ray tracing (bisection on the ray parameter, JS); 0° plane-wave transmit (vertical, exact through flat layers). Beamformer: homogeneous \(\hat c\), pixel tied to the two-way apex time. Each element contributes a Gaussian pulse (60 % bandwidth) at its residual delay; CF (32.3) of those signals. Same physics as m32_global_sos_focus.m without noise or the ROI average. No attenuation, no element directivity.

CF(ĉ)V_h = z/TV_rmslocal speed at the target

Fat thickness 10 mm

Fat speed 1450 m/s

Liver speed 1580 m/s

Target depth 40 mm

Receive F-number 1.0

Assumed speed ĉ (phase panel) 1540 m/s

Try: with the defaults the CF peaks within about a metre per second of V_rms, 33 m/s below the liver it is meant to measure. Set the fat to 0: the peak moves onto the liver speed. Push the target to 90 mm: the average approaches the liver but never reaches it, and the peak broadens (the aperture is capped at 38.4 mm, so F grows). Raise the F-number to 3: the peak widens as F² predicts (32.1). In the phase panel, move ĉ off the optimum: the residual is a parabola; at the optimum only a small quartic term remains, the non-hyperbolic moveout of the layered medium.

32.3From average to local: layered media, Dix and the price of differentiation

A focus maximizer applied at many depths returns an average profile \(\bar c(z)\). The local speed is hidden in how that average changes with depth. In the vertical straight-ray model the one-way time to depth \(z\) is the integral of slowness, and the measured average is \(\bar c(z) = z/T(z)\). Written in slowness, the relation is linear, and inverting it is a derivative:

\[ \begin{gathered} z\,\bar s(z) = \int_0^{z} s(z')\,dz' \;\Longrightarrow\; s(z) = \frac{d}{dz}\bigl[z\,\bar s(z)\bigr] \\ s_k \approx \frac{z_k \bar s_k - z_{k-1}\bar s_{k-1}}{\Delta z}, \qquad \sigma_{c,\mathrm{loc}} \approx \frac{\sqrt{z_k^{2} + z_{k-1}^{2}}}{\Delta z}\,\sigma_{\bar c} \approx \sqrt2\,\frac{z}{\Delta z}\,\sigma_{\bar c} \end{gathered} \]

(32.5)

\(\bar s(z) = 1/\bar c(z)\): average slowness from the surface to depth \(z\), s/m
\(s(z)\): local slowness at depth \(z\)
\(\Delta z\): spacing of two average measurements, m
\(\sigma_{\bar c},\ \sigma_{c,\mathrm{loc}}\): standard deviation of each average and of the local estimate, m/s (independent errors, \(c_{\mathrm{loc}} \approx \bar c\))
\(z\prime\): integration variable (depth), m
\(k\): index of the measurement depth \(z_k\)
\(z_k,\ \bar s_k\): depth and average slowness at measurement \(k\)
\(\bar c\): average speed \(1/\bar s\), m/s

Differencing two averages taken 1 mm apart at 30 mm multiplies their noise by 42: 2 m/s on each average becomes about 85 m/s on the local speed (83.5 from the formula; m32_layered_inversion.m predicts 85.8 because the slowness-to-speed conversion adds the factor \((c_{\mathrm{loc}}/\bar c)^{2}\), and finds 87.1 by Monte Carlo; 142 and 148 at 50 mm). Locality costs \(z/\Delta z\) in precision, and no estimator escapes that factor without a prior.

Seismologists met the same problem with stacking velocities, which are RMS averages (32.4). Dix's interval formula inverts them exactly for flat layers:

\[ c_{\mathrm{int},k}^{2} = \frac{V_{\mathrm{rms},k}^{2}\,T_k - V_{\mathrm{rms},k-1}^{2}\,T_{k-1}}{T_k - T_{k-1}}, \qquad \text{straight-ray:}\quad c_{\mathrm{int},k} = \frac{z_k - z_{k-1}}{T_k - T_{k-1}},\ \ z_k = \hat c_k T_k \]

(32.6)

\(c_{\mathrm{int},k}\): interval (local) speed between reflectors \(k-1\) and \(k\), m/s
\(T_k\): one-way vertical time to reflector \(k\), from the echo time
\(V_{\mathrm{rms},k},\ \hat c_k\): focus-optimal (RMS) speed measured at reflector \(k\)
\(z_k\): depth of reflector \(k\), m
\(k\): reflector index

On the CF-optimal speeds of m32_global_sos_focus.m the liver interval speed comes out at 1580.1–1580.6 m/s by Dix and 1580.5–1582.5 m/s by the straight-ray rule (true 1580): the right formula matters at the 1–2 m/s level, the noise amplification of (32.5) at the 10–100 m/s level.

Two ways out of the noise. Regularization (Module 13) trades resolution for precision: penalize the first difference of the local slowness and pick \(\lambda\) by the discrepancy principle, \(\lVert \mathbf{W}(\mathbf{A}\mathbf{s} - \mathbf{b})\rVert^{2} \approx m\). A layer prior uses the interfaces visible in the B-mode image and estimates one speed per layer, which turns the ill-posed derivative into a well-posed least-squares problem with three unknowns:

\[ \hat{\mathbf{s}}_\lambda = \arg\min_{\mathbf{s}} \lVert \mathbf{W}(\mathbf{A}\mathbf{s} - \mathbf{b})\rVert^{2} + \lambda^{2}\lVert \mathbf{D}\mathbf{s}\rVert^{2}, \qquad \hat{\mathbf{s}}_{\mathrm{layer}} = \mathbf{P}\,(\mathbf{B}^{T}\mathbf{B})^{-1}\mathbf{B}^{T}\mathbf{W}\mathbf{b},\ \ \mathbf{B} = \mathbf{W}\mathbf{A}\mathbf{P} \]

(32.7)

\(\mathbf{A}\): cumulative-sum operator: \(A_{kj} = \Delta z\) for cells \(j\) above measurement depth \(z_k\)
\(\mathbf{b}\): data \(b_k = z_k/\bar c_k\), s
\(\mathbf{W}\): whitening, \(1/(z_k\sigma_{\bar s})\) on the diagonal
\(\mathbf{D}\): first-difference matrix
\(\mathbf{P}\): cell-to-layer indicator matrix (interfaces from the image)
\(\mathbf{s}\): slowness per cell, s/m
\(\hat{\mathbf{s}}_\lambda,\ \hat{\mathbf{s}}_{\mathrm{layer}}\): regularized and layer-constrained slowness estimates
\(\lambda\): regularization weight (balance of data fit and smoothness)
\(\lVert\cdot\rVert,\ \arg\min\): Euclidean norm, minimizing argument
\(\mathbf{B}\): \(\mathbf{W}\mathbf{A}\mathbf{P}\), the weighted forward operator on layers
\((\cdot)^{T}\): transpose

m32_layered_inversion.m, 2 m/s noise per average: Tikhonov (λ = 1.5×10⁵) gives the liver at 1548.6 m/s with 6.8 m/s RMSE but rounds the interfaces; the layer fit gives 1449.3 ± 0.8, 1587.6 ± 2.5 and 1548.9 ± 0.7 m/s for fat, muscle and liver (true 1450, 1585, 1550). The average at 60 mm is 1539.2 m/s, 10.8 m/s below the liver.

This is the logic of Jakovljevic et al. (J. Acoust. Soc. Am., 2018): Anderson-type averages first, then a model relating local speeds to path averages, solved by gradient descent. In two-layer phantoms and simulations their local estimates had biases below 6.2 m/s where the matched averages were biased by up to 66 m/s, at the cost of a larger standard deviation that a 5 × 5 mm region brought below 7 m/s. Ali and Dahl later moved the same idea to travel-time tomography from average speeds (IEEE IUS 2019), and Brevett, Sanabria, Ali and Dahl to common-midpoint gathers of non-beamformed data (IEEE IUS 2022).

Pitfall

Reporting the focus-optimal speed of a deep liver ROI as "the liver speed". Under 10 mm of fat (1450 over 1580 m/s) it reads about 45, 34, 22 and 17 m/s low (\(V_{\mathrm{rms}}\)) at 30, 40, 60 and 80 mm, more than the 14 m/s between the two Imbault cut-offs. A global estimate is a property of the path, not of the tissue at the end of it.

Demo 32.2

Average-to-local inversion: naive differencing, Tikhonov and a layer prior on the same noisy data

SimulationVertical straight-ray model (32.5) on 60 cells of 1 mm; averages \(\bar c(z_k)\) from 4 mm down, every 1–4 mm, with Gaussian noise of the set standard deviation (seeded, JS). Linear least-squares solutions of (32.7) by Cholesky. Same model as m32_layered_inversion.m (different random numbers, so the readouts differ slightly). No refraction, no lateral variation, and the averages are taken as the harmonic \(z/T\), not \(V_{\mathrm{rms}}\).

true localmeasured averagenaive differenceTikhonovlayer fit

Medium

Noise on each average 2 m/s

Spacing of the averages 1 mm

Tikhonov λ (log₁₀) 5.2

pick λ by the discrepancy principle

Try: the grey naive curve swings by ±100 m/s at depth even though every average is within 2 m/s: that is the z/Δz factor of (32.5). Untick the discrepancy box and slide λ: small λ follows the noise, large λ flattens the fat–muscle step and drags the liver toward the average. The layer fit is the best of the three whenever the interfaces are right; pick the lesion case and it misses the lesion entirely, because the prior has no layer there, and the liver value absorbs it (about 1564 instead of 1580 m/s with an error bar of ±0.7), while Tikhonov shows a blurred dip. Widen the spacing to 4 mm: the naive noise drops fourfold, the resolution with it.

32.4Echo-mode tomography I: what CUTE measures

Averages along vertical paths cannot see lateral structure. Jaeger and colleagues (Ultrasound Med. Biol., 2015) added angle diversity with ordinary pulse-echo hardware: beamform the same region with \(c_0\) from several steered plane-wave transmits, and compare the phase of the complex images pixel by pixel. If the tissue is homogeneous at \(c_0\), every steered image of a given speckle pixel has the same phase. If not, each transmit direction accumulates a different travel-time error on its way to the pixel, and the echo phase shifts. In the straight-ray, first-order model (Section 13.8 built its linear algebra):

\[ L_\theta(x,z) = \frac{1}{\cos\theta}\int_0^{z} \Delta s\bigl(x - (z - z')\tan\theta,\ z'\bigr)\,dz', \qquad \Delta\tau_{\theta}(x,z) = L_\theta - L_0, \qquad \Delta\phi = 2\pi f_0\,\Delta\tau \]

(32.8)

\(L_\theta(x,z)\): travel-time error accumulated along the straight path from pixel \((x,z)\) to the array at angle \(\theta\) from vertical, s
\(\Delta s\): slowness deviation (32.2)
\(\Delta\tau_\theta\): echo delay shift between the image transmitted at \(\theta\) and at 0°, with the receive beamforming kept fixed (it cancels), s
\(\Delta\phi\): measured phase shift at the centre frequency \(f_0\), rad
\(x,\ z,\ z\prime\): lateral position, depth of the pixel and integration depth, m
\(\theta\): transmit steering angle, rad
\(\Delta\tau\): delay shift, s
\(f_0\): centre frequency, Hz

A +40 m/s inclusion of radius 5 mm (\(\Delta s\) = −16.4 ns/mm) gives up to 169 ns of shift at 15° below its centre, where the 0° path crosses the whole inclusion and the 15° path misses it (m32_cute_phase_shift.m; 164 ns from the inclusion chord, 5.7 ns from the 4 mm fat layer of the script), 5.3 rad at 5 MHz. The measurement is linear in \(\Delta s\) and differential: only differences between paths are seen.

Two practical limits follow at once. Wrapping: a phase is known modulo 2π, so \(\lvert\Delta\tau\rvert\) must stay below \(1/(2f_0)\) = 100 ns at 5 MHz. CUTE therefore measures between small angle steps and accumulates: a 2° step reaches only 56 ns (1.8 rad) for the same inclusion. Decorrelation: the phase of a speckle pixel is random, and only its change is meaningful if the two images contain the same speckle. Steering the transmit alone changes the speckle pattern. Stähli et al. (Ultrasonics, 2020) therefore paired transmit angle \(\varphi\) and receive angle \(\psi\) and compared only pairs that share the mid-angle \(\gamma = (\varphi + \psi)/2\), whose isochrones (and speckle) stay in place:

\[ \begin{gathered} d_{\gamma}(x,z) = \bigl[L_{\varphi_n} + L_{\psi_m}\bigr] - \bigl[L_{\varphi_{n'}} + L_{\psi_{m'}}\bigr], \qquad \tfrac12(\varphi_n + \psi_m) = \tfrac12(\varphi_{n'} + \psi_{m'}) = \gamma \\ \text{flat layer, } (\theta,-\theta)\ \text{vs}\ (0,0)\text{:}\quad d = 2(\sec\theta - 1)\,S \end{gathered} \]

(32.9)

\(\varphi,\ \psi\): transmit and receive steering angles
\(\gamma\): common mid-angle
\(S\): vertical one-way delay error of a laterally uniform layer, s
\(d_{\gamma}\): difference in travel-time error between two transmit–receive angle pairs with the same mid-angle \(\gamma\), s
\(n,\ n\prime,\ m,\ m\prime\): indices of the transmit and receive angles of the two pairs
\(L\): travel-time error along one ray, as in (32.8), s
\(\theta\): steering angle, rad
\(\sec\): secant, \(1/\cos\)

A laterally uniform layer is nearly invisible: 4 mm of fat at 1450 m/s delays everything below it by 161 ns, yet changes a 15° transmit-only difference by 5.7 ns and the common-mid-angle difference by 11.4 ns, the \((\sec\theta - 1)\) and \(2(\sec\theta - 1)\) factors (3.5 % and 7.1 %). Layers, the thing liver measurements care about most, are the weakest part of the CUTE signal.

The noise of a phase-shift datum is the phase of a sample correlation. For \(N\) independent speckle samples in the kernel and a correlation coefficient \(\rho\) between the two images, the large-\(N\) result (the Cramér–Rao bound of interferometric phase, Module 12) is

\[ \hat\phi = \arg\sum_{k=1}^{N} a_{2,k}\,a_{1,k}^{*}, \qquad \sigma_{\phi} \approx \frac{\sqrt{1-\rho^{2}}}{\rho\sqrt{2N}}, \qquad \sigma_\tau = \frac{\sigma_\phi}{2\pi f_0} \]

(32.10)

\(a_{1,k},\ a_{2,k}\): complex pixel values of the two beamformed images in the kernel
\(N\): number of independent speckle cells in the kernel
\(\rho\): correlation coefficient between the two images (speckle decorrelation plus noise)
\(\sigma_\phi,\ \sigma_\tau\): phase and delay jitter of one datum
\(\hat\phi\): estimated phase shift between the two images, rad
\(k\): pixel index in the kernel
\(\arg\): phase angle, rad
\((\cdot)^{*}\): complex conjugate
\(f_0\): centre frequency, Hz

Monte Carlo in m32_cute_phase_shift.m: N = 64, ρ = 0.95 gives 0.0295 rad (formula 0.0291), 0.94 ns at 5 MHz; ρ = 0.99 gives 0.41 ns. The formula fails at low ρ and small N (0.355 versus 0.306 rad at ρ = 0.5, N = 16) where the phase distribution grows tails. Correlation is worth more than kernel size: going from ρ = 0.9 to 0.99 is worth a factor 3.4 in σ, the same as 12 times more samples.

Put the pieces together for the hardest case: the speed of a 5 mm layer from common-mid-angle data at ±20° with 0.92 ns jitter (ρ = 0.95, N = 64) is uncertain by 2.0 ns/mm, i.e. 5 m/s, from a single pair of pixel data (the layer delay \(S\) is the datum just below the layer minus the datum just above it, so \(\sqrt2 \times 0.92\) ns divided by \(2(\sec 20° - 1) \times 5\) mm). Averaging laterally helps; the axial blur of Section 32.5 does not. Stähli et al. also showed that a quantitative forward model needs a second correction: the misfocused echo is reconstructed at a shifted position, which adds a phase shift of its own. With both features their liver-phantom reconstructions came within a few m/s of the through-transmission values (standard error 2.5–10 m/s) and a volunteer liver read 1564 ± 4 m/s.

Schematic Figure 32.A · CUTE geometry: transmit-only versus common-mid-angle pairs

Not to scale. Angles in Stähli et al. (2020): plane-wave transmits from −27.5° to 27.5°, target transmit and receive angles −25° to 25° in 2° steps, L7-4 probe at 5 MHz.

Demo 32.3

CUTE phase-shift maps: inclusion, fat layer, steering angle, wrapping and speckle jitter

Approx. modelStraight-ray, first-order delay model (32.8)–(32.9) on a 0.25 mm grid (38.4 × 40 mm), marched row by row in JS with linear interpolation, smooth inclusion edge (0.3 mm). Optional jitter: independent Gaussian delay noise per pixel with \(\sigma_\tau\) from (32.10). No refraction, no echo-position shift, no beamforming; same model as m32_cute_phase_shift.m.

negative shiftpositive shiftcolumn profile±σ_τ bandwrap limit

Datum

Display

Steering angle θ 15°

Inclusion contrast +40 m/s

Inclusion radius 5 mm

Inclusion depth 15 mm

Fat layer (1450 m/s) 4 mm

Centre frequency 5 MHz

Speckle correlation ρ 0.95

Kernel samples N 64

Profile column x 0 mm

add the jitter of (32.10) to the map

Try: at 15° the shift below the inclusion reaches the full chord delay, about 170 ns: switch to wrapped phase and count the wraps. Choose the 2° step: no wrap, and this is why CUTE differences small steps and accumulates. Set the inclusion contrast to 0: what is left is the fat layer, a few ns even though the layer delays every echo by 161 ns. Common mid-angle doubles both the inclusion and the layer signal. Raise f₀ to 8 MHz: more phase per ns, more wraps, and smaller σ_τ at the same ρ. Drop ρ to 0.6 with the jitter on and watch the layer signal disappear.

32.5Echo-mode tomography II: inversion, limited angle and the forward model

Stack every phase-shift datum of every pixel into one linear system. Each row of \(\mathbf{A}\) is a difference of straight-ray weight vectors (path lengths through the cells of the slowness grid), and the inversion is the general-form Tikhonov problem of Module 13 with a gradient penalty:

\[ \begin{gathered} \mathbf{d} = \mathbf{A}\,\Delta\mathbf{s} + \mathbf{n}, \qquad \Delta\hat{\mathbf{s}} = \bigl(\mathbf{A}^{T}\mathbf{A} + \lambda^{2}\mathbf{G}^{T}\mathbf{G}\bigr)^{-1}\mathbf{A}^{T}\mathbf{d} \\ \mathbb{E}\bigl[\Delta\hat{\mathbf{s}}\bigr] = \mathbf{R}\,\Delta\mathbf{s}, \qquad \mathbf{R} = \bigl(\mathbf{A}^{T}\mathbf{A} + \lambda^{2}\mathbf{G}^{T}\mathbf{G}\bigr)^{-1}\mathbf{A}^{T}\mathbf{A} \end{gathered} \]

(32.11)

\(\mathbf{d}\): all delay-shift data (32.8)–(32.9), ns
\(\mathbf{A}\): differences of ray-weight rows, mm (path length per cell)
\(\Delta\mathbf{s}\): slowness deviation per cell, ns/mm
\(\mathbf{G}\): lateral and axial first differences between neighbouring cells
\(\mathbf{R}\): resolution matrix: row \(k\) is the point-spread function of the reconstruction at cell \(k\); \(\operatorname{tr}\mathbf{R}\) counts the resolved degrees of freedom
\(\mathbf{n}\): measurement noise, ns
\(\Delta\hat{\mathbf{s}}\): regularized estimate of the slowness deviation
\(\lambda\): regularization weight
\(\mathbb{E}\): expectation over the noise
\((\cdot)^{T},\ (\cdot)^{-1}\): transpose, matrix inverse

\(\mathbf{R}\) does not depend on the data, only on the acquisition design and \(\lambda\). Compute it before acquiring anything: it is the design tool of SoS tomography, as the beam pattern is of an array.

Limited angle. Rays within ±θ_max of vertical measure, by the Fourier-slice theorem, only a wedge of spatial frequencies around the \(k_x\) axis (Section 13.8). Everything near \(k_z\), horizontal edges and the axial profile, sits in the missing cone. The resolution kernels are therefore elongated in depth and get worse with depth: m32_sos_tomography.m (θ_max = 25°, 7 angles, 2.4 mm cells) finds axial/lateral FWHM of 2.5/2.4, 2.6/2.4 and 3.4/2.4 mm at 9, 19 and 31 mm for transmit-only data. Layers and the absolute level are in the same blind spot: a constant \(\Delta s_0\) shifts every datum by only \(\Delta s_0 z(\sec\theta - 1)\), so the mean speed must come from somewhere else. That is why CUTE reconstructions are run with an a priori speed and why Stähli et al. (2020, Figs. 4–6 of the preprint) and Jaeger et al. (2022, Figs. 9–10) show their maps for several a priori values.

Forward-model error dominates noise. The same script separates the error sources by running each design twice. With data generated by the cell model itself plus 2 ns noise (noise and regularization only) the RMSE is 3.8 m/s for transmit-only and 1.9 m/s for common mid-angle. With noise-free data from fine straight rays through a smooth truth (discretization and regularization only) it is 22.0 and 30.8 m/s. The model, not the noise, is the error budget, and adding more data rows does not fix it. In the ray model a common-mid-angle row \((W_i + W_j) - (W_{i-1} + W_{j+1})\) is exactly a combination of transmit-only rows of the same pixel (the script verifies a residual of 5×10⁻¹⁵): CMA adds no ray information. Its benefit lives outside this model, in the speckle correlation of (32.10) and the echo-position term that Stähli et al. added. Spatial-domain formulations of the same reconstruction (Sanabria et al., Phys. Med. Biol., 2018) and learned forward models are attempts to shrink exactly this model error.

Intuition

Echo-mode SoS tomography is transmission tomography with the sources and detectors on the same side, so every ray is nearly vertical. It sees lateral changes well, axial changes poorly and uniform offsets hardly at all, which is the opposite of what a liver protocol needs (a layered abdominal wall over a uniform organ). Every successful method adds what the rays cannot give: an a priori mean, interfaces from the image, or a better forward model.

Demo 32.4

Limited-angle SoS tomography: reconstruction and resolution-matrix kernels for two acquisition designs

Simulation16 × 12 cells of 2.4 mm (x ±19.2 mm, z 2–30.8 mm) under a 38.4 mm array. Data at every cell centre from straight-ray delays at K angles in ±θ_max (rays must reach the aperture); transmit-only = adjacent-angle differences, CMA = adjacent pairs with a common mid-angle. "Cell-consistent" data use the inversion's own ray weights (0.3 mm steps, nearest cell); "fine-ray" data integrate the smooth truth with 0.1 mm steps (model error). Tikhonov with a gradient penalty (32.11), Cholesky in JS. Smaller grid than m32_sos_tomography.m (16 × 14), same construction.

colour: 1440 (blue) to 1640 m/s (pink), white = 1540axial profilelateral profile

Phantom

Design

Maximum angle θ_max 25°

Number of angles K 7

Delay noise 2 ns

λ (log₁₀) 0

Data

Try: click a cell in either panel to see its resolution kernel (a row of R): narrow laterally, long axially, and longer with depth. Shrink θ_max to 10°: the kernels become vertical streaks, the missing cone made visible. Choose the layered phantom: neither design recovers the 5 mm fat layer, and the reconstruction drifts in mean level. Switch the data to fine-ray and press the oracle button for each design: the error now comes from the model, not the 2 ns noise; set the noise to 0 and it barely changes.

32.6Autofocus as optimization: differentiable beamforming, learned and full-wave inversion

Module 14 climbed a single \(\hat c\) by gradient ascent on a focus metric. Differentiable beamforming does the same for a whole map. Simson et al. (MICCAI 2023) wrote delay-and-sum in an automatic-differentiation framework so that a loss computed on the beamformed data could be differentiated with respect to every cell of a sound-speed map, and optimized phase-shift, speckle-brightness and coherence losses by gradient descent. With straight rays, the chain rule shows what that gradient is:

\[ \tau_i(\mathbf{s}) = \sum_j \ell_{ij}\, s_j \;\Longrightarrow\; \frac{\partial \mathcal{L}}{\partial s_j} = \sum_i \frac{\partial \mathcal{L}}{\partial \tau_i}\,\ell_{ij} = \bigl[\mathbf{A}^{T}\,\nabla_{\tau}\mathcal{L}\bigr]_j \]

(32.12)

\(\tau_i\): beamforming delay of channel/path \(i\) (transmit plus receive) for a pixel, s
\(\ell_{ij}\): length of path \(i\) in cell \(j\), m (the ray weights of (32.11))
\(\mathcal{L}\): scalar loss on the beamformed data (phase error, brightness, coherence)
\(\nabla_{\tau}\mathcal{L}\): sensitivity of the loss to each delay, computed by automatic differentiation through the interpolator
\(s_j\): slowness of cell \(j\), s/m
\(\mathbf{A}\): ray-weight matrix with entries \(\ell_{ij}\), as in (32.11)
\(\partial\): partial derivative
\(i,\ j\): path index and cell index

The gradient is a back-projection of delay residuals along the rays: one step of differentiable autofocus is one step of a nonlinear, iterative tomography with the same straight-ray \(\mathbf{A}\) as CUTE, recomputed at the current map. It inherits the missing cone, and the iteration count acts as the regularization parameter (Section 13.5).

The later formulation (Simson et al., arXiv:2410.03008) replaced the loss with a common-midpoint phase error, the mean absolute phase of correlations between lag-1 common-midpoint subaperture pairs, keeping only pairs with correlation above 0.7, plus total-variation regularization of the map:

\[ \mathcal{L}_{\mathrm{CMPE}} = \frac{1}{N_x N_{\mathrm{CM}}} \sum_{\mathbf{x}_f \in \mathrm{FOV}}\ \sum_{(p,q)\in \mathrm{CM}} \Bigl\lvert \angle \hat\Gamma^{\mathrm{CM}}_{pq}(\mathbf{x}_f) \Bigr\rvert, \qquad \mathcal{L} = \mathcal{L}_{\mathrm{CMPE}} + \alpha\,\mathrm{TV}\bigl(c(\mathbf{x})\bigr) \]

(32.13)

\(\hat\Gamma^{\mathrm{CM}}_{pq}(\mathbf{x}_f)\): complex correlation at pixel \(\mathbf{x}_f\) between the images of subaperture pairs \(p\) and \(q\) that share a common midpoint
\(N_x,\ N_{\mathrm{CM}}\): number of pixels and of common-midpoint pairs
\(\alpha\): weight of the total-variation penalty
\(\angle\): phase angle, rad
\(\mathrm{TV}\): total variation of the sound-speed map
\(c(\mathbf{x})\): sound-speed map being optimized, m/s
\(\mathrm{FOV}\): field of view (set of pixels \(\mathbf{x}_f\))
\(\mathcal{L}_{\mathrm{CMPE}}\): common-midpoint phase-error loss, rad

The common-midpoint pairing is the aperture-domain analogue of CUTE's common mid-angle: it keeps the speckle correlated so that the phase means delay, not speckle. Reported errors (arXiv version 3): −11 m/s in a simulated inclusion, 8–14 m/s below calibrated values in layered and commercial phantoms, from 12 m/s below to 42 m/s above in graphite–gelatin phantoms.

Why phase losses have local minima. A phase is periodic in delay with period \(1/f_0\). Once a residual delay difference exceeds half a period the wrapped phase points the wrong way, and gradient descent locks onto a wrong cycle, the cycle skipping of Section 14.4. The ways out are the ones seismology uses: start from a good global estimate (Section 32.2), compare only nearby (lag-1) and well-correlated signals, go from low to high frequency, or use envelope and coherence losses that are smooth but less sharp. Demo 32.5 maps the loss over a two-layer model and lets you start the descent anywhere.

Learned inversion. Feigin, Freedman and Anthony (IEEE Trans. Biomed. Eng., 2020; preprint 2018) trained a fully convolutional network to map raw channel data of three plane waves directly to a sound-speed map. It was trained on 6,026 k-Wave simulations of random elliptical inclusions (1300–1800 m/s) in speckle and reached 20.5 m/s RMSE on simulated test data; on a phantom and in vivo the results were qualitative. The forward model is now implicit in the training set, so the model error of Section 32.5 becomes a domain gap that no regularizer can see.

Full-wave inversion. Straight rays ignore refraction and diffraction, which are not small at fat–muscle interfaces. Geophysics solved the same problem with wave-equation migration velocity analysis; Ali's thesis (Stanford, 2021) laid the groundwork for it in pulse-echo ultrasound. On the IEEE IUS 2026 QUASOS challenge, which scored bias, variance, lesion detectability and robustness to near-field aberration on phantom and simulated RF data, the published final sound-speed scores (lower is better) put R. Ali (Rochester) first at 22.6 and the Bern CUTE group (M. Jaeger) ninth at 72.5 among 11 entries.

Demo 32.5

Differentiable autofocus on a two-layer model: loss landscape, gradient path and local minima

SimulationFour point targets on axis at 20, 30, 40 and 50 mm under a flat fat layer of known thickness; true receive times by Snell ray tracing (JS). The beamformer's model is two layers with unknown speeds \((\hat c_1, \hat c_2)\) and straight rays, pixel tied to the two-way apex time, as in a straight-ray differentiable beamformer. Losses: 1 − mean CF (32.3), or the mean absolute wrapped phase difference between elements \(\ell\) apart, normalized by π (a simplified stand-in for phase-error losses such as (32.13), not CMPE itself). Gradient by central differences; normalized gradient descent with a decaying step.

truthdescent pathgrid minimumclick the map to start from there

Loss

Element lag ℓ (phase loss) 8

Centre frequency 5 MHz

Fat thickness (known) 12 mm

True fat speed 1450 m/s

True liver speed 1570 m/s

Initial step 12 m/s

Try: with the CF loss the valley is long and tilted: fat and liver speed trade against each other (a slower fat and a faster liver give nearly the same average), which is the layer ambiguity of Sections 32.3–32.5 seen as an optimization landscape. Switch to the phase loss: at lag 8 and 5 MHz it still converges from every corner. Raise the lag to 16, or the frequency to 10 MHz, and click near the lower-left corner: the descent stalls on a wrapped plateau, the cycle skipping of Module 14. Lag 1 is smooth but shallow. The grid minimum sits about 2 m/s from the true fat speed: the straight-ray model ignores refraction, a bias no optimizer removes.

32.7Error budget, validation and clinical numbers

How precise can a sound-speed estimate be? Take the cleanest case: a wavefront-curvature fit (Anderson and Trahey's idea) to receive arrival times \(t_i = t_a + (\sqrt{z^{2}+x_i^{2}} - z)/c\) with \(z = c\,t_a\), independent per-element jitter \(\sigma_t\) (whose own bound, set by SNR, bandwidth and decorrelation, is the delay CRLB of Modules 11 and 12), and an unknown common offset. The Fisher information gives, in the parabolic limit,

\[ \sigma_c \;\ge\; \frac{c^{2} z\,\sigma_t}{\sqrt{N}\,\operatorname{sd}(x_i^{2})} \;=\; c^{2} z\,\sigma_t \sqrt{\frac{45}{4N}}\,\frac{1}{X^{2}} \;=\; \frac{4c^{2}F^{2}\sigma_t}{z}\sqrt{\frac{45}{4N}}, \qquad b_{\mathrm{cal}} \approx -\frac{c^{2} z\,\delta}{X^{2}} \]

(32.14)

\(\sigma_t\): per-element arrival-time jitter, s
\(N,\ X\): number of elements and half-aperture, uniform spacing; \(\operatorname{sd}(x^{2}) = X^{2}\sqrt{4/45}\) for a continuous uniform aperture
\(F\): F-number \(z/(2X)\)
\(b_{\mathrm{cal}}\): bias from a systematic delay error that grows quadratically to \(\delta\) at the aperture edge (lens, element delay calibration), m/s
\(c,\ z\): sound speed (m/s) and depth (m)
\(x_i\): lateral position of element \(i\) in the aperture, m
\(\operatorname{sd}\): standard deviation over the aperture
\(\delta\): systematic delay error at the aperture edge, s
\(\sigma_c\): standard deviation of the speed estimate, m/s

m32_sos_error_budget.m: 64 elements (18.9 mm), z = 30 mm, \(\sigma_t\) = 5 ns: exact CRLB 1.69 m/s, closed form 1.62 (1.67 for the continuous aperture), Monte Carlo 1.68 m/s. The closed form grows linearly with depth at a fixed aperture (the exact bound only once F ≳ 1; at F < 1 it grows more slowly) and falls about as \(N^{-2.5}\) at fixed pitch. A 1 ns edge calibration error biases \(c\) by −0.83 m/s, comparable to the CRLB, and it never averages out.

The precision floor is not the problem. Averaging over a region of a few hundred independent kernels pushes the random part well below 1 m/s, and the systematic terms take over. The complete budget for a liver measurement adds, roughly in decreasing order:

\[ \hat c - c_{\mathrm{liver}} = \underbrace{\bigl(V_{\mathrm{rms}} - c_{\mathrm{liver}}\bigr)}_{\text{layer bias, if an average}} + \underbrace{b_{\mathrm{model}}}_{\text{rays, refraction}} + \underbrace{b_{\mathrm{cal}}}_{\text{lens, delays}} + \underbrace{b_{\mathrm{reg}}}_{\text{prior, }\lambda} + \underbrace{\varepsilon}_{\sigma_c/\sqrt{M}} \]

(32.15)

\(b_{\mathrm{model}}\): forward-model error (straight rays, no echo-position term, discretization), m/s
\(b_{\mathrm{reg}}\): regularization and prior bias (a priori mean, smoothing across interfaces)
\(\varepsilon,\ M\): random error and number of independent kernels averaged
\(\hat c,\ c_{\mathrm{liver}}\): estimated and true liver sound speed, m/s
\(V_{\mathrm{rms}}\): as in (32.4), m/s
\(b_{\mathrm{cal}}\): as in (32.14), m/s
\(\lambda\): regularization weight, as in (32.7)
\(\sigma_c\): standard deviation of one kernel’s estimate, m/s

Typical sizes from this module: layer bias 17–45 m/s under 10 mm of fat (Section 32.3), model error tens of m/s in the tomography toy (Section 32.5), calibration about 1 m/s per ns, random error well under the rest after averaging. The targets: 14 m/s between the two Imbault cut-offs, and the roughly 40 m/s drop for grade-1 steatosis that Stähli et al. cite.

Validation. Phantom references come from through-transmission time of flight, quoted at ±5 m/s by Stähli et al. (2020) and Jaeger et al. (2022), so a phantom study cannot resolve errors much below that. Jaeger et al. (Phys. Med. Biol., 2022) carried CUTE to a convex probe (C6-1X, 3 MHz) for liver depths: in a four-layer abdominal phantom the liver agreed with its nominal value, axial resolution was about 10 mm at 30 mm depth, and in a volunteer the superficial liver read about 1580 m/s, while below about 50 mm the estimate deteriorated and, in vivo, varied by up to 100 m/s with the a priori speed below 60 mm. QUASOS (IEEE IUS 2026) added the first shared benchmark with phantom and simulated RF data, with and without near-field aberrators.

Ultrasound link

Sound speed is upstream of everything quantitative. Attenuation and backscatter estimates (Module 27) are computed on images formed with \(c_0\); shear-wave tracking (Module 30) and ULM localization (Module 31) inherit the depth scaling \(\hat c/c\) and the defocus; and aberration correction (Module 33) is, in its distributed form, a sound-speed map turned back into delays. A 1 % speed error is a 0.3 mm depth error at 30 mm, larger than the localization precision ULM claims.

Demo 32.6

Error budget of a liver sound-speed measurement

Approx. modelRandom term: parabolic CRLB (32.14) with the discrete-aperture spread, divided by √M for M independent kernels; local estimates by differencing over Δz multiply it by √2 z/Δz (32.5). Calibration: quadratic edge error \(\delta\), bias \(c^{2}z\delta/X^{2}\). Layer bias of an average: \(V_{\mathrm{rms}}\) of (32.4) minus the liver speed. Model error: entered by hand (no formula exists; Section 32.5 suggests tens of m/s for straight-ray tomography). Biases are added by magnitude (worst case), the random term in quadrature; the calibration bias is applied unchanged to both estimates.

random (1σ)biasestotaldashed: 14 m/s cut-off gap and 40 m/s grade-1 shift

Reported quantity

ROI depth 40 mm

Receive elements 64

Per-element jitter σ_t 5 ns

Independent kernels averaged M 100

Differencing span Δz (local) 10 mm

Edge calibration error δ 2 ns

Fat thickness (1450 m/s) over liver 1570 12 mm

Model error (entered) 5 m/s

Try: with the defaults the average is dominated by the layer bias (−37 m/s), more than twice the gap between the healthy and steatotic cut-offs. Switch to the local estimate: the layer bias disappears, but differencing over 10 mm at 40 mm multiplies the random term by 5 (32.5); average more kernels or widen Δz to win it back. Then raise the calibration error: a 10 ns lens-delay error costs several m/s at any M. The lesson of the capstone in one panel: precision is cheap, locality and calibration are not.

32.8Visual summary

Thirteen questions a sound-speed expert answers in five seconds. Each figure is computed from the formula or model named in its badge.

32.9How experts think about it

A sound-speed expert keeps three questions apart. What does the datum measure (a curvature, a phase shift between two paths, a focus metric), and along which paths? Which average, or which linear combination of the local slowness, does that make it? And what blind spots does the path geometry leave (layers, the mean, axial detail), so what prior or extra physics must fill them? Precision is rarely the limit; locality, calibration and the forward model are.

A speed error is a curvature error first and a depth error second.

|Δc| < c²F²/(2f₀z)7.9 m/s at F/1, 30 mm, 5 MHz1 % → 0.3 mm at 30 mm

With the pixel tied to two-way time the on-axis delay is always right; the edge element is off by \((\Delta c/c)X^{2}/(zc)\) (32.1). Setting the edge phase below π/4 gives the tolerance, which scales as F². The depth error \((\Delta c/c)z\) is linear and only matters for geometry and for registration with other modalities.

Exact formula

Receive only, point-like pixel, π/4 criterion; plane-wave transmits add their own error, and real tissue adds non-quadratic aberration.

Work in slowness: travel times are linear in 1/c, not in c.

1 ns/mm ↔ 2.37 m/s+40 m/s = −16.4 ns/mm

\(\tau = \int s\,dl\) is linear in \(s\), so every inversion in this module (layered, tomographic, gradient) is a linear or linearized problem in \(\Delta s\) (32.2). Converting back with \(\Delta c \approx -c_0^{2}\Delta s\) costs \(\Delta c^{2}/c\), 1.0 m/s (2.5 % of the correction) at 40 m/s; averaging speeds instead of slownesses biases means by \(\operatorname{Var}(c)/c\).

Exact formula

Linear only for fixed ray paths; bent rays make the problem nonlinear in \(s\) as well.

A focus maximizer measures the RMS speed of the path, not the tissue at the end of it.

V_rms − V_h = 1.1 m/sliver − ĉ* = 33 m/s10 mm fat, 40 mm

Small-offset moveout under flat layers is hyperbolic with \(V_{\mathrm{rms}}\) (Dix), so CF, coherence, brightness and curvature fits all return it (32.4). In m32_global_sos_focus.m the CF optimum sits 0.6–1.0 m/s above \(V_{\mathrm{rms}}\) at every depth; Demo 32.1 finds 0.1–0.3 m/s for a single point.

SimulationDemo 32.1 model: Snell-traced point target at 40 mm under 10 mm of fat, 128 elements.

Flat layers, small offsets; large apertures and strong contrasts add non-hyperbolic moveout, and lateral variation breaks the 1D picture.

The CF peak is broad: read its centre, never its width.

CF = 0.5 at Φ = 2.73 rad±67 m/s: 19.2 mm aperture at 30 mm (F/1.6)

A quadratic phase across a uniform aperture gives \(\mathrm{CF} = \lvert\langle e^{j\Phi u^{2}}\rangle\rvert^{2}\) (32.3). The half-maximum sits at an edge phase of 2.73 rad, tens of m/s from the peak, while the peak position can be read to a fraction of a m/s with parabolic interpolation. Width tracks F² and frequency, not tissue.

Exact formula

Point-like target and uniform weighting; speckle, apodization and noise flatten and bias the curve.

Locality costs z/Δz in precision.

σ_loc ≈ √2 (z/Δz) σ_avg2 → 85 m/s at 30 mm, Δz = 1 mm

The local slowness is the derivative of \(z\bar s(z)\) (32.5); differencing two averages of independent error multiplies it by \(\sqrt{z_k^{2} + z_{k-1}^{2}}/\Delta z\). m32_layered_inversion.m: 85.8 predicted, 87.1 by Monte Carlo at 30 mm. Every local estimator pays this, explicitly or through its regularizer.

Exact formula

Independent errors on the averages; correlated errors (common calibration, shared data) amplify less but bias more.

RMS averages need Dix; harmonic averages need straight-ray differences.

Dix: 1580.1–1580.6straight-ray: 1580.5–1582.5true 1580 m/s

Applying the straight-ray rule \(\Delta z/\Delta T\) to focus-optimal (RMS) speeds mixes two definitions; Dix's \(\sqrt{\Delta(V^{2}T)/\Delta T}\) (32.6) is the consistent inverse. On noise-free CF optima the difference is 0.4–1.9 m/s: negligible against noise, not against a claimed 1 m/s accuracy.

SimulationInterval speeds printed by m32_global_sos_focus.m (Octave), liver true 1580 m/s.

Flat layers and accurate echo times; with noise, both are dominated by the z/Δz amplification.

A layer prior is the strongest regularizer, and the most dangerous one.

±0.7–2.5 m/s per layera missed lesion biases its layer

Interfaces from the B-mode image turn the derivative into a three-unknown least-squares fit (32.7): 1449.3 ± 0.8, 1587.6 ± 2.5, 1548.9 ± 0.7 m/s in m32_layered_inversion.m. The error bars describe noise only; an unmodelled lesion moves the layer value while the error bar stays small (figure: noise-free fit of a liver containing a 10 mm lesion).

SimulationLayer fit (32.7) on noise-free averages, fat 8 mm | muscle 6 mm | liver 1580 with a 32–42 mm lesion of the plotted speed; JS.

Assumes flat, correctly segmented interfaces; tilted or curved boundaries and in-layer heterogeneity violate it silently.

Phase wraps at half a period: steer in small steps and accumulate.

1/(2f₀) = 100 ns at 5 MHz15°: 169 ns2°: 56 ns

The CUTE datum is a phase (32.8), so only \(\lvert\Delta\tau\rvert < 1/(2f_0)\) is unambiguous. A +40 m/s, 5 mm inclusion produces 169 ns between 0° and 15° but 56 ns for a 2° step (m32_cute_phase_shift.m). Stähli et al. used 2° steps and accumulated over 10°.

Exact formula

Accumulation adds the errors of each step; strong contrasts or high frequencies may need unwrapping across space as well.

In phase tracking, correlation is worth more than kernel size.

σ_φ ≈ √(1−ρ²)/(ρ√(2N))0.93 ns at ρ 0.95, N 64, 5 MHz

The bound of (32.10) falls as \(1/\sqrt N\) but as \(\sqrt{1-\rho^{2}}\) near ρ = 1: from 0.9 to 0.99 is a factor 3.4, worth 12 times more samples. That is the reason for common-mid-angle and common-midpoint pairings: they keep the speckle in place.

Exact formula

Large N, circular Gaussian speckle; at low ρ and small N the phase distribution grows tails (0.355 versus 0.306 rad at ρ 0.5, N 16).

Layers and the mean speed are in echo-mode tomography's blind spot.

tx-only: (sec θ − 1)SCMA: 2(sec θ − 1)S15°: 3.5 % / 7.1 %

A laterally uniform layer, or a constant offset of \(\Delta s\), shifts every steered path almost equally, so the differences that CUTE measures see only the \(\sec\theta - 1\) excess (32.9). 4 mm of fat delays everything by 161 ns and moves the 15° datum by 5.7 ns. Reconstructions need an a priori mean or an independent global estimate.

Exact formula

Straight rays, flat layer; refraction at strong interfaces adds a small extra signal.

Compute the resolution matrix before you acquire anything.

axial 2.5 → 3.4 mm with depthlateral 2.4 mm (cell)

\(\mathbf{R} = (\mathbf{A}^{T}\mathbf{A} + \lambda^{2}\mathbf{G}^{T}\mathbf{G})^{-1}\mathbf{A}^{T}\mathbf{A}\) (32.11) needs only the geometry. Its rows show the missing cone: m32_sos_tomography.m finds axial FWHM of 2.5, 2.6, 3.4 mm (transmit-only) and 2.6, 2.9, 4.2 mm (CMA) at 9, 19, 31 mm, lateral fixed at the cell size.

SimulationAxial FWHM of resolution-matrix rows printed by m32_sos_tomography.m (θ_max 25°, 7 angles, oracle λ).

Depends on λ and on the penalty; FWHM below the cell size cannot be read on this grid.

In SoS tomography the forward model, not the noise, is the error budget.

noise only: 3.8 / 1.9 m/smodel only: 22.0 / 30.8 m/s

Running the same inversion on data from its own cell model and on fine-ray data through a smooth truth separates the two (m32_sos_tomography.m, transmit-only / CMA). Removing the noise altogether moves the fine-ray RMSE by 0.2 m/s (22.2 to 22.0); refining the forward model (echo position, refraction, discretization) is what moves the number. That is the argument of Stähli et al. (2020).

SimulationOracle-λ RMSE printed by m32_sos_tomography.m.

Toy grid and straight rays on both sides; real data add refraction and speckle effects that make the model term larger, not smaller.

Differentiable beamforming is tomography with a back-projected gradient.

∂L/∂s = Aᵀ ∂L/∂τiterations act as λ

With straight rays the delay of each path is \(\sum_j \ell_{ij}s_j\), so the gradient of any loss is a back-projection of delay sensitivities along those rays (32.12). Each step is a nonlinear Landweber step on the same \(\mathbf{A}\) as CUTE; it inherits the missing cone and needs the same priors (TV in the CMPE formulation).

Schematic

Straight-ray beamformers only; a full-wave forward model replaces \(\ell_{ij}\) by wave-equation sensitivity kernels.

Phase losses skip cycles past half a period: start from the global estimate.

wrap at 1/(2f₀)lag 8 safe, lag 16 not (Demo 32.5)

A wrapped phase is periodic in the delay residual, so the loss has secondary minima spaced by one period (Section 14.4). The figure is the edge lag-ℓ phase loss for a homogeneous error at 30 mm, 5 MHz, 30 mm aperture: longer lags sharpen the main minimum and bring the side minima closer. Initialize with a coherence-based global speed and use lag-1, well-correlated pairs.

Exact formulaMean |wrapped phase difference| between elements ℓ apart for the quadratic residual of (32.1), normalized by π.

Homogeneous error, point target; speckle adds random phase that raises the floor and the correlation threshold.

Precision is cheap: σ_c ∝ F²σ_t/(z√N).

1.69 m/s per kernel64 el., 30 mm, 5 ns÷√M by averaging

The curvature CRLB (32.14) is 1.69 m/s for one kernel at 30 mm (Monte Carlo 1.68); a few hundred independent kernels bring it to 0.1 m/s. At fixed pitch it falls about as \(N^{-2.5}\) because the curvature grows as \(X^{2}\). Quote a CRLB only next to the biases of (32.15).

Exact formula

Independent per-element jitter, known apex time; speckle-induced jitter is correlated across neighbouring elements, so the true bound is worse.

Calibrate the lens and the element delays before blaming the algorithm.

b ≈ −c²zδ/X²−0.83 m/s per ns (64 el., 30 mm)

A systematic delay error that grows to δ at the aperture edge is indistinguishable from a curvature change (32.14). A lens whose speed or thickness differs from the nominal one produces exactly such an error, and it does not average out. Lens speed and thickness belong in the estimated parameters or in a calibration step.

Exact formula

Quadratic error pattern; linear (tilt) patterns bias position instead, and random element errors add jitter, not bias.

A learned sound-speed network inherits its simulator.

6,026 k-Wave samples20.5 m/s RMSE (simulated)

Feigin et al. (2020) mapped three plane waves of raw channel data to a speed map and reached 20.5 m/s RMSE on simulated test data; in vivo the results were qualitative. The forward model moved from the equations into the training distribution, so model error becomes domain gap: random ellipses in speckle are not an abdominal wall.

A statement about generalization, not capacity; physics-anchored hybrids (unrolled or model-based networks) reduce but do not remove it.

32.10Cheat sheet

Concept	Formula	Ultrasound use
Receive defocus from a speed error	\(\delta\tau(x) \approx (\Delta c/c)\,x^{2}/(zc)\)	104 ns at the edge of 38.4 mm, 30 mm, 20 m/s
Focusing tolerance	\(\lvert\Delta c\rvert < c^{2}F^{2}/(2f_0 z)\)	7.9 m/s at F/1, 30 mm, 5 MHz; scales as F²
Depth scaling	\(\hat z = (\hat c/c)\,z\)	1 % → 0.3 mm at 30 mm
Slowness linearization	\(\Delta c \approx -c_0^{2}\Delta s\)	1 ns/mm = 2.37 m/s; error \(\Delta c^{2}/c\)
Coherence factor	\(\lvert\sum s_i\rvert^{2}/(N\sum\lvert s_i\rvert^{2})\)	CF = 0.5 at an edge phase of 2.73 rad
Harmonic and RMS averages	\(V_h = \sum w_i c_i\), \(V_{\mathrm{rms}}^{2} = \sum w_i c_i^{2}\)	focus returns \(V_{\mathrm{rms}}\); difference ≈ Var/(2V)
Average to local (straight ray)	\(s(z) = d[z\bar s(z)]/dz\)	noise × \(\sqrt2\,z/\Delta z\): 85 m/s from 2 m/s at 30 mm, 1 mm
Dix interval speed	\(c_{\mathrm{int}}^{2} = \Delta(V_{\mathrm{rms}}^{2}T)/\Delta T\)	consistent inverse of focus-optimal speeds
CUTE datum	\(\Delta\tau = L_\theta - L_0\), \(L_\theta = \sec\theta\int \Delta s\,dz'\)	wrap at 1/(2f₀) = 100 ns at 5 MHz
Layer visibility	\((\sec\theta - 1)S\), CMA \(2(\sec\theta - 1)S\)	3.5 % / 7.1 % of the layer delay at 15°
Phase jitter of a correlation	\(\sigma_\phi \approx \sqrt{1-\rho^{2}}/(\rho\sqrt{2N})\)	0.93 ns at ρ 0.95, N 64, 5 MHz
Resolution matrix	\(\mathbf{R} = (\mathbf{A}^{T}\mathbf{A} + \lambda^{2}\mathbf{G}^{T}\mathbf{G})^{-1}\mathbf{A}^{T}\mathbf{A}\)	kernels elongated axially (missing cone)
Autofocus gradient	\(\partial\mathcal{L}/\partial\mathbf{s} = \mathbf{A}^{T}\nabla_\tau\mathcal{L}\)	differentiable beamforming = back-projection
Curvature CRLB	\(\sigma_c = 4c^{2}F^{2}\sigma_t\sqrt{45/(4N)}/z\)	1.7 m/s: 64 el., 30 mm, 5 ns
Calibration bias	\(b \approx -c^{2}z\delta/X^{2}\)	−0.83 m/s per ns of edge error (64 el., 30 mm)
Liver numbers (Imbault 2017)	1492–1604 m/s; cut-offs 1541 and 1555 m/s	the 14 m/s a liver estimator must resolve

32.11Interview and qualifying-exam questions

I1Why does a 1540 m/s assumption hurt a high-frequency, low-F-number image more than a low-frequency, high-F-number one? · screen

The speed error appears as a curvature error, \(\delta\tau(X) = (\Delta c/c)X^{2}/(zc)\); in phase that is \(2\pi f_0\delta\tau\). Requiring less than π/4 at the edge gives \(\lvert\Delta c\rvert < c^{2}F^{2}/(2f_0z)\): inversely proportional to frequency and proportional to F². At 5 MHz and 30 mm: 7.9 m/s at F/1, 71 m/s at F/3. Follow-up: and the depth error? \((\Delta c/c)z\), independent of aperture: 0.3 mm per percent at 30 mm.

I2You sweep the assumed speed and maximize coherence in a liver ROI at 50 mm under 10 mm of subcutaneous fat. What number did you measure? · screen

An average of the path, specifically the RMS speed \(V_{\mathrm{rms}} = \sqrt{\sum w_i c_i^{2}}\) (time weights), because small-offset moveout under flat layers is hyperbolic with \(V_{\mathrm{rms}}\) (Dix). With 1450 over 1580 m/s it reads about 1553, 27 m/s below the liver, and only about 1 m/s above \(z/T\). Follow-up: how do you get the liver itself? Measure at several depths and invert (Dix or regularized), or use a local method (tomography, layer fit with the B-mode interfaces).

I3Derive the noise of a local speed obtained by differencing two average speeds. · onsite

\(z\bar s(z) = \int_0^{z}s\), so \(s_k = (z_k\bar s_k - z_{k-1}\bar s_{k-1})/\Delta z\). With independent errors \(\sigma_{\bar s} = \sigma_{\bar c}/c^{2}\) on each average, \(\sigma_{s} = \sigma_{\bar s}\sqrt{z_k^{2} + z_{k-1}^{2}}/\Delta z\), i.e. \(\sigma_{c,\mathrm{loc}} \approx \sqrt2(z/\Delta z)\sigma_{\bar c}\). 2 m/s averages 1 mm apart at 30 mm give 85 m/s. Follow-up: what fixes it? Larger Δz (resolution), Tikhonov with a discrepancy-chosen λ, or a layer prior from the image, each a different bias.

I4Explain CUTE to a beamforming engineer in two minutes, including what limits it. · onsite

Beamform the same region with \(c_0\) from several steered plane waves (and receive angles); for each pixel the echo phase changes between steering angles in proportion to the difference of the travel-time errors along the two paths, \(\Delta\tau = L_\theta - L_{\theta'}\), \(L_\theta = \sec\theta\int\Delta s\). Stack all pixels and angles into \(\mathbf{d} = \mathbf{A}\Delta\mathbf{s}\) and invert with regularization. Limits: phase wraps at \(1/(2f_0)\) (small steps, accumulate), speckle decorrelation between angles (common mid-angle pairing), limited angle (axial blur, layers and the mean nearly invisible), and the forward model (straight rays, echo-position shift). Follow-up: what did Stähli et al. change? Common mid-angle pairs and the echo-position term.

I5Why is a laterally uniform fat layer hard for echo-mode tomography? Quantify. · onsite

All steered paths cross it, with lengths \(d/\cos\theta\); the difference to 0° is \((\sec\theta - 1)S\), where \(S\) is the vertical delay. For 4 mm at 1450 m/s, \(S\) = 161 ns, but the 15° datum moves by 5.7 ns (11.4 ns for a symmetric common-mid-angle pair), comparable to a few times the 1 ns jitter of one datum. A constant offset of the mean behaves the same way. Follow-up: so how do you anchor the mean? A global coherence estimate, a calibrated a priori value, or physics beyond straight rays.

I6What is the phase jitter of a CUTE datum, and what is the cheapest way to halve it? · onsite

\(\sigma_\phi \approx \sqrt{1-\rho^{2}}/(\rho\sqrt{2N})\), the interferometric-phase CRLB; at ρ = 0.95, N = 64: 0.029 rad, 0.93 ns at 5 MHz. Halving it by kernel size needs 4N (resolution loss); raising ρ from 0.95 to 0.987 does the same. Hence pairing schemes that keep speckle correlated (common mid-angle, common midpoint), small angle steps and good SNR. Follow-up: when does the formula fail? Low ρ, small N: the distribution grows tails and wraps.

I7Your SoS tomography gives 25 m/s RMSE on a phantom. More angles or a better forward model? · research

Check which term dominates first. Rerun the inversion on data synthesized by your own forward model with your measured noise level: if that RMSE is small (a few m/s, as in m32_sos_tomography.m: 3.8 and 1.9 m/s), the 25 m/s is model error (straight rays, refraction, echo-position shift, discretization, lens) and more angles will not help. Then look at the resolution matrix to see whether the phantom features are resolvable at all. Follow-up: how would you validate the forward model alone? Full-wave simulation of the phantom geometry, compared datum by datum.

I8What is the relationship between differentiable-beamforming autofocus and CUTE? · research

With straight rays, each beamforming delay is \(\sum_j \ell_{ij}s_j\); the chain rule makes the gradient of any delay-based loss \(\mathbf{A}^{T}\nabla_\tau\mathcal{L}\), a back-projection along the same rays CUTE uses. Gradient descent is a nonlinear iterative tomography that recomputes the delays (and could recompute bent rays) at each step; the iteration count and TV play the role of λ. It shares CUTE's missing cone and benefits from the same correlated pairings (common midpoint). Follow-up: what does it buy? Losses evaluated on the actual beamformed data, any differentiable forward model, and easy priors.

I9Where do phase-based autofocus losses get stuck, and how would you initialize? · research

Wherever a residual delay difference exceeds half a period, the wrapped phase reverses and the loss has a secondary minimum: cycle skipping (Module 14). Long subaperture lags and high frequency make it worse. Initialize with a coherence- or CF-based global speed (smooth, wide basin), use lag-1 pairs with a correlation threshold (as in CMPE), and go from low to high frequency. Follow-up: can coherence alone be the loss? Yes, it is smooth but flat near the optimum and biased toward the path average.

I10Write the error budget for a liver sound-speed number reported to a hepatologist. · research

Random: curvature or phase CRLB divided by √M (sub-m/s after averaging). Biases: layer bias if an average is reported (17–45 m/s under 10 mm of fat at 30–80 mm), forward-model error (refraction, straight rays, echo position: tens of m/s in tomography), calibration (≈ 0.8 m/s per ns of edge delay error for 64 elements at 30 mm; lens speed and thickness), regularization and the a priori mean. Compare with 14 m/s between the Imbault cut-offs. Follow-up: how do you validate? Phantoms with through-transmission references (±5 m/s), then patients against MRI-PDFF.

I11Would you trust a deep network that outputs a sound-speed map from channel data? · research

Only inside its training distribution. Feigin et al. (2020) trained on 6,026 k-Wave simulations of random ellipses and reported 20.5 m/s RMSE on simulated test data, with qualitative results in vivo: the forward model is now the simulator, and anatomy, lens, attenuation and noise mismatches become an unseen domain gap. I would want physics anchoring (unrolled or model-based inversion), uncertainty estimates and benchmark results on independent data such as QUASOS. Follow-up: what is a good sanity test? Feed a homogeneous phantom with known speed and a layered phantom; check the mean and the layer bias.

32.12MATLAB lab

Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run. The numbers below are the Octave output.

Download the scripts: m32_global_sos_focus.m · m32_layered_inversion.m · m32_cute_phase_shift.m · m32_sos_tomography.m · m32_sos_error_budget.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m32_global_sos_focus.mSnell-traced channel data under a fat layer, CF versus assumed speed at 9 depths, V_h versus V_rms, Dix and straight-ray interval speeds · change c_fat, d_fat, c_liv, Nel · runtime ≈ 2 s

%% m32_global_sos_focus.m
% Global sound speed by focus-quality maximization in a two-layer medium.
% Purpose : Simulate channel data (0 deg plane-wave transmit, 64-element
%           receive, 5 MHz) from wire targets under a 10 mm fat layer
%           (1450 m/s) over liver (1580 m/s). Receive times are ray-traced
%           with Snell's law (refraction included). For each target, sweep
%           the assumed speed c_hat, beamform a small ROI tied to two-way
%           time, and keep the c_hat that maximizes the coherence factor.
%           Compare with two averages of the true profile:
%             harmonic (time-average) V_h = z / T,   T = sum d_i/c_i
%             RMS (Dix)                V_rms = sqrt( sum c_i^2 dT_i / T )
%           Then invert the measured averages to interval (local) speeds
%           with the straight-ray rule (dz/dT) and with Dix (d(V^2 T)/dT).
% Change  : c_fat, d_fat, c_liv, f0, the receive aperture Nel.
% Expect  : best c_hat follows V_rms within about 1 m/s (0.6-1.0 m/s
%           high) and sits 1.5-2.3 m/s above V_h; V_rms exceeds V_h by
%           w1 w2 (c2 - c1)^2 / (2 V_h), about 1.1 m/s at 40 mm; the interval
%           speed below the fat is recovered to < 1 m/s by Dix and to
%           0.5-2.5 m/s by the straight-ray rule; the shallow fat target
%           gives 1450.
% Runtime : about 3 s in Octave.
clear; close all; rng(1);

c_fat = 1450; d_fat = 10e-3; c_liv = 1580;   % layer speeds (m/s), fat thickness (m)
f0 = 5e6; fs = 40e6; sig_t = 0.12e-6;        % centre frequency, sampling, Gaussian envelope std (s)
Nel = 64; pitch = 0.3e-3;                    % receive elements, pitch (m)
xe = ((0:Nel-1) - (Nel-1)/2)*pitch;          % element positions (m)
zt = [6 15 20 25 30 35 40 45 50]*1e-3;       % wire depths on axis (m)
nT = numel(zt);

% ---- ray tracing: one-way time from element (xe,0) to wire (0,z), Snell refraction
Trec = zeros(nT, Nel); Ttx = zeros(nT, 1);
for k = 1:nT
    z = zt(k);
    if z <= d_fat
        Trec(k, :) = sqrt(z^2 + xe.^2)/c_fat; Ttx(k) = z/c_fat; continue;
    end
    d2 = z - d_fat; X = abs(xe);
    lo = zeros(size(X)); hi = (1/c_liv)*(1 - 1e-12)*ones(size(X));   % ray parameter p (s/m)
    for it = 1:60
        p = (lo + hi)/2;
        f = d_fat*p*c_fat./sqrt(1 - (p*c_fat).^2) + d2*p*c_liv./sqrt(1 - (p*c_liv).^2);
        hi(f > X) = p(f > X); lo(f <= X) = p(f <= X);
    end
    p = (lo + hi)/2;
    Trec(k, :) = d_fat./(c_fat*sqrt(1 - (p*c_fat).^2)) + d2./(c_liv*sqrt(1 - (p*c_liv).^2));
    Ttx(k) = d_fat/c_fat + d2/c_liv;         % 0 deg plane wave: vertical path through flat layers
end

% ---- channel data as baseband IQ (analytic RF times exp(-j w0 t))
t = (0:round(80e-6*fs))/fs;                  % fast time (s)
IQ = zeros(numel(t), Nel);
for k = 1:nT
    tau = Ttx(k) + Trec(k, :);               % two-way arrival per element (s)
    D = bsxfun(@minus, t(:), tau);           % time relative to arrival
    IQ = IQ + exp(-D.^2/(2*sig_t^2)).*exp(-1j*2*pi*f0*repmat(tau, numel(t), 1));
end
IQ = IQ + 0.01*(randn(size(IQ)) + 1j*randn(size(IQ)))/sqrt(2);   % -40 dB noise

% ---- sweep c_hat, coherence factor over a time-tied ROI around each wire
chat = 1400:1:1650;
CF = zeros(nT, numel(chat));
Nt = numel(t);
xr = (-0.1:0.05:0.1)*1e-3;                   % ROI lateral positions (m), inside the PSF main lobe
for k = 1:nT
    t0 = Ttx(k) + min(Trec(k, :));           % two-way apex time of the wire (s)
    tr = t0 + (-0.05:0.01:0.05)*1e-6;         % ROI in two-way time (s), +-39 um
    [TR, XR] = ndgrid(tr, xr); TR = TR(:); XR = XR(:);   % pixels, column vectors
    for q = 1:numel(chat)
        ch = chat(q); zh = ch*TR/2;          % pixel depth tied to time
        tau = zh/ch + sqrt(zh.^2 + bsxfun(@minus, XR, xe).^2)/ch;   % 0 deg tx + rx delay, Npix x Nel
        u = tau*fs + 1; i0 = floor(u); fr = u - i0; i0 = min(max(i0, 1), Nt - 1);
        lin = bsxfun(@plus, i0, (0:Nel-1)*Nt);                    % linear indices into IQ
        s = (1 - fr).*IQ(lin) + fr.*IQ(lin + 1);                    % linear interpolation of baseband IQ
        s = s.*exp(1j*2*pi*f0*tau);          % remodulate to the focused phase
        CF(k, q) = sum(abs(sum(s, 2)).^2)/(Nel*sum(abs(s(:)).^2));  % energy-weighted CF over the ROI
    end
end

% ---- results versus the two averages
Vh = zeros(nT, 1); Vr = zeros(nT, 1); best = zeros(nT, 1); T1 = zeros(nT, 1);
fprintf('  z (mm)   best c_hat   V_h (harm.)   V_rms (Dix)   best - V_h   best - V_rms   CF max\n');
for k = 1:nT
    z = zt(k); d1 = min(z, d_fat); d2 = max(0, z - d_fat);
    dT = [d1/c_fat, d2/c_liv]; T1(k) = sum(dT);
    Vh(k) = z/T1(k); Vr(k) = sqrt((c_fat^2*dT(1) + c_liv^2*dT(2))/T1(k));
    [m, i] = max(CF(k, :)); i = min(max(i, 2), numel(chat) - 1);
    y = CF(k, i-1:i+1); best(k) = chat(i) + 0.5*(y(1) - y(3))/(y(1) - 2*y(2) + y(3));   % parabolic refinement
    fprintf('  %5.0f   %9.1f   %10.1f   %10.1f   %9.1f   %11.1f   %6.3f\n', z*1e3, best(k), Vh(k), Vr(k), best(k) - Vh(k), best(k) - Vr(k), m);
end

% ---- interval speeds from the measured averages (one-way times T1 known from the echo times)
fprintf('\nInterval speed between successive wires (true liver %d m/s):\n', c_liv);
fprintf('  interval (mm)   straight-ray dz/dT   Dix sqrt(d(V^2 T)/dT)\n');
for k = 3:nT
    zs_k = best(k)*T1(k); zs_j = best(k-1)*T1(k-1);
    c_str = (zs_k - zs_j)/(T1(k) - T1(k-1));
    c_dix = sqrt((best(k)^2*T1(k) - best(k-1)^2*T1(k-1))/(T1(k) - T1(k-1)));
    fprintf('  %4.0f-%-4.0f      %10.1f          %10.1f\n', zt(k-1)*1e3, zt(k)*1e3, c_str, c_dix);
end
k40 = find(zt == 40e-3); w1 = (d_fat/c_fat)/T1(k40); w2 = 1 - w1;
fprintf('Prediction V_rms - V_h ~ w1 w2 (c_liv - c_fat)^2/(2 V_h) at 40 mm: %.2f m/s (exact %.2f)\n', ...
    w1*w2*(c_liv - c_fat)^2/(2*Vh(k40)), Vr(k40) - Vh(k40));

figure('Position', [50 50 1000 420]);
subplot(1, 2, 1);
plot(chat, CF(zt == 15e-3, :)/max(CF(zt == 15e-3, :)), 'b', chat, CF(zt == 30e-3, :)/max(CF(zt == 30e-3, :)), 'r', ...
     chat, CF(zt == 50e-3, :)/max(CF(zt == 50e-3, :)), 'k', 'LineWidth', 1.4);
xlabel('assumed sound speed c_{hat} (m/s)'); ylabel('CF / max'); legend('z = 15 mm', 'z = 30 mm', 'z = 50 mm', 'Location', 'southwest');
title('Coherence factor versus assumed speed'); grid on;
subplot(1, 2, 2);
plot(best, zt*1e3, 'ko', Vh, zt*1e3, 'b-', Vr, zt*1e3, 'r--', 'LineWidth', 1.4); set(gca, 'YDir', 'reverse');
xlabel('speed (m/s)'); ylabel('depth (mm)'); legend('CF-optimal c_{hat}', 'V_h = z/T', 'V_{rms} (Dix)', 'Location', 'southwest');
title('The focus-optimal speed is an average'); grid on;

What to look for: the CF-optimal speed is 1450.4 m/s at 6 mm (inside the fat), then 1493.1, 1514.5, … 1553.7 m/s from 15 to 50 mm, always 0.6–1.0 m/s above \(V_{\mathrm{rms}}\) and 1.5–2.3 m/s above \(V_h\); at 40 mm \(V_{\mathrm{rms}} - V_h\) = 1.07 m/s, exactly the prediction of (32.4). Interval speeds of the liver (true 1580): 1580.1–1580.6 m/s by Dix, 1580.5–1582.5 m/s by straight-ray differences.

m32_layered_inversion.mNoisy averages to local speeds: naive differencing (with Monte Carlo of its noise), Tikhonov with the discrepancy principle, layer-constrained least squares · change sig_c, dz, the profile · runtime ≈ 1 s

%% m32_layered_inversion.m
% From average to local sound speed in a layered medium (straight-ray model).
% Purpose : A global estimator (Section 32.2) returns the average speed down
%           to each depth, c_avg(z) = z / T(z), T(z) = int_0^z s(z') dz'.
%           In slowness this is linear: z * s_avg(z) = int_0^z s dz', so the
%           local slowness is the derivative of z * s_avg. Differencing two
%           noisy averages amplifies their noise by about sqrt(2) z / dz.
%           Compare three inversions of the same noisy c_avg(z):
%             (1) naive differencing on the measurement grid,
%             (2) Tikhonov with a first-difference penalty, lambda chosen by
%                 the discrepancy principle (residual = expected noise),
%             (3) layer-constrained least squares (interfaces known from the
%                 B-mode image, one unknown speed per layer).
% Change  : sig_c (m/s noise on each average), dz (depth step), the profile.
% Expect  : naive local error ~ 1.41 z/dz * sig_c (about 85 m/s at 30 mm for
%           sig_c = 2 m/s, dz = 1 mm); Tikhonov recovers the liver to a few
%           m/s but rounds the interfaces; the layer fit gives each layer to
%           +-1 to 3 m/s, worst for the thin muscle layer. Monte Carlo exceeds
%           the linear noise rule where errors reach ~10 % (1/s is nonlinear).
% Runtime : about 1 s in Octave.
clear; close all; rng(1);

dz = 1e-3; zg = (dz:dz:60e-3)';              % cell bottoms (m), 60 cells of 1 mm
n = numel(zg); zc = zg - dz/2;               % cell centres
iface = [8e-3 14e-3];                        % fat | muscle | liver interfaces (m)
cl = [1450 1585 1550];                       % fat, muscle, liver speeds (m/s)
lay = 1 + (zc > iface(1)) + (zc > iface(2)); % layer index of every cell
c_true = cl(lay)'; s_true = 1./c_true;       % local speed and slowness

km = (4:60)'; zk = zg(km); m = numel(km);    % averages measured at 4..60 mm
A = double(bsxfun(@le, (1:n), km))*dz;       % b_k = z_k * s_avg(z_k) = sum_{j<=k} s_j dz
cavg_true = zk./(A*s_true);
sig_c = 2;                                   % noise on each c_avg (m/s)
cavg = cavg_true + sig_c*randn(m, 1);
b = zk./cavg;                                % measured z * s_avg (s)

% (1) naive differencing
s_naive = diff(b)/dz; z_naive = zk(2:end) - dz/2; c_naive = 1./s_naive;

% predicted noise amplification of (1): sigma_s,k = sigma_savg sqrt(z_k^2 + z_{k-1}^2)/dz
sig_savg = sig_c./cavg_true.^2;
pred = sqrt(zk(2:end).^2 + zk(1:end-1).^2)/dz .* sig_savg(2:end) .* c_true(km(2:end)).^2;   % in m/s

% Monte Carlo check of (1)
R = 400; E = zeros(R, m-1);
for r = 1:R
    bb = zk./(cavg_true + sig_c*randn(m, 1));
    E(r, :) = (1./(diff(bb)/dz))' - c_true(km(2:end))';
end
mc = std(E)';
k30 = find(abs(z_naive - 29.5e-3) < 1e-6); k50 = find(abs(z_naive - 49.5e-3) < 1e-6);
fprintf('Naive differencing, std of local c: predicted %.1f / MC %.1f m/s at 30 mm; %.1f / %.1f at 50 mm (rule 1.41 z/dz sig_c: %.0f, %.0f)\n', ...
    pred(k30), mc(k30), pred(k50), mc(k50), sqrt(2)*30*sig_c, sqrt(2)*50*sig_c);

% (2) Tikhonov with first-difference penalty, discrepancy principle
D = diff(eye(n));
W = diag(1./(zk.*sig_savg));                 % whiten: std of b_k = z_k * sigma_savg
lams = logspace(1, 9, 120); res = zeros(size(lams));
for q = 1:numel(lams)
    s_q = (A'*W'*W*A + lams(q)^2*(D'*D)) \ (A'*W'*W*b);
    res(q) = norm(W*(A*s_q - b))^2;
end
qd = find(res <= m, 1, 'last'); lam = lams(qd);   % largest lambda with chi^2 <= number of data
s_tik = (A'*W'*W*A + lam^2*(D'*D)) \ (A'*W'*W*b); c_tik = 1./s_tik;
liv = zc > 20e-3;
fprintf('Tikhonov (discrepancy, lambda = %.2g): liver (20-60 mm) mean %.1f m/s, RMSE %.1f m/s; whole-profile RMSE %.1f m/s\n', ...
    lam, mean(c_tik(liv)), sqrt(mean((c_tik(liv) - c_true(liv)).^2)), sqrt(mean((c_tik - c_true).^2)));

% (3) layer-constrained least squares: s = P * s_layer, P(j,l) = 1 if cell j is in layer l
P = double(bsxfun(@eq, lay, 1:3));
B = W*A*P; s_lay = B \ (W*b); C = inv(B'*B);  % covariance of the layer slownesses
c_lay = 1./s_lay; sd_lay = sqrt(diag(C)).*c_lay.^2;
fprintf('Layer fit: fat %.1f +- %.1f, muscle %.1f +- %.1f, liver %.1f +- %.1f m/s (true %d, %d, %d)\n', ...
    c_lay(1), sd_lay(1), c_lay(2), sd_lay(2), c_lay(3), sd_lay(3), cl);
fprintf('Average at 60 mm: %.1f m/s, i.e. %.1f m/s below the liver value it is often mistaken for\n', ...
    cavg_true(end), cl(3) - cavg_true(end));

figure('Position', [40 40 1050 430]);
subplot(1, 2, 1);
plot(c_true, zc*1e3, 'k-', 'LineWidth', 2); hold on;
plot(cavg, zk*1e3, 'b.', 'MarkerSize', 10);
plot(c_naive, z_naive*1e3, 'Color', [0.8 0.5 0.5]);
plot(c_tik, zc*1e3, 'r-', 'LineWidth', 1.6);
stairs([c_lay(lay)' c_lay(3)], [0; zg]*1e3, 'g--', 'LineWidth', 1.4);
set(gca, 'YDir', 'reverse'); xlim([1300 1800]); ylim([0 60]); grid on;
xlabel('sound speed (m/s)'); ylabel('depth (mm)');
legend('true local', 'measured average', 'naive difference', 'Tikhonov', 'layer fit', 'Location', 'northeast');
title(sprintf('Average-to-local inversion, sigma_c = %g m/s', sig_c));
subplot(1, 2, 2);
plot(z_naive*1e3, mc, 'k.', z_naive*1e3, pred, 'r-', 'LineWidth', 1.4); grid on;
xlabel('depth (mm)'); ylabel('std of local c (m/s)'); legend('Monte Carlo', 'sqrt(z_k^2 + z_{k-1}^2)/dz rule', 'Location', 'northwest');
title('Noise amplification of naive differencing');

What to look for: naive local noise 85.8 m/s predicted and 87.1 m/s by Monte Carlo at 30 mm (142.4 and 147.6 at 50 mm) from 2 m/s per average; Tikhonov at λ = 1.5×10⁵ gives the liver at 1548.6 m/s (RMSE 6.8 m/s, 10.5 m/s over the whole profile); the layer fit gives 1449.3 ± 0.8, 1587.6 ± 2.5 and 1548.9 ± 0.7 m/s (true 1450, 1585, 1550); the average at 60 mm is 1539.2 m/s.

m32_cute_phase_shift.mStraight-ray CUTE delay maps (transmit-only, common mid-angle, 2° step) for an inclusion under fat, and a Monte Carlo of the phase jitter of correlated speckle · change th_deg, dc_inc, rho_list, N_list · runtime ≈ 3 s

%% m32_cute_phase_shift.m
% CUTE-style echo phase shifts between steered acquisitions (ray model) and
% the phase jitter that limits them (speckle Monte Carlo).
% Purpose : (A) Straight-ray, first-order model of the beamforming delay
%           error of every pixel r for a steering angle th:
%              L_th(r) = sec(th) * int_0^z ds(x - (z - z') tan th, z') dz',
%           ds = 1/c - 1/c0 (c0 = 1540 m/s assumed by the beamformer).
%           Two data types are mapped for a +40 m/s inclusion under a fat layer:
%              transmit-only pair   d = L_th - L_0             (receive fixed)
%              common mid-angle     d = L_th + L_-th - 2 L_0   (tx th, rx -th)
%           Echo-position shifts and refraction are neglected.
%           (B) The phase of the complex correlation of two speckle kernels of
%           N independent samples and correlation rho has std
%              sigma_phi ~ sqrt(1 - rho^2) / (rho sqrt(2N))   (large N),
%           checked by Monte Carlo; sigma_tau = sigma_phi / (2 pi f0).
% Change  : th_deg, dc_inc, r_inc, the fat layer, rho_list, N_list.
% Expect  : at 15 deg the shift below the inclusion centre (0 deg path through
%           the full chord, 15 deg path missing it) reaches the chord delay
%           (164 ns) plus 5.7 ns from the fat layer, ~169 ns, several phase wraps at 5 MHz; a 2 deg step
%           stays below 100 ns (no wrap), zero above the inclusion; the layer adds only (sec th - 1) of its
%           vertical delay; the jitter formula holds within ~10 % for N = 64
%           and fails at low rho and small N (phase wrapping, bias).
% Runtime : about 2 s in Octave.
clear; close all; rng(1);

c0 = 1540; f0 = 5e6;
dx = 0.1; x = -19.2:dx:19.2; z = 0:dx:40;    % grid (mm)
[X, Z] = meshgrid(x, z);
dc_inc = 40; r_inc = 5; zi = 15;             % inclusion: +40 m/s, radius 5 mm, centre depth 15 mm
c = c0 + dc_inc./(1 + exp((hypot(X, Z - zi) - r_inc)/0.3));   % smooth edge, 0.3 mm
d_fat = 4; c(Z < d_fat) = 1450;              % 4 mm fat layer at the top
ds = (1./c - 1/c0)*1e6;                      % slowness deviation (ns/mm)

th_deg = 15; ths = [-th_deg 0 th_deg 2]*pi/180; L = cell(1, 4);
for a = 1:4
    th = ths(a); La = zeros(size(ds));
    La(1, :) = ds(1, :)*dx/cos(th)/2;        % half cell at the surface
    for k = 2:numel(z)                       % march down: shift the previous row by dz tan(th)
        prev = interp1(x, La(k-1, :), x - dx*tan(th), 'linear', 'extrap');
        La(k, :) = prev + 0.5*(ds(k, :) + ds(k-1, :))*dx/cos(th);
    end
    L{a} = La;                               % ns, delay error of each pixel for angle th
end
d_tx  = L{3} - L{2};                         % transmit-only pair (th vs 0)
d_cma = L{3} + L{1} - 2*L{2};                % common mid-angle pair (th,-th) vs (0,0)
in_ap = abs(X) < 12;                         % keep away from the lateral edges
fprintf('Inclusion %+d m/s: ds = %.2f ns/mm, vertical delay through its centre %.1f ns\n', dc_inc, ...
    (1/(c0 + dc_inc) - 1/c0)*1e6, (1/(c0 + dc_inc) - 1/c0)*1e6*2*r_inc);
fprintf('Peak |d| at th = %d deg: transmit-only %.1f ns (%.2f rad at %.0f MHz), common mid-angle %.1f ns (%.2f rad)\n', th_deg, ...
    max(abs(d_tx(in_ap))), 2*pi*f0*max(abs(d_tx(in_ap)))*1e-9, f0/1e6, max(abs(d_cma(in_ap))), 2*pi*f0*max(abs(d_cma(in_ap)))*1e-9);
d2 = L{4} - L{2};
fprintf('Adjacent-angle step of 2 deg (tx-only): peak |d| %.1f ns (%.2f rad): no wrapping; wrapping starts at |d| = 1/(2 f0) = %.0f ns\n', ...
    max(abs(d2(in_ap))), 2*pi*f0*max(abs(d2(in_ap)))*1e-9, 1e9/(2*f0));
k30 = find(abs(z - 30) < 1e-9); i0 = find(abs(x + 15) < 1e-9);   % a column far from the inclusion
S = (1/1450 - 1/c0)*1e6*d_fat;
fprintf('Fat layer only (x = -15 mm, z = 30 mm): vertical %.1f ns, tx-only %.2f ns (pred (sec-1)S = %.2f), CMA %.2f ns (pred 2(sec-1)S = %.2f)\n', ...
    L{2}(k30, i0), d_tx(k30, i0), (sec(ths(3)) - 1)*S, d_cma(k30, i0), 2*(sec(ths(3)) - 1)*S);

% (B) phase jitter of the complex correlation of two speckle kernels
rho_list = [0.5 0.7 0.9 0.95 0.99]; N_list = [16 64]; R = 4000;
fprintf('\n  rho     N   sigma_phi MC (rad)   formula   sigma_tau at 5 MHz (ns)\n');
SPH = zeros(numel(N_list), numel(rho_list));
for iN = 1:numel(N_list)
    N = N_list(iN);
    for ir = 1:numel(rho_list)
        rho = rho_list(ir); phi0 = 0.4;
        a1 = (randn(N, R) + 1j*randn(N, R))/sqrt(2);
        w  = (randn(N, R) + 1j*randn(N, R))/sqrt(2);
        a2 = (rho*a1 + sqrt(1 - rho^2)*w)*exp(1j*phi0);
        ph = angle(sum(a2.*conj(a1), 1)) - phi0;
        SPH(iN, ir) = std(ph);
        fo = sqrt(1 - rho^2)/(rho*sqrt(2*N));
        fprintf('  %.2f  %3d   %10.4f        %8.4f   %8.2f\n', rho, N, SPH(iN, ir), fo, SPH(iN, ir)/(2*pi*f0)*1e9);
    end
end
th = 20*pi/180; dzl = 5; sig_tau = sqrt(1 - 0.95^2)/(0.95*sqrt(128))/(2*pi*f0)*1e9;
% layer delay S = datum below minus datum above the layer: two independent data, sqrt(2) sigma_tau
sig_s = sqrt(2)*sig_tau/(2*(sec(th) - 1)*dzl);
fprintf('\nLayer slowness from CMA data below minus above a %d mm layer at +-20 deg, sigma_tau = %.2f ns per datum (rho 0.95, N 64): sigma_s = %.1f ns/mm = %.0f m/s\n', ...
    dzl, sig_tau, sig_s, sig_s*1e-6*c0^2);

figure('Position', [40 40 1150 380]); colormap(jet);
subplot(1, 3, 1); imagesc(x, z, c); axis image; colorbar; title('true c (m/s)'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(1, 3, 2); imagesc(x, z, d_tx); axis image; colorbar; caxis([-180 180]);
title(sprintf('tx-only: L_{%d} - L_0 (ns)', th_deg)); xlabel('x (mm)');
subplot(1, 3, 3); imagesc(x, z, d_cma); axis image; colorbar; caxis([-180 180]);
title(sprintf('CMA: L_{%d} + L_{-%d} - 2L_0 (ns)', th_deg, th_deg)); xlabel('x (mm)');
figure('Position', [60 60 560 380]);
semilogy(rho_list, SPH(1, :), 'bo', rho_list, SPH(2, :), 'ro', 'LineWidth', 1.4); hold on;
rr = linspace(0.45, 0.995, 200);
semilogy(rr, sqrt(1 - rr.^2)./(rr*sqrt(2*16)), 'b-', rr, sqrt(1 - rr.^2)./(rr*sqrt(2*64)), 'r-');
grid on; xlabel('correlation coefficient rho'); ylabel('phase std (rad)');
legend('MC N = 16', 'MC N = 64', 'formula N = 16', 'formula N = 64', 'Location', 'southwest');
title('Phase jitter of the complex correlation');

What to look for: inclusion slowness −16.44 ns/mm; peak shift at 15° of 168.9 ns (5.31 rad) transmit-only and 335.4 ns (10.54 rad) common mid-angle; 56.2 ns for a 2° step (no wrap below 100 ns); the fat layer alone gives 5.69 and 11.37 ns, matching \((\sec\theta - 1)S\) and \(2(\sec\theta - 1)S\). Phase jitter: 0.0295 rad by Monte Carlo against 0.0291 from (32.10) at ρ = 0.95, N = 64; at ρ = 0.5, N = 16 the formula underestimates (0.306 versus 0.355 rad). A 5 mm layer from CMA data at ±20° with 0.92 ns jitter on the data above and below it (√2 × 0.92 ns on their difference): σ_s = 2.0 ns/mm = 5 m/s.

m32_sos_tomography.mLimited-angle echo-mode tomography with transmit-only and CMA data, Tikhonov with gradient penalty, resolution-matrix kernels, singular spectra, and an error budget that separates noise from model error · change thmax, K, noise_ns, the truth · runtime ≈ 6 s

%% m32_sos_tomography.m
% Limited-angle echo-mode SoS tomography: two acquisition designs, the
% resolution matrix versus depth, and the singular spectrum.
% Purpose : Unknown: slowness deviation ds = 1/c - 1/c0 on 16 x 14 cells of
%           2.4 mm (x = +-19.2 mm, z = 2 to 35.6 mm) under a 38.4 mm array.
%           Data at every cell centre, built from straight-ray delay errors
%           L_phi(r) (line integral of ds from r to the array at angle phi):
%             'tx-only'  transmit steered, receive fixed at 0 deg; adjacent
%                        transmit angles differenced: d = L_phi(k) - L_phi(k-1)
%             'CMA'      transmit phi_i and receive phi_j from the same angle set;
%                        pairs grouped by mid-angle (phi_i + phi_j)/2 and
%                        adjacent pairs in a group differenced
%           Tikhonov with a gradient penalty; oracle lambda (minimum RMSE
%           against the truth), the best case for each design. R = (A'A + lam^2 G'G)^-1 A'A: row k is the point-spread
%           function of the reconstruction at cell k.
%           Data come from 0.05 mm ray steps through a smooth truth; the
%           inversion uses 0.3 mm steps with nearest-cell weights.
% Change  : thmax, K, noise_ns, the truth.
% Expect  : CMA gives about 2.2x more rows, yet in this straight-ray model
%           every CMA row (W_i + W_j) - (W_i-1 + W_j+1) is a combination of
%           the tx-only rows of the same pixel (printed residual ~1e-15): CMA
%           adds no ray information. Its real-data advantage (higher speckle
%           correlation, no echo-position bias) is outside this model.
%           The error budget is printed per design: with data consistent with
%           the cell model the error comes from noise + regularization; with
%           fine-ray data it is dominated by model (discretization) error,
%           which is larger than the 2 ns noise effect. Both designs smear
%           axially more than laterally and worse with depth; the 5 mm fat
%           layer comes back near 1506 m/s instead of 1450.
% Runtime : about 6 s in Octave.
clear; close all; rng(2);

c0 = 1540; Nx = 16; Nz = 14; h = 2.4;        % reference speed, grid, cell size (mm)
xc = ((0:Nx-1) - (Nx-1)/2)*h; zc = 2 + ((0:Nz-1) + 0.5)*h;
np = Nx*Nz; Xh = 19.2;                       % unknowns, array half-width (mm)
sg = @(r, r0) 1./(1 + exp((r - r0)/0.4));
truth = @(x, z) 1540 + 40*sg(hypot(x - 2, z - 20), 4.5) - 90*(z < 5);   % inclusion + 5 mm fat (1450)
dsf = @(c) (1./c - 1/c0)*1e6;                % ns/mm
[XX, ZZ] = meshgrid(xc, zc); ct = reshape(truth(XX, ZZ)', [], 1);  % index = iz*Nx + ix + 1
cellof = @(x, z) min(max(floor((z - 2)/h), 0), Nz-1)*Nx + min(max(floor((x + Xh)/h), 0), Nx-1) + 1;

% straight-ray integral from point (x0,z0) to the surface at angle phi
% (weights over cells for the inversion; fine numerical value for the data)
ray_w = @(x0, z0, phi, st) accumarray(cellof(x0 - (0:st:z0/cos(phi))'*sin(phi), z0 - (0:st:z0/cos(phi))'*cos(phi)), st, [np 1])';
ray_d = @(x0, z0, phi) sum(dsf(truth(x0 - (0:0.05:z0/cos(phi))*sin(phi), z0 - (0:0.05:z0/cos(phi))*cos(phi))))*0.05;

% gradient penalty
I = speye(np); G = [];
for iz = 0:Nz-1, for ix = 0:Nx-2, k = iz*Nx + ix + 1; G = [G; I(k+1, :) - I(k, :)]; end, end %#ok<AGROW>
for iz = 0:Nz-2, for ix = 0:Nx-1, k = iz*Nx + ix + 1; G = [G; I(k+Nx, :) - I(k, :)]; end, end %#ok<AGROW>
GG = full(G'*G);

thmax = 25; K = 7; noise_ns = 2; phis = linspace(-thmax, thmax, K)*pi/180;
designs = {'tx-only', 'CMA'}; fw = @(p, dd) dd*0.01*sum(interp1(1:numel(p), p, 1:0.01:numel(p)) >= max(p)/2);
fprintf('theta_max = %d deg, K = %d angles, delay noise %.1f ns\n', thmax, K, noise_ns);
fprintf('  design    rows   lambda    trace(R)  RMSE (m/s)  incl. mean (true %.0f)  fat mean (true 1450)   axial / lateral FWHM (mm) at z = 9, 19, 31 mm\n', truth(2, 20));
res = struct(); rowres = 0;
for dsg = 1:2
    Arows = []; d = [];
    for iz = 1:Nz
        for ix = 1:Nx
            x0 = xc(ix); z0 = zc(iz);
            ok = abs(x0 - z0*tan(phis)) <= Xh;            % ray must reach the aperture
            Wp = zeros(K, np); Dp = zeros(K, 1);
            for a = 1:K
                if ok(a), Wp(a, :) = ray_w(x0, z0, phis(a), 0.3); Dp(a) = ray_d(x0, z0, phis(a)); end
            end
            if dsg == 1                                     % tx-only: receive fixed at 0 deg (cancels)
                for a = 2:K
                    if ok(a) && ok(a-1), Arows = [Arows; Wp(a, :) - Wp(a-1, :)]; d = [d; Dp(a) - Dp(a-1)]; end %#ok<AGROW>
                end
            else                                            % CMA: tx i, rx j, group by mid-angle
                Tp = []; for a = 2:K, if ok(a) && ok(a-1), Tp = [Tp; Wp(a, :) - Wp(a-1, :)]; end, end %#ok<AGROW>
                Cp = [];
                for gsum = 2:2*K
                    pr = [];
                    for i = max(1, gsum-K):min(K, gsum-1), j = gsum - i; if i <= j && ok(i) && ok(j), pr = [pr; i j]; end, end %#ok<AGROW>
                    for q = 2:size(pr, 1)
                        w1 = Wp(pr(q, 1), :) + Wp(pr(q, 2), :); w0 = Wp(pr(q-1, 1), :) + Wp(pr(q-1, 2), :);
                        Arows = [Arows; w1 - w0]; d = [d; Dp(pr(q, 1)) + Dp(pr(q, 2)) - Dp(pr(q-1, 1)) - Dp(pr(q-1, 2))]; %#ok<AGROW>
                        Cp = [Cp; w1 - w0]; %#ok<AGROW>
                    end
                end
                if ~isempty(Cp) && ~isempty(Tp), rowres = max(rowres, norm(Cp - (Cp*pinv(Tp))*Tp, 'fro')/norm(Cp, 'fro')); end
            end
        end
    end
    m = size(Arows, 1); nse = noise_ns*randn(m, 1);
    AA = Arows'*Arows; lams = logspace(-1, 4, 50);
    orac = @(dd) min(arrayfun(@(L) sqrt(mean((1./(1/c0 + ((AA + L^2*GG) \ (Arows'*dd))*1e-6) - ct).^2)), lams));
    rmse_cons = orac(Arows*dsf(ct) + nse);          % data consistent with the cell model: noise + regularization only
    rmse_nf = orac(d);                              % fine-ray data without noise: model error + regularization only
    dn = d + nse; Ab = Arows'*dn; rr = zeros(size(lams));
    for q = 1:numel(lams), sq = (AA + lams(q)^2*GG) \ Ab; rr(q) = sqrt(mean((1./(1/c0 + sq*1e-6) - ct).^2)); end
    [~, qd] = min(rr); lam = lams(qd);              % oracle lambda (minimum RMSE against the truth)
    ds_hat = (AA + lam^2*GG) \ Ab; c_hat = 1./(1/c0 + ds_hat*1e-6);
    Rm = (AA + lam^2*GG) \ AA;
    inc = hypot(reshape(XX', [], 1) - 2, reshape(ZZ', [], 1) - 20) < 3;
    fat = reshape(ZZ', [], 1) < 5;
    ixc = find(abs(xc - 1.2) < 1e-9);
    fwa = zeros(1, 3); fwl = zeros(1, 3); zq = [9 19 31];
    for t = 1:3
        iz = find(abs(zc - zq(t)) < 1.3, 1); k = (iz - 1)*Nx + ixc;
        row = reshape(Rm(k, :), Nx, Nz)';           % Nz x Nx kernel
        fwa(t) = fw(max(row(:, ixc), 0)', h); fwl(t) = fw(max(row(iz, :), 0), h);
    end
    fprintf('  %-8s %5d   %6.2f   %7.1f    %7.1f      %7.1f                %7.1f            %4.1f/%3.1f  %4.1f/%3.1f  %4.1f/%3.1f\n', designs{dsg}, m, lam, trace(Rm), ...
        sqrt(mean((c_hat - ct).^2)), mean(c_hat(inc)), mean(c_hat(fat)), fwa(1), fwl(1), fwa(2), fwl(2), fwa(3), fwl(3));
    fprintf('            RMSE with cell-consistent data + noise: %.1f m/s; fine-ray data without noise: %.1f m/s\n', rmse_cons, rmse_nf);
    sv = svd(Arows); res(dsg).sv = sv/sv(1); res(dsg).c = c_hat; res(dsg).R = Rm;
end
fprintf('Singular values above 1e-2 of the largest: tx-only %d, CMA %d (of %d unknowns)\n', sum(res(1).sv > 1e-2), sum(res(2).sv > 1e-2), np);
fprintf('CMA rows of each pixel lie in the span of its tx-only rows: max relative residual %.1e\n', rowres);

figure('Position', [30 30 1150 640]); colormap(jet);
subplot(2, 3, 1); imagesc(xc, zc, reshape(ct, Nx, Nz)'); axis image; caxis([1440 1590]); colorbar; title('truth (m/s)'); ylabel('z (mm)');
subplot(2, 3, 2); imagesc(xc, zc, reshape(res(1).c, Nx, Nz)'); axis image; caxis([1440 1590]); colorbar; title('tx-only, Tikhonov');
subplot(2, 3, 3); imagesc(xc, zc, reshape(res(2).c, Nx, Nz)'); axis image; caxis([1440 1590]); colorbar; title('CMA, Tikhonov');
k = (find(abs(zc - 19) < 1.3, 1) - 1)*Nx + find(abs(xc - 1.2) < 1e-9);
subplot(2, 3, 4); semilogy(res(1).sv, 'b', 'LineWidth', 1.4); hold on; semilogy(res(2).sv, 'r', 'LineWidth', 1.4); grid on;
xlabel('index'); ylabel('\sigma_i / \sigma_1'); legend('tx-only', 'CMA', 'Location', 'southwest'); title('singular spectrum of A');
subplot(2, 3, 5); imagesc(xc, zc, reshape(res(1).R(k, :), Nx, Nz)'); axis image; colorbar; title('R row at (1.2, 19.4) mm, tx-only'); xlabel('x (mm)');
subplot(2, 3, 6); imagesc(xc, zc, reshape(res(2).R(k, :), Nx, Nz)'); axis image; colorbar; title('R row at (1.2, 19.4) mm, CMA'); xlabel('x (mm)');

What to look for: 1146 rows (transmit-only) and 2498 rows (CMA) for 224 unknowns; oracle-λ RMSE 22.2 and 30.6 m/s, inclusion mean 1558.6 and 1535.5 m/s (true 1580), fat row about 1506 m/s (true 1450) for both; axial/lateral FWHM 2.5/2.4, 2.6/2.4, 3.4/2.4 mm (transmit-only) and 2.6/2.4, 2.9/2.4, 4.2/2.4 mm (CMA) at 9, 19, 31 mm. The error budget: 3.8 and 1.9 m/s with cell-consistent data plus noise, 22.0 and 30.8 m/s with noise-free fine-ray data. Each pixel's CMA rows lie in the span of its transmit-only rows (residual 5.1×10⁻¹⁵). An earlier version of this script's header claimed CMA would give a larger trace(R) and a smaller axial kernel; the printed numbers show the opposite, and the header now explains why.

m32_sos_error_budget.mFisher information and CRLB of a wavefront-curvature speed estimate, checked by Monte Carlo; CRLB table versus depth and aperture; bias of a quadratic calibration error · change sigma_t, N, z_list · runtime ≈ 2 s

%% m32_sos_error_budget.m
% Error budget of a global sound-speed estimate from wavefront curvature:
% how per-element delay jitter becomes m/s (Fisher information / CRLB).
% Purpose : A point (or a coherent speckle patch) at depth z returns echoes
%           whose receive arrival times across the aperture follow
%              t_i = t_a + ( sqrt(z^2 + x_i^2) - z ) / c,   z = c t_a,
%           with t_a the apex (on-axis) one-way time, measured from the echo.
%           Unknowns: c and an arrival offset t0 (common delay). The
%           parabolic (small-aperture) limit gives a closed form
%              sigma_c = c^2 z sigma_t sqrt(45/(4N)) / X^2
%                      = 4 c^2 F^2 sigma_t sqrt(45/(4N)) / z,   F = z/(2X),
%           for N elements uniformly spread over [-X, X] and independent
%           per-element jitter sigma_t. (A) checks it against the exact
%           Fisher matrix of the hyperbolic model and a Monte Carlo
%           least-squares fit; (B) tabulates sigma_c against depth and
%           aperture; (C) shows the bias from ignoring a 1 ns/element
%           systematic delay error (e.g. lens or element delay calibration).
% Change  : sigma_t, N, pitch, z_list, the systematic error pattern.
% Expect  : closed form, exact CRLB and Monte Carlo agree within ~5 % (1.62,
%           1.69, 1.68 m/s for 64 elements, 5 ns jitter, z = 30 mm);
%           sigma_c grows linearly with z at a fixed aperture and falls as
%           N^-2.5 at fixed pitch (aperture^2 times sqrt(N)); a 1 ns
%           quadratic calibration error at the aperture edge biases c by
%           about -0.8 m/s (0.05 % of c per ns of the ~970 ns curvature),
%           comparable to the CRLB and, unlike jitter, it never averages out.
% Runtime : about 2 s in Octave.
clear; close all; rng(3);

c = 1540; pitch = 0.3e-3;                    % true speed (m/s), pitch (m)
sigma_t = 5e-9;                              % per-element arrival-time jitter (s)
model = @(cc, t0, x, ta) t0 + ta + (sqrt((cc*ta)^2 + x.^2) - cc*ta)/cc;   % one-way receive times (s)

% (A) one configuration: 64 elements, z = 30 mm
N = 64; x = ((0:N-1) - (N-1)/2)*pitch; X = max(abs(x)); z = 30e-3; ta = z/c;
h = 1e-3;                                     % finite-difference step for the Jacobian (m/s and s)
J = [(model(c + h, 0, x, ta) - model(c - h, 0, x, ta))'/(2*h), ones(N, 1)];   % d t / d(c, t0)
F = (J'*J)/sigma_t^2; Ci = inv(F); crlb = sqrt(Ci(1, 1));
Xrms = sqrt(mean(x.^4) - mean(x.^2)^2);      % discrete spread of x^2 (continuous: X^2 sqrt(4/45))
closed = c^2*z*sigma_t/(sqrt(N)*Xrms);
closed_u = c^2*z*sigma_t*sqrt(45/(4*N))/X^2;
R = 3000; est = zeros(R, 1);
for r = 1:R
    tm = model(c, 0, x, ta) + sigma_t*randn(1, N);
    p = [c; 0];                              % Gauss-Newton on (c, t0)
    for it = 1:6
        res = tm' - model(p(1), p(2), x, ta)';
        Jg = [(model(p(1) + h, p(2), x, ta) - model(p(1) - h, p(2), x, ta))'/(2*h), ones(N, 1)];
        p = p + Jg \ res;
    end
    est(r) = p(1);
end
fprintf('(A) N = %d, aperture %.1f mm, z = %.0f mm (F = %.2f), sigma_t = %.0f ns\n', N, 2*X*1e3, z*1e3, z/(2*X), sigma_t*1e9);
fprintf('    exact CRLB %.2f m/s | parabolic closed form %.2f (uniform-aperture form %.2f) | Monte Carlo std %.2f, bias %+.2f m/s\n', ...
    crlb, closed, closed_u, std(est), mean(est) - c);

% (B) table: sigma_c versus depth and aperture (exact CRLB)
z_list = [10 20 30 40 60 80]*1e-3; N_list = [32 64 128];
fprintf('\n(B) CRLB of c (m/s), sigma_t = %.0f ns\n   z (mm) ', sigma_t*1e9); fprintf('   N = %3d (F)  ', N_list); fprintf('\n');
S = zeros(numel(z_list), numel(N_list));
for iz = 1:numel(z_list)
    fprintf('   %5.0f  ', z_list(iz)*1e3);
    for in = 1:numel(N_list)
        Nn = N_list(in); xx = ((0:Nn-1) - (Nn-1)/2)*pitch; tz = z_list(iz)/c;
        Jn = [(model(c + h, 0, xx, tz) - model(c - h, 0, xx, tz))'/(2*h), ones(Nn, 1)];
        Ci = inv(Jn'*Jn)*sigma_t^2; S(iz, in) = sqrt(Ci(1, 1));
        fprintf('   %7.2f (%4.2f)', S(iz, in), z_list(iz)/(2*max(abs(xx))));
    end
    fprintf('\n');
end

% (C) a systematic error: quadratic delay error of 1 ns at the aperture edge
e_sys = 1e-9*(x/X).^2;                       % s, e.g. lens/element delay calibration
p = [c; 0]; tm = model(c, 0, x, ta) + e_sys;
for it = 1:8
    res = tm' - model(p(1), p(2), x, ta)';
    Jg = [(model(p(1) + h, p(2), x, ta) - model(p(1) - h, p(2), x, ta))'/(2*h), ones(N, 1)];
    p = p + Jg \ res;
end
fprintf('\n(C) 1 ns quadratic calibration error at the aperture edge: bias %+.2f m/s (CRLB at 5 ns jitter: %.2f m/s)\n', p(1) - c, crlb);
fprintf('    rule: the curvature of a 1540 m/s wavefront across this aperture is X^2/(2 z c) = %.0f ns; 1 ns of it is %.2f %% of c\n', ...
    X^2/(2*z*c)*1e9, 100*1e-9/(X^2/(2*z*c))/2);

figure('Position', [60 60 980 380]);
subplot(1, 2, 1); hist(est, 40); xlabel('estimated c (m/s)'); ylabel('count');
title(sprintf('Monte Carlo, std %.2f m/s, CRLB %.2f m/s', std(est), crlb)); grid on;
subplot(1, 2, 2); loglog(z_list*1e3, S, 'o-', 'LineWidth', 1.4); grid on;
xlabel('depth z (mm)'); ylabel('CRLB of c (m/s)'); legend('N = 32', 'N = 64', 'N = 128', 'Location', 'northwest');
title(sprintf('Wavefront-curvature CRLB, sigma_t = %.0f ns', sigma_t*1e9));

What to look for: 64 elements, 30 mm, 5 ns: exact CRLB 1.69 m/s, closed form 1.62 m/s (1.67 continuous), Monte Carlo 1.68 m/s with +0.02 m/s bias. The table runs from 0.20 m/s (128 elements, 10 mm) to 24.5 m/s (32 elements, 80 mm), linear in depth and about \(N^{-2.5}\) in aperture. A 1 ns quadratic edge error biases the estimate by −0.83 m/s.

32.13Self-check

Q1At 7.5 MHz, F/1.5 and 25 mm, what speed error puts the edge phase at π/4?

\(\lvert\Delta c\rvert = c^{2}F^{2}/(2f_0z)\) = 2.37×10⁶ × 2.25 / (2 × 7.5×10⁶ × 0.025) = 14.2 m/s.

Q2A region at 1600 m/s is beamformed at 1540 m/s. Slowness deviation, and the linearized speed?

\(\Delta s\) = 1/1600 − 1/1540 = −24.35 ns/mm; \(-c_0^{2}\Delta s\) = 57.7 m/s against the true 60 m/s (error \(\Delta c^{2}/c\) ≈ 2.3 m/s).

Q315 mm of fat (1450 m/s) over liver (1560 m/s). V_h and V_rms at 60 mm?

\(\Delta T\) = 10.34 and 28.85 µs, T = 39.19 µs, \(w_1\) = 0.264. \(V_h\) = 60/39.19 = 1531.0 m/s; \(V_{\mathrm{rms}}\) = √(0.264 × 1450² + 0.736 × 1560²) = 1531.7 m/s. Both are about 29 m/s below the liver.

Q4Averages at 40 and 45 mm, each with 1 m/s noise. Noise on the local speed between them?

\(\sqrt{45^{2} + 40^{2}}/5\) × 1 = 12.0 m/s.

Q5A 3 mm-radius region at −30 m/s (relative to 1540). Maximum vertical chord delay, and does a 10° transmit-only step wrap at 6 MHz?

\(\Delta s\) = 1/1510 − 1/1540 = 12.9 ns/mm. The longest chord through a circle is its diameter in any direction, 6 mm, so no datum can exceed 77 ns, even if one path crosses the full chord and the other misses the region. Wrapping starts at 1/(2f₀) = 83 ns at 6 MHz: no wrap, whatever the step.

Q6CUTE kernel with N = 100 independent samples, ρ = 0.9, 5 MHz. Delay jitter?

\(\sigma_\phi\) = √0.19 / (0.9 × √200) = 0.0342 rad; \(\sigma_\tau\) = 0.0342/(2π × 5×10⁶) = 1.09 ns.

Q7Curvature fit with 128 elements (pitch 0.3 mm), 60 mm, 10 ns jitter. CRLB, and how many independent kernels to reach 0.5 m/s?

From the table of m32_sos_error_budget.m, 0.60 m/s at 5 ns, so 1.2 m/s at 10 ns (linear in σ_t). M = (1.2/0.5)² ≈ 6 kernels.

32.14Figures from the literature

Measured data behind the models of this module. Figure numbers were checked against the open preprints named in each entry (the published versions may number them differently); links only, no images reproduced.

Common mid-angle and the echo-position term: phantom and volunteer SoS maps

P. Stähli, M. Kuriakose, M. Frenz and M. Jaeger, "Improved forward model for quantitative pulse-echo speed-of-sound imaging," Ultrasonics 108, 106168 (2020); preprint arXiv:1902.10639. doi:10.1016/j.ultras.2020.106168

Look at Fig. 1 (isochrones for different transmit/receive angle pairs: only pairs sharing a mid-angle keep their relative position), Figs. 4–6 (two-layer, four-layer and laterally varying phantoms with liver 1555 and fat 1420 m/s; rows: old model, CMA model, new model; columns: four a priori speeds) and Fig. 7 (volunteer abdominal wall and liver, old versus new model). Why: it is the measured version of Sections 32.4–32.5: the gain from the CMA pairing, the extra gain from the echo-position term, and the dependence on the a priori mean.

CUTE with a convex probe: liver phantom, volunteer and the a priori speed

M. Jaeger, P. Stähli, N. Korta Martiartu, P. Salemi Yolgunlu, T. Frappart, C. Fraschini and M. Frenz, "Pulse-echo speed-of-sound imaging using convex probes," Phys. Med. Biol. (2022); preprint arXiv:2202.10633. doi:10.1088/1361-6560/ac96c6

Look at Fig. 5 (four-layer liver phantom: B-mode, SoS maps for two beamforming grids and two regularizations, depth profiles), Fig. 6 (volunteer liver, about 1580 m/s superficially), Fig. 8 (images for maximum angle spreads of 22, 6 and 3 steps) and Figs. 9–10 (phantom and in vivo maps for a priori speeds from 1510 to 1590 m/s). Why: Fig. 8 is Demo 32.4's θ_max slider on real data, and Figs. 9–10 show the blind spot for the mean: robust in the phantom, up to about 100 m/s of spread below 60 mm in vivo.

Learned sound-speed inversion from raw channel data

M. Feigin, D. Freedman and B. W. Anthony, "A deep learning framework for single-sided sound speed inversion in medical ultrasound," IEEE Trans. Biomed. Eng. 67(4), 1142–1151 (2020); preprint arXiv:1810.00322 (figure numbers below follow the preprint). doi:10.1109/TBME.2019.2931195

Look at Fig. 6 (16 simulated test cases: reference versus three-plane-wave network output), Fig. 7 (absolute error for one and three plane waves), Fig. 8 (polyurethane phantom with an inclusion: B-mode, ground truth and reconstructions) and Fig. 10 (neck cross-section in vivo). Why: compare the clean simulated recoveries with the artefacts on the phantom: the domain gap of Section 32.6 in one page.

Differentiable-beamforming autofocus with a common-midpoint phase loss

W. Simson, L. Zhuang, B. N. Frey, S. J. Sanabria, J. J. Dahl and D. Hyun, "Ultrasound autofocusing: common midpoint phase error optimization via differentiable beamforming," arXiv:2410.03008 (version 3). Earlier work: W. Simson et al., "Differentiable beamforming for ultrasound autofocusing," MICCAI 2023, LNCS 14229, 428–437. doi:10.1007/978-3-031-43999-5_41

Look at Fig. 1 (straight-ray model of two common-midpoint subaperture pairs at lag 4 and their regions of support), Fig. 3 (simulated inclusion: lag-1 midpoint-offset phase error versus CMPE, at constant and correct speed), Fig. 5 (in vivo rat liver before and after autofocusing, with CMPE maps and coherence) and Fig. 7 (estimated versus calibrated phantom speeds). Why: Fig. 1 is (32.12) drawn by the authors, and Fig. 7 is the honest error budget: 8–14 m/s below calibrated values in layered and commercial phantoms.

32.15Go deeper

M. E. Anderson and G. E. Trahey, "The direct estimation of sound speed using pulse-echo ultrasound," J. Acoust. Soc. Am. 104(5), 3099–3106 (1998). doi:10.1121/1.423889. Average speed from the arrival-time profile.
C. H. Dix, "Seismic velocities from surface measurements," Geophysics (1955). The RMS-velocity moveout and the interval formula of (32.6).
M. Jakovljevic, S. Hsieh, R. Ali, G. Chau, D. Hyun and J. J. Dahl, "Local speed of sound estimation in tissue using pulse-echo ultrasound: model-based approach," J. Acoust. Soc. Am. 144, 254–266 (2018). Average-to-local in two-layer media.
M. Imbault, A. Faccinetto, B.-F. Osmanski, A. Tissier, T. Deffieux, J.-L. Gennisson, V. Vilgrain and M. Tanter, "Robust sound speed estimation for ultrasound-based hepatic steatosis assessment," Phys. Med. Biol. (2017). doi:10.1088/1361-6560/aa6226. Liver speed against MRI-PDFF and biopsy (open access).
M. Jaeger, K. G. Held, S. E. Peeters, S. Preisser, M. Grünig and M. Frenz, "Computed ultrasound tomography in echo mode for imaging speed of sound using pulse-echo sonography: proof of principle," Ultrasound Med. Biol. 41(1), 235–250 (2015). doi:10.1016/j.ultrasmedbio.2014.05.019. The original CUTE paper: 0.6 % contrast resolution, 1 mm lateral by 3 mm axial.
S. J. Sanabria, E. Ozkan, M. Rominger and O. Goksel, "Spatial domain reconstruction for imaging speed-of-sound with pulse-echo ultrasound: simulation and in vivo study," Phys. Med. Biol. 63(21) (2018). doi:10.1088/1361-6560/aae2fb.
R. Ali, T. Brevett, D. Hyun, L. L. Brickson and J. J. Dahl, "Distributed aberration correction techniques based on tomographic sound speed estimates," IEEE Trans. Ultrason. Ferroelectr. Freq. Control (2022); R. Ali et al., "Sound speed estimation for distributed aberration correction in laterally varying media," IEEE Trans. Comput. Imaging (2023). From sound-speed maps back to delays (Module 33).
R. Ali, "Sound speed estimation and phase-aberration correction for medical ultrasound imaging," Ph.D. thesis, Stanford University (2021). Includes the groundwork for wave-equation migration velocity analysis.
IEEE IUS 2026 QUASOS challenge, final scores: UW Quantitative Ultrasound Lab page.
Lab cross-links: coherence and the van Cittert–Zernike theorem (Module 07); delay estimation (Module 11); CRLB (Module 12); Tikhonov, resolution matrices and limited-angle SoS tomography (Module 13); autofocus and cycle skipping (Module 14); QUS (Module 27); CF and SLSC (Module 25); elastography and ULM, which inherit the speed (Modules 30, 31); aberration correction (Module 33).

← 31 ULM end to end 33 Aberration correction →