Report an errorErrata
Part V · Module 33

Aberration correction (capstone)

Aberration is what is left of the sound-speed problem after the beamformer has done its best with one speed: a delay error that changes across the aperture, and usually across the image too. Correcting it is delay estimation (Module 11) on signals whose coherence is set by the van Cittert–Zernike theorem (Module 07), turned into an eigenvector problem by the SVD (Module 06) and into a matrix change of basis by the reflection matrix (Module 05). This capstone builds the phase-screen model, the classical estimators (nearest-neighbour correlation, speckle brightness, time reversal), the isoplanatic patch that limits all of them, correction through a sound-speed map (Module 32), and the reflection and distortion matrices that make the correction local, and closes with transcranial, iterative and learned methods and the metrics that judge them.

You will be able to

  • Model an aberrator as a random delay screen and predict its Strehl ratio, PSF halo and coherence loss from two numbers, the rms delay and the correlation length
  • Derive why nearest-neighbour cross-correlation drifts as a random walk, when multi-lag least squares fixes it, and why a blurred transmit focus biases every speckle-based estimate until you iterate
  • Show that speckle-brightness maximization, time reversal and DORT all return the dominant eigenvector of an aperture covariance, and say what that eigenvector focuses on
  • Compute the isoplanatic patch of a screen at depth \(d\), including the distance beyond which a borrowed correction is worse than none
  • Turn a Module 32 sound-speed map into a space-variant correction and budget its errors (straight rays, blur, speed bias)
  • Build the focused reflection matrix and the distortion matrix of a speckle medium, extract a local aberration law by SVD, iterate it, and read the singular spectrum as a count of isoplanatic modes

Assumes

Exact formulaclosed form plotted as is, or published numbers re-plottedSimulationsolvers that run here in JS: Fraunhofer focal fields, reflection matrices of simulated speckle, power iteration, Monte CarloApprox. modelsimplified physics; assumptions under the badgeSchematicdrawing for intuition, not to scaleDefaults: \(c_0\) = 1540 m/s, \(f_0\) = 5 MHz (\(\lambda\) = 0.308 mm), pitch 0.3 mm, 128-element linear array (38.4 mm) unless a demo says otherwise.

33.1The phase-screen model: two numbers that predict the damage

Fat, muscle and connective tissue in the body wall have speeds from about 1450 to 1600 m/s and irregular boundaries. After the beamformer has removed the geometric delays for \(c_0\), each element still sees an excess delay \(\tau(x)\) that depends on the path. The simplest model that captures the damage puts all of it in a thin screen against the array: a random delay function of the element position, Gaussian with rms \(\sigma_\tau\) and a correlation length \(\ell\). In matrix language (Module 05) it is a diagonal matrix of phasors applied on transmit and again on receive, and its correction is the conjugate.

\[ \phi(x) = 2\pi f_0\,\tau(x), \qquad \langle \tau(x)\,\tau(x+\Delta x)\rangle = \sigma_\tau^{2}\,\rho(\Delta x), \qquad \rho(\Delta x) = e^{-\Delta x^{2}/\ell^{2}}, \qquad \sigma_\phi = 2\pi f_0\,\sigma_\tau \]
(33.1)
\(\tau(x)\)
residual (aberration) delay at aperture position \(x\) after geometric focusing with \(c_0\), s
\(\phi(x)\)
the same delay as a phase at the centre frequency, rad
\(\sigma_\tau,\ \sigma_\phi\)
rms delay (ns) and rms phase (rad)
\(\rho,\ \ell\)
correlation coefficient of the screen and its correlation length (Gaussian model), mm
\(f_0\)
centre frequency, MHz
\(\langle\cdot\rangle\)
ensemble average over screen realizations
\(\Delta x\)
separation between two aperture positions, mm
At 5 MHz, 10 ns of delay is 0.314 rad. A 40 ns rms screen is 1.26 rad: a fifth of the 200 ns period (72°), enough to remove two thirds of the focal intensity. The model is narrowband (one phase per element) and puts the aberrator at the array, which makes it exactly isoplanatic; Section 33.4 moves it away.

Strehl ratio. The on-axis focal intensity relative to the unaberrated focus is the squared mean phasor across the aperture. For a Gaussian screen its expectation is exact, and it has two simple limits:

\[ S = \Bigl\lvert \frac{1}{N}\sum_{i=1}^{N} e^{j\phi_i} \Bigr\rvert^{2}, \qquad \mathbb{E}[S] = \frac{1}{N^{2}}\sum_{i,j} e^{-\sigma_\phi^{2}\,[1-\rho(x_i - x_j)]} \;\approx\; \begin{cases} e^{-\sigma_\phi^{2}} & \sigma_\phi \lesssim 1,\ \ell \ll D \\[4pt] \dfrac{\sqrt{\pi}\,\ell}{\sigma_\phi\,D} & \sigma_\phi \gg 1 \end{cases} \]
(33.2)
\(S\)
Strehl ratio: on-axis intensity with the screen over intensity without, one way
\(N,\ x_i\)
number of elements and their positions
\(D\)
aperture length \(N p\), mm
\(\mathbb{E}[S]\)
ensemble mean over screen realizations
\(\phi_i\)
aberration phase of element \(i\), rad (\(\phi(x_i)\) of (33.1))
\(\sigma_\phi,\ \rho,\ \ell\)
as in (33.1)
\(j\)
imaginary unit
\(i,\ j\)
element indices (in the double sum, \(j\) is an index, not the imaginary unit)
The double sum uses \(\mathbb{E}[e^{j(\phi_i-\phi_j)}] = e^{-\frac12\mathrm{Var}(\phi_i - \phi_j)}\) for Gaussian phases. m33_phase_screen_psf.m (128 elements, \(\ell\) = 4 mm, 400 screens): at 40 ns the Monte Carlo Strehl is 0.311, the double sum 0.315, \(e^{-\sigma_\phi^2}\) only 0.206; at 80 ns 0.078, 0.079 and 0.002, where the large-\(\sigma\) form gives 0.073. Two way (same screen on transmit and receive) the Strehl of a realization is exactly \(S^{2}\); its mean \(\mathbb{E}[S^{2}]\) = 0.131 at 40 ns is above \(\mathbb{E}[S]^{2}\) = 0.099 because \(S\) fluctuates between realizations.

Why the Maréchal-type \(e^{-\sigma_\phi^{2}}\) underestimates: elements closer than \(\ell\) share their phase and keep adding coherently, so a large aperture behaves like \(D/\ell\) partly independent cells, and a fraction of order \(\ell/(\sigma_\phi D)\) of the intensity survives in the focus. The intensity that leaves the focus forms a halo. Its width follows from the slope of the screen: the local wavefront tilt has rms \(\sqrt2\,c\,\sigma_\tau/\ell\), so for a strong screen the mean intensity is a Gaussian of rms width \(s\). And on receive the same structure function multiplies the van Cittert–Zernike coherence of speckle (Module 07):

\[ s = \frac{\sqrt2\,c\,\sigma_\tau\,z}{\ell}\ \ (\sigma_\phi \gg 1), \qquad R(m) = R_0(m)\; e^{-\sigma_\phi^{2}\,[1-\rho(m p)]}, \qquad R_0(m) = 1 - \frac{m}{N} \]
(33.3)
\(s\)
rms width of the Gaussian halo of the mean focal intensity, mm (−6 dB width 3.33 s)
\(z\)
focal depth, mm
\(R(m),\ R_0(m)\)
receive spatial coherence at lag \(m\) elements with and without the screen (uniform transmit aperture, ideal transmit focus)
\(p\)
element pitch
\(c\)
sound speed, m/s
\(\sigma_\tau,\ \sigma_\phi,\ \rho,\ \ell\)
as in (33.1)
\(m\)
element lag (separation in elements)
\(N\)
number of elements
The −6 dB width of the mean PSF stays at the diffraction value (0.38 mm at F/1, 40 mm) up to about 30 ns, then jumps to 1.4 mm at 40 ns and 3.5–5.6 mm at 60–80 ns, approaching the halo prediction (4.4 and 5.8 mm). Coherence at 40 ns: lag 10 drops from 0.92 to 0.465 (formula 0.467), lag 40 from 0.69 to 0.138 (0.142). A broadened transmit spot lowers \(R_0\) further, which is why coherence is both a symptom and, later, a metric.
Demo 33.1

Phase screen generator: the screen, the PSF it leaves, and the coherence it removes

SimulationOne Gaussian screen realization (33.1) on 128 elements (white noise convolved with a Gaussian kernel, seeded, JS). Focal-plane (Fraunhofer) field of the focused aperture at z = 40 mm, monochromatic; two way = same screen on transmit and receive. The dashed mean-intensity curve and the coherence curves are the exact Gaussian-screen expressions (33.2), (33.3). Same model as m33_phase_screen_psf.m.
no screenthis realizationexact ensemble mean (one way)
Try: at 40 ns the realization's Strehl scatters around the ensemble value of 0.31, well above e−σ² = 0.21. Push στ to 100 ns: the core collapses into a halo whose width grows as στ/ℓ (33.3). Halve ℓ at fixed στ: the Strehl falls and the halo widens, because the slopes double. Lower the frequency to 2.5 MHz: the same delays cost a quarter of the phase variance. Two way squares the loss.
FamilyWhat it estimatesGuideValid overWhere in this module
Nearest-neighbour cross-correlation (Flax and O'Donnell, 1988)element-to-element delay differences, integratedpoint target or speckleone isoplanatic patch33.2, Demo 33.2
Speckle brightness maximization (Nock, Trahey and Smith, 1989)the law that maximizes the energy of a regionspeckleone patch, one beam33.3, Demo 33.3
Time reversal and DORT (Prada and Fink)eigenvectors of the time-reversal operatorpoint scatterers, or a speckle covarianceone patch per eigenvector33.3
Sound-speed map (Module 32)delays integrated through an estimated mapthe mapeverywhere the map is right33.5, Demo 33.5
Reflection and distortion matrices (Lambert, Aubry and coworkers)local laws from the SVD of re-centred focused datavirtual guide stars in specklemany patches, transmit and receive separately33.6, Demo 33.6

33.2Estimating the law from echoes: nearest-neighbour correlation, drift and least squares

Flax and O'Donnell (IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 1988) made the screen measurable without a guide star. Adjacent elements look at the same focal region through almost the same patch of screen, so their echoes are highly correlated, for a point reflector and, by van Cittert–Zernike, for speckle (Module 07: lag-one coherence \(1 - 1/N\) for an ideal focus). Cross-correlate each element with its neighbour (Module 11), read the delay difference, and integrate:

\[ \hat d_i = \tau_{i+1} - \tau_i + \varepsilon_i, \qquad \hat\tau_k = \sum_{i=1}^{k-1}\hat d_i, \qquad \mathrm{Var}(\hat\tau_k - \tau_k) = (k-1)\,\sigma_\varepsilon^{2}, \qquad \overline{e^{2}}\,\Big\rvert_{\text{mean, ramp removed}} \approx \frac{N}{15}\,\sigma_\varepsilon^{2} \]
(33.4)
\(\hat d_i\)
measured delay difference between elements \(i+1\) and \(i\), s
\(\varepsilon_i,\ \sigma_\varepsilon\)
its error and the error standard deviation (independent pairs)
\(\hat\tau_k\)
integrated estimate at element \(k\), defined up to a constant
\(\overline{e^{2}}\)
mean square error across the aperture after the mean and the linear ramp are removed
\(\tau_i,\ \tau_k\)
true aberration delay of element \(i\) or \(k\), s
\(i,\ k\)
element indices
\(\mathrm{Var}\)
variance
\(N\)
number of elements
The integration turns independent errors into a random walk: the raw error grows as \(\sqrt{k}\), and even after removing the best line (which only steers the beam, and which speckle cannot determine anyway) the rms error is \(\sigma_\varepsilon\sqrt{N/15}\) = 2.9 \(\sigma_\varepsilon\) for 128 elements (checked by Monte Carlo: \(N/15\) to 0.1 %). If instead the error is element noise, each noise sample enters two adjacent differences with opposite signs and the sum telescopes: no drift at all. Real data sit between the two.

Using more pairs breaks the random walk. Every pair \((i, j)\) up to a lag \(L\) gives a measured difference \(\hat d_{ij}\) of known reliability, and the profile is the weighted least-squares solution of an overdetermined system (Module 13) whose unknowns sit on the nodes of a graph and whose data sit on its edges, with one extra row to fix the constant:

\[ \hat{\mathbf{t}} = \arg\min_{\mathbf{t}} \sum_{1 \le j - i \le L} w_{ij}\,\bigl(\tau_j - \tau_i - \hat d_{ij}\bigr)^{2} + \bigl(\mathbf{1}^{T}\mathbf{t}\bigr)^{2} = \bigl(\mathbf{A}^{T}\mathbf{W}\mathbf{A} + \mathbf{1}\mathbf{1}^{T}\bigr)^{-1}\mathbf{A}^{T}\mathbf{W}\,\hat{\mathbf{d}}, \qquad w_{ij} = \frac{\lvert\rho_{ij}\rvert^{2}}{1 - \lvert\rho_{ij}\rvert^{2}} \]
(33.5)
\(\mathbf{t}\)
vector of element delays \(\tau_i\), s
\(\mathbf{A}\)
difference operator: one row per pair, −1 at \(i\), +1 at \(j\)
\(\mathbf{W}\)
diagonal weights; \(w_{ij}\) is proportional to the inverse variance of a correlation phase (Module 11, Section 32.4)
\(\rho_{ij}\)
measured correlation coefficient of the pair
\(L\)
maximum lag used
\(\hat{\mathbf{t}},\ \hat{\mathbf{d}}\)
estimated delay vector and vector of measured pair differences \(\hat d_{ij}\), s
\(\tau_i,\ \tau_j\)
delays of elements \(i\) and \(j\), s
\(w_{ij}\)
weight of pair \((i,j)\), dimensionless
\(\mathbf{1}\)
vector of ones (the term pins the mean delay to zero)
\((\cdot)^{T},\ (\cdot)^{-1}\)
transpose, matrix inverse
\(\arg\min\)
minimizing argument
\(\mathbf{A}^{T}\mathbf{W}\mathbf{A}\) is a weighted graph Laplacian whose only null vector is the constant; the \(\mathbf{1}\mathbf{1}^{T}\) term removes it. With \(L = 1\) the solution is the integrated nearest-neighbour estimate; each added lag adds a path between distant elements that does not accumulate the intermediate errors.

m33_nn_crosscorr_correction.m synthesizes speckle channel data (128 elements, 60 % bandwidth, dynamic receive focus) through a 40 ns, 4 mm screen applied on transmit and receive, and finds two regimes. Noiseless speckle: the adjacent-element correlation is 0.992, the pair jitter only 0.68 ns, yet the residual after correction is 5.7 ns rms, three times the random-walk prediction, and least squares does not reduce it. That residual is a bias: the aberrated transmit beam insonifies a broad, irregular spot, and the echoes estimate the law that focuses on that spot, not on the intended point. Transmitting again with the pass-1 estimate narrows the spot and halves the residual (2.6 ns, Strehl 0.992). 0 dB channel SNR: correlation 0.50, jitter 4.5 ns, and now the drift dominates: after the transmit iteration the nearest-neighbour residual is 8.3 ns (Strehl 0.93) and least squares over 8 lags gives 4.0 ns (Strehl 0.98). The random-walk bound \(\sigma_\varepsilon\sqrt{N/15}\) = 13.1 ns overestimates because part of the error is element noise, which telescopes.

Pitfall

Judging an aberration estimator by its pair jitter. Sub-nanosecond jitter is easy in bright speckle; the error that matters is the bias from the blurred transmit focus, which no amount of averaging or least squares removes. Iterate the transmit (Måsøy, Varslot and Angelsen reported near-ideal correction after two to three transmit iterations, J. Acoust. Soc. Am., 2005), or estimate where the focus is already good.

Demo 33.2

Nearest-neighbour integration drifts; multi-lag least squares does not

Approx. modelError model only, no channel data: each pair estimate \(\hat d_{ij}\) = true difference + Gaussian error with the interferometric standard deviation \(\sqrt{1-\rho_m^{2}}/(\rho_m\sqrt{2N_k})/(2\pi f_0)\), \(N_k\) = 32 independent samples per kernel, \(\rho_m = \rho_1^{\,m}\) (exponential decorrelation with lag, assumed), independent across pairs; or, with the second error type, element phase noise that enters both members of a pair. 128 elements, 5 MHz, 150 Monte Carlo trials, LS by Cholesky (33.5). Bias from the transmit focus is not modelled (see m33_nn_crosscorr_correction.m).
nearest neighbour, integratedleast squares, lags 1…Ltrue screen / σε√(N/15)
Try: with independent pair errors the nearest-neighbour error (L = 1) is largest at the two ends of the aperture (about 4 σε) with a hump in the middle, and averages σε√(N/15) over the elements. Raise L: least squares flattens it, quickly at first and then slowly, because the long lags are noisier (ρm = ρ1m). Lower ρ1 to 0.6: long lags become useless and the optimum L shrinks. Switch to element noise: the integrated estimate no longer drifts, and LS gains little; the hump was a property of the error correlations, not of integration as such.

33.3Speckle brightness, time reversal and DORT: one beam, one eigenvector

Nock, Trahey and Smith (J. Acoust. Soc. Am., 1989) replaced the correlation by a figure of merit: adjust the element delays to maximize the brightness of a speckle region, which is largest when the beam is focused. Zhao and Trahey (1991) compared such quality factors for diffuse and point targets. Brightness has a structure worth seeing. Let \(\mathbf{s}\) be the vector of focused channel signals of a pixel and \(\mathbf{w}\) the receive weights; the beamformed output is \(\mathbf{w}^{H}\mathbf{s}\) and its mean energy is a quadratic form in the aperture covariance:

\[ B(\mathbf{w}) = \mathbb{E}\bigl\lvert \mathbf{w}^{H}\mathbf{s}\bigr\rvert^{2} = \mathbf{w}^{H}\mathbf{C}\,\mathbf{w}, \qquad \mathbf{C} = \mathbb{E}\bigl[\mathbf{s}\mathbf{s}^{H}\bigr] = \mathbf{\Phi}\,\mathbf{C}_0\,\mathbf{\Phi}^{H}, \qquad \mathbf{\Phi} = \mathrm{diag}\bigl(e^{j\phi_i}\bigr), \qquad \max_{\lVert\mathbf{w}\rVert^{2} = N} B \;\Rightarrow\; \mathbf{w} \propto \mathbf{u}_1(\mathbf{C}) = \mathbf{\Phi}\,\mathbf{u}_1(\mathbf{C}_0) \]
(33.6)
\(B\)
mean speckle brightness of the beamformed pixel
\(\mathbf{s}\)
focused (geometrically delayed) channel signals, one per element
\(\mathbf{C},\ \mathbf{C}_0\)
aperture covariance with and without the screen; \(\mathbf{C}_0\) is the van Cittert–Zernike matrix of the transmit focal intensity
\(\mathbf{\Phi}\)
diagonal screen matrix at the array (Module 05)
\(\mathbf{u}_1\)
eigenvector of the largest eigenvalue
\(\mathbf{w}\)
receive weight vector, \(N\) elements
\(\mathbb{E}\)
expectation over speckle
\((\cdot)^{H}\)
conjugate transpose
\(\mathrm{diag}\)
diagonal matrix of the listed entries
\(\phi_i,\ j\)
screen phase of element \(i\) (rad) and imaginary unit
\(\propto\)
proportional to
\(\max,\ \lVert\cdot\rVert\)
maximize, subject to fixed weight norm
The unconstrained maximizer is the dominant eigenvector (Varslot, Mo, Angelsen and Krogstad, J. Acoust. Soc. Am., 2004, showed that the brightness-maximizing filter focuses on the regions of high insonification of the aberrated transmit beam). Its phase is the screen phase plus the phase of \(\mathbf{u}_1(\mathbf{C}_0)\), which is flat (up to a ramp from the spot's centroid) only when the insonified spot is compact and symmetric. A phase-only search (\(\lvert w_i\rvert = 1\)) is solved by coordinate ascent, \(w_i \leftarrow e^{j\arg(\mathbf{C}\mathbf{w})_i}\), which never decreases \(B\).

Mallart and Fink (J. Acoust. Soc. Am., 1991; the aberration and focusing-criterion paper followed in 1994) gave the van Cittert–Zernike reading: the speckle behind a focused transmit radiates like an incoherent source shaped like the transmit focal intensity, so \(\mathbf{C}\) is the covariance a real source of that shape would produce through the screen. A point-like source gives a rank-one \(\mathbf{C}\) and a perfect estimate; an aberrated, broad source spreads the eigenvalues, and an asymmetric one biases the eigenvector beyond a removable ramp. That is the bias of Section 33.2 again, and the same remedy applies: transmit with the estimate, re-measure, repeat.

Time reversal and DORT. With discrete scatterers the inter-element response matrix is \(\mathbf{K} = \mathbf{H}\,\mathrm{diag}(a)\,\mathbf{H}^{T}\) (monochromatic, single scattering, \(\mathbf{H}\) the element-to-scatterer Green's matrix including the aberrator). One time-reversal round trip (emit \(\mathbf{e}\), receive \(\mathbf{K}\mathbf{e}\), phase-conjugate, re-emit) returns \(\mathbf{K}\mathbf{K}^{*}\mathbf{e}^{*}\); two round trips apply \(\mathbf{K}^{*}\mathbf{K}\) (the time-reversal operator, equal to \(\mathbf{K}^{H}\mathbf{K}\) since \(\mathbf{K}\) is symmetric), and iterating is power iteration, which converges to the brightest scatterer. Prada and Fink (Wave Motion, 1994) decomposed the time-reversal operator instead:

\[ \mathbf{K}^{H}\mathbf{K}\,\mathbf{v}_s = \lambda_s\,\mathbf{v}_s, \qquad \mathbf{v}_s \propto \mathbf{h}_s^{*}, \qquad \lambda_s = \lvert a_s\rvert^{2}\,\lVert\mathbf{h}_s\rVert^{4} \quad \text{(well-resolved scatterers: } \mathbf{h}_s^{H}\mathbf{h}_{s'} \approx 0) \]
(33.7)
\(\mathbf{K}\)
inter-element response matrix at one frequency (transmit on column, receive on row)
\(\mathbf{h}_s\)
Green's vector from scatterer \(s\) to the elements, aberrator included
\(a_s\)
reflectivity of scatterer \(s\)
\(\mathbf{v}_s,\ \lambda_s\)
eigenvectors and eigenvalues of the time-reversal operator
\(\lambda_s\)
eigenvalue belonging to scatterer \(s\) (listed with \(\mathbf{v}_s\))
\((\cdot)^{H},\ (\cdot)^{*}\)
conjugate transpose, complex conjugate
\(\lVert\cdot\rVert,\ \lvert\cdot\rvert\)
vector norm, magnitude
\(\propto\)
proportional to
One significant eigenvalue per resolved scatterer (Module 06), and each eigenvector is the conjugate of that scatterer's Green's vector, aberration included: emitting it focuses through the screen. In speckle there are no resolved scatterers, and \(\mathbf{K}^{H}\mathbf{K}\) has a slowly decaying spectrum; the distortion matrix of Section 33.6 is the trick that makes the speckle look like one virtual point.
Demo 33.3

Speckle brightness and the dominant eigenvector: the estimate follows the transmit spot

Simulation64 elements (19.2 mm), focus 30 mm, 5 MHz, monochromatic Fraunhofer model; Gaussian screen at the array on transmit and receive. Expected speckle covariance \(\mathbf{C} = \sum_s I_{\mathrm{tx}}(x_s)\,\mathbf{h}_s\mathbf{h}_s^{H}\) over a dense grid of focal-plane positions (infinite speckle averaging), \(I_{\mathrm{tx}}\) the transmit intensity with the current correction. Dominant eigenvector by power iteration; phase-only brightness by coordinate ascent. Mean and ramp removed before comparing laws.
ideal transmit / true lawtransmit used in this passeigenvector estimatephase-only brightness
Try: set στ = 60 ns (ℓ = 3 mm, first screen). Pass 0 transmits through the screen: the spot is broad and asymmetric, λ2/λ1 = 0.74, and the eigenvector misses the law by 0.28 rad rms; the brightness is only 0.39 of the ideal because the transmit half of the focus is still aberrated. Press "Transmit with estimate": 0.017 rad and B/Bideal = 0.93; once more: 0.001 rad, 1.000, and λ2/λ1 settles at 0.30, the value for an unaberrated sinc² spot. At 40 ns the first pass is already within 0.04 rad: in this idealized model (infinite speckle averaging, one depth, one frequency) only asymmetric spots bias the eigenvector, while the broadband, finite-speckle simulation of Section 33.2 kept a 0.18 rad (5.7 ns) first-pass bias. The phase-only brightness search tracks the eigenvector throughout.

33.4Distributed aberration and the isoplanatic patch

A screen against the array is a diagonal matrix: one law corrects every point of the image. Real aberrators are thick and sit millimetres to centimetres from the probe, so rays to different image points cross different parts of them. Keep a thin screen but move it to depth \(d\). The straight ray from element \(x_u\) to a point \((x_p, z)\) crosses the screen at \(x_u(1 - d/z) + x_p\,d/z\), so moving the point laterally by \(\Delta x\) slides the footprint of the whole aperture across the screen by \(\delta = \Delta x\,d/z\). A law measured on one point and borrowed for another leaves the residual \(\phi(a_u + \delta) - \phi(a_u)\):

\[ \delta = \Delta x\,\frac{d}{z}, \qquad S(\Delta x) \approx e^{-2\sigma_\phi^{2}\,[1 - \rho(\delta)]}, \qquad \Delta x_{\times} = \frac{\ell\,z\sqrt{\ln 2}}{d}, \qquad \Delta x_{1/2} = \frac{\ell\,z}{d}\sqrt{-\ln\!\left(1 - \frac{\ln 2}{2\sigma_\phi^{2}}\right)} \]
(33.8)
\(d,\ z\)
depth of the screen and of the image points, mm
\(\Delta x,\ \delta\)
lateral distance from the guide point, and the resulting slide of the aperture footprint on the screen, mm
\(S(\Delta x)\)
Strehl ratio at \(\Delta x\) with the law of the guide point (large aperture, \(\ell \ll D\))
\(\Delta x_\times\)
distance at which the borrowed correction is no better than no correction (\(\rho = 1/2\))
\(\Delta x_{1/2}\)
isoplanatic half-width at \(S = 1/2\) (defined for \(\sigma_\phi^{2} > \ln 2/2\))
\(\sigma_\phi,\ \rho,\ \ell\)
as in (33.1)
\(\ln\)
natural logarithm
\(e\)
Euler’s number
The residual of a borrowed law has variance \(2\sigma_\phi^{2}(1-\rho)\), which exceeds the uncorrected \(\sigma_\phi^{2}\) once \(\rho < 1/2\): outside the patch, a correction is worse than none. m33_isoplanatic_patch.m (40 ns, \(\ell\) = 4 mm, z = 40 mm, 128 elements, finite-aperture exact formula and Monte Carlo): half-widths at \(S = 1/2\) of 16.5, 8.3 and 4.3 mm for d = 5, 10 and 20 mm; the large-aperture form gives 8.0 mm at 10 mm. The borrowed law falls to the uncorrected level at 22.5, 11.0 and 5.1 mm, against \(\Delta x_\times\) = 26.6, 13.3 and 6.7 mm (the finite aperture keeps the uncorrected Strehl above \(e^{-\sigma_\phi^{2}}\)). The exponential in (33.8) is not a small-phase expansion: like \(e^{-\sigma_\phi^{2}}\) in (33.2), it is the exact Gaussian average in the limit \(\ell \ll D\). At finite aperture, elements closer than \(\ell\) add the floor of Section 33.1, which lifts the uncorrected Strehl (0.349 at d = 10 mm, against 0.206) more than the borrowed one (0.349 against 0.304 at 11 mm) and so moves the crossing inward.

Three consequences shape every modern method. First, the patch scales as \(z\ell/d\): a shallow, coarse aberrator is benign, a deep or fine one is not, which is why the skull (thick, so its far side sits at a sizeable fraction of the target depth, finely structured, and with a phase error of several radians, which shrinks \(\Delta x_{1/2}\) further as about \(1/\sigma_\phi\)) is the hardest case. Second, in depth the footprint also changes, so patches are finite axially too. Third, the patch is a property of the screen statistics, not of the estimator: no single-law method, however accurate at its guide point, beats \(\Delta x_\times\). Methods that work over the whole field must either estimate many local laws (Section 33.6) or model the medium between the probe and the image (Section 33.5).

Schematic A screen at depth d: rays to two image points cross it at footprints shifted by δ = Δx·d/z. Not to scale.
array, aperture D screen at depth dφ(x): rms σ, corr. length ℓ guide pointx + Δx footprints shift by δ = Δx d/z d z
Demo 33.4

How far does one law reach? The isoplanatic patch versus screen depth

Exact formulaExact mean Strehl of a Gaussian screen at depth d for a 128-element aperture (the finite-aperture double sum written in m33_isoplanatic_patch.m, evaluated by lags), with its large-aperture limit (33.8) dashed. Straight rays, thin screen, monochromatic, image points at one depth. The Monte Carlo check is in the script.
borrowed law, exactlarge-aperture limitno correctionΔx×, Δx1/2
Try: the defaults reproduce the script's 10 mm row: S = 0.95 at 2 mm, 0.52 at 8 mm, the uncorrected level at 11.0 mm. Slide d to 1 mm: the screen is almost at the array and one law covers the whole image (S = 0.99 at 8 mm). Slide d to 30 mm: the half-width shrinks to 3.4 mm and the borrowed law loses to no correction beyond 2.9 mm; the uncorrected Strehl rises to 0.58 because the aperture's footprint on a deep screen is smaller and spans fewer correlation lengths. Go deeper (z = 120 mm, d = 10 mm): the half-width grows to 25 mm, as z/d predicts; shallow targets under a near-field aberrator are the patch-limited ones.

33.5Correction from a sound-speed map: the Module 32 route

Module 32 estimated the slowness field \(s(\mathbf r) = 1/c(\mathbf r)\). Once a map exists, aberration correction is a forward computation: integrate the slowness along the path from each element to each pixel and beamform with those times. The correction is space-variant by construction, so the isoplanatic patch disappears as a concept; what remains is the accuracy of the map and of the path model:

\[ T_u(\mathbf r) = \int_{\Gamma^{\star}_u(\mathbf r)} s\,dl, \qquad \hat T_u(\mathbf r) = \int_{\Gamma_u(\mathbf r)} \hat s\,dl, \qquad T_u - \hat T_u = \underbrace{\int_{\Gamma_u} (s - \hat s)\,dl}_{\text{map error}} \; + \; \underbrace{\Bigl(\int_{\Gamma^{\star}_u} s\,dl - \int_{\Gamma_u} s\,dl\Bigr)}_{\text{path error}\ \le\ 0} \]
(33.9)
\(T_u(\mathbf r)\)
true one-way travel time from element \(u\) to pixel \(\mathbf r\), s
\(\Gamma^{\star}_u\)
the true (Fermat, refracted) ray; \(\Gamma_u\) the path used by the beamformer, usually straight
\(\hat s\)
estimated slowness map, s/m
\(\mathbf r\)
pixel position, mm
\(u\)
element index
\(s\)
true slowness map, s/m
\(dl\)
path length element, m
\(\hat T_u(\mathbf r)\)
travel time predicted from \(\hat s\) along \(\Gamma_u\), s
The path error is never positive: Fermat's ray is the fastest, so straight rays always overestimate the time. Its size grows with the speed contrast and the interface slope. A Module 32 map has limited resolution (its resolution matrix, Section 32.5), so the map term is a low-pass version of the truth minus the truth, plus any bias in the mean (CUTE's blind spot).

m33_sos_map_correction.m puts a wavy fat layer (1450 m/s, interface at 12 ± 1.5 mm, period 15 mm) over 1540 m/s tissue and compares, at 25 points 40 mm deep across ±12 mm: no correction (45 ns rms, mean Strehl 0.08), one law measured at the centre (30 ns, 0.43; it keeps S ≥ 0.8 only within ±3 mm), straight rays through the true map (5.2 ns, 0.97), straight rays through a map blurred to 3 mm resolution (8.2 ns, 0.93) and through a map whose fat speed is 10 m/s too high (7.3 ns, 0.95). The straight-ray residual (3.4 ns at the centre, 11 ns at −12 mm) is the path error alone. Ali et al. (IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 2022) developed this tomographic route to distributed aberration correction on measured data.

Ultrasound link

Map-based correction and matrix correction answer different questions. The map predicts delays everywhere but inherits its own resolution, bias and path model; the matrix measures the delays that actually focus, but only where there is enough signal and only per patch. Using the map as the starting law and the matrix to measure the residual combines their strengths, and the residual itself is a datum for refining the map (Module 32's tomography with a better forward model).

Demo 33.5

One measured law versus a sound-speed map: straight rays, blur and bias

Approx. modelFat layer from the skin to h(x) = h₀ + A sin(2πx/P) over 1540 m/s. True times by Fermat with one refracting interface (coarse search plus golden-section refinement per element and point, JS). Map corrections by straight rays (analytic fat length); the blurred map is the true map convolved with a 2D Gaussian of the set FWHM (vertical part in closed form, lateral part by 7-point Gauss–Hermite quadrature), as in m33_sos_map_correction.m. 128 elements, 5 MHz, points at 40 mm; residuals with mean and ramp removed.
noneone law (x = 0)straight-ray mapblurred mapbiased fat speed
Try: with the defaults the one-law curve is a narrow bump around x = 0 while all three map curves stay at or above 0.85 across ±12 mm. Raise the amplitude to 2 mm: the straight-ray Strehl at −12 mm falls to 0.72 (19.5 ns), the path error growing with the interface slope, while at the centre it stays at 0.97. Shorten the period to 12 mm and widen the blur to 6 mm: the blurred map misses the wiggles (28 ns mean) while the sharp straight-ray map still gives 8.7 ns. Set the bias to −30 m/s: a smooth, nearly quadratic error (18 ns mean), a defocus that a global sound-speed tweak (Module 32) would absorb. The sliders stop at slopes where a straight ray could cross the interface twice, which this model does not handle. The geometry sliders recompute when released.

33.6The reflection matrix and the distortion matrix: local laws from speckle

Everything so far used one beam, or one guide point, at a time. Matrix imaging records the full set of responses between every transmit and every receive channel, the reflection matrix, and then changes basis in software (Module 05). The acquisition can be synthetic aperture (one element at a time, \(\mathbf{K}\) directly) or plane waves (Module 23), which are related to the element basis by a unitary-like transform at each frequency. Lambert, Cobus, Couade, Fink and Aubry (Phys. Rev. X, 2020) projected it onto a focused basis at each depth, focusing in post-processing on transmit at \(x_{\mathrm{in}}\) and on receive at \(x_{\mathrm{out}}\):

\[ \mathbf{R}_{xx}(z) = \mathbf{G}_0^{H}\,\mathbf{K}\,\mathbf{G}_0^{*}, \qquad G_0(u, x) = e^{-jk\sqrt{(x_u - x)^{2} + z^{2}}}, \qquad I(\Delta x) = \bigl\langle \lvert R_{xx}(x_{\mathrm{in}} + \Delta x,\, x_{\mathrm{in}})\rvert^{2} \bigr\rangle_{x_{\mathrm{in}}} \;\propto\; \int \lvert P_{\mathrm{out}}(y - \Delta x)\rvert^{2}\,\lvert P_{\mathrm{in}}(y)\rvert^{2}\,dy \]
(33.10)
\(\mathbf{K}\)
reflection matrix in the element basis at one frequency (or one time gate), \(N_e \times N_e\)
\(\mathbf{G}_0\)
free-space propagator from the elements to the focal points at depth \(z\) for the assumed \(c_0\), \(N_e \times N_x\)
\(\mathbf{R}_{xx}\)
focused reflection matrix: responses between virtual transducers at \(x_{\mathrm{in}}\) and \(x_{\mathrm{out}}\)
\(I(\Delta x)\)
common-midpoint intensity profile, averaged over a region
\(P_{\mathrm{in}},\ P_{\mathrm{out}}\)
transmit and receive point-spread functions
\(j,\ k\)
imaginary unit and wavenumber \(\omega/c_0\), rad/m
\(u,\ x,\ x_u\)
element index, focal-point position and element position, m
\(x_{\mathrm{in}},\ x_{\mathrm{out}}\)
transmit-side and receive-side virtual transducer positions, m
\(G_0(u, x)\)
entry of \(\mathbf{G}_0\)
\((\cdot)^{H},\ (\cdot)^{*}\)
conjugate transpose, complex conjugate
\(\langle\cdot\rangle\)
average over \(x_{\mathrm{in}}\)
\(y\)
integration variable (lateral position), m
The diagonal of \(\mathbf{R}_{xx}\) is the confocal image. In speckle (single scattering), the off-diagonal intensity is the cross-correlation of the transmit and receive PSF intensities, so the width of \(I(\Delta x)\) measures the local focusing quality without any guide star. Lambert et al. defined the focusing criterion \(F = w_0/w\), the ratio of the diffraction-limited width to the measured width of this profile, and used it to choose the sound speed that maximizes focusing, layer by layer (their Fig. 2).

The distortion matrix (Lambert, Cobus, Frappart, Fink and Aubry, PNAS, 2020; generalized to many patches in their two-part "Ultrasound matrix imaging" papers, IEEE Trans. Med. Imaging, 2022) keeps the focused transmit and expresses the receive side in an aperture basis: the plane-wave (far-field) basis in the PNAS paper, the transducer basis in Part II and here, either one for a screen near the probe. Each column is the wavefront received from one virtual source at \(x_{\mathrm{in}}\). Removing the geometric wavefront the source would produce in a homogeneous medium (multiplying by its conjugate) leaves only the distortion:

\[ \mathbf{R}_{ux} = \mathbf{K}\,\mathbf{G}_0^{*}, \qquad \mathbf{D} = \mathbf{R}_{ux} \circ \mathbf{G}_0^{*}, \qquad \mathbf{C} = \mathbf{D}\mathbf{D}^{H} = \mathbf{U}\,\mathbf{\Sigma}^{2}\,\mathbf{U}^{H} \;\approx\; \mathbf{\Phi}\,\mathbf{C}_{\mathrm{vr}}\,\mathbf{\Phi}^{H}, \qquad \hat\phi_u = \arg U_{u1} \]
(33.11)
\(\mathbf{R}_{ux}\)
reflection matrix with focused transmit (columns \(x_{\mathrm{in}}\)) and element receive (rows \(u\))
\(\circ\)
Hadamard (element-wise) product
\(\mathbf{D}\)
distortion matrix: each column is a received wavefront re-centred on its own focal point
\(\mathbf{C},\ \mathbf{U},\ \mathbf{\Sigma}\)
its correlation matrix, left singular vectors and singular values
\(\mathbf{C}_{\mathrm{vr}}\)
van Cittert–Zernike covariance of a virtual reflector shaped like the transmit focal intensity, common to all columns within a patch
\(\hat\phi_u\)
estimated aberration law (phase of the first singular vector)
\(\mathbf{K},\ \mathbf{G}_0\)
as in (33.10)
\(\mathbf{\Phi}\)
diagonal screen matrix, as in (33.6)
\(U_{u1}\)
entry \(u\) of the first left singular vector \(\mathbf{U}_{\cdot 1}\)
\(\arg\)
phase angle, rad
\((\cdot)^{H}\)
conjugate transpose
Re-centring makes every focal point in an isoplanatic patch look like the same source at the origin, seen through the same screen; averaging the columns over the patch averages over independent speckle realizations (one per resolution cell). This is (33.6) with the speckle average done by the matrix. \(\hat\phi\) is the receive-side law; by reciprocity (the same screen on both paths) the script and Demo 33.6 apply it to both sides of \(\mathbf{K}\) and recompute: the corrected transmit shrinks the virtual reflector toward a point and the estimate converges. Lambert et al. instead alternate the correction between the transmit and receive sides and shrink the window at each step.

m33_distortion_matrix_svd.m runs this on a simulated speckle layer at 20 mm behind a 40 ns, 2 mm screen (48 elements, 61 focal points over 12 mm, monochromatic). Screen at the array: the residual falls from 0.83 rad (Strehl 0.50) to 0.31, 0.07 and 0.018 rad in three iterations, and the focusing quality \(\Gamma\) (defined below) from 0.48 to 0.87, 0.99 and 1.00. Screen at 8 mm (anisoplanatic over the field): one global law raises \(\Gamma\) only from 0.48 to 0.69 and dims a bright point at the centre by 7.6 dB, while three 4 mm windows reach 0.75, 0.93 and 0.88 (mean 0.86) and brighten the point by 3.1 dB. Two cautions from the same script: the singular spectrum of \(\mathbf{D}\) does not collapse to one value even after perfect correction (\(\sigma_2/\sigma_1\) = 0.52), because the virtual reflector keeps the size of the diffraction-limited spot; and an entrywise-normalized \(\mathbf{C}/\lvert\mathbf{C}\rvert\), which makes the reflector point-like, is no longer positive semidefinite, so its eigenvalues need care when counted.

\[ \Gamma = \frac{I(0)\big/\sum_{\Delta x} I(\Delta x)}{\bigl[I(0)\big/\sum_{\Delta x} I(\Delta x)\bigr]_{\text{no aberration}}}, \qquad F = \frac{w_0}{w}, \qquad \frac{\lvert \mathbf{u}_1^{H}\mathbf{v}^{(n)}\rvert}{\lVert\mathbf{v}^{(n)}\rVert} \to 1\ \text{ with error } \propto \Bigl(\frac{\lambda_2}{\lambda_1}\Bigr)^{n} \]
(33.12)
\(\Gamma\)
confocal fraction of the common-midpoint energy, normalized by its aberration-free value; a Strehl-like focusing measure used in this module's demos and scripts
\(F\)
focusing criterion of Lambert et al.: diffraction-limited over measured width of \(I(\Delta x)\)
\(\mathbf{v}^{(n)}\)
power-iteration vector after \(n\) multiplications by \(\mathbf{C}\) (Module 06)
\(\lambda_1,\ \lambda_2\)
the two largest eigenvalues of \(\mathbf{C}\)
\(I(\Delta x)\)
common-midpoint intensity profile, as in (33.10)
\(w_0,\ w\)
diffraction-limited and measured width of \(I(\Delta x)\), mm
\(\mathbf{u}_1\)
dominant eigenvector of \(\mathbf{C}\)
\(\propto\)
proportional to
\(n\)
iteration count
Width-based measures saturate for moderate aberration because the core of the PSF keeps its width while energy moves into a halo (Section 33.1); an energy fraction does not. Neither sees a pure shift of both foci: a residual tilt warps the image without blurring it (Module 31's bias field).
Schematic One recording, three bases. Element basis (what is measured), focused basis (what is imaged and how well), distortion basis (what is corrected).
K: element × element R_xx: x_out × x_in D: element × x_in measured; no visible structure diagonal = image; band width = focus rows share exp(jφ_u):rank ≈ 1 per patch G0^H · K · G0* (K G0*) ∘ G0*
Demo 33.6

Reflection matrix, distortion matrix and SVD, computed live on simulated speckle

Simulation48 elements (pitch 0.3 mm), 5 MHz, monochromatic. 400 complex Gaussian scatterers plus one bright point (amplitude 6) in a thin layer at z = 20 mm (stands in for a time gate). Gaussian screen at depth d; the straight ray from element to scatterer crosses it at \(x_u + (x_s - x_u)d/z\). \(\mathbf{K} = \mathbf{G}\,\mathrm{diag}(a)\,\mathbf{G}^{T}\); 61 focal points over ±6 mm; \(\mathbf{R}_{xx}\), \(\mathbf{D}\), \(\mathbf{C}\) as in (33.10)–(33.11); first singular vector by power iteration, the next nine by deflation; correction applied to both sides of \(\mathbf{K}\) in post-processing. Same model as m33_distortion_matrix_svd.m (different random numbers). No attenuation, no multiple scattering, no noise.
before / true lawcurrent estimate / afterwindow
Try: with the screen at the array (d = 0), press Iterate three times: the band around the diagonal of Rxx collapses to a line, Γ goes from 0.40 to 1.00, the residual from 0.80 to 0.03 rad, and the bright point at x = 0 rises by 4.3 dB. Then set d = 8 mm: four iterations on the whole field stall at Γ = 0.63, while four in the centre window reach 0.91 there (the whole field only 0.60) and lift the bright point by 4 dB. Compare the laws estimated in the left, centre and right windows: those are the isoplanatic modes. At 80 ns one iteration already removes most of the error (0.39 rad, Γ 0.84) and four reach 0.02 rad. Watch σ2/σ1: it falls but stops near 0.5, the size of the diffraction-limited virtual reflector, not zero. The screen sliders rebuild the medium when released.

33.7Frontiers and judgement: iterative matrix methods, the skull, learned correction and the metrics

Iterative, alternating, shrinking. The full matrix pipeline of Lambert et al. (IEEE Trans. Med. Imaging, 2022, Part II) alternates the correction between the transmit and receive sides, works in either the plane-wave or the transducer basis, and shrinks the local window at each step (from 12 × 16 mm to 6.2 × 8.1 mm in their in vivo example, each window holding of the order of 5000 independent virtual sources). On a gallbladder in a volunteer they report the transverse −3 dB resolution going from 1.69 to 0.50 mm and their focusing factor rising from 0.31 to 0.48 over four steps. The structure echoes CLASS from optics (closed-loop accumulation of single scattering; Kang et al., Nat. Commun., 2017): in a reflection matrix expressed in the angular (pupil) basis, find the input and the output phase corrections that maximize the total confocal (single-scattering) intensity, updating one side while holding the other, and iterate. "CLASS-like" in ultrasound usually means exactly that alternating, energy-maximizing update; CLASS itself searches the phase channel by channel and uses no SVD. The distortion-matrix SVD can stand in for each half-step: with a receive law \(\mathbf{w}\) the confocal energy of a patch is \(\mathbf{w}^{H}\mathbf{C}\,\mathbf{w}\), CLASS maximizes it over unit-modulus laws, and \(\mathbf{u}_1\) maximizes it when the amplitudes are left free (this is the module's reading, not a statement from either paper).

Through the skull. Bone is roughly twice as fast as soft tissue, strongly attenuating and variable in thickness, and it sits against the probe. The thin-screen picture then fails three ways: amplitude varies as much as phase, the aberrator is thick (the patch is finite in depth too), and for shallow targets \(d/z\) is large, so patches are small (33.8). Three strategies coexist. Priors: CT-derived speed and density maps give the delays, as in transskull therapy focusing (Aubry and coworkers, J. Acoust. Soc. Am., 2003); this is the Section 33.5 route with a better map. Guide stars: in transcranial ULM, Demené et al. (Nat. Biomed. Eng., 2021) used the microbubbles themselves as guide stars to estimate the aberration law before localizing, patch by patch, which Module 31 showed is necessary because an uncorrected tilt displaces every bubble. Matrices: 2D-array matrix imaging (Bureau et al., Nat. Commun., 2023) estimates the transmission matrix voxel by voxel in 3D, with transcranial imaging as the target.

Learned correction. Networks trained on simulated channel data can output aberration laws, corrected delays or corrected images directly. They inherit the Module 32 caveat: the forward model moves into the training distribution, so skull, abdominal-wall and probe-lens mismatches become domain gaps. Physics-anchored designs (unrolled matrix or autofocus iterations, differentiable beamformers, Module 14) keep the measured data in the loop and are the safer path to deployment.

MetricWhat it measuresBlind spotLink
Coherence factor, lag-one coherence, SLSCphase alignment across the aperture; for unit amplitudes CF of a point is the receive Strehllowered by noise and clutter as well as aberration; blind to a pure shiftModules 07, 25
Strehl, Γ (energy fraction)energy in the focus relative to idealneeds a reference; blind to a common shift of both foci33.1, 33.6
FWHM of a point, \(F = w_0/w\)core widthsaturates: the core keeps its width while energy moves to the halo (0.38 → 0.40 mm at 20 ns, Strehl already 0.73)33.1
Contrast, CNR, gCNRlesion detectabilitydynamic range and compression (contrast, CNR); gCNR biased upward for small regionsModule 26
Speckle brightnessfocus quality of a regionconfounded by gain, attenuation and TGC changes33.3
Geometric fidelity (target positions)tilts and image warpinginvisible to every focus metric above; needs known targets or registrationModule 31
Pitfall

Reporting an aberration-correction gain with a coherence-weighted image (CF or SLSC) and a dB scale chosen after the fact. CF and SLSC improve when the law improves, but also when noise drops, so they double-count; and contrast in dB depends on the display range. Report Strehl-like or Γ-like focus measures, gCNR with the region sizes, and target positions, each before and after.

33.8Visual summary

Twelve questions an aberration-correction expert answers in five seconds. Each figure is computed from the formula, model or simulation named in its badge.

33.9How experts think about it

An aberration expert asks four questions in order. How big is the delay error in fractions of a period, and on what spatial scale does it vary (\(\sigma_\phi\), \(\ell\))? Where is the aberrator relative to the probe and the target, which sets the isoplanatic patch (\(z\ell/d\))? What guide does the estimator use (a point, a speckle spot, a virtual reflector, a map), and how big and how symmetric is that guide, which sets the bias? And which metric would reveal a failure, given that width metrics saturate and no focus metric sees a shift? Estimation precision is rarely the limit; bias, locality and the forward model are.

Aberration is a delay measured in fractions of a period: 10 ns is 0.31 rad at 5 MHz.
σφ = 2πf₀στ40 ns → 1.26 rad at 5 MHz0.63 rad at 2.5 MHz

What destroys the coherent sum is the phase spread, and the same delay is twice the phase at twice the frequency. With \(S \approx e^{-\sigma_\phi^{2}}\), halving \(f_0\) quarters the exponent: a 40 ns screen leaves 0.67 at 2.5 MHz and 0.21 at 5 MHz (large aperture). High-frequency probes need correction most.

Exact formula
Narrowband; a broadband pulse averages the phase over its band, and delays beyond half a period make phase-only reasoning ambiguous.
Strehl is e−σ² for a weak screen, never much below √π ℓ/(σD) for a strong one, and squared for two-way paths.
40 ns: 0.315 exact vs 0.20680 ns: 0.079 vs 0.073two way E[S²] = 0.131 at 40 ns

\(\mathbb{E}[S] = N^{-2}\sum_{ij} e^{-\sigma_\phi^{2}(1-\rho_{ij})}\) (33.2): element pairs closer than \(\ell\) keep adding in phase, so a fraction of order \(\ell/(\sigma_\phi D)\) survives. 128 elements, \(\ell\) = 4 mm (m33_phase_screen_psf.m). Transmit and receive through the same screen multiply their focal amplitudes, so the two-way Strehl of a realization is \(S^{2}\) (mean 0.131 at 40 ns, against \(\mathbb{E}[S]^{2}\) = 0.099).

Exact formula
Gaussian statistics and correlation; real body walls are layered, non-Gaussian and also modulate amplitude.
The PSF keeps its width and loses its energy, until the halo takes over.
20 ns: −6 dB width 0.40 mm, S 0.7340 ns: 1.4 mmhalo rms s = √2 c στ z/ℓ

For \(\sigma_\phi < 1\) the lost energy spreads into a broad, low halo while the core keeps the diffraction width (0.38 mm at F/1, 40 mm). For \(\sigma_\phi \gg 1\) the core is gone and the mean intensity is the Gaussian halo set by the rms wavefront tilt \(\sqrt2\,c\sigma_\tau/\ell\) (33.3): −6 dB width 3.33 s = 5.8 mm at 80 ns against 5.5 mm exact.

SimulationRe-plotted from the m33_phase_screen_psf.m output (Monte Carlo, 400 screens).
Mean PSF; a single realization is a lumpy speckle-like spot, and its peak can sit off axis.
Aberration multiplies the van Cittert–Zernike curve by e−σ²(1−ρ).
lag 10: 0.92 → 0.47long lags → e−σ² R₀

The channel signal is \(e^{j\phi_i}u_i\) with \(u\) the speckle of Module 07; for a Gaussian screen \(\mathbb{E}[e^{j(\phi_i-\phi_j)}] = e^{-\sigma_\phi^{2}(1-\rho)}\) multiplies the triangle (33.3). At 40 ns, \(\ell\) = 4 mm: lag 10 is 0.465 by Monte Carlo against 0.467, lag 40 is 0.138 against 0.142. This is why coherence is a focus metric (Module 32) and an aberration detector.

Exact formula
Screen on receive with an ideal transmit; an aberrated transmit narrows \(R_0\) as well, and noise lowers every lag above zero.
Nearest-neighbour integration is a random walk: σε√(N/15) after detrending.
raw error ∝ √k128 elements: 2.9 σεelement noise telescopes

Summing \(N-1\) independent differences (33.4) gives a walk; Brownian motion minus its least-squares line has mean square \(T/15\) (Monte Carlo: \(N/15\) to 0.1 %). At 0 dB channel SNR the pair jitter was 4.5 ns, bound 13 ns, measured 8.3 ns: element noise enters adjacent pairs with opposite signs and partly cancels.

Exact formula
Independent pair errors; a smooth bias does not drift but does not average out either.
Least squares over many lags breaks the walk; nothing breaks a bias.
0 dB: NN 8.3 → LS 4.0 nsnoiseless: 5.7 vs 6.3 ns

Extra lags give paths between distant elements that skip the intermediate errors (33.5). With noise (m33_nn_crosscorr_correction.m, pass 2) LS halves the residual and lifts the Strehl from 0.93 to 0.98. Without noise the residual is a bias common to all lags (the blurred transmit spot), and LS changes nothing.

SimulationRe-plotted from the m33_nn_crosscorr_correction.m output (16 screens, 128 elements).
Weights need the pair correlations; very long lags decorrelate and add little.
Every speckle estimator measures the law of the spot it insonifies, so iterate the transmit.
broadband sim: 5.7 → 2.6 nsDemo 33.3: 0.28 → 0.017 → 0.001 rad2–3 iterations (Måsøy 2005)

\(\mathbf{C} = \mathbf{\Phi}\mathbf{C}_{\mathrm{vr}}\mathbf{\Phi}^{H}\), so the eigenvector phase is the screen plus the phase of the virtual reflector's own eigenvector (33.6), flat only for a compact, symmetric spot. Each transmit with a better law shrinks the spot. Måsøy, Varslot and Angelsen reported near-ideal correction after two to three transmit iterations.

SimulationDemo 33.3 readouts, 60 ns, ℓ = 3 mm, first screen.
Strong aberration can lock the first pass onto a side lump of the spot; the transmit must see the same patch.
Brightness, time reversal and DORT all return the top eigenvector of a covariance.
B = wHCwone eigenvalue per resolved scatterer

Mean brightness is a quadratic form in the aperture covariance (33.6); iterative time reversal is power iteration on \(\mathbf{K}^{H}\mathbf{K}\); DORT diagonalizes it (33.7). Prada and Fink found four significant eigenvalues out of 96 for four wires. The figure: three points behind a 40 ns screen, 64 elements, three eigenvalues.

SimulationMonochromatic K of three point scatterers (30 mm, 3 mm apart, reflectivities 1, 0.7, 0.4) behind a screen; eigenvalues by power iteration with deflation, JS.
Resolved point scatterers; in speckle the spectrum decays slowly and the size of the virtual reflector decides the bias.
Outside the isoplanatic patch a correction is worse than none.
residual variance 2σ²(1 − ρ)crossover at ρ = 1/2Δx× = ℓz√ln2/d

A borrowed law leaves \(\phi(a+\delta) - \phi(a)\), of variance \(2\sigma_\phi^{2}(1-\rho(\delta))\), which exceeds \(\sigma_\phi^{2}\) once \(\rho < 1/2\) (33.8). For d = 10 mm, \(\ell\) = 4 mm, z = 40 mm: 13.3 mm in the large-aperture limit, 11.0 mm with the finite aperture (m33_isoplanatic_patch.m).

Exact formula
Thin Gaussian screen; thick aberrators decorrelate in depth as well.
The isoplanatic patch scales as zℓ/d.
Δx½ = 16.5, 8.3, 4.3 mmd = 5, 10, 20 mm40 ns, ℓ 4 mm, z 40 mm

Moving the point by \(\Delta x\) slides the aperture footprint on the screen by \(\Delta x\,d/z\) (33.8), so the patch is the correlation length magnified by \(z/d\). Shallow targets under a thick near-field aberrator (the skull, for cortical vessels) are the worst case; deep targets under a thin wall are the best.

Exact formula
Straight rays and one target depth; refraction and finite screen thickness change the factor.
Straight rays always overestimate the travel time.
path error ≤ 03.4 ns at the centre, 11 ns at −12 mm

Fermat's ray is the fastest path, so any other path, including the straight one, takes longer (33.9). With the exact map, straight-ray delays left 3.4 ns rms at the centre and 11 ns at −12 mm behind the wavy fat layer (m33_sos_map_correction.m): the forward-model error of Module 32 reappears as residual aberration.

Approx. modelDemo 33.5 geometry at xp = −12 mm: straight-ray minus Fermat delay across the aperture.
First arrivals in smooth media; caustics and multiple arrivals need full-wave models.
A map corrects everywhere it is right; a measured law corrects exactly, but only in its patch.
mean Strehl: one law 0.43straight-ray map 0.97blurred map 0.93

Across ±12 mm behind the wavy layer: no correction 0.08, one law 0.43 (S ≥ 0.8 only within ±3 mm), the true map with straight rays 0.97, a 3 mm-resolution map 0.93, a 10 m/s speed bias 0.95. Each route fails differently: the law by locality, the map by resolution, bias and path model.

SimulationRe-plotted from the m33_sos_map_correction.m output.
One geometry; the balance shifts with aberrator depth, contrast and map quality.
The diagonal of Rxx is the image; the width of its off-diagonal band is the focus.
I(Δx) = |Pin|² ⋆ |Pout|²F = w₀/w (Lambert 2020)

In speckle, the mean intensity at offset \(\Delta x\) from the diagonal is the cross-correlation of the transmit and receive PSF intensities (33.10), so the focus can be measured anywhere without a guide star. The figure uses the exact mean PSF of a 30 ns and a 60 ns screen.

Exact formula
Single scattering; multiple scattering adds a broad off-diagonal background, which Lambert et al. used to map it.
Re-centre, correlate, take the first singular vector: the distortion matrix turns speckle into a guide star.
0.83 → 0.31 → 0.07 → 0.018 radΓ 0.48 → 1.00

Removing each focal point's geometric wavefront makes all points of a patch look like one virtual source (33.11); the correlation over columns averages speckle; the first singular vector carries the law; correcting the input shrinks the virtual source; repeat. Numbers from m33_distortion_matrix_svd.m, screen at the array.

SimulationRe-plotted from the m33_distortion_matrix_svd.m output (d = 0).
Needs many independent resolution cells per window and one patch per window; noise and multiple scattering bias the correlation.
Once the field spans several patches, local windows beat one global law.
Γ: global 0.69, local mean 0.86bright point −7.6 dB vs +3.1 dB

With the screen at 8 mm (patch about 5 mm at 20 mm depth), one law for the 12 mm field is a compromise that defocuses part of it; three 4 mm windows each find their own mode (m33_distortion_matrix_svd.m). The window size is a trade: smaller is more local, larger averages more speckle.

SimulationRe-plotted from the m33_distortion_matrix_svd.m output (d = 8 mm).
Patch boundaries are not known in advance; the singular spectrum and the focus maps guide the choice.
The singular spectrum of D counts patches only after you account for the virtual reflector.
σ₂/σ₁: 0.68 → 0.52 after correctionC/|C| is not positive semidefinite

Even a perfectly corrected isoplanatic window has a virtual reflector the size of the diffraction-limited spot, whose van Cittert–Zernike covariance has several significant eigenvalues (Demo 33.3 settles at \(\lambda_2/\lambda_1\) = 0.30 for a sinc² spot). Lambert et al. normalized \(\mathbf{C}\) entrywise to make the reflector point-like and counted three patches from the entropy of its eigenvalues (\(H \approx 2.85\)); that normalized matrix can have negative eigenvalues.

SimulationRe-plotted from the m33_distortion_matrix_svd.m output.
Counting also depends on window size, noise and the frequency band used.
Width metrics saturate, energy metrics do not, and no focus metric sees a common shift.
w₀/w = 0.93 at 20 ns, 0.81 at 30 nsE[S] = 0.73 and 0.50

The figure shows, for the exact mean PSF of a 128-element aperture, the width ratio and the Strehl versus \(\sigma_\tau\): at 30 ns the width ratio is still 0.81 while half the focal energy is gone; at 40 ns it drops to 0.27 as the core breaks up. A residual linear phase moves both foci together, warping the image without blurring it (Module 31's bias field); only known targets or registration reveal it.

Exact formula
Mean-PSF statements; single realizations and point targets fluctuate.
A learned aberration corrector inherits its simulator.
forward model → training setkeep the data in the loop

A network trained on simulated screens learns the statistics of those screens: Gaussian thin screens are not an abdominal wall, a skull or a lens mismatch. Unrolled or differentiable versions of the iterations in this module (Module 14) keep the measured matrix in the loop, so a domain gap shows up as slower convergence instead of a confident wrong law.

A statement about generalization, not capacity.

33.10Cheat sheet

ConceptFormulaUltrasound use
Screen phase\(\sigma_\phi = 2\pi f_0\sigma_\tau\)0.314 rad per 10 ns at 5 MHz
Mean Strehl (Gaussian screen)\(N^{-2}\sum_{ij}e^{-\sigma_\phi^{2}(1-\rho_{ij})}\)0.315 at 40 ns, \(\ell\) 4 mm, 128 el.
Weak and strong limits\(e^{-\sigma_\phi^{2}}\); \(\sqrt\pi\,\ell/(\sigma_\phi D)\)two way ≈ \(S^{2}\)
Halo of a strong screen\(s = \sqrt2\,c\sigma_\tau z/\ell\)−6 dB width 3.33 s
Coherence under aberration\(R(m) = R_0(m)\,e^{-\sigma_\phi^{2}(1-\rho(mp))}\)lag 10: 0.92 → 0.47 at 40 ns
Nearest-neighbour drift\(\sigma_\varepsilon\sqrt{N/15}\) (detrended)2.9 \(\sigma_\varepsilon\) for 128 elements
Multi-lag least squares\((\mathbf{A}^{T}\mathbf{W}\mathbf{A} + \mathbf{1}\mathbf{1}^{T})^{-1}\mathbf{A}^{T}\mathbf{W}\hat{\mathbf{d}}\)halves the noisy residual (8.3 → 4.0 ns)
Brightness\(B = \mathbf{w}^{H}\mathbf{C}\mathbf{w}\), max at \(\mathbf{u}_1\)speckle-brightness and eigenvector laws coincide
DORT\(\mathbf{K}^{H}\mathbf{K}\mathbf{v} = \lambda\mathbf{v}\), \(\mathbf{v}_s \propto \mathbf{h}_s^{*}\)one eigenvalue per resolved scatterer
Footprint shift\(\delta = \Delta x\,d/z\)patch ∝ \(z\ell/d\)
Borrowed-law Strehl\(e^{-2\sigma_\phi^{2}(1-\rho(\delta))}\)worse than none when \(\rho < 1/2\)
Crossover distance\(\Delta x_\times = \ell z\sqrt{\ln 2}/d\)13.3 mm (d 10, \(\ell\) 4, z 40 mm)
Map-based delays\(T = \int_{\Gamma}\hat s\,dl\)straight rays overestimate (path error ≤ 0)
Focused reflection matrix\(\mathbf{R}_{xx} = \mathbf{G}_0^{H}\mathbf{K}\mathbf{G}_0^{*}\)diagonal = image, band = focus
Distortion matrix\(\mathbf{D} = (\mathbf{K}\mathbf{G}_0^{*}) \circ \mathbf{G}_0^{*}\), law = \(\arg\mathbf{u}_1\)0.83 → 0.018 rad in 3 iterations
Power iterationerror \(\propto (\lambda_2/\lambda_1)^{n}\)10 iterations for \(10^{-3}\) at ratio 0.5

33.11Interview and qualifying-exam questions

I1Model an abdominal-wall aberrator in two numbers and predict the focal loss at 5 MHz for 40 ns rms. · screen

A thin delay screen at the array: rms \(\sigma_\tau\) and correlation length \(\ell\) (Gaussian). \(\sigma_\phi = 2\pi f_0\sigma_\tau\) = 1.26 rad. Infinite-aperture Strehl \(e^{-1.58}\) = 0.21; with 128 elements and \(\ell\) = 4 mm the exact mean is 0.315 because nearby elements stay in phase; two way about 0.1. The −6 dB width hardly changes until about 30 ns, then the halo takes over. Follow-up: what if the frequency halves? \(\sigma_\phi\) halves and the exponent quarters: about 0.67.

I2Why does nearest-neighbour cross-correlation drift, and how do you fix it? · screen

Integrating \(N-1\) independent differences is a random walk: error \(\propto\sqrt{k}\), \(\sigma_\varepsilon\sqrt{N/15}\) after removing the mean and ramp. Fixes: weighted least squares over many lags (graph Laplacian plus a constant constraint), better SNR, longer kernels. Element noise telescopes and does not drift. Follow-up: in bright speckle the jitter is sub-ns but the residual is 5.7 ns; why? A bias from the blurred transmit spot, removed by transmitting with the estimate and re-measuring.

I3Show that speckle-brightness maximization and DORT give the same aberration law. · onsite

Brightness of the beamformed pixel is \(\mathbf{w}^{H}\mathbf{C}\mathbf{w}\) with \(\mathbf{C}\) the aperture covariance, maximized over \(\lVert\mathbf{w}\rVert^{2} = N\) by the dominant eigenvector. Iterative time reversal is power iteration on \(\mathbf{K}^{H}\mathbf{K}\), and DORT diagonalizes it; for speckle behind a focused transmit, \(\mathbf{C} = \mathbf{\Phi}\mathbf{C}_{\mathrm{vr}}\mathbf{\Phi}^{H}\), so the eigenvector is the screen times the eigenvector of the virtual reflector. Follow-up: when is that second factor flat? When the insonified spot is compact and symmetric, which is why you iterate.

I4Derive the isoplanatic patch for a thin screen at depth d. When is a correction worse than none? · onsite

The straight ray from \(x_u\) to \((x_p, z)\) crosses the screen at \(x_u(1-d/z) + x_p d/z\); moving the point by \(\Delta x\) shifts the footprint by \(\delta = \Delta x\,d/z\). The borrowed law leaves a residual of variance \(2\sigma_\phi^{2}(1-\rho(\delta))\): Strehl \(e^{-2\sigma_\phi^{2}(1-\rho)}\), worse than uncorrected when \(\rho < 1/2\), at \(\Delta x_\times = \ell z\sqrt{\ln 2}/d\) = 13.3 mm for \(\ell\) 4, z 40, d 10 mm. Follow-up: finite aperture? The crossover moves in (11.0 mm) because the uncorrected Strehl is higher than \(e^{-\sigma_\phi^{2}}\).

I5You have a CUTE sound-speed map from Module 32. Turn it into an aberration correction and write its error budget. · onsite

Integrate slowness from every element to every pixel (straight rays or eikonal) and beamform with those times. Errors: map resolution (a low-pass truth: 3 mm blur cost 8.2 against 5.2 ns rms in m33_sos_map_correction.m), map bias (a 10 m/s error in the fat speed: 7.3 ns, mostly defocus), path model (straight rays always overestimate times: 3–11 ns behind a wavy fat layer), and the mean speed CUTE cannot see. Follow-up: how would you check it in vivo? Measure the residual law with a matrix method; a large smooth residual points to the mean or the path model.

I6What is the focused reflection matrix, and how does it measure focusing without a point target? · onsite

\(\mathbf{R}_{xx} = \mathbf{G}_0^{H}\mathbf{K}\mathbf{G}_0^{*}\): responses between virtual transducers at \(x_{\mathrm{in}}\) and \(x_{\mathrm{out}}\) after focusing in software. The diagonal is the confocal image; in speckle the mean intensity versus offset is the cross-correlation of transmit and receive PSF intensities, so its width (Lambert et al.'s \(F = w_0/w\)) or its confocal energy fraction measures the local focus anywhere. Follow-up: what spoils it? Multiple scattering adds a broad background (useful in itself) and noise adds a floor.

I7Walk through distortion-matrix correction. Why re-centre, why SVD, why iterate, and how big should the window be? · research

Focus on transmit, keep element receive, multiply each column by the conjugate geometric wavefront of its focal point: every point of a patch now looks like one virtual source at the origin through the same screen. \(\mathbf{D}\mathbf{D}^{H}\) averages the speckle over the columns; its first singular vector is screen × virtual-reflector eigenvector. Correcting the input shrinks the reflector, so iterate (0.83 → 0.018 rad in three steps in the script). Window: many resolution cells for averaging, smaller than a patch; Lambert et al. shrank it from 12 × 16 to 6.2 × 8.1 mm. Follow-up: transmit vs receive laws with plane-wave transmits? Estimate both, alternating sides.

I8Your speckle-based estimator converges, but the image is shifted relative to CT. What happened? · research

A residual linear phase (tilt) is unobservable from speckle: it moves both foci together, so focus metrics (Γ, F, CF) look perfect while every point is displaced, a bias field as in Module 31's transcranial ULM. Remedies: anchor the tilt with known targets, a sound-speed map, or bubbles as guide stars with known geometry; in matrix methods watch for linear ramps removed during the estimation. Follow-up: how big? 10 ns across 10 mm at 50 mm depth moves a point by 77 µm (Module 31).

I9What breaks when you move from abdominal wall to transcranial imaging, and what do you do instead? · research

Bone is about twice as fast, attenuating and variable in thickness, so amplitude matters as much as phase, the aberrator is thick (finite patches in depth) and \(d/z\) is large for shallow targets (small patches). Options: CT-based priors (map route), guide stars (Demené et al. 2021 derived one law per patch from microbubbles before ULM), 3D matrix imaging with 2D arrays (Bureau et al. 2023). Follow-up: what is the cheapest first step? Lower frequency: the same delays cost a quarter of the phase variance at half the frequency.

I10Which metrics would you report for a new aberration-correction method, and which would you avoid? · research

A Strehl-like or energy-fraction focus measure (Γ or the confocal energy), point-target FWHM where available, gCNR with region sizes (Module 26), and target positions for geometric fidelity; all before and after, at several depths and lateral positions to show the patch. Avoid CF or SLSC images as the main evidence (they reward noise suppression too) and contrast in dB on a display range chosen afterwards. Follow-up: why not FWHM alone? It saturates: 0.38 → 0.40 mm at 20 ns while the Strehl is already 0.73.

I11Propose a method that combines a sound-speed map and a matrix estimate. What does each contribute? · research

Use the map (Module 32) to compute space-variant starting delays, so the residual aberration is small and the patch large; then estimate the residual laws with local distortion matrices, which converge fast because the virtual reflector is already compact. Feed the residual delays back as tomographic data with the correct (bent-ray) forward model to update the map, and iterate. The map supplies locality and the tilt/mean anchor; the matrix supplies the measured truth where the map's resolution, bias or path model fail. Follow-up: what would you validate first? On a phantom with a known aberrator, compare map-only, matrix-only and combined Γ maps and target positions.

33.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: m33_phase_screen_psf.m · m33_nn_crosscorr_correction.m · m33_distortion_matrix_svd.m · m33_isoplanatic_patch.m · m33_sos_map_correction.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m33_phase_screen_psf.mGaussian screen on 128 elements: Monte Carlo Strehl against the exact mean and both limits, mean PSF width against the halo, receive coherence against (33.3) · change sig_list, ell, f0, z · runtime ≈ 5 s 
%% m33_phase_screen_psf.m
% Thin near-field phase screen: Strehl ratio, PSF broadening and the loss of
% spatial coherence, by Monte Carlo against the closed forms.
% Purpose : A Gaussian random delay screen tau(x) (rms sigma_tau, correlation
%           rho(dx) = exp(-dx^2/l^2)) sits on a 128-element array (pitch
%           0.3 mm, f0 = 5 MHz). Focal-plane (Fraunhofer) PSF of the focused
%           aperture, one-way p(x) = sum_i exp(j phi_i) exp(-j k x_i x / z).
%           (A) On-axis Strehl S = |sum exp(j phi_i)|^2 / N^2 against
%                exp(-sigma_phi^2) (infinite aperture) and the exact mean
%                E[S] = (1/N^2) sum_ij exp(-sigma_phi^2 (1 - rho_ij)).
%           (B) Mean PSF intensity: Monte Carlo against the exact
%                E|p(x)|^2 = sum_ij exp(-j k x (x_i - x_j)/z) exp(-sigma^2(1-rho_ij)),
%                and its -6 dB width (intensity 1/4) against the halo of the
%                large-sigma limit, a Gaussian of rms s = sqrt(2) c sigma_tau z / l
%                (-6 dB width 3.33 s).
%           (C) Receive coherence of speckle (ideal transmit, screen on
%                receive): R(m) = (1 - m/N) exp(-sigma_phi^2 (1 - rho(m p))).
% Change  : sig_list (ns), ell (mm), f0, z.
% Expect  : Strehl follows exp(-sigma_phi^2) for small sigma and then the
%           slower sqrt(pi) l/(sigma_phi D); the mean PSF keeps its core until
%           sigma_phi ~ 1 rad and then becomes the halo; coherence falls as predicted.
% Runtime : about 5 s in Octave.
clear; close all; rng(1);

c = 1540; f0 = 5e6; lam = c/f0; k = 2*pi/lam;
N = 128; p = 0.3e-3; xe = ((0:N-1) - (N-1)/2)*p; D = N*p;   % element positions (m), aperture 38.4 mm
z = 40e-3; ell = 4e-3;                                       % focal depth, screen correlation length (m)
sig_list = [10 20 30 40 60 80]*1e-9;                         % rms delay (s)
R = 400;                                                     % realizations

% screen generator: white noise on a fine grid convolved with a Gaussian kernel
dxs = p/4; xs = (-60e-3:dxs:60e-3);                          % fine grid (m), much longer than the aperture
g = exp(-2*(((-3*ell):dxs:(3*ell))/ell).^2); g = g/sqrt(sum(g.^2));   % autocorrelation exp(-dx^2/l^2), unit variance
idx = round((xe - xs(1))/dxs) + 1;                           % element samples on the fine grid
screen = @(sig) sig*conv(randn(1, numel(xs)), g, 'same');    % one realization tau(xs), rms sig

DX = xe' - xe; RHO = exp(-(DX/ell).^2);                      % rho_ij
fprintf('Aperture %.1f mm, z = %.0f mm (F/%.2f), l = %.0f mm, lambda = %.3f mm\n', D*1e3, z*1e3, z/D, ell*1e3, lam*1e3);
fprintf('Large-sigma asymptote of E[S]: sqrt(pi) l / (sigma_phi D) (one correlation cell of N)\n\n');
fprintf(' sig_tau  sig_phi   S MC     E[S] exact  exp(-s^2)  asympt.   S2way MC  E[S]^2\n');
xf = linspace(-4e-3, 4e-3, 801);                             % lateral focal coordinate (m)
E = exp(-1j*k*xf'*xe/z);                                     % Fraunhofer kernel (801 x N)
w6 = zeros(size(sig_list)); w6x = w6; blur = w6; S_mc = w6; S_ex = w6;
for a = 1:numel(sig_list)
    sig = sig_list(a); sp = 2*pi*f0*sig;
    S = zeros(R, 1); S2 = S; I = zeros(numel(xf), 1);
    chk = 0;
    for r = 1:R
        t = screen(sig); ph = 2*pi*f0*t(idx);
        chk = chk + mean(t(idx).^2);
        S(r) = abs(sum(exp(1j*ph)))^2/N^2; S2(r) = S(r)^2;
        I = I + abs(E*exp(1j*ph(:))).^2/N^2;
    end
    I = I/R; Sex = sum(sum(exp(-sp^2*(1 - RHO))))/N^2;
    S_mc(a) = mean(S); S_ex(a) = Sex;
    fprintf('%5.0f ns  %5.3f    %.3f    %.3f       %.3f      %.3f     %.3f     %.3f\n', sig*1e9, sp, mean(S), Sex, exp(-sp^2), sqrt(pi)*ell/(sp*D), mean(S2), Sex^2);
    % exact mean PSF intensity
    Iex = real(sum((E*exp(-sp^2*(1 - RHO))).*conj(E), 2))/N^2;
    above = xf(I >= max(I)/4); w6(a) = max(above) - min(above);           % -6 dB (intensity 1/4) width, MC
    above = xf(Iex >= max(Iex)/4); w6x(a) = max(above) - min(above);      % exact
    blur(a) = sqrt(2)*c*sig*z/ell;                               % rms width s of the Gaussian intensity halo
    if a == 1, fprintf('   (rms of generated screen at the elements: %.2f ns)\n', sqrt(chk/R)*1e9); end
    if a == 4, I40 = I; Iex40 = Iex; end
end
Pid = abs(E*ones(N, 1)).^2/N^2; above = xf(Pid >= 1/4); w0 = max(above) - min(above);
fprintf('\n-6 dB width of the mean one-way PSF intensity (mm); unaberrated %.3f mm (1.21 lambda F = %.3f)\n', w0*1e3, 1.207*lam*z/D*1e3);
fprintf(' sig_tau   MC     exact   halo 3.33 s, s = sqrt(2) c sig z / l (-6 dB width of the Gaussian halo)\n');
for a = 1:numel(sig_list)
    fprintf('%5.0f ns  %.3f   %.3f   %.3f\n', sig_list(a)*1e9, w6(a)*1e3, w6x(a)*1e3, 3.330*blur(a)*1e3);
end

% (C) receive coherence of speckle with an ideal transmit and the screen on receive
sig = 40e-9; sp = 2*pi*f0*sig; Ns = 1500; Rc = 300; lags = 0:40; Rm = zeros(size(lags)); Pa = Rm; Pb = Rm;
for r = 1:Rc
    t = screen(sig); ph = 2*pi*f0*t(idx);
    xsc = (rand(Ns, 1) - 0.5)*30*lam*z/D; b = sin(pi*xsc*D/(lam*z) + eps)./(pi*xsc*D/(lam*z) + eps);   % ideal transmit beam amplitude
    a_s = (randn(Ns, 1) + 1j*randn(Ns, 1)).*b;
    s = exp(1j*ph(:)).*(exp(-1j*k*xe'*xsc'/z)*a_s);          % received aperture vector
    for m = lags
        Rm(m+1) = Rm(m+1) + real(sum(s(1:N-m).*conj(s(1+m:N))));   % accumulate numerator over realizations
        Pa(m+1) = Pa(m+1) + sum(abs(s(1:N-m)).^2); Pb(m+1) = Pb(m+1) + sum(abs(s(1+m:N)).^2);
    end
end
Rm = Rm./sqrt(Pa.*Pb);                                        % ensemble (depth-averaged) coherence estimate
Rth = (1 - lags/N).*exp(-sp^2*(1 - exp(-(lags*p/ell).^2)));
R0 = 1 - lags/N;
fprintf('\nReceive coherence, 40 ns screen (sigma_phi = %.2f rad), l = 4 mm\n lag   MC     formula   no screen\n', sp);
for m = [1 2 5 10 13 20 40]
    fprintf('%3d   %.3f   %.3f     %.3f\n', m, Rm(m+1), Rth(m+1), R0(m+1));
end

figure('Position', [50 50 1000 330]);
subplot(1, 3, 1); sl = linspace(0, 90, 200)*1e-9; spl = 2*pi*f0*sl; Sl = zeros(size(sl));
for q = 1:numel(sl), Sl(q) = sum(sum(exp(-spl(q)^2*(1 - RHO))))/N^2; end
plot(sl*1e9, exp(-spl.^2), 'k--', sl*1e9, Sl, 'b-', sig_list*1e9, S_mc, 'ro'); xlabel('\sigma_\tau (ns)'); ylabel('Strehl S'); title('On-axis Strehl, one way'); legend('exp(-\sigma_\phi^2)', 'exact mean', 'Monte Carlo');
subplot(1, 3, 2); plot(xf*1e3, 10*log10(Pid), 'k', xf*1e3, 10*log10(I40), 'r', xf*1e3, 10*log10(Iex40), 'b--'); ylim([-40 0]); xlabel('x (mm)'); ylabel('dB'); title('Mean PSF intensity, 40 ns'); legend('no screen', 'MC', 'exact');
subplot(1, 3, 3); plot(lags, R0, 'k', lags, Rth, 'b-', lags, Rm, 'ro'); xlabel('lag m (elements)'); ylabel('coherence'); title('Receive coherence, 40 ns'); legend('no screen', 'formula', 'MC');

What to look for: Strehl (Monte Carlo / exact / \(e^{-\sigma_\phi^{2}}\)) 0.920/0.922/0.906 at 10 ns, 0.730/0.727/0.674 at 20 ns, 0.311/0.315/0.206 at 40 ns, 0.078/0.079/0.002 at 80 ns, where the asymptote gives 0.073; two way 0.131 at 40 ns. The −6 dB width of the mean PSF: 0.380 mm without the screen, 0.400 at 20 ns, 1.38 at 40 ns, 5.58 at 80 ns (halo prediction 5.80). Coherence at 40 ns: lag 10 0.465 (formula 0.467), lag 40 0.138 (0.142).

m33_nn_crosscorr_correction.mBroadband speckle channel data through a 40 ns screen: nearest-neighbour integration versus multi-lag least squares, one transmit iteration, noiseless and 0 dB SNR · change sig_tau, ell, Lmax, snr_list · runtime ≈ 10 s 
%% m33_nn_crosscorr_correction.m
% Aberration law from speckle: nearest-neighbour cross-correlation (Flax and
% O'Donnell), its random-walk drift, multi-lag least squares, and one
% transmit iteration.
% Purpose : 128 elements, pitch 0.3 mm, 5 MHz, 60 % bandwidth, focus 40 mm.
%           A Gaussian delay screen (rms sig_tau, correlation length l) sits on
%           the array, on transmit and on receive. Speckle IQ data are
%           synthesized per element after the geometric focusing delays:
%             s_i(t) = sum_s b_s a_s env(t - t_is) exp(-j 2 pi f0 t_is),
%           t_is = residual geometric delay of scatterer s after dynamic receive
%           focusing (lateral offset only) + its two-way time + tau_i, and b_s the
%           (aberrated or corrected) transmit beam at the scatterer.
%           NN estimate : d_i = angle(sum_t s_i conj(s_i+1)) / (2 pi f0),
%                         tau_hat = cumulative sum of d_i (a random walk).
%           LS estimate : all pairs up to lag L, weighted by |rho_ij|^2/(1-|rho_ij|^2),
%                         solved for tau with sum(tau) = 0.
%           Mean and linear ramp are removed from truth and estimates before
%           comparison (a ramp only steers the beam).
%           Prediction for the detrended NN error: sigma_eps sqrt(N/15)
%           (Brownian motion minus its least-squares line).
% Change  : sig_tau, ell, Lmax, R, Ns.
% Expect  : noiseless speckle: sub-ns pair jitter, but a smooth bias of a few
%           ns from the blurred (aberrated) transmit focus; LS cannot remove a
%           bias, a second pass with the corrected transmit halves it.
%           0 dB channel SNR: jitter dominates, the NN error drifts as a random
%           walk (sigma_eps sqrt(N/15) after detrending is the bound for
%           independent pair errors; element noise enters adjacent pairs with
%           opposite signs and partly cancels) and LS over 8 lags cuts it.
% Runtime : about 12 s in Octave.
clear; close all; rng(2);

c = 1540; f0 = 5e6; lam = c/f0; k = 2*pi/lam;
N = 128; p = 0.3e-3; xe = ((0:N-1) - (N-1)/2)*p; D = N*p; zf = 40e-3;
sig_tau = 40e-9; ell = 4e-3; Lmax = 8; R = 16; Ns = 400;
snr_list = [Inf 0];                                          % channel SNR (dB): noiseless, and 0 dB per element
st = 1/(2*pi*0.6*f0/2.355);                                 % envelope std for 60 % (-6 dB) bandwidth
t = (-2e-6:0.05e-6:2e-6)';                                  % fast time around the focus (s), IQ at 20 MHz

dxs = p/4; xs = -60e-3:dxs:60e-3; g = exp(-2*(((-3*ell):dxs:(3*ell))/ell).^2); g = g/sqrt(sum(g.^2));
idx = round((xe - xs(1))/dxs) + 1;
detr = @(v) v - [ones(N, 1) xe(:)]*([ones(N, 1) xe(:)]\v(:));   % remove mean and linear ramp

% LS design matrix for all pairs with lag 1..Lmax (+ one row for sum(tau) = 0)
pairs = []; for L = 1:Lmax, pairs = [pairs; (1:N-L)' ((1:N-L) + L)']; end
A = zeros(size(pairs, 1), N); for q = 1:size(pairs, 1), A(q, pairs(q, 1)) = -1; A(q, pairs(q, 2)) = 1; end

TAU = zeros(N, R); XS = zeros(Ns, R); ZS = XS; AS = XS;      % the same screens and scatterers for every SNR
for r = 1:R
    tmp = sig_tau*conv(randn(1, numel(xs)), g, 'same'); TAU(:, r) = tmp(idx)';
    XS(:, r) = (rand(Ns, 1) - 0.5)*6e-3; ZS(:, r) = zf + (rand(Ns, 1) - 0.5)*4e-3;
    AS(:, r) = (randn(Ns, 1) + 1j*randn(Ns, 1))/sqrt(2);
end
for isnr = 1:numel(snr_list)
snr = snr_list(isnr);
drift = zeros(N, 1); res = zeros(R, 2, 3); eps_nn = zeros(R, 2); S_before = zeros(R, 1); S_after = zeros(R, 2, 2);
drift = zeros(N, 1);
for r = 1:R
    tau = TAU(:, r); phi = 2*pi*f0*tau; xsc = XS(:, r); zsc = ZS(:, r); a_s = AS(:, r);
    tx_law = zeros(N, 1);                                    % pass 1: transmit with the geometric law only
    S_before(r) = abs(mean(exp(1j*phi)))^2;
    for pass = 1:2
        b = exp(-1j*k*xsc*xe/zf)*exp(1j*(phi - tx_law))/N;   % transmit beam at each scatterer (Fraunhofer)
        b = b.*exp(-((zsc - zf)/2e-3).^2);                  % soft depth weighting of the kernel
        S = zeros(numel(t), N);
        for i = 1:N
            tis = (sqrt((xsc - xe(i)).^2 + zsc.^2) - sqrt(xe(i)^2 + zsc.^2) + 2*(zsc - zf))/c + tau(i);   % dynamic receive focus
            S(:, i) = (exp(-(t - tis').^2/(2*st^2)).*exp(-1j*2*pi*f0*tis'))*(b.*a_s);
        end
        if isfinite(snr), S = S + sqrt(mean(abs(S(:)).^2)/10^(snr/10))*(randn(size(S)) + 1j*randn(size(S)))/sqrt(2); end
        % NN estimate
        d = zeros(N-1, 1); rho1 = d;
        for i = 1:N-1
            q = sum(S(:, i).*conj(S(:, i+1))); d(i) = angle(q)/(2*pi*f0);    % estimate of tau_{i+1} - tau_i
            rho1(i) = abs(q)/sqrt(sum(abs(S(:, i)).^2)*sum(abs(S(:, i+1)).^2));
        end
        eps_nn(r, pass) = std(d - diff(tau));
        tnn = [0; cumsum(d)];
        % LS multi-lag estimate
        dd = zeros(size(pairs, 1), 1); w = dd;
        for q = 1:size(pairs, 1)
            i = pairs(q, 1); j = pairs(q, 2); cc = sum(S(:, i).*conj(S(:, j)));
            dd(q) = angle(cc)/(2*pi*f0); rh = min(abs(cc)/sqrt(sum(abs(S(:, i)).^2)*sum(abs(S(:, j)).^2)), 0.999);
            w(q) = rh^2/(1 - rh^2);
        end
        W = sqrt(w); tls = [A.*W; ones(1, N)]\[dd.*W; 0];
        et = detr(tau); e1 = detr(tnn) - et; e2 = detr(tls) - et;
        res(r, pass, :) = [std(et) std(e1) std(e2)];
        S_after(r, pass, :) = [abs(mean(exp(1j*2*pi*f0*e1)))^2 abs(mean(exp(1j*2*pi*f0*e2)))^2];
        if pass == 1, drift = drift + e1.^2; mr1(r) = mean(rho1); else, mr2(r) = mean(rho1); end
        tx_law = 2*pi*f0*tls;                                % pass 2: transmit with the LS estimate
    end
end
fprintf('\n===== channel SNR = %g dB =====\n', snr);
fprintf('Screen: rms %.0f ns, l = %.0f mm; %d elements; %d realizations\n', sig_tau*1e9, ell*1e3, N, R);
fprintf('rms of the detrended true screen: %.1f ns; Strehl of the raw screen %.3f\n', mean(res(:, 1, 1))*1e9, mean(S_before));
for pass = 1:2
    if pass == 1, lab = 'pass 1 (aberrated transmit)'; mr = mean(mr1); else, lab = 'pass 2 (transmit corrected by pass-1 LS)'; mr = mean(mr2); end
    fprintf('%s: mean adjacent-element |rho| = %.3f\n', lab, mr);
    fprintf('  NN pair jitter sigma_eps = %.2f ns\n', mean(eps_nn(:, pass))*1e9);
    fprintf('  NN residual (detrended) %.2f ns, predicted sigma_eps*sqrt(N/15) = %.2f ns\n', mean(res(:, pass, 2))*1e9, mean(eps_nn(:, pass))*sqrt(N/15)*1e9);
    fprintf('  LS residual, lags 1..%d: %.2f ns (%.1f times smaller than NN)\n', Lmax, mean(res(:, pass, 3))*1e9, mean(res(:, pass, 2))/mean(res(:, pass, 3)));
    fprintf('  Strehl after correction: NN %.3f, LS %.3f\n', mean(S_after(:, pass, 1)), mean(S_after(:, pass, 2)));
end

if isnr == numel(snr_list), drift_noisy = drift/R; end
end
figure('Position', [50 50 900 330]);
subplot(1, 2, 1); plot(1:N, sqrt(drift_noisy)*1e9, 'r'); xlabel('element'); ylabel('rms error (ns)');
title('NN error after detrending, 0 dB SNR, pass 1');
subplot(1, 2, 2); plot(xe*1e3, detr(tau)*1e9, 'k', xe*1e3, detr(tnn)*1e9, 'r', xe*1e3, detr(tls)*1e9, 'b--');
xlabel('element position (mm)'); ylabel('delay (ns)'); title('Last realization, pass 2'); legend('true', 'NN', 'LS');

What to look for: detrended screen 33.3 ns, Strehl 0.304. Noiseless: adjacent correlation 0.992, pair jitter 0.68 ns, residual 5.74 ns (NN) and 6.34 ns (LS), Strehl 0.965/0.956; after one transmit iteration 2.62/2.67 ns, Strehl 0.992. At 0 dB SNR: correlation 0.50, jitter 4.5–4.6 ns; pass 1 NN 10.5 and LS 7.9 ns; pass 2 NN 8.31 and LS 3.98 ns (Strehl 0.932 and 0.984), against a random-walk bound of 13.1 ns.

m33_distortion_matrix_svd.mReflection, focused and distortion matrices of a simulated speckle layer behind a screen; first singular vector, iteration, Γ, global versus local windows · change sig_tau, ell, d_scr, n_iter · runtime ≈ 2 s 
%% m33_distortion_matrix_svd.m
% Reflection matrix, focused basis and distortion matrix of a speckle medium
% behind a phase screen; aberration law from the first eigenvector, iteration,
% and the eigenvalue spectrum when the screen is away from the array.
% Purpose : Monochromatic (5 MHz) model, 48 elements, pitch 0.3 mm. Scatterers
%           (complex Gaussian, 400 of them, plus one bright point at x = 0) lie
%           in a thin layer at z = 20 mm, standing in for a time gate. A random
%           delay screen (rms sig_tau, correlation length l) sits at depth d;
%           the straight ray from element x_e to (x_p, z) crosses it at
%           x_e + (x_p - x_e) d/z. Green's function G = exp(-j k r + j phi).
%             K    = G diag(a) G.'                    (element-basis reflection matrix)
%             R_ux = K conj(G0) diag-corrected        (focused on transmit at x_in)
%             D    = R_ux .* conj(G0)                 (remove the geometric wavefront)
%             C    = D D'                             (= U S^2 U', SVD of D)
%             phi_hat = angle(first eigenvector of C), by power iteration.
%           Iteration: correct the input side with phi_hat (post-processing,
%           no new acquisition) and re-estimate. Focusing quality Gamma: the
%           confocal (zero-offset) fraction of the common-midpoint intensity
%           profile of R_xx, normalized by its value without the screen
%           (a Strehl-like ratio; not the FWHM-based F of Lambert et al.).
% Change  : sig_tau, ell, d_scr, Ns, n_iter.
% Expect  : d = 0: residual phase falls from ~1 rad to a few hundredths in
%           3 iterations, Gamma from ~0.5 to 1, the singular spectrum of D
%           collapses onto sigma_1; d = 8 mm (anisoplanatic): a slower singular
%           spectrum, and local windows reach a higher Gamma than one global law.
%           For d > 0 Gamma, not the residual against one point's law, is the
%           metric: a window law may focus well on slightly displaced points.
% Runtime : about 3 s in Octave.
clear; close all; rng(3);

c = 1540; f0 = 5e6; lam = c/f0; k = 2*pi/lam;
Ne = 48; p = 0.3e-3; xe = ((0:Ne-1) - (Ne-1)/2)*p; z = 20e-3;
xin = (-6e-3:0.2e-3:6e-3); Nx = numel(xin);
Ns = 400; sig_tau = 40e-9; ell = 2e-3; n_iter = 4;
xs_sc = [(rand(Ns, 1) - 0.5)*18e-3; 0]; a_sc = [(randn(Ns, 1) + 1j*randn(Ns, 1))/sqrt(2); 6];   % speckle + point at x = 0

dxs = 0.05e-3; xg = -30e-3:dxs:30e-3; gk = exp(-2*(((-3*ell):dxs:(3*ell))/ell).^2); gk = gk/sqrt(sum(gk.^2));
scr = 2*pi*f0*sig_tau*conv(randn(1, numel(xg)), gk, 'same');                    % phase screen phi(x) on the grid (rad)
phiAt = @(x) interp1(xg, scr, x, 'linear');
detr = @(v) v - [ones(numel(v), 1) xe(:)]*([ones(numel(v), 1) xe(:)]\v(:));
resid = @(v) detr(unwrap(angle(exp(1j*v(:)))));                                % wrapped residual -> unwrapped, mean and ramp removed
G0 = exp(-1j*k*sqrt((xe' - xin).^2 + z^2));                                     % Ne x Nx free-space
Ss = @(v) abs(mean(exp(1j*v)))^2;                                               % Strehl of a residual phase vector

for d_scr = [0 8e-3]
    fprintf('\n===== screen at d = %.0f mm (sigma_tau = %.0f ns -> %.2f rad rms, l = %.0f mm) =====\n', d_scr*1e3, sig_tau*1e9, 2*pi*f0*sig_tau, ell*1e3);
    G = exp(-1j*k*sqrt((xe' - xs_sc').^2 + z^2)).*exp(1j*phiAt(xe' + (xs_sc' - xe')*d_scr/z));   % Ne x (Ns+1)
    K = G*diag(a_sc)*G.';                                                        % reflection matrix, element basis
    phiTrue = @(xp) phiAt(xe' + (xp - xe')*d_scr/z);                             % true law for a focal point xp
    % ideal reference width w0 from the screen-free medium
    Gf = exp(-1j*k*sqrt((xe' - xs_sc').^2 + z^2)); Kf = Gf*diag(a_sc)*Gf.';
    dl = -2e-3:0.05e-3:2e-3;                                                     % output offsets for the CMP profile (m)
    cmp = @(Kc, cols) arrayfun(@(dd) mean(abs(sum(conj(exp(-1j*k*sqrt((xe' - xin(cols) - dd).^2 + z^2))).*(Kc*conj(G0(:, cols))), 1)).^2), dl);
    cfr = @(I) I(dl == 0)/sum(I);                                                % confocal fraction of the CMP energy
    w0 = cfr(cmp(Kf, 1:Nx));
    windows = {1:Nx}; Gfin = []; gain = [];                                                            % global window
    if d_scr > 0, windows = [windows, {1:21, 21:41, 41:61}]; end
    for wi = 1:numel(windows)
        cols = windows{wi}; xc = mean(xin(cols));
        psi = zeros(Ne, 1);                                                      % current correction (input and output)
        if wi == 1, fprintf(' window: all %d focal points\n', numel(cols)); else, fprintf(' window %d: x = %.1f to %.1f mm\n', wi - 1, xin(cols(1))*1e3, xin(cols(end))*1e3); end
        for it = 0:n_iter
            Kc = diag(exp(-1j*psi))*K*diag(exp(-1j*psi));                       % correction applied to both sides
            Rux = Kc*conj(G0(:, cols));                                          % focused on transmit
            Dm = Rux.*conj(G0(:, cols));                                         % distortion matrix
            C = Dm*Dm'; Ch = C;                                               % (an entrywise-normalized C./|C| is not positive semidefinite)
            v = ones(Ne, 1)/sqrt(Ne);                                            % power iteration
            for q = 1:60, v = Ch*v; v = v/norm(v); end
            lamC = sort(real(eig(Ch)), 'descend');
            Rxx = G0(:, cols)'*Kc*conj(G0(:, cols)); F = cfr(cmp(Kc, cols))/w0;
            sv = sqrt(max(lamC, 0)); sv = sv/sv(1);                              % singular values of D, normalized
            if d_scr == 0
                err = resid(phiTrue(xc) - psi);
                fprintf('  iter %d: residual phase %.3f rad rms (Strehl %.3f), Gamma = %.2f, sigma2/sigma1 = %.3f\n', it, std(err), Ss(err), F, sv(2));
            elseif it == 0 || it == n_iter
                fprintf('  iter %d: Gamma = %.2f, sigma2/sigma1 = %.3f\n', it, F, sv(2));
            end
            if it == 0 && wi == 1
                pn = lamC/sum(lamC); pn = pn(pn > 0); H = -sum(pn.*log(pn));
                fprintf('   singular values of D (sigma_i/sigma_1): %s ... effective rank exp(H) = %.2f\n', sprintf('%.3f ', sv(1:6)), exp(H));
                conf0 = abs(diag(Rxx)); Rxx0 = Rxx;
            end
            if it < n_iter, psi = psi + detr(unwrap(angle(v))); end               % update (mean and ramp removed: a ramp only shifts the image)
        end
        Gfin(wi) = F; i0 = find(abs(xin(cols)) < 1e-9);
        if ~isempty(i0), cf = abs(diag(Rxx)); gain(wi) = cf(i0); end
        if wi == 1, RxxN = Rxx; g0 = conf0(abs(xin) < 1e-9); end
    end
    fprintf(' bright point at x = 0: confocal amplitude gain, global law %.2f (%.1f dB)\n', gain(1)/g0, 20*log10(gain(1)/g0));
    if d_scr > 0
        fprintf(' bright point at x = 0: confocal amplitude gain, law of window 2 %.2f (%.1f dB)\n', gain(3)/g0, 20*log10(gain(3)/g0));
        fprintf(' final Gamma: global law %.2f, mean of the three local laws %.2f\n', Gfin(1), mean(Gfin(2:4)));
    end
    if d_scr == 0
        figure('Position', [50 50 1100 330]);
        subplot(1, 3, 1); imagesc(xin*1e3, xin*1e3, abs(Rxx0)); axis image; colormap(gray); xlabel('x_{in} (mm)'); ylabel('x_{out} (mm)'); title('|R_{xx}| before');
        subplot(1, 3, 2); imagesc(xin*1e3, xin*1e3, abs(RxxN)); axis image; xlabel('x_{in} (mm)'); title(sprintf('|R_{xx}| after %d iterations', n_iter));
        subplot(1, 3, 3); plot(xe*1e3, detr(phiTrue(0)), 'k', xe*1e3, detr(phiTrue(0)) - resid(phiTrue(0) - psi), 'r--'); xlabel('element (mm)'); ylabel('phase (rad)'); title('aberration law: true vs estimate'); legend('true', 'estimate');
    end
end

What to look for: screen at the array: residual 0.827, 0.314, 0.070, 0.018 rad (Strehl 0.495 → 1.000) and Γ 0.48, 0.87, 0.99, 1.00 over iterations 0–3; σ₂/σ₁ from 0.677 to 0.515; the bright point gains 2.8 dB. Screen at 8 mm: global Γ 0.48 → 0.69, windows 0.41 → 0.75, 0.65 → 0.93, 0.42 → 0.88 (mean 0.86); the bright point loses 7.6 dB with the global law and gains 3.1 dB with its local law.

m33_isoplanatic_patch.mOne law applied across the field for a screen at 5, 10 and 20 mm: Monte Carlo against the exact finite-aperture formula, patch widths and the crossover · change sig_tau, ell, d_list, z · runtime ≈ 25 s 
%% m33_isoplanatic_patch.m
% Isoplanatic patch of a phase screen away from the array: how far from the
% guide point one aberration law stays useful, by Monte Carlo against the
% exact Gaussian-screen formula and its large-aperture limit.
% Purpose : 128 elements (pitch 0.3 mm), 5 MHz, focal depth z = 40 mm. A thin
%           Gaussian delay screen (rms sig_tau, correlation rho(dx) =
%           exp(-dx^2/l^2)) lies at depth d. The straight ray from element x_u
%           to the point (x_p, z) crosses it at a_u + x_p d/z, a_u = x_u(1 - d/z).
%           The law measured at x_p = 0 is applied at x_p = dx: the residual is
%           r_u = phi(a_u + delta) - phi(a_u), delta = dx d / z.
%           Exact mean Strehl (Gaussian statistics):
%             E[S] = (1/N^2) sum_uv exp(-sig^2 [2 - 2 rho(delta) - 2 rho(da)
%                                        + rho(da + delta) + rho(da - delta)])
%           Large aperture: S ~ exp(-2 sig^2 (1 - rho(delta))).
%           Correction beats no correction while rho(delta) > 1/2, i.e.
%             |dx| < dx_x = l z sqrt(ln 2) / d.
% Change  : sig_tau, ell, d_list, z.
% Expect  : Monte Carlo on the exact curve; the patch shrinks as z/d; beyond
%           dx_x the borrowed law is worse than none.
% Runtime : about 25 s in Octave.
clear; close all; rng(4);

c = 1540; f0 = 5e6; N = 128; p = 0.3e-3; xe = ((0:N-1) - (N-1)/2)*p; z = 40e-3;
sig_tau = 40e-9; ell = 4e-3; sp = 2*pi*f0*sig_tau; d_list = [5 10 20]*1e-3; R = 200;
rho = @(x) exp(-(x/ell).^2);
dxs = 0.05e-3; xg = -40e-3:dxs:40e-3; g = exp(-2*(((-3*ell):dxs:(3*ell))/ell).^2); g = g/sqrt(sum(g.^2));
dx = [0:0.5:12, 13:1:30]*1e-3;                                                           % lateral distance from the guide point
fprintf('sigma_tau = %.0f ns (%.2f rad), l = %.0f mm, z = %.0f mm; uncorrected E[S] = %.3f\n', sig_tau*1e9, sp, ell*1e3, z*1e3, ...
    sum(sum(exp(-sp^2*(1 - rho(xe' - xe)))))/N^2);
fprintf('\n  d     dx_x = l z sqrt(ln2)/d (exact crossing)   dx(S=1/2) exact   MC S at 0, 2, 4, 8 mm         exact S at 0, 2, 4, 8 mm\n');
Smc = zeros(numel(d_list), numel(dx)); Sex = Smc;
for a = 1:numel(d_list)
    d = d_list(a); au = xe*(1 - d/z); DA = au' - au;
    for q = 1:numel(dx)
        del = dx(q)*d/z;
        Sex(a, q) = sum(sum(exp(-sp^2*(2 - 2*rho(del) - 2*rho(DA) + rho(DA + del) + rho(DA - del)))))/N^2;
    end
    acc = zeros(1, numel(dx));
    for r = 1:R
        phi = 2*pi*f0*sig_tau*conv(randn(1, numel(xg)), g, 'same');
        p0 = interp1(xg, phi, au);
        for q = 1:numel(dx)
            acc(q) = acc(q) + abs(mean(exp(1j*(interp1(xg, phi, au + dx(q)*d/z) - p0))))^2;
        end
    end
    Smc(a, :) = acc/R;
    S0 = sum(sum(exp(-sp^2*(1 - rho(DA)))))/N^2;                              % uncorrected at this d (screen compressed by 1 - d/z)
    half = interp1(Sex(a, :), dx, 0.5); ii = find(Sex(a, :) < S0, 1); cross = interp1(Sex(a, ii-1:ii), dx(ii-1:ii), S0);
    k = arrayfun(@(v) find(abs(dx - v) < 1e-9), [0 2 4 8]*1e-3);
    fprintf('%3.0f mm    %6.2f mm (%5.2f mm)                %6.2f mm         %s   %s  (uncorrected %.3f)\n', d*1e3, ell*z*sqrt(log(2))/d*1e3, cross*1e3, half*1e3, ...
        sprintf('%.3f ', Smc(a, k)), sprintf('%.3f ', Sex(a, k)), S0);
end
fprintf('\nLarge-aperture limit exp(-2 sig^2 (1 - rho)) at delta = l: %.3f; correction crosses the uncorrected level at rho = 1/2\n', exp(-2*sp^2*(1 - exp(-1))));
fprintf('Patch half-width at S = 1/2 (large aperture): dx = (z/d) l sqrt(-ln(1 - ln2/(2 sig^2))) = %.2f mm at d = 10 mm\n', ...
    z/10e-3*ell*sqrt(-log(1 - log(2)/(2*sp^2)))*1e3);

figure('Position', [50 50 600 380]); cols = 'brk';
hl = zeros(1, numel(d_list) + 1);
for a = 1:numel(d_list)
    hl(a) = plot(dx*1e3, Sex(a, :), [cols(a) '-']); hold on; plot(dx*1e3, Smc(a, :), [cols(a) 'o']);
end
hl(end) = plot(dx*1e3, exp(-sp^2)*ones(size(dx)), 'k:'); xlim([0 16]); xlabel('distance from the guide point (mm)'); ylabel('Strehl S');
title('One law, many points (lines exact, circles MC)'); legend(hl, 'd = 5 mm', 'd = 10 mm', 'd = 20 mm', 'exp(-\sigma^2)');

What to look for: uncorrected Strehl 0.315 (screen at the array). At d = 10 mm the Monte Carlo Strehl at 2, 4 and 8 mm is 0.952, 0.827, 0.517 against 0.953, 0.831, 0.522 exact; half-widths 16.52, 8.32 and 4.31 mm for d = 5, 10, 20 mm; crossovers 22.5, 11.0 and 5.1 mm against \(\Delta x_\times\) = 26.6, 13.3 and 6.7 mm; large-aperture half-width 7.96 mm at d = 10 mm.

m33_sos_map_correction.mWavy fat layer: Fermat travel times against one measured law, straight rays through the true, a blurred and a biased map · change c_fat, amp, per, blur_mm, dc_bias · runtime ≈ 4 s 
%% m33_sos_map_correction.m
% Correction from a sound-speed map (the Module 32 route): delays computed
% through the map versus one measured law, and the price of straight rays,
% a blurred map and a biased layer speed.
% Purpose : 128 elements (pitch 0.3 mm), 5 MHz, beamforming speed c0 = 1540
%           m/s. A fat layer (1450 m/s) runs from the skin to a wavy interface
%           h(x) = 12 + 1.5 sin(2 pi x / 15 mm) mm over 1540 m/s tissue.
%           True one-way times from element x_u to a point (x_p, z): Fermat's
%           principle with one refracting interface (1D search over the
%           crossing point). The aberration law is the excess over c0:
%             tau_u(x_p) = T_true(u, p) - r_up / c0.
%           Corrections compared at points x_p = -12 ... 12 mm, z = 40 mm:
%             none; one law measured at x_p = 0 and used everywhere;
%             straight-ray integration of the true map;
%             straight rays through a blurred map (Gaussian, 3 mm FWHM);
%             straight rays with the fat speed 10 m/s too high.
%           Residual = law - correction, mean and ramp removed; Strehl at 5 MHz.
% Change  : c_fat, amp, period, blur_mm, dc_bias, z.
% Expect  : the single law fails a few mm off the guide point; the
%           straight-ray map leaves a few ns (refraction ignored), the blurred
%           map more, the biased speed a smooth curvature error.
% Runtime : about 5 s in Octave.
clear; close all;

c0 = 1540; c_fat = 1450; f0 = 5e6; N = 128; p = 0.3e-3; xe = ((0:N-1) - (N-1)/2)*p;
amp = 1.5e-3; per = 15e-3; h = @(x) 12e-3 + amp*sin(2*pi*x/per);
z = 40e-3; xp_list = (-12:1:12)*1e-3; blur_mm = 3; dc_bias = 10;
xi = (-30e-3:0.02e-3:30e-3);                                                    % candidate crossing points
detr = @(v) v(:) - [ones(numel(v), 1) xe(:)]*([ones(numel(v), 1) xe(:)]\v(:));
Ss = @(v) abs(mean(exp(1j*2*pi*f0*v)))^2;

% gridded map for straight-ray integration
dg = 0.1e-3; xg = -30e-3:dg:30e-3; zg = 0:dg:z; [XG, ZG] = meshgrid(xg, zg);
s_true = 1/c0*ones(size(XG)); s_true(ZG < h(XG)) = 1/c_fat;
sig = blur_mm*1e-3/2.355; kx = -4*sig:dg:4*sig; gk = exp(-kx.^2/(2*sig^2)); gk = gk/sum(gk);
s_blur = conv2(conv2(s_true, gk, 'same'), gk', 'same');                         % blurred map (tomography-like resolution)
s_blur(:, [1:40, end-39:end]) = s_true(:, [1:40, end-39:end]);                  % avoid zero-padding at the lateral edges
s_blur(1:40, :) = s_true(1:40, :); s_blur(end-39:end, :) = s_true(end-39:end, :);
s_bias = 1/c0*ones(size(XG)); s_bias(ZG < h(XG)) = 1/(c_fat + dc_bias);
tt = linspace(0, 1, 801)';                                                      % samples along each straight ray
ray = @(S, xp) sum(interp2(xg, zg, S, xe + (xp - xe).*tt, z*tt*ones(1, N), 'linear'))/numel(tt).*hypot(xp - xe, z);   % straight-ray times, all elements

law = @(xp) arrayfun(@(xu) min(hypot(xi - xu, h(xi))/c_fat + hypot(xp - xi, z - h(xi))/c0), xe) - hypot(xp - xe, z)/c0;   % Fermat
L0 = law(0);
res = zeros(numel(xp_list), 5);
for q = 1:numel(xp_list)
    xp = xp_list(q); Lt = law(xp);
    ts = ray(s_true, xp) - hypot(xp - xe, z)/c0;
    tb = ray(s_blur, xp) - hypot(xp - xe, z)/c0;
    tbi = ray(s_bias, xp) - hypot(xp - xe, z)/c0;
    R = [detr(Lt), detr(Lt - L0), detr(Lt - ts), detr(Lt - tb), detr(Lt - tbi)];
    res(q, :) = std(R);
    for m = 1:5, St(q, m) = Ss(R(:, m)); end
end
names = {'none', 'one law (x=0)', 'straight-ray map', 'blurred map', 'biased fat speed'};
fprintf('Fat 1450 m/s to h(x) = 12 + 1.5 sin(2 pi x/15) mm, z = 40 mm, 5 MHz, 128 elements\n');
fprintf('Residual rms (ns, mean and ramp removed) / Strehl\n   x_p   ');
fprintf('%-20s', names{:}); fprintf('\n');
for q = 1:3:numel(xp_list)
    fprintf('%5.0f mm ', xp_list(q)*1e3);
    for m = 1:5, fprintf('%6.1f / %.2f        ', res(q, m)*1e9, St(q, m)); end
    fprintf('\n');
end
fprintf('\nMean over x_p = -12..12 mm:  '); for m = 1:5, fprintf('%s %.1f ns (S %.2f);  ', names{m}, mean(res(:, m))*1e9, mean(St(:, m))); end; fprintf('\n');
ii = find(St(:, 2) >= 0.8); fprintf('One law keeps S >= 0.8 for |x_p| <= %.0f mm\n', max(abs(xp_list(ii)))*1e3);

figure('Position', [50 50 900 330]);
subplot(1, 2, 1); imagesc(xg*1e3, zg*1e3, 1./s_true); axis image; colorbar; hold on; plot([xe(1) 0 xe(end)]*1e3, [0 z 0]*1e3, 'w-'); xlabel('x (mm)'); ylabel('z (mm)'); title('sound-speed map (m/s)');
subplot(1, 2, 2); plot(xp_list*1e3, St(:, 1), 'k', xp_list*1e3, St(:, 2), 'r', xp_list*1e3, St(:, 3), 'b', xp_list*1e3, St(:, 4), 'b--', xp_list*1e3, St(:, 5), 'g');
xlabel('lateral position of the point (mm)'); ylabel('Strehl after correction'); legend(names, 'Location', 'south'); ylim([0 1.05]);

What to look for: mean over ±12 mm: none 45.0 ns (Strehl 0.08), one law 30.3 ns (0.43; S ≥ 0.8 only within ±3 mm), straight-ray map 5.2 ns (0.97), blurred map 8.2 ns (0.93), biased fat speed 7.3 ns (0.95). The straight-ray residual is 3.4 ns at the centre and 11.1 ns at −12 mm: the path error alone.

33.13Self-check

Q1A 25 ns rms screen at 7.5 MHz. σφ and the large-aperture Strehl?

\(\sigma_\phi = 2\pi \times 7.5\times10^{6} \times 25\times10^{-9}\) = 1.18 rad; \(e^{-1.388}\) = 0.25 (higher with a finite aperture).

Q2Nearest-neighbour pair jitter 2 ns on a 192-element array, independent pairs. rms error after removing mean and ramp?

\(2\sqrt{192/15}\) = 7.2 ns.

Q3Screen at 15 mm, ℓ = 5 mm, target depth 60 mm. Beyond what lateral distance is a borrowed law worse than none (large aperture)?

\(\Delta x_\times = \ell z\sqrt{\ln 2}/d\) = 5 × 60 × 0.833 / 15 = 16.7 mm.

Q4Same geometry, σφ = 1.5 rad. Isoplanatic half-width at S = 1/2?

\(\ln 2/(2 \times 2.25)\) = 0.154; \(\sqrt{-\ln(1 - 0.154)}\) = 0.409; × \(z\ell/d\) = 20 mm gives 8.2 mm.

Q5Power iteration on a correlation matrix with λ2/λ1 = 0.5. Iterations to reach an error of 10−3?

\(n = \ln 10^{-3}/\ln 0.5\) = 10.

Q630 ns rms at 5 MHz, the same screen on transmit and receive. Large-aperture two-way Strehl?

One way \(e^{-0.888}\) = 0.41; two way about \(0.41^{2}\) = 0.17.

Q760 ns, ℓ = 4 mm, focus at 40 mm. Predicted −6 dB width of the mean PSF?

\(s = \sqrt2 \times 1540 \times 60\times10^{-9} \times 0.04/0.004\) = 1.31 mm; 3.33 s = 4.35 mm (the script measures 3.50 mm by Monte Carlo and 3.62 mm exact: at σφ = 1.9 rad a remnant of the core still narrows the −6 dB width).

33.14Figures from the literature

Measured data behind the models of this module. Figure numbers follow the published versions (checked against arXiv 1912.05417v2, 1911.03147 and 2103.02036v3, which match them); links only, no images reproduced.

The distortion matrix on a phantom behind a plexiglass aberrator

W. Lambert, L. A. Cobus, T. Frappart, M. Fink and A. Aubry, "Distortion matrix approach for ultrasound imaging of human soft tissues," Proc. Natl. Acad. Sci. U.S.A. 117(26), 14645–14656 (2020); preprint arXiv:1912.05417 (under the title "… of random scattering media"). doi:10.1073/pnas.1921533117

Look at Fig. 4 (time-reversal reading of D: each focused input gives a tilted wavefront in the far field, removing the tilt gives D, its correlation matrix acts as a virtual reflector shaped like the focal spot, and normalization makes it point-like), Fig. 5 (normalized eigenvalues of the correlation matrix, phases of the first three eigenvectors, entropy versus model speed, and confocal and Strehl maps corrected with each of the three eigenvectors) and Supplementary Fig. S2 (phases of the reflection matrix, of the free-space propagation matrix and of the distortion matrix at 30 mm: the hidden correlations appear only in D), with Fig. 2 for the images and Strehl maps before and after. Why: it is Demo 33.6 on measured data, in the plane-wave receive basis: 15 mm of plexiglass (about 2750 m/s) on a 1542 m/s phantom, entropy H ≈ 2.85 hence three isoplanatic modes, a full-field Strehl of about 0.4 and a 14 dB contrast gain; Fig. 6 repeats it in vivo on a calf.

The focused reflection matrix: focusing criterion, sound speed and multiple scattering

W. Lambert, L. A. Cobus, M. Couade, M. Fink and A. Aubry, "Reflection matrix approach for quantitative imaging of scattering media," Phys. Rev. X 10, 021048 (2020); preprint arXiv:1911.03147. doi:10.1103/PhysRevX.10.021048

Look at Fig. 2 (Rxx at 18 and 30 mm under bovine tissue, the focusing criterion F versus the assumed speed, and the optimized speeds of 1573 and 1546 m/s), Fig. 3 (maps of F for the homogeneous and the two-layer models) and Fig. 6 (in vivo liver with a four-layer speed model and a map of the multiple-scattering ratio). Why: Fig. 2 is (33.10) on real data and the link to Module 32: the same matrix measures focus, speed and multiple scattering.

Local laws over many patches, in vivo

W. Lambert, L. A. Cobus, J. Robin, M. Fink and A. Aubry, "Ultrasound matrix imaging. Part II: The distortion matrix for aberration correction over multiple isoplanatic patches," IEEE Trans. Med. Imaging (2022). doi:10.1109/TMI.2022.3199483; preprint arXiv:2103.02036

Look at Fig. 1 (Rxx of a gallbladder at 39 mm: original, after the confocal filter, after the first and the fourth correction), Fig. 6 (the correlation matrix over one area, its normalized eigenvalues, the first two eigenvectors and the phase law before and after ramp removal), Fig. 7 (focusing factor versus depth for several correction schemes) and Fig. 8 (numerical validation with true and estimated laws). Why: Fig. 6 shows the ramp ambiguity of Section 33.6 explicitly, and Fig. 7 shows where the correction saturates with depth.

DORT: one eigenvalue per scatterer, one focusing law per eigenvector

C. Prada and M. Fink, "Eigenmodes of the time reversal operator: a solution to selective focusing in multiple-target media," Wave Motion 20, 151–163 (1994). author PDF (Institut Langevin)

Look at Fig. 10 (the 96 eigenvalues of the time-reversal operator for four wires: four significant) and Fig. 11 (phase laws of the first four eigenvectors and the fields they produce, each focused on one wire); Figs. 14–16 repeat it for two spheres. Why: the measured version of (33.7) and of Module 06's DORT rule.

Per-patch aberration laws from microbubbles, through the human skull

C. Demené, J. Robin, A. Dizeux, B. Heiles, M. Pernot, M. Tanter and F. Perren, "Transcranial ultrafast ultrasound localization microscopy of brain vasculature in patients," Nat. Biomed. Eng. 5, 219–228 (2021). doi:10.1038/s41551-021-00697-x

Look at the aberration-correction comparison in the Methods and Supplementary Information (ULM with and without correction, with the microbubble echoes used to derive the laws); the supplementary figure number could not be verified for this entry, so none is quoted. Why: the isoplanatic patch of Section 33.4 and the guide-star logic of Section 33.3, applied in patients, and the reason Module 31 corrects before localizing.

33.15Go deeper

  • S. W. Flax and M. O'Donnell, "Phase-aberration correction using signals from point reflectors and diffuse scatterers: basic principles," IEEE Trans. Ultrason. Ferroelectr. Freq. Control (1988). doi:10.1109/58.9333; and the companion "…: measurements" by M. O'Donnell and S. W. Flax (1988), doi:10.1109/58.9334. Nearest-neighbour correlation (Section 33.2).
  • L. Nock, G. E. Trahey and S. W. Smith, "Phase aberration correction in medical ultrasound using speckle brightness as a quality factor," J. Acoust. Soc. Am. 85(5), 1819–1833 (1989). doi:10.1121/1.397889. D. Zhao and G. E. Trahey, "Comparisons of image quality factors for phase aberration correction with diffuse and point targets: theory and experiments," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 38(2), 125–132 (1991). doi:10.1109/58.68469.
  • R. Mallart and M. Fink, "The van Cittert–Zernike theorem in pulse echo measurements," J. Acoust. Soc. Am. 90(5), 2718–2727 (1991), and "Adaptive focusing in scattering media through sound-speed inhomogeneities: the van Cittert Zernike approach and focusing criterion," J. Acoust. Soc. Am. 96(6), 3721–3732 (1994). The virtual-source reading of speckle covariance and the coherence-based focusing criterion (Section 33.3).
  • T. Varslot, E. Mo, B. Angelsen and H. Krogstad, "Eigenfunction analysis of stochastic backscatter for characterization of acoustic aberration in medical ultrasound imaging," J. Acoust. Soc. Am. 115(6), 3068–3076 (2004); S.-E. Måsøy, T. Varslot and B. Angelsen, "Iteration of transmit-beam aberration correction in medical ultrasound imaging," J. Acoust. Soc. Am. 117(1), 450–461 (2005): close to ideal correction after two to three transmit-beam iterations for all eight aberrators tested.
  • S. Kang et al., "High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering," Nat. Commun. 8 (2017). doi:10.1038/s41467-017-02117-8. CLASS, the optical ancestor of alternating input/output correction.
  • F. Bureau, J. Robin, A. Le Ber, W. Lambert, M. Fink and A. Aubry, "Three-dimensional ultrasound matrix imaging," Nat. Commun. 14 (2023). doi:10.1038/s41467-023-42338-8.
  • 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 69(5), 1714–1726 (2022). doi:10.1109/TUFFC.2022.3162836. The map route of Section 33.5.
  • Lab cross-links: phase screen as a diagonal matrix, focused basis (Module 05); SVD, power iteration, DORT (Module 06); van Cittert–Zernike (Module 07); delay estimation (Module 11); CRLB (Module 12); least squares and regularization (Module 13); autofocus, Wirtinger gradients, cycle skipping (Module 14); gCNR (Module 26); plane-wave compounding (Module 23); CF, SLSC (Module 25); full-wave simulators for validation (Module 21); refraction (Module 15); ULM bias field and transcranial ULM (Module 31); sound-speed maps (Module 32).